ProblemStat.inc.hpp 17.9 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
9
10
11
12
namespace AMDiS {

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

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

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

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

38
39
40
41
42
  int globalRefinements = 0;
  Parameters::get(meshName + "->global refinements", globalRefinements);
  if (globalRefinements > 0) {
    mesh->globalRefine(globalRefinements);
  }
43
44


45
46
47
48
49
50
51
52
  // create fespace
  if (feSpaces) {
    warning("feSpaces already created");
  }
  else {
    if (initFlag.isSet(INIT_FE_SPACE) ||
        (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
      createFeSpaces();
53
    }
54

55
56
57
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
      feSpaces = adoptProblem->getFeSpaces();
58
    }
59
  }
60

61
62
  if (!feSpaces)
    warning("no feSpaces created\n");
63
64


65
66
67
  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();
68

69
70
71
72
73
  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
    solution = adoptProblem->getSolution();
    rhs = adoptProblem->getRhs();
    systemMatrix = adoptProblem->getSystemMatrix();
  }
74

75

76
77
78
79
80
81
82
  // create solver
  if (linearSolver) {
    warning("solver already created\n");
  }
  else {
    if (initFlag.isSet(INIT_SOLVER))
      createSolver();
83

84
85
86
    if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
      test_exit(!linearSolver, "solver already created\n");
      linearSolver = adoptProblem->getSolver();
87
    }
88
89
  }

90
91
  if (!linearSolver) {
    warning("no solver created\n");
92
93
  }

94

95
96
97
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
98

99
100
  solution->compress();
}
101

102

103
// add matrix/vector operator terms
104

105
106
107
108
109
110
111
112
template <class Traits>
  template <class OperatorType>
void ProblemStat<Traits>::
addMatrixOperator(OperatorType const& op, int i, int j, double* factor, double* estFactor)
{
  matrixOperators[i][j].push_back({std::make_shared<OperatorType>(op), factor, estFactor});
  matrixChanging[i][j] = true;
}
113

114
115
116
117
118
119
120
template <class Traits>
void ProblemStat<Traits>::
addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j, double* factor, double* estFactor)
{
  matrixBoundaryOperators[i][j].push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor, b});
  matrixChanging[i][j] = true;
}
121

122
123
template <class Traits>
void ProblemStat<Traits>::
124
addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops)
125
126
127
128
{
  for (auto op : ops){
    int r = op.first.first;
    int c = op.first.second;
129
    matrixOperators[r][c].push_back({std::make_shared<ElementOperator>(op.second)});
130
    matrixChanging[r][c] = true;
131
  }
132
}
133

134
135
template <class Traits>
void ProblemStat<Traits>::
136
addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> > ops)
137
138
139
140
{
  for (auto op : ops) {
    int r = op.first.first;
    int c = op.first.second;
141
    matrixOperators[r][c].push_back({std::make_shared<ElementOperator>(op.second.first)});
142
    matrixChanging[r][c] = matrixChanging[r][c] || op.second.second;
143
  }
144
}
145
146


147
148
149
150
151
152
153
154
template <class Traits>
  template <class OperatorType>
void ProblemStat<Traits>::
addVectorOperator(OperatorType const& op, int i, double* factor, double* estFactor)
{
  vectorOperators[i].push_back({std::make_shared<OperatorType>(op), factor, estFactor});
  vectorChanging[i] = true;
}
155

156
157
158
159
160
161
162
template <class Traits>
void ProblemStat<Traits>::
addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i, double* factor, double* estFactor)
{
  vectorBoundaryOperators[i].push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor, b});
  vectorChanging[i] = true;
}
163

164
165
template <class Traits>
void ProblemStat<Traits>::
166
addVectorOperator(std::map< int, ElementOperator> ops)
167
168
{
  for (auto op : ops) {
169
    vectorOperators[op.first].push_back({std::make_shared<ElementOperator>(op.second)});
170
171
172
    vectorChanging[op.first] = true;
  }
}
173

174
175
template <class Traits>
void ProblemStat<Traits>::
176
addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops)
177
178
{
  for (auto op : ops) {
179
    vectorOperators[op.first].push_back({std::make_shared<ElementOperator>(op.second.first)});
180
    vectorChanging[op.first] = vectorChanging[op.first] || op.second.second;
181
  }
182
}
183
184


185
186
187
188
189
190
191
192
// Adds a Dirichlet boundary condition
template <class Traits>
  template <class Predicate, class Values>
void ProblemStat<Traits>::
addDirichletBC(Predicate const& predicate, int row, int col, Values const& values)
{
  static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
    "Function passed to addDirichletBC for predicate does not model the Functor<bool(WorldVector)> concept");
193

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

197
198
  test_exit(row >= 0 && row < nComponents, "row number out of range: ", row);
  test_exit(col >= 0 && col < nComponents, "col number out of range: ", col);
199

200
//     boundaryConditionSet = true;
201

202
203
204
  using BcType = DirichletBC<WorldVector>;
  dirichletBc[row][col].emplace_back( new BcType{predicate, values} );
}
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
template <class Traits>
void ProblemStat<Traits>::
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
  Timer t;

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

  solution->compress();
  linearSolver->solve(systemMatrix->getMatrix(),
                      solution->getVector(), rhs->getVector(),
                      solverInfo);

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

    if (solverInfo.getAbsResidual() >= 0.0) {
      if (solverInfo.getRelResidual() >= 0.0)
        msg("Residual norm: ||b-Ax|| = ", solverInfo.getAbsResidual(),
                          ", ||b-Ax||/||b|| = ", solverInfo.getRelResidual());
      else
        msg("Residual norm: ||b-Ax|| = ", solverInfo.getAbsResidual());
    }
  }
233

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

246

247
template <class Traits>
248
249
250
  template <int R, int C>
bool ProblemStat<Traits>::
assembleMatrix(bool asmMatrix_, const index_t<R> _r = {}, const index_t<C> _c = {}) const
251
{
252
253
  return asmMatrix_ && (!matrixAssembled[R][C] || matrixChanging[R][C]);
}
254

255
256
257
258
259
260
261
template <class Traits>
  template <int R>
bool ProblemStat<Traits>::
assembleVector(bool asmVector_, const index_t<R> _r = {}) const
{
  return asmVector_ && (!vectorAssembled[R] || vectorChanging[R]);
}
262

263
264
265
266
template <class Traits>
void ProblemStat<Traits>::
initMatrixVector(bool asmMatrix_, bool asmVector_)
{
267
268
269
270
  forEach(range_<0, nComponents>, [&,this](auto const _r)
  {
    static const int R = decltype(_r)::value;
    msg(this->getFeSpace(_r).size(), " DOFs for FeSpace[", R, "]");
271

272
    if (assembleVector(asmVector_, _r)) {
273
274
      rhs->compress(_r);
      rhs->getVector(_r) = 0.0;
275
276
277
278
279
280
281
282

      // init vector operators
      for (auto& op : vectorOperators[R])
        op.init(this->getFeSpace(_r));
      for (auto& op : vectorBoundaryOperators[R])
        op.init(this->getFeSpace(_r));
      for (auto& op : vectorIntersectionOperators[R])
        op.init(this->getFeSpace(_r));
283
    }
284

285
286
287
    forEach(range_<0, nComponents>, [&,this](auto const _c)
    {
      static const int C = decltype(_c)::value;
288

289
290
      bool asmMatrix = assembleMatrix(asmMatrix_, _r, _c);
      (*systemMatrix)(_r, _c).init(asmMatrix);
291

292
293
294
295
296
297
298
299
      if (asmMatrix) {
        // init matrix operators
        for (auto& op : matrixOperators[R][C])
          op.init(this->getFeSpace(_r), this->getFeSpace(_c));
        for (auto& op : matrixBoundaryOperators[R][C])
          op.init(this->getFeSpace(_r), this->getFeSpace(_c));
        for (auto& op : matrixIntersectionOperators[R][C])
          op.init(this->getFeSpace(_r), this->getFeSpace(_c));
300

301
302
303
        // init boundary condition
        for (int c = 0; c < nComponents; ++c)
          for (auto bc : dirichletBc[R][c])
304
            bc->init(c == C, (*systemMatrix)(_r, _c), (*solution)[_c], (*rhs)[_r]);
305
      }
306
307
308
    });
  });
}
309

310

311
312
313
314
315
316
317
318
319
320
321
322
323
template <class Traits>
std::size_t finishMatrixVector(bool asmMatrix_, bool asmVector_)
{
  std::size_t nnz = 0;
  forEach(range_<0, nComponents>, [&,this](auto const _r)
  {
    static const int R = decltype(_r)::value;
    if (assembleVector(asmVector_, _r))
      vectorAssembled[R] = true;
    forEach(range_<0, nComponents>, [&,this](auto const _c)
    {
      static const int C = decltype(_c)::value;
      bool asmMatrix = assembleMatrix(asmMatrix_, _r, _c);
324

325
      (*systemMatrix)(_r, _c).finish();
326

327
      if (asmMatrix)
328
329
        matrixAssembled[R][C] = true;

330
      if (asmMatrix) {
331
332
333
334
335
336
        // finish boundary condition
        for (int c = 0; c < nComponents; ++c) {
          for (int r = 0; r < nComponents; ++r) {
            if (r != R && c != C)
              continue;
            for (auto bc : dirichletBc[r][c])
337
              bc->finish(r == R, c == C, (*systemMatrix)(_r, _c), (*solution)[_c], (*rhs)[_r]);
338
          }
339
        }
340

341
342
343
        nnz += dofmatrix.getNnz();
    });
  });
344

345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
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
414
415
  return nnz;
}


template <class Traits>
  template <class Element>
void ProblemStat<Traits>::
LocalFiniteElement<typename Traits::FeSpaces> initElement(Element const& element)
{
  LocalFiniteElement<FeSpaces> localFiniteElem(*feSpaces);

  forEach(range_<0, nComponents>, [&,this](auto const _i)
  {
    auto& localView = localFiniteElem.localView(_i);
    auto& localIndexSet  = localFiniteElem.localIndexSet(_i);

    localView.bind(element);
    localIndexSet.bind(element);
  });

  return localFiniteElem;
}


template <class Traits>
void ProblemStat<Traits>::
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag, bool asmMatrix_, bool asmVector_)
{
  Timer t;

  auto gridView = this->leafGridView();

  // 1. update global feSpace. This is necessary after mesh change
  forEach(range_<0, nComponents>, [&,this](auto const _r) {
    this->getFeSpace(_r).update(gridView);
  });

  // 2. init matrix and rhs vector and initialize dirichlet boundary conditions
  initMatrixVector(asmMatrix_, asmVector_);

  // 3. assemble operators
  for (auto const& element : elements(gridView))
  {
    auto localFiniteElem = initElement(element);

    forEach(range_<0, nComponents>, [&,this](auto const _r)
    {
      forEach(range_<0, nComponents>, [&,this](auto const _c)
      {
        static const int R = decltype(_r)::value;
        static const int C = decltype(_c)::value;

        using MatrixData = typename ProblemStat<Traits>::template MatrixData<R, C>;
        using VectorData = typename ProblemStat<Traits>::template VectorData<R>;

        auto mat = MatrixData{(*systemMatrix)(_r, _c),
          matrixOperators[R][C], matrixBoundaryOperators[R][C], matrixIntersectionOperators[R][C],
          assembleMatrix(asmMatrix_, _r, _c)};
        auto vec = VectorData{(*rhs)[_r],
          vectorOperators[R], vectorBoundaryOperators[R], vectorIntersectionOperators[R],
          assembleVector(asmVector_, _r) && R==C};

        // assemble the DOFMatrix block and the corresponding rhs vector if r==c
        this->assemble(mat, vec, element);
      });
    });
  }

  // 4. finish matrix insertion and apply dirichlet boundary conditions
  std::size_t nnz = finishMatrixVector(asmMatrix_, asmVector_);

416
417
418
  msg("fillin of assembled matrix: ", nnz);
  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
419

420

421
422
423
424
425
426
427
428
template <class Traits>
void ProblemStat<Traits>::
writeFiles(AdaptInfo& adaptInfo, bool force)
{
  Timer t;
  filewriter->write(adaptInfo.getTime(), *solution);
  msg("writeFiles needed ", t.elapsed(), " seconds");
}
429
430


431
template <class Traits>
432
  template <int R, int C, class LocalFE>
433
void ProblemStat<Traits>::
434
assemble(LocalFe& localFiniteElem, MatrixData<R,C> lhs, VectorData<R> rhs)
435
{
436
437
438
439
440
441
  if (((lhs.operators.empty() &&
        lhs.boundary_operators.empty() &&
        lhs.intersection_operators.empty()) || !lhs.assemble) &&
      ((rhs.operators.empty() &&
        rhs.boundary_operators.empty() &&
        rhs.intersection_operators.empty()) || !rhs.assemble))
442
443
    return; // nothing to do

444
445
  const index_t<R> _r{};
  const index_t<C> _c{};
446

447
448
449
450
451
  auto& rowLocalView = localFiniteElem.localView(_r);
  auto& rowIndexSet  = localFiniteElem.localIndexSet(_r);

  auto& colLocalView = localFiniteElem.localView(_c);
  auto& colIndexSet  = localFiniteElem.localIndexSet(_c);
452

453
454
455
456
457
458
459
460
461
462
463
464
465
466
467

  if (lhs.assemble) {
    ElementMatrix elementMatrix;
    bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix,
                  lhs.operators, lhs.boundary_operators, lhs.intersection_operators);
    if (add)
      addElementMatrix(lhs.matrix, rowIndexSet, colIndexSet, elementMatrix);
  }

  if (rhs.assemble) {
    ElementVector elementVector;
    bool add = getElementVector(rowLocalView, elementVector,
                  rhs.operators, rhs.boundary_operators, rhs.intersection_operators);
    if (add)
      addElementVector(rhs.vector, rowIndexSet, elementVector);
468
  }
469

470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
}


template <class Traits>
  template <class RowView, class ColView>
bool ProblemStat<Traits>::
getElementMatrix(RowView const& rowLocalView,
                 ColView const& colLocalView,
                 ElementMatrix& elementMatrix,
                 std::list<Scaled<ElementOperator>>& operators,
                 std::list<Scaled<IntersectionOperator>>& boundary_operators,
                 std::list<Scaled<IntersectionOperator>>& intersection_operators)
{
  auto const nRows = rowLocalView.tree().finiteElement().size();
  auto const nCols = colLocalView.tree().finiteElement().size();
485

486
487
  auto const& element = rowLocalView.element();
  auto const& gridView = rowLocalView.globalBasis().gridView();
488

489
490
491
  // fills the entire matrix with zeroes
  elementMatrix.change_dim(nRows, nCols);
  set_to_zero(elementMatrix);
492

493
  bool add = false;
494

495
496
497
498
499
500
  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;
    }
  };
501

502
503
  // assemble element operators
  assemble_operators(element, operators);
504

505
506
507
508
509
510
511
512
513
  // assemble intersection operators
  if (!intersection_operators.empty()
      || (!boundary_operators.empty() && element.hasBoundaryIntersections()))
  {
    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?
514
    }
515
  }
516

517
518
  return add;
}
519

520

521
522
523
524
525
526
527
528
529
530
template <class Traits>
  template <class RowView>
bool ProblemStat<Traits>::
getElementVector(RowView const& rowLocalView,
                 ElementVector& elementVector,
                 std::list<Scaled<ElementOperator>>& operators,
                 std::list<Scaled<IntersectionOperator>>& boundary_operators,
                 std::list<Scaled<IntersectionOperator>>& intersection_operators)
{
  auto const nRows = rowLocalView.tree().finiteElement().size();
531

532
533
  auto const& element = rowLocalView.element();
  auto const& gridView = rowLocalView.globalBasis().gridView();
534

535
536
537
  // Set all entries to zero
  elementVector.change_dim(nRows);
  set_to_zero(elementVector);
538

539
  bool add = false;
540

541
542
543
544
545
546
  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;
    }
  };
547

548
549
  // assemble element operators
  assemble_operators(element, operators);
550

551
552
553
  // assemble intersection operators
  if (!intersection_operators.empty()
      || (!boundary_operators.empty() && element.hasBoundaryIntersections()))
554
  {
555
556
557
558
559
    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?
560
    }
561
  }
562

563
564
  return add;
}
565

566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581

template <class Traits>
  template <class Matrix, class RowIndexSet, class ColIndexSet>
void ProblemStat<Traits>::
addElementMatrix(Matrix& dofmatrix,
                 RowIndexSet const& rowIndexSet,
                 ColIndexSet const& colIndexSet,
                 ElementMatrix const& elementMatrix)
{
  for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
    // The global index of the i−th vertex of the element
    auto const row = rowIndexSet.index(i);
    for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
      // The global index of the j−th vertex of the element
      auto const col = colIndexSet.index(j);
      dofmatrix(row, col) += elementMatrix[i][j];
582
583
    }
  }
584
}
585
586


587
588
589
590
591
592
593
594
595
596
597
template <class Traits>
  template <class Vector, class IndexSet>
void ProblemStat<Traits>::
addElementVector(Vector& dofvector,
                 IndexSet const& indexSet,
                 ElementVector const& elementVector)
{
  for (std::size_t i = 0; i < size(elementVector); ++i) {
    // The global index of the i-th vertex
    auto const row = indexSet.index(i);
    dofvector[row] += elementVector[i];
598
  }
599
}
600

601
} // end namespace AMDiS