ProblemStat.inc.hpp 16.7 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>
110
    template <class OperatorType>
111
112
  void ProblemStatSeq<Traits>::
  addMatrixOperator(OperatorType const& op, int i, int j, double* factor, double* estFactor)
113
  {
114
115
    matrixOperators[i][j].emplace_back(std::make_shared<OperatorType>(op), factor, estFactor);
    matrixChanging[i][j] = true;
116
117
  }

118
119

  template <class Traits>
120
  void ProblemStatSeq<Traits>::
121
  addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j, double* factor, double* estFactor)
122
  {
123
    matrixBoundaryOperators[i][j].emplace_back(std::make_shared<IntersectionOperator>(op), factor, estFactor, b);
124
    matrixChanging[i][j] = true;
125
126
127
  }

  template <class Traits>
128
  void ProblemStatSeq<Traits>::
129
  addMatrixOperator(std::map< std::pair<int,int>, std::shared_ptr<ElementOperator> > ops)
130
  {
131
132
133
134
135
136
    for (auto op : ops){
      int r = op.first.first;
      int c = op.first.second;
      matrixOperators[r][c].emplace_back(op.second);
      matrixChanging[r][c] = true;
    }
137
  }
138
139

  template <class Traits>
140
  void ProblemStatSeq<Traits>::
141
  addMatrixOperator(std::map< std::pair<int,int>, std::pair<std::shared_ptr<ElementOperator>,bool> > ops)
142
143
  {
    for (auto op : ops) {
144
145
146
147
      int r = op.first.first;
      int c = op.first.second;
      matrixOperators[r][c].emplace_back(op.second.first);
      matrixChanging[r][c] = matrixChanging[r][c] || op.second.second;
148
149
    }
  }
150

151
  template <class Traits>
152
    template <class OperatorType>
153
154
  void ProblemStatSeq<Traits>::
  addVectorOperator(OperatorType const& op, int i, double* factor, double* estFactor)
155
  {
156
    vectorOperators[i].emplace_back(std::make_shared<OperatorType>(op), factor, estFactor);
157
    vectorChanging[i] = true;
158
  }
159

160
  template <class Traits>
161
  void ProblemStatSeq<Traits>::
162
  addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i, double* factor, double* estFactor)
163
  {
164
    vectorBoundaryOperators[i].emplace_back(std::make_shared<IntersectionOperator>(op), factor, estFactor, b);
165
    vectorChanging[i] = true;
166
  }
167

168

169
  template <class Traits>
170
  void ProblemStatSeq<Traits>::
171
  addVectorOperator(std::map< int, std::shared_ptr<ElementOperator> > ops)
172
  {
173
174
175
176
    for (auto op : ops) {
      vectorOperators[op.first].emplace_back(op.second);
      vectorChanging[op.first] = true;
    }
177
  }
178

179
  template <class Traits>
180
  void ProblemStatSeq<Traits>::
181
  addVectorOperator(std::map< int, std::pair<std::shared_ptr<ElementOperator>, bool> > ops)
182
183
  {
    for (auto op : ops) {
184
      vectorOperators[op.first].emplace_back(op.second.first);
185
186
187
      vectorChanging[op.first] = vectorChanging[op.first] || op.second.second;
    }
  }
188

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

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

201
202
203
    test_exit(row >= 0 && row < nComponents, "row number out of range: ", row);
    test_exit(col >= 0 && col < nComponents, "col number out of range: ", col);

204
//     boundaryConditionSet = true;
205

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


211
  template <class Traits>
212
213
  void ProblemStatSeq<Traits>::
  solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
214
215
  {
    Timer t;
216

217
218
219
    SolverInfo solverInfo(name + "->solver");
    solverInfo.setCreateMatrixData(createMatrixData);
    solverInfo.setStoreMatrixData(storeMatrixData);
220

221
    solution->compress();
222
    linearSolver->solve(systemMatrix->getMatrix(),
223
224
                        solution->getVector(), rhs->getVector(),
                        solverInfo);
225

226
    if (solverInfo.getInfo() > 0) {
227
228
      msg("solution of discrete system needed ", t.elapsed(), " seconds");

229
230
      if (solverInfo.getAbsResidual() >= 0.0) {
        if (solverInfo.getRelResidual() >= 0.0)
231
232
          msg("Residual norm: ||b-Ax|| = ", solverInfo.getAbsResidual(),
                           ", ||b-Ax||/||b|| = ", solverInfo.getRelResidual());
233
        else
234
          msg("Residual norm: ||b-Ax|| = ", solverInfo.getAbsResidual());
235
236
      }
    }
237

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


251
  template <class Traits>
252
253
  void ProblemStatSeq<Traits>::
  buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag, bool asmMatrix_, bool asmVector_)
254
255
  {
    Timer t;
256

257
258
    forEach(range_<0, nComponents>, [this](auto const _r) {
      this->getFeSpace(_r).update(this->leafGridView());
259
    });
260

261

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

268
      bool asmVector = asmVector_ && (!vectorAssembled[R] || vectorChanging[R]);
269

270
      if (asmVector) {
Praetorius, Simon's avatar
Praetorius, Simon committed
271
272
	      rhs->compress(_r);
        rhs->getVector(_r) = 0.0;
273
274
      }

275
      forEach(range_<0, nComponents>, [&,this](auto const _c)
276
      {
Praetorius, Simon's avatar
Praetorius, Simon committed
277
        static const int C = decltype(_c)::value;
278

Praetorius, Simon's avatar
Praetorius, Simon committed
279
280
        using MatrixData = typename ProblemStatSeq<Traits>::template MatrixData<R, C>;
        using VectorData = typename ProblemStatSeq<Traits>::template VectorData<R>;
281

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

287
        bool assembleMatrix = asmMatrix_ && (!matrixAssembled[R][C] || matrixChanging[R][C]);
Praetorius, Simon's avatar
Praetorius, Simon committed
288
        bool assembleVector = asmVector  && R == C;
289

290
291
        if (assembleMatrix) {
          // init boundary condition
292
293
294
          for (int c = 0; c < nComponents; ++c)
            for (auto bc : dirichletBc[R][c])
              bc->init(c == C, dofmatrix, solution_vec, rhs_vec);
295
        }
296

297
298
        auto mat = MatrixData{dofmatrix, matrixOperators[R][C], matrixBoundaryOperators[R][C], matrixIntersectionOperators[R][C], assembleMatrix};
        auto vec = VectorData{rhs_vec, vectorOperators[R], vectorBoundaryOperators[R], vectorIntersectionOperators[R], assembleVector};
299

300
        // assemble the DOFMatrix block and the corresponding rhs vector, of r==c
301

302
        dofmatrix.init(assembleMatrix);
303
        this->assemble(mat, vec, this->leafGridView());
304
        dofmatrix.finish();
305

306
        if (assembleMatrix)
307
          matrixAssembled[R][C] = true;
308
        if (assembleVector)
Praetorius, Simon's avatar
Praetorius, Simon committed
309
          vectorAssembled[R] = true;
310

311
        if (assembleMatrix) {
312
          // finish boundary condition
313
314
315
          for (int c = 0; c < nComponents; ++c)
            for (auto bc : dirichletBc[R][c])
              bc->finish(c == C, dofmatrix, solution_vec, rhs_vec);
316
317
318

          nnz += dofmatrix.getNnz();
        }
319
320
      });
    });
321
322
323
324

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

    msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
325
  }
326

327

328
329
  template <class Traits>
  void ProblemStatSeq<Traits>::writeFiles(AdaptInfo& adaptInfo, bool force)
330
  {
331
    Timer t;
332
    filewriter->write(adaptInfo.getTime(), *solution);
333
    msg("writeFiles needed ", t.elapsed(), " seconds");
334
335
336
  }


337
  template <class Traits>
338
    template <class LhsData, class RhsData, class GV>
339
  void ProblemStatSeq<Traits>::assemble(LhsData lhs, RhsData rhs, GV const& gridView)
340
  {
341
342
    auto const& rowFeSpace = lhs.matrix.getRowFeSpace();
    auto const& colFeSpace = lhs.matrix.getColFeSpace();
343

344
    if ((lhs.operators.empty() || !lhs.assemble) &&
345
346
        (rhs.operators.empty() || !rhs.assemble))
      return; // nothing to do
347

348
349
350
351
    for (auto scaledOp : lhs.operators)
      scaledOp.op->init(rowFeSpace, colFeSpace);
    for (auto scaledOp : rhs.operators)
      scaledOp.op->init(rowFeSpace, colFeSpace);
352

353
    auto rowLocalView = rowFeSpace.localView();
354
    auto rowIndexSet  = rowFeSpace.localIndexSet();
355

356
    auto colLocalView = colFeSpace.localView();
357
    auto colIndexSet  = colFeSpace.localIndexSet();
358

359
    for (auto const& element : elements(gridView)) {
360
361
      rowLocalView.bind(element);
      rowIndexSet.bind(rowLocalView);
362

363
364
      colLocalView.bind(element);
      colIndexSet.bind(colLocalView);
365

366
      if (lhs.assemble) {
367
        ElementMatrix elementMatrix;
368
        bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix,
369
                                    lhs.operators, lhs.boundary_operators, lhs.intersection_operators);
370
        if (add)
371
          addElementMatrix(lhs.matrix, rowIndexSet, colIndexSet, elementMatrix);
372
      }
373

374
      if (rhs.assemble) {
375
        ElementVector elementVector;
376
        bool add = getElementVector(rowLocalView, elementVector,
377
                                    rhs.operators, rhs.boundary_operators, rhs.intersection_operators);
378
        if (add)
379
          addElementVector(rhs.vector, rowIndexSet, elementVector);
380
      }
381
    }
382
  }
383
384


385
  template <class Traits>
386
    template <class RowView, class ColView>
387
388
389
  bool ProblemStatSeq<Traits>::getElementMatrix(RowView const& rowLocalView,
                                                ColView const& colLocalView,
                                                ElementMatrix& elementMatrix,
390
391
392
                                                std::list<Scaled<ElementOperator>>& operators,
                                                std::list<Scaled<IntersectionOperator>>& boundary_operators,
                                                std::list<Scaled<IntersectionOperator>>& intersection_operators)
393
  {
394
395
    auto const nRows = rowLocalView.tree().finiteElement().size();
    auto const nCols = colLocalView.tree().finiteElement().size();
396

397
398
399
    auto const& element = rowLocalView.element();
    auto const& gridView = rowLocalView.globalBasis().gridView();

400
    // fills the entire matrix with zeroes
401
402
    elementMatrix.change_dim(nRows, nCols);
    set_to_zero(elementMatrix);
403

404
    bool add = false;
405

406
407
408
409
410
411
412
    auto assemble_operators = [&](auto& e, auto& operator_list) {
      for (auto scaled : operator_list) {
        bool add_op = scaled.op->getElementMatrix(e, rowLocalView, colLocalView, elementMatrix, scaled.factor);
        add = add || add_op;
      }
    };

413
    // assemble element operators
414
    assemble_operators(element, operators);
415

416
417
418
419
420
421
422
    // assemble intersection operators
    if (!boundary_operators.empty() || !intersection_operators.empty()) {
      for (auto const& intersection : intersections(gridView, element)) {
        if (intersection.boundary())
          assemble_operators(intersection, boundary_operators);
        else
          assemble_operators(intersection, intersection_operators); // TODO: assemble intersection operators also on the boundary?
423
424
425
      }
    }

426
    return add;
427
  }
428
429


430
  template <class Traits>
431
    template <class RowView>
432
433
  bool ProblemStatSeq<Traits>::getElementVector(RowView const& rowLocalView,
                                                ElementVector& elementVector,
434
435
436
                                                std::list<Scaled<ElementOperator>>& operators,
                                                std::list<Scaled<IntersectionOperator>>& boundary_operators,
                                                std::list<Scaled<IntersectionOperator>>& intersection_operators)
437
  {
438
    auto const nRows = rowLocalView.tree().finiteElement().size();
439

440
441
442
    auto const& element = rowLocalView.element();
    auto const& gridView = rowLocalView.globalBasis().gridView();

443
    // Set all entries to zero
444
445
    elementVector.change_dim(nRows);
    set_to_zero(elementVector);
446

447
    bool add = false;
448

449
450
451
452
453
454
455
    auto assemble_operators = [&](auto& e, auto& operator_list) {
      for (auto scaled : operator_list) {
        bool add_op = scaled.op->getElementVector(e, rowLocalView, elementVector, scaled.factor);
        add = add || add_op;
      }
    };

456
    // assemble element operators
457
    assemble_operators(element, operators);
458

459
460
461
462
463
464
465
    // assemble intersection operators
    if (!boundary_operators.empty() || !intersection_operators.empty()) {
      for (auto const& intersection : intersections(gridView, element)) {
        if (intersection.boundary())
          assemble_operators(intersection, boundary_operators);
        else
          assemble_operators(intersection, intersection_operators); // TODO: assemble intersection operators also on the boundary?
466
467
468
      }
    }

469
    return add;
470
  }
471
472


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


492
  template <class Traits>
493
    template <class Vector, class IndexSet>
494
495
496
  void ProblemStatSeq<Traits>::addElementVector(Vector& dofvector,
                                                IndexSet const& indexSet,
                                                ElementVector const& elementVector)
497
  {
498
    for (std::size_t i = 0; i < size(elementVector); ++i) {
499
      // The global index of the i-th vertex
500
      auto const row = indexSet.index(i);
501
      dofvector[row] += elementVector[i];
502
503
    }
  }
504

505
} // end namespace AMDiS