ProblemStat.inc.hpp 16.1 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
  void ProblemStatSeq<Traits>::addMatrixOperator(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)
  {
130
    addMatrixOperator(make_shared<OperatorType>(op), i, j, factor, estFactor);
131
132
133
  }

  template <class Traits>
134
  void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, shared_ptr<OperatorType> > ops)
135
136
137
138
  {
    for (auto op : ops)
      addMatrixOperator(op.second, op.first.first, op.first.second);
  }
139
140
141
142
143
144
145
146
147
148
149

  template <class Traits>
  void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, std::pair<shared_ptr<OperatorType>,bool> > ops)
  {
    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>
152
  void ProblemStatSeq<Traits>::addVectorOperator(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)
  {
171
    addVectorOperator(make_shared<OperatorType>(op), i, factor, estFactor);
172
  }
173

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

181
182
183
184
185
186
187
188
189
  template <class Traits>
  void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, std::pair<shared_ptr<OperatorType>, bool> > ops)
  {
    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
    For<0, nComponents>::loop([this](auto const _r) {
262
	    this->getNodeFactory(_r).initializeIndices();
263
    });
264

265

266
    size_t nnz = 0;
267
    For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector_](auto const r)
268
    {
269
      msg(this->getFeSpace(r).size(), " DOFs for FeSpace[", r, "]");
270

271
      bool asmVector = asmVector_ && (!vectorAssembled[int(r)] || vectorChanging[int(r)]);
272

273
      if (asmVector) {
274
275
	      rhs->compress(r);
        rhs->getVector(r) = 0.0;
276
277
      }

278
      For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector](auto const c)
279
      {
280
281
282
283
        static const int _c = decltype(c)::value;
        static const int _r = decltype(r)::value;
        index_<_r> r{};

284
285
        using MatrixData = typename ProblemStatSeq<Traits>::template MatrixData<_r, _c>;
        using VectorData = typename ProblemStatSeq<Traits>::template VectorData<_r>;
286

287
        // The DOFMatrix which should be assembled
288
289
290
        auto& dofmatrix    = (*systemMatrix)(r, c);
        auto& solution_vec = (*solution)[c];
        auto& rhs_vec      = (*rhs)[r];
291

292
        auto row_col = std::make_pair(_r, _c);
293
294
        bool assembleMatrix = asmMatrix_ && (!matrixAssembled[row_col] || matrixChanging[row_col]);
        bool assembleVector = asmVector  && _r == _c;
295

296
297
298
299
300
301
302
303
304
305
        int r_ = 0, c_ = 0;
        if (assembleMatrix) {
          // init boundary condition
          for (auto bc_list : dirichletBc) {
            std::tie(r_, c_) = bc_list.first;
            if (r_ == _r) {
                    for (auto bc : bc_list.second)
                      bc->init(c_ == _c, dofmatrix, solution_vec, rhs_vec);
            }
          }
306
        }
307

308
        auto mat = MatrixData{dofmatrix, matrixOperators[row_col], matrixFactors[row_col], assembleMatrix};
309
        auto vec = VectorData{rhs_vec,   vectorOperators[_r], vectorFactors[_r], assembleVector};
310

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

313
314
315
        dofmatrix.init(assembleMatrix);
        this->assemble(mat, vec, elements(*meshView));
        dofmatrix.finish();
316

317
318
319
        if (assembleMatrix)
          matrixAssembled[row_col] = true;
        if (assembleVector)
320
          vectorAssembled[_r] = true;
321

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

          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
354
355
    template <class LhsData, class RhsData, class Elements>
  void ProblemStatSeq<Traits>::assemble(LhsData lhs, RhsData rhs,
                                        Elements const& elements)
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) {
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,
407
                                                std::list<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
    return add;
430
  }
431
432


433
  template <class Traits>
434
    template <class RowView>
435
436
  bool ProblemStatSeq<Traits>::getElementVector(RowView const& rowLocalView,
                                                ElementVector& elementVector,
437
                                                std::list<shared_ptr<OperatorType>>& operators,
438
                                                std::list<double*> const& factors)
439
  {
440
    auto const nRows = rowLocalView.tree().finiteElement().size();
441

442
443
    test_exit_dbg(operators.size() == factors.size(),
                  "Nr. of operators and factors must be consistent!");
444

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

449
    bool add = false;
450

451
452
453
454
455
456
    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;
    }
457

458
    return add;
459
  }
460
461


462
  template <class Traits>
463
    template <class Matrix, class RowIndexSet, class ColIndexSet>
464
465
466
467
  void ProblemStatSeq<Traits>::addElementMatrix(Matrix& dofmatrix,
                                                RowIndexSet const& rowIndexSet,
                                                ColIndexSet const& colIndexSet,
                                                ElementMatrix const& elementMatrix)
468
  {
469
    for (size_t i = 0; i < num_rows(elementMatrix); ++i) {
470
      // The global index of the i−th vertex of the element
471
472
      auto const row = rowIndexSet.index(i);
      for (size_t j = 0; j < num_cols(elementMatrix); ++j) {
473
        // The global index of the j−th vertex of the element
474
        auto const col = colIndexSet.index(j);
475
476
        dofmatrix(row, col) += elementMatrix[i][j];
      }
477
478
    }
  }
479
480


481
  template <class Traits>
482
    template <class Vector, class IndexSet>
483
484
485
  void ProblemStatSeq<Traits>::addElementVector(Vector& dofvector,
                                                IndexSet const& indexSet,
                                                ElementVector const& elementVector)
486
  {
487
    for (size_t i = 0; i < size(elementVector); ++i) {
488
      // The global index of the i-th vertex
489
      auto const row = indexSet.index(i);
490
      dofvector[row] += elementVector[i];
491
492
    }
  }
493

494
} // end namespace AMDiS