ProblemStat.inc.hpp 16.8 KB
Newer Older
1
2
#pragma once

3
4
5
#include <dune/amdis/AdaptInfo.hpp>
#include <dune/amdis/common/Loops.hpp>
#include <dune/amdis/common/Timer.hpp>
6
7
8

namespace AMDiS
{
9
10
11
  template <class Traits>
  void ProblemStatSeq<Traits>::initialize(Flag initFlag,
                                          Self* adoptProblem,
12
                                          Flag adoptFlag)
13
14
15
  {
    // === create meshes ===
    if (mesh) {
16
      warning("mesh already created");
17
18
19
20
    }
    else {
      if (initFlag.isSet(CREATE_MESH) ||
          (!adoptFlag.isSet(INIT_MESH) &&
21
          (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) {
22
23
        createMesh();
      }
24

25
26
27
      if (adoptProblem &&
          (adoptFlag.isSet(INIT_MESH) ||
          adoptFlag.isSet(INIT_SYSTEM) ||
28
          adoptFlag.isSet(INIT_FE_SPACE))) {
29
30
        mesh = adoptProblem->getMesh();
      }
31

32
33
      componentMeshes.resize(nComponents, mesh.get());
    }
34

35
    if (!mesh)
36
37
38
      warning("no mesh created");


39
40
    int globalRefinements = 0;
    Parameters::get(meshName + "->global refinements", globalRefinements);
41
42
    if (globalRefinements > 0) {
      mesh->globalRefine(globalRefinements);
43
44
45
46
47
      msg("After refinement:");
      msg("#elements = "   , mesh->size(0));
      msg("#faces/edges = ", mesh->size(1));
      msg("#vertices = "   , mesh->size(dim));
      msg("");
48
    }
49
50
51



52
53
    // === create fespace ===
    if (feSpaces) {
54
      warning("feSpaces already created");
55
56
57
    }
    else {
      if (initFlag.isSet(INIT_FE_SPACE) ||
58
          (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
59
60
        createFeSpaces();
      }
61

62
      if (adoptProblem &&
63
          (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
64
65
        feSpaces = adoptProblem->getFeSpaces();
      }
66
    }
67
68

    if (!feSpaces)
69
70
71
72
      warning("no feSpaces created\n");

    msg("#DOFs = ", this->getFeSpace<0>().size());

73
74
75
76
    // === create system ===
    if (initFlag.isSet(INIT_SYSTEM))
      createMatricesAndVectors();

77
    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
78
79
80
      solution = adoptProblem->getSolution();
      rhs = adoptProblem->getRhs();
      systemMatrix = adoptProblem->getSystemMatrix();
81
82
    }

83
    // === create solver ===
84
    if (linearSolver) {
85
      warning("solver already created\n");
86
    }
87
    else {
88
89
90
      if (initFlag.isSet(INIT_SOLVER))
        createSolver();

91
      if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
92
        test_exit(!linearSolver, "solver already created\n");
93
94
        linearSolver = adoptProblem->getSolver();
      }
95
96
    }

97
    if (!linearSolver) {
98
      warning("no solver created\n");
99
    }
100
101
102
103

    // === create file writer ===
    if (initFlag.isSet(INIT_FILEWRITER))
      createFileWriter();
104

105
106
107
108
109
    solution->compress();
  }


  template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
110
  void ProblemStatSeq<Traits>::addMatrixOperator(std::shared_ptr<OperatorType> op,
111
                                                 int i, int j,
112
                                                 double* factor,
113
114
115
                                                 double* estFactor)
  {
    // TODO: currently the factors are ignored
116
    assert( estFactor == NULL );
117

118
    matrixOperators[{i,j}].push_back(op);
119
    matrixFactors[{i,j}].push_back(factor);
120
    matrixChanging[{i,j}] = true;
121
122
  }

123
124

  template <class Traits>
125
  void ProblemStatSeq<Traits>::addMatrixOperator(OperatorType& op,
126
                                                 int i, int j,
127
                                                 double* factor,
128
129
                                                 double* estFactor)
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
130
    addMatrixOperator(std::make_shared<OperatorType>(op), i, j, factor, estFactor);
131
132
133
  }

  template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
134
  void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, std::shared_ptr<OperatorType> > ops)
135
136
137
138
  {
    for (auto op : ops)
      addMatrixOperator(op.second, op.first.first, op.first.second);
  }
139
140

  template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
141
  void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, std::pair<std::shared_ptr<OperatorType>,bool> > ops)
142
143
144
145
146
147
148
149
  {
    for (auto op : ops) {
      auto row_col = op.first;
      matrixOperators[row_col].push_back(op.second.first);
      matrixFactors[row_col].push_back(NULL);
      matrixChanging[row_col] = matrixChanging[row_col] || op.second.second;
    }
  }
150

151
  template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
152
  void ProblemStatSeq<Traits>::addVectorOperator(std::shared_ptr<OperatorType> op,
153
                                                 int i,
154
                                                 double* factor,
155
156
                                                 double* estFactor)
  {
157
158
    // TODO: currently the est-factors are ignored
    assert( estFactor == NULL );
159

160
    vectorOperators[i].push_back(op);
161
    vectorFactors[i].push_back(factor);
162
    vectorChanging[i] = true;
163
  }
164

165
  template <class Traits>
166
  void ProblemStatSeq<Traits>::addVectorOperator(OperatorType& op,
167
                                                 int i,
168
                                                 double* factor,
169
170
                                                 double* estFactor)
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
171
    addVectorOperator(std::make_shared<OperatorType>(op), i, factor, estFactor);
172
  }
173

174
  template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
175
  void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, std::shared_ptr<OperatorType> > ops)
176
177
178
179
  {
    for (auto op : ops)
      addVectorOperator(op.second, op.first);
  }
180

181
  template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
182
  void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, std::pair<std::shared_ptr<OperatorType>, bool> > ops)
183
184
185
186
187
188
189
  {
    for (auto op : ops) {
      vectorOperators[op.first].push_back(op.second.first);
      vectorFactors[op.first].push_back(NULL);
      vectorChanging[op.first] = vectorChanging[op.first] || op.second.second;
    }
  }
190

191
  /// Adds a Dirichlet boundary condition
192
  template <class Traits>
193
    template <class Predicate, class Values>
194
195
  void ProblemStatSeq<Traits>::addDirichletBC(Predicate const& predicate,
                                              int row, int col,
196
                                              Values const& values)
197
  {
198
    static_assert( Concepts::Callable<Predicate, WorldVector>,
199
      "Function passed to addDirichletBC for predicate does not model the Callable<WorldVector> concept");
200
201

    static_assert( Concepts::Callable<Values, WorldVector>,
202
203
      "Function passed to addDirichletBC for values does not model the Callable<WorldVector> concept");

204
205
206
    test_exit(row >= 0 && row < nComponents, "row number out of range: ", row);
    test_exit(col >= 0 && col < nComponents, "col number out of range: ", col);

207
//     boundaryConditionSet = true;
208

209
    using BcType = DirichletBC<WorldVector>;
210
    dirichletBc[{row, col}].emplace_back( new BcType{predicate, values} );
211
  }
212
213


214
  template <class Traits>
215
216
  void ProblemStatSeq<Traits>::solve(AdaptInfo& adaptInfo,
                                     bool createMatrixData,
217
                                     bool storeMatrixData)
218
219
  {
    Timer t;
220

221
222
223
    SolverInfo solverInfo(name + "->solver");
    solverInfo.setCreateMatrixData(createMatrixData);
    solverInfo.setStoreMatrixData(storeMatrixData);
224

225
    solution->compress();
226
    linearSolver->solve(systemMatrix->getMatrix(),
227
228
                        solution->getVector(), rhs->getVector(),
                        solverInfo);
229

230
    if (solverInfo.getInfo() > 0) {
231
232
      msg("solution of discrete system needed ", t.elapsed(), " seconds");

233
234
      if (solverInfo.getAbsResidual() >= 0.0) {
        if (solverInfo.getRelResidual() >= 0.0)
235
236
          msg("Residual norm: ||b-Ax|| = ", solverInfo.getAbsResidual(),
                           ", ||b-Ax||/||b|| = ", solverInfo.getRelResidual());
237
        else
238
          msg("Residual norm: ||b-Ax|| = ", solverInfo.getAbsResidual());
239
240
      }
    }
241

242
243
244
    if (solverInfo.doBreak()) {
      std::stringstream tol_str;
      if (solverInfo.getAbsTolerance() > 0
245
          && solverInfo.getAbsResidual() > solverInfo.getAbsTolerance())
246
247
        tol_str << "absTol = " << solverInfo.getAbsTolerance() << " ";
      if (solverInfo.getRelTolerance() > 0
248
          && solverInfo.getRelResidual() > solverInfo.getRelTolerance())
249
        tol_str << "relTol = " << solverInfo.getRelTolerance() << " ";
250
      error_exit("Tolerance ", tol_str.str(), " could not be reached!");
251
    }
252
  }
253
254


255
256
  template <class Traits>
  void ProblemStatSeq<Traits>::buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag,
257
                                                 bool asmMatrix_, bool asmVector_)
258
259
  {
    Timer t;
260

261
262
    forEach(range_<0, nComponents>, [this](auto const _r) {
      this->getFeSpace(_r).update(this->leafGridView());
263
    });
264

265

266
267
    std::size_t nnz = 0;
    forEach(range_<0, nComponents>, [&,this](auto const _r)
268
    {
269
270
      static const int R = decltype(_r)::value;
      msg(this->getFeSpace(_r).size(), " DOFs for FeSpace[", R, "]");
271

272
      bool asmVector = asmVector_ && (!vectorAssembled[R] || vectorChanging[R]);
273

274
      if (asmVector) {
Praetorius, Simon's avatar
Praetorius, Simon committed
275
276
	      rhs->compress(_r);
        rhs->getVector(_r) = 0.0;
277
278
      }

279
      forEach(range_<0, nComponents>, [&,this](auto const _c)
280
      {
Praetorius, Simon's avatar
Praetorius, Simon committed
281
        static const int C = decltype(_c)::value;
282

Praetorius, Simon's avatar
Praetorius, Simon committed
283
284
        using MatrixData = typename ProblemStatSeq<Traits>::template MatrixData<R, C>;
        using VectorData = typename ProblemStatSeq<Traits>::template VectorData<R>;
285

286
        // The DOFMatrix which should be assembled
Praetorius, Simon's avatar
Praetorius, Simon committed
287
288
289
        auto& dofmatrix    = (*systemMatrix)(_r, _c);
        auto& solution_vec = (*solution)[_c];
        auto& rhs_vec      = (*rhs)[_r];
290

Praetorius, Simon's avatar
Praetorius, Simon committed
291
        auto row_col = std::make_pair(R, C);
292
        bool assembleMatrix = asmMatrix_ && (!matrixAssembled[row_col] || matrixChanging[row_col]);
Praetorius, Simon's avatar
Praetorius, Simon committed
293
        bool assembleVector = asmVector  && R == C;
294

295
296
297
298
299
        int r_ = 0, c_ = 0;
        if (assembleMatrix) {
          // init boundary condition
          for (auto bc_list : dirichletBc) {
            std::tie(r_, c_) = bc_list.first;
Praetorius, Simon's avatar
Praetorius, Simon committed
300
301
302
            if (r_ == R) {
              for (auto bc : bc_list.second)
                bc->init(c_ == C, dofmatrix, solution_vec, rhs_vec);
303
304
            }
          }
305
        }
306

307
        auto mat = MatrixData{dofmatrix, matrixOperators[row_col], matrixFactors[row_col], assembleMatrix};
Praetorius, Simon's avatar
Praetorius, Simon committed
308
        auto vec = VectorData{rhs_vec,   vectorOperators[R], vectorFactors[R], assembleVector};
309

310
        // assemble the DOFMatrix block and the corresponding rhs vector, of r==c
311

312
        dofmatrix.init(assembleMatrix);
313
314
        this->assemble(mat, vec, this->leafGridView());
        // TODO: assemble boundary conditions
315
        dofmatrix.finish();
316

317
318
319
        if (assembleMatrix)
          matrixAssembled[row_col] = true;
        if (assembleVector)
Praetorius, Simon's avatar
Praetorius, Simon committed
320
          vectorAssembled[R] = true;
321

322
        if (assembleMatrix) {
323
324
325
          // finish boundary condition
          for (auto bc_list : dirichletBc) {
            std::tie(r_, c_) = bc_list.first;
Praetorius, Simon's avatar
Praetorius, Simon committed
326
327
328
            if (r_ == R) {
              for (auto bc : bc_list.second)
                bc->finish(c_ == C, dofmatrix, solution_vec, rhs_vec);
329
330
331
332
333
            }
          }

          nnz += dofmatrix.getNnz();
        }
334
335
      });
    });
336
337
338
339

    msg("fillin of assembled matrix: ", nnz);

    msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
340
  }
341

342

343
344
  template <class Traits>
  void ProblemStatSeq<Traits>::writeFiles(AdaptInfo& adaptInfo, bool force)
345
  {
346
    Timer t;
347
    filewriter->write(adaptInfo.getTime(), *solution);
348
    msg("writeFiles needed ", t.elapsed(), " seconds");
349
350
351
  }


352
  template <class Traits>
353
    template <class LhsData, class RhsData, class GV>
354
  void ProblemStatSeq<Traits>::assemble(LhsData lhs, RhsData rhs,
355
                                        GV const& gridView)
356
  {
357
358
    auto const& rowFeSpace = lhs.matrix.getRowFeSpace();
    auto const& colFeSpace = lhs.matrix.getColFeSpace();
359

360
    if ((lhs.operators.empty() || !lhs.assemble) &&
361
362
        (rhs.operators.empty() || !rhs.assemble))
      return; // nothing to do
363
364


365
    for (auto op : lhs.operators)
366
      op->init(rowFeSpace, colFeSpace);
367
    for (auto op : rhs.operators)
368
      op->init(rowFeSpace, colFeSpace);
369

370
    auto rowLocalView = rowFeSpace.localView();
371
    auto rowIndexSet  = rowFeSpace.localIndexSet();
372

373
    auto colLocalView = colFeSpace.localView();
374
    auto colIndexSet  = colFeSpace.localIndexSet();
375

376
    for (auto const& element : elements(gridView)) {
377
378
      rowLocalView.bind(element);
      rowIndexSet.bind(rowLocalView);
379

380
381
      colLocalView.bind(element);
      colIndexSet.bind(colLocalView);
382

383
      if (lhs.assemble) {
384
        ElementMatrix elementMatrix;
385
        bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix,
386
                                    lhs.operators, lhs.factors);
387
        if (add)
388
          addElementMatrix(lhs.matrix, rowIndexSet, colIndexSet, elementMatrix);
389
      }
390

391
      if (rhs.assemble) {
392
        ElementVector elementVector;
393
        bool add = getElementVector(rowLocalView, elementVector,
394
                                    rhs.operators, rhs.factors);
395
        if (add)
396
          addElementVector(rhs.vector, rowIndexSet, elementVector);
397
      }
398
    }
399
  }
400
401


402
  template <class Traits>
403
    template <class RowView, class ColView>
404
405
406
  bool ProblemStatSeq<Traits>::getElementMatrix(RowView const& rowLocalView,
                                                ColView const& colLocalView,
                                                ElementMatrix& elementMatrix,
Praetorius, Simon's avatar
Praetorius, Simon committed
407
                                                std::list<std::shared_ptr<OperatorType>>& operators,
408
                                                std::list<double*> const& factors)
409
  {
410
411
    auto const nRows = rowLocalView.tree().finiteElement().size();
    auto const nCols = colLocalView.tree().finiteElement().size();
412

413
414
    test_exit_dbg(operators.size() == factors.size(),
                  "Nr. of operators and factors must be consistent!");
415
416

    // fills the entire matrix with zeroes
417
418
    elementMatrix.change_dim(nRows, nCols);
    set_to_zero(elementMatrix);
419

420
    bool add = false;
421

422
423
424
    auto op_it = operators.begin();
    auto fac_it = factors.begin();
    for (; op_it != operators.end(); ++op_it, ++fac_it) {
425
      bool add_op = (*op_it)->getElementMatrix(rowLocalView, colLocalView, elementMatrix, *fac_it);
426
427
      add = add || add_op;
    }
428

429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
#if 0
    auto const& element = rowLocalView.element();
    if (element.hasBoundaryIntersections()) {
      auto const& gridView = rowLocalView.globalBasis().gridView();
      for (auto const& intersection = intersections(gridView, element)) {
        if (!intersection.boundary())
          continue;

        auto op_it = boundary_operators.begin();
        auto fac_it = boundary_factors.begin();
        for (; op_it != boundary_operators.end(); ++op_it, ++fac_it) {
          bool add_op = (*op_it)->getBoundaryElementMatrix(intersection, rowLocalView, colLocalView, elementMatrix, *fac_it);
          add = add || add_op;
        }
      }
    }
#endif

447
    return add;
448
  }
449
450


451
  template <class Traits>
452
    template <class RowView>
453
454
  bool ProblemStatSeq<Traits>::getElementVector(RowView const& rowLocalView,
                                                ElementVector& elementVector,
Praetorius, Simon's avatar
Praetorius, Simon committed
455
                                                std::list<std::shared_ptr<OperatorType>>& operators,
456
                                                std::list<double*> const& factors)
457
  {
458
    auto const nRows = rowLocalView.tree().finiteElement().size();
459

460
461
    test_exit_dbg(operators.size() == factors.size(),
                  "Nr. of operators and factors must be consistent!");
462

463
    // Set all entries to zero
464
465
    elementVector.change_dim(nRows);
    set_to_zero(elementVector);
466

467
    bool add = false;
468

469
470
471
472
473
474
    auto op_it = operators.begin();
    auto fac_it = factors.begin();
    for (; op_it != operators.end(); ++op_it, ++fac_it) {
      bool add_op = (*op_it)->getElementVector(rowLocalView, elementVector, *fac_it);
      add = add || add_op;
    }
475

476
    return add;
477
  }
478
479


480
  template <class Traits>
481
    template <class Matrix, class RowIndexSet, class ColIndexSet>
482
483
484
485
  void ProblemStatSeq<Traits>::addElementMatrix(Matrix& dofmatrix,
                                                RowIndexSet const& rowIndexSet,
                                                ColIndexSet const& colIndexSet,
                                                ElementMatrix const& elementMatrix)
486
  {
487
    for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
488
      // The global index of the i−th vertex of the element
489
      auto const row = rowIndexSet.index(i);
490
      for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
491
        // The global index of the j−th vertex of the element
492
        auto const col = colIndexSet.index(j);
493
494
        dofmatrix(row, col) += elementMatrix[i][j];
      }
495
496
    }
  }
497
498


499
  template <class Traits>
500
    template <class Vector, class IndexSet>
501
502
503
  void ProblemStatSeq<Traits>::addElementVector(Vector& dofvector,
                                                IndexSet const& indexSet,
                                                ElementVector const& elementVector)
504
  {
505
    for (std::size_t i = 0; i < size(elementVector); ++i) {
506
      // The global index of the i-th vertex
507
      auto const row = indexSet.index(i);
508
      dofvector[row] += elementVector[i];
509
510
    }
  }
511

512
} // end namespace AMDiS