common.H 15.4 KB
Newer Older
greole's avatar
Ir (#1)  
greole committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
/*---------------------------------------------------------------------------*\
License
    This file is part of OGL.

    OGL is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Class
    Foam::GKOCG

Author: Gregor Olenik <go@hpsim.de>

SourceFiles
    GKOCG.C

\*---------------------------------------------------------------------------*/
#ifndef OGL_COMMON_H
#define OGL_COMMON_H

#include "fvCFD.H"
#include "regIOobject.H"

greole's avatar
greole committed
33
#include "../IOExecutorHandler/IOExecutorHandler.H"
greole's avatar
greole committed
34
#include "../IOGKOMatrixHandler/IOGKOMatrixHandler.H"
greole's avatar
greole committed
35
36
37
#include "../IOHandler/IOPreconditioner/IOPreconditioner.H"
#include "../IOHandler/IOSortingIdxHandler/IOSortingIdxHandler.H"
#include "../common/StoppingCriterion.H"
greole's avatar
greole committed
38

greole's avatar
Ir (#1)  
greole committed
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include <ginkgo/ginkgo.hpp>
#include <map>

#define SIMPLE_TIME(NAME, F)                                            \
    auto start_##NAME = std::chrono::steady_clock::now();               \
    F auto end_##NAME = std::chrono::steady_clock::now();               \
    std::cout << "Gingko " #NAME " : "                                  \
              << std::chrono::duration_cast<std::chrono::milliseconds>( \
                     end_##NAME - start_##NAME)                         \
                     .count()                                           \
              << " ms\n";

namespace Foam {

greole's avatar
greole committed
53
54
template <class MatrixType>
class HostMatrix : public MatrixType::solver {
55
private:
greole's avatar
Ir (#1)  
greole committed
56
57
58
59
60
61
    const label nCells_;

    const label nNeighbours_;

    const label nElems_;

greole's avatar
greole committed
62
    mutable std::vector<scalar> values_;
greole's avatar
Ir (#1)  
greole committed
63

greole's avatar
greole committed
64
    mutable std::vector<label> col_idxs_;
greole's avatar
Ir (#1)  
greole committed
65

greole's avatar
greole committed
66
    mutable std::vector<label> row_idxs_;
greole's avatar
Ir (#1)  
greole committed
67
68

public:
greole's avatar
greole committed
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
    HostMatrix(const word &fieldName, const MatrixType &matrix,
               const FieldField<Field, scalar> &interfaceBouCoeffs,
               const FieldField<Field, scalar> &interfaceIntCoeffs,
               const lduInterfaceFieldPtrsList &interfaces,
               const dictionary &solverControls)
        : MatrixType::solver(fieldName, matrix, interfaceBouCoeffs,
                             interfaceIntCoeffs, interfaces, solverControls),
          nCells_(matrix.diag().size()),
          nNeighbours_(matrix.lduAddr().upperAddr().size()),
          nElems_(nCells_ + 2 * nNeighbours_){};

    HostMatrix(const word &fieldName, const MatrixType &matrix,
               const dictionary &solverControls)
        : MatrixType::solver(fieldName, matrix, solverControls),
          nCells_(matrix.diag().size()),
          nNeighbours_(matrix.lduAddr().upperAddr().size()),
          nElems_(nCells_ + 2 * nNeighbours_){};

87
    void init_host_sparsity_pattern() const
greole's avatar
Ir (#1)  
greole committed
88
    {
greole's avatar
greole committed
89
90
91
92
93
94
95
96
97
98
        col_idxs_.resize(0);
        row_idxs_.resize(0);

        col_idxs_.reserve(nElems());
        row_idxs_.reserve(nElems());

        // fill vectors unsorted
        for (label i = 0; i < nNeighbours(); ++i) {
            row_idxs_.push_back(this->matrix().lduAddr().lowerAddr()[i]);
            col_idxs_.push_back(this->matrix().lduAddr().upperAddr()[i]);
greole's avatar
Ir (#1)  
greole committed
99
100
        }

greole's avatar
greole committed
101
102
103
        for (label i = 0; i < nCells(); ++i) {
            col_idxs_.push_back(i);
            row_idxs_.push_back(i);
greole's avatar
Ir (#1)  
greole committed
104
105
        }

greole's avatar
greole committed
106
107
108
        for (label i = 0; i < nNeighbours(); ++i) {
            row_idxs_.push_back(this->matrix().lduAddr().upperAddr()[i]);
            col_idxs_.push_back(this->matrix().lduAddr().lowerAddr()[i]);
greole's avatar
Ir (#1)  
greole committed
109
        }
110
111
    }

112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
    void sort_host_matrix_sparsity_pattern(
        const IOField<label> *sorting_idxs) const
    {
        std::vector<label> tmp_col_idxs(nElems());
        std::vector<label> tmp_row_idxs(nElems());

        for (label i = 0; i < nElems(); i++) tmp_col_idxs[i] = col_idxs_[i];
        for (label i = 0; i < nElems(); i++) tmp_row_idxs[i] = row_idxs_[i];

        for (label i = 0; i < nElems(); i++) {
            label j = sorting_idxs->operator[](i);
            col_idxs_[i] = tmp_col_idxs[j];
            row_idxs_[i] = tmp_row_idxs[j];
        }
    };


129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
    void update_host_matrix_data() const
    {
        // reset vectors
        values_.resize(0);

        values_.reserve(nElems());

        // fill vectors unsorted
        for (label i = 0; i < nNeighbours(); ++i) {
            values_.push_back(this->matrix().lower()[i]);
        }

        for (label i = 0; i < nCells(); ++i) {
            values_.push_back(this->matrix().diag()[i]);
        }

        for (label i = 0; i < nNeighbours(); ++i) {
            values_.push_back(this->matrix().upper()[i]);
        }
greole's avatar
greole committed
148
    };
greole's avatar
Ir (#1)  
greole committed
149

150
151
152
153
154
155
156
157
158
159
160
161
    void sort_host_matrix_data(const IOField<label> *sorting_idxs) const
    {
        std::vector<scalar> tmp_values(nElems());

        for (label i = 0; i < nElems(); i++) tmp_values[i] = values_[i];

        for (label i = 0; i < nElems(); i++) {
            label j = sorting_idxs->operator[](i);
            values_[i] = tmp_values[j];
        }
    };

162
    label nCells() const { return nCells_; };
greole's avatar
Ir (#1)  
greole committed
163

164
    label nElems() const { return nElems_; };
greole's avatar
Ir (#1)  
greole committed
165

166
    label nNeighbours() const { return nNeighbours_; };
greole's avatar
Ir (#1)  
greole committed
167

greole's avatar
greole committed
168
    std::vector<scalar> &values() const { return values_; };
greole's avatar
Ir (#1)  
greole committed
169

greole's avatar
greole committed
170
    std::vector<label> &col_idxs() const { return col_idxs_; };
greole's avatar
Ir (#1)  
greole committed
171

greole's avatar
greole committed
172
173
    std::vector<label> &row_idxs() const { return row_idxs_; };
};
greole's avatar
Ir (#1)  
greole committed
174

greole's avatar
greole committed
175
void export_system(const word fieldName, const mtx *A, const vec *x,
176
                   const vec *b, const word time);
greole's avatar
Ir (#1)  
greole committed
177
178


greole's avatar
greole committed
179
180
181
182
183
184
185
// Logs the number of iteration executed
struct IterationLogger : gko::log::Logger {
    void on_iteration_complete(const gko::LinOp *,
                               const gko::size_type &num_iterations,
                               const gko::LinOp *residual,
                               const gko::LinOp *res_norm,
                               const gko::LinOp *) const override
greole's avatar
Ir (#1)  
greole committed
186
    {
greole's avatar
greole committed
187
188
        this->num_iters = num_iterations;
    }
greole's avatar
Ir (#1)  
greole committed
189

greole's avatar
greole committed
190
191
192
193
    IterationLogger(std::shared_ptr<const gko::Executor> exec)
        : exec(exec),
          gko::log::Logger(exec, gko::log::Logger::iteration_complete_mask)
    {}
greole's avatar
Ir (#1)  
greole committed
194

greole's avatar
greole committed
195
    gko::size_type get_iters() { return num_iters; }
greole's avatar
Ir (#1)  
greole committed
196

greole's avatar
greole committed
197
198
199
200
private:
    std::shared_ptr<const gko::Executor> exec;
    mutable gko::size_type num_iters{0};
};
greole's avatar
Ir (#1)  
greole committed
201
202


greole's avatar
greole committed
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
template <class MatrixType, class SolverFactory>
class lduLduBase : public HostMatrix<MatrixType>,
                   public SolverFactory,
                   public StoppingCriterion,
                   public IOSortingIdxHandler,
                   public IOGKOMatrixHandler,
                   public IOPreconditioner {
public:
    lduLduBase(const word &fieldName, const lduMatrix &matrix,
               const FieldField<Field, scalar> &interfaceBouCoeffs,
               const FieldField<Field, scalar> &interfaceIntCoeffs,
               const lduInterfaceFieldPtrsList &interfaces,
               const dictionary &solverControls)
        : HostMatrix<lduMatrix>(fieldName, matrix, interfaceBouCoeffs,
                                interfaceIntCoeffs, interfaces, solverControls),
          SolverFactory(solverControls),
          StoppingCriterion(solverControls),
          IOSortingIdxHandler(
              matrix.mesh().thisDb(), this->nElems(),
              solverControls.lookupOrDefault<Switch>("sort", true)),
          IOGKOMatrixHandler(matrix.mesh().thisDb(), solverControls, fieldName),
          IOPreconditioner(matrix.mesh().thisDb(), solverControls, fieldName)
    {
        init_base();
    }


    //- Construct from matrix components and solver controls
    lduLduBase(const word &fieldName, const MatrixType &matrix,
               const dictionary &solverControls)
        : HostMatrix<MatrixType>(fieldName, matrix, solverControls),
          SolverFactory(solverControls),
          StoppingCriterion(solverControls),
          IOSortingIdxHandler(
              matrix.mesh().thisDb(), this->nElems(),
              solverControls.lookupOrDefault<Switch>("sort", true)),
          IOGKOMatrixHandler(matrix.mesh().thisDb(), solverControls, fieldName),
          IOPreconditioner(matrix.mesh().thisDb(), solverControls, fieldName)
    {
        init_base();
    }

    void init_base()
    {
        // if sys_matrix is not stored updating is neccesary
        // initially
        bool stored = get_sys_matrix_stored();
        if (!stored) {
            SIMPLE_TIME(init_host_sparsity_pattern,
                        this->init_host_sparsity_pattern();)
            std::cout << " matrix not stored update host matrix " << std::endl;
            SIMPLE_TIME(update_host_mtx, this->update_host_matrix_data();)
        } else {
            // if sys_matrix is  stored updating is only neccesary
            // when requested explictly
            if (get_update_sys_matrix()) {
                std::cout
                    << " matrix is stored and  update of host matrix requested"
                    << std::endl;
                SIMPLE_TIME(exp_update_host_mtx,
                            this->update_host_matrix_data();)
            }
        }

        // after updating the host matrix the host matrix needs to be sorted
        if (!get_is_sorted()) {
            std::cout << "matrix is not yet sorted, compute sorting idxs "
                      << std::endl;
            SIMPLE_TIME(compute_sort, this->compute_sorting_idxs(
                                          this->row_idxs(), this->col_idxs(),
                                          this->nCells());)
        }

        if (!stored && get_sort()) {
            std::cout << "matrix is not yet sorted,  sorting host matrix "
                      << std::endl;
            SIMPLE_TIME(
                sort_host_mtx_sparsity_pattern,
                this->sort_host_matrix_sparsity_pattern(get_sorting_idxs());)
            SIMPLE_TIME(sort_host_mtx,
                        this->sort_host_matrix_data(this->get_sorting_idxs());)
        }

        init_device_matrix(this->matrix().mesh().thisDb(), this->values(),
                           this->col_idxs(), this->row_idxs(), this->nElems(),
                           this->nCells(), !get_update_sys_matrix());
    }


    template <class OFField>
    void init_vector_views(std::vector<val_array> &target,
                           OFField &source) const
    {
        for (int i = 0; i < 3; i++) {
            target.push_back(
                val_array::view(ref_exec(), 3 * this->nCells(),
                                const_cast<scalar *>(&source[0][i])));
        }
    }

    template <class Type>
    SolverPerformance<Type> solve_impl_(Field<Type> &psi) const
    {
        std::shared_ptr<gko::Executor> device_exec =
            this->get_device_executor();
        // --- Setup class containing solver performance data
        // Implement
        word preconditionerName(this->controlDict_.lookup("preconditioner"));

312
313
314
        SolverPerformance<Type> solverPerf(
            preconditionerName + this->get_device_executor_name(),
            this->fieldName_);
greole's avatar
greole committed
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361

        std::vector<val_array> source_views{};

        init_vector_views(source_views, this->matrix().source());

        std::vector<std::shared_ptr<vec>> b{};
        for (int i = 0; i < 3; i++) {
            b.push_back(vec::create(ref_exec(), gko::dim<2>(this->nCells(), 1),
                                    source_views[i], 3));
        }

        init_initial_guess_vector(psi, this->matrix().mesh().thisDb(),
                                  this->nCells());
        std::vector<std::shared_ptr<vec>> x{};
        this->get_initial_guess(x);

        std::vector<std::shared_ptr<const gko::stop::CriterionFactory>>
            criterion_vec{};

        criterion_vec.push_back(build_stopping_criterion(device_exec, 1.0));

        // Generate solver
        auto solver_gen = this->create_solver(device_exec, criterion_vec);

        std::shared_ptr<mtx> gkomatrix = get_gkomatrix();

        auto solver = solver_gen->generate(gko::share(gkomatrix));

        for (int i = 0; i < 3; i++) {
            SIMPLE_TIME(solve, solver->apply(gko::lend(b[i]), gko::lend(x[i]));)
        }

        // copy back
        copy_result_back(psi, this->nCells());

        return solverPerf;
    };

    solverPerformance solve_impl_(word typeName, scalarField &psi,
                                  const scalarField &source,
                                  const direction cmpt = 0) const
    {
        std::shared_ptr<gko::Executor> device_exec =
            this->get_device_executor();

        // --- Setup class containing solver performance data
        solverPerformance solverPerf(
362
363
            lduMatrix::preconditioner::getName(this->controlDict_) +
                this->get_device_executor_name() + typeName,
greole's avatar
greole committed
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
            this->fieldName());


        scalarField pA(this->nCells());
        scalarField wA(this->nCells());
        this->matrix().Amul(wA, psi, this->interfaceBouCoeffs_,
                            this->interfaces_, cmpt);
        scalar norm_factor = this->normFactor(psi, source, wA, pA);

        auto source_view = val_array::view(ref_exec(), this->nCells(),
                                           const_cast<scalar *>(&source[0]));

        std::vector<val_array> source_views{source_view};
        auto b = vec::create(ref_exec(), gko::dim<2>(this->nCells(), 1),
                             source_views[0], 1);

        init_initial_guess_vector(psi, this->matrix().mesh().thisDb(),
                                  this->nCells());
        std::vector<std::shared_ptr<vec>> x{};
        this->get_initial_guess(x);

        std::vector<std::shared_ptr<const gko::stop::CriterionFactory>>
            criterion_vec{};

        criterion_vec.push_back(
            build_stopping_criterion(device_exec, norm_factor));

        // Generate solver
        std::shared_ptr<mtx> gkomatrix = get_gkomatrix();
        SIMPLE_TIME(init_preconditioner,
                    this->init_preconditioner(this->matrix().mesh().thisDb(),
                                              gkomatrix, device_exec);)

        auto bj_generated = this->get_preconditioner();
        auto solver_gen =
            this->create_solver(device_exec, criterion_vec, bj_generated);

        // Instantiate a ResidualLogger logger.
        auto logger = std::make_shared<IterationLogger>(device_exec);

        // Add the previously created logger to the solver factory. The
        // logger will be automatically propagated to all solvers created
        // from this factory.
        solver_gen->add_logger(logger);

        // auto b_clone = gko::clone(b);


        if (get_export())
            export_system(this->fieldName(), gko::lend(gkomatrix),
414
415
                          gko::lend(x[0]), gko::lend(b),
                          this->matrix().mesh().thisDb().time().timeName());
greole's avatar
greole committed
416
417


418
        solverPerf.initialResidual() = this->get_init_res_norm();
greole's avatar
greole committed
419
420
421
422
423
424
425
426
427
428
429

        auto solver = solver_gen->generate(gko::share(gkomatrix));

        // Solve system
        SIMPLE_TIME(solve, solver->apply(gko::lend(b), gko::lend(x[0]));)
        std::cout << " copy back" << std::endl;

        this->copy_result_back(psi, this->nCells());

        // b_clone->copy_from(gko::lend(b));

430
        solverPerf.finalResidual() = this->get_res_norm();
greole's avatar
greole committed
431
432
433
434
435
436
437
438
439

        std::cout << " get iters " << std::endl;
        solverPerf.nIterations() = logger->get_iters();

        return solverPerf;
    };
};


greole's avatar
Ir (#1)  
greole committed
440
441
442
}  // namespace Foam

#endif