ProblemStat.inc.hpp 15.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
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
248
249
250
251
template <class Traits>
void ProblemStat<Traits>::
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag, bool asmMatrix_, bool asmVector_)
{
  Timer t;
252

253
254
255
256
  // update global feSpace, i.e. necessary after mesh change
  forEach(range_<0, nComponents>, [this](auto const _r) {
    this->getFeSpace(_r).update(this->leafGridView());
  });
257

258

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

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

    if (asmVector) {
      rhs->compress(_r);
      rhs->getVector(_r) = 0.0;
    }
271

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

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

279
280
281
282
      // The DOFMatrix which should be assembled
      auto& dofmatrix    = (*systemMatrix)(_r, _c);
      auto& solution_vec = (*solution)[_c];
      auto& rhs_vec      = (*rhs)[_r];
283

284
285
      bool assembleMatrix = asmMatrix_ && (!matrixAssembled[R][C] || matrixChanging[R][C]);
      bool assembleVector = asmVector  && R == C;
286

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

294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
      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};

      // assemble the DOFMatrix block and the corresponding rhs vector if r==c

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

      if (assembleMatrix)
        matrixAssembled[R][C] = true;
      if (assembleVector)
        vectorAssembled[R] = true;

      if (assembleMatrix) {
        // 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])
              bc->finish(r == R, c == C, dofmatrix, solution_vec, rhs_vec);
          }
317
        }
318

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

324
325
326
  msg("fillin of assembled matrix: ", nnz);
  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
327

328

329
330
331
332
333
334
335
336
template <class Traits>
void ProblemStat<Traits>::
writeFiles(AdaptInfo& adaptInfo, bool force)
{
  Timer t;
  filewriter->write(adaptInfo.getTime(), *solution);
  msg("writeFiles needed ", t.elapsed(), " seconds");
}
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
362
363
364
365
366
367
368
369
370
371
372
373
374
375
template <class Traits>
  template <class LhsData, class RhsData, class GV>
void ProblemStat<Traits>::
assemble(LhsData lhs, RhsData rhs, GV const& gridView)
{
  auto const& rowFeSpace = lhs.matrix.getRowFeSpace();
  auto const& colFeSpace = lhs.matrix.getColFeSpace();

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

  for (auto scaledOp : lhs.operators)
    scaledOp.op->init(rowFeSpace, colFeSpace);
  for (auto scaledOp : rhs.operators)
    scaledOp.op->init(rowFeSpace, colFeSpace);

  auto rowLocalView = rowFeSpace.localView();
  auto rowIndexSet  = rowFeSpace.localIndexSet();

  auto colLocalView = colFeSpace.localView();
  auto colIndexSet  = colFeSpace.localIndexSet();

  for (auto const& element : elements(gridView)) {
    // TODO: use only one localView and localIndexSet if feSpaces are equal
    rowLocalView.bind(element);
    colLocalView.bind(element);
    rowIndexSet.bind(rowLocalView); // NOTE: expensive operation!
    colIndexSet.bind(colLocalView);

    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);
    }
376

377
378
379
380
381
382
383
    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);
    }
384

385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
    rowIndexSet.unbind();
    rowLocalView.unbind();
    colLocalView.unbind();
    colIndexSet.unbind();
  }
}


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();
405

406
407
  auto const& element = rowLocalView.element();
  auto const& gridView = rowLocalView.globalBasis().gridView();
408

409
410
411
  // fills the entire matrix with zeroes
  elementMatrix.change_dim(nRows, nCols);
  set_to_zero(elementMatrix);
412

413
  bool add = false;
414

415
416
417
418
419
420
  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;
    }
  };
421

422
423
  // assemble element operators
  assemble_operators(element, operators);
424

425
426
427
428
429
430
431
432
433
  // 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?
434
    }
435
  }
436

437
438
  return add;
}
439

440

441
442
443
444
445
446
447
448
449
450
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();
451

452
453
  auto const& element = rowLocalView.element();
  auto const& gridView = rowLocalView.globalBasis().gridView();
454

455
456
457
  // Set all entries to zero
  elementVector.change_dim(nRows);
  set_to_zero(elementVector);
458

459
  bool add = false;
460

461
462
463
464
465
466
  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;
    }
  };
467

468
469
  // assemble element operators
  assemble_operators(element, operators);
470

471
472
473
  // assemble intersection operators
  if (!intersection_operators.empty()
      || (!boundary_operators.empty() && element.hasBoundaryIntersections()))
474
  {
475
476
477
478
479
    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?
480
    }
481
  }
482

483
484
  return add;
}
485

486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501

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];
502
503
    }
  }
504
}
505
506


507
508
509
510
511
512
513
514
515
516
517
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];
518
  }
519
}
520

521
} // end namespace AMDiS