ProblemStat.inc.hpp 14.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
12
13
14
15
16
17
18
19
20
  template <class Traits>
  void ProblemStatSeq<Traits>::initialize(Flag initFlag,
                                          Self* adoptProblem,
                                          Flag adoptFlag) 
  {
    // === create meshes ===
    if (mesh) {
      AMDIS_WARNING("mesh already created");
    }
    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
31
32
33
        mesh = adoptProblem->getMesh();
      }
      
      componentMeshes.resize(nComponents, mesh.get());
    }
34

35
36
37
38
39
40
    if (!mesh)
      AMDIS_WARNING("no mesh created");
    
    
    int globalRefinements = 0;
    Parameters::get(meshName + "->global refinements", globalRefinements);
41
42
43
44
45
46
47
48
49
50
    if (globalRefinements > 0) {
      mesh->globalRefine(globalRefinements);
      AMDIS_MSG("After refinement:");
      AMDIS_MSG("#elements = "    << mesh->size(0));
      AMDIS_MSG("#faces/edges = " << mesh->size(1));
      AMDIS_MSG("#vertices = "    << mesh->size(dim));
      AMDIS_MSG("");
    }
    
      
51
52
53
54
55
56
57
    
    // === create fespace ===
    if (feSpaces) {
      AMDIS_WARNING("feSpaces already created");
    }
    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
69
70

    if (!feSpaces)
      AMDIS_WARNING("no feSpaces created\n");
    
71
    AMDIS_MSG("#DOFs = "     << this->getFeSpace<0>().size());
72
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
86
      AMDIS_WARNING("solver already created\n");
    }
87
    else {
88
89
90
      if (initFlag.isSet(INIT_SOLVER))
        createSolver();

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

97
    if (!linearSolver) {
98
      AMDIS_WARNING("no solver created\n");
99
    }
100
101
102
103
104
105
106
107
108
109

    // === create file writer ===
    if (initFlag.isSet(INIT_FILEWRITER))
      createFileWriter();
    
    solution->compress();
  }


  template <class Traits>
110
  void ProblemStatSeq<Traits>::addMatrixOperator(shared_ptr<OperatorType> op, 
111
112
113
114
115
                                                 int i, int j,
                                                 double* factor, 
                                                 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
121
  }

122
123
124
125
126
127
128

  template <class Traits>
  void ProblemStatSeq<Traits>::addMatrixOperator(OperatorType& op, 
                                                 int i, int j,
                                                 double* factor, 
                                                 double* estFactor)
  {
129
    addMatrixOperator(make_shared<OperatorType>(op), i, j, factor, estFactor);
130
131
132
  }

  template <class Traits>
133
  void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, shared_ptr<OperatorType> > ops)
134
135
136
137
  {
    for (auto op : ops)
      addMatrixOperator(op.second, op.first.first, op.first.second);
  }
138
139
  
  template <class Traits>
140
  void ProblemStatSeq<Traits>::addVectorOperator(shared_ptr<OperatorType> op, 
141
142
143
144
                                                 int i,
                                                 double* factor, 
                                                 double* estFactor)
  {
145
146
    // TODO: currently the est-factors are ignored
    assert( estFactor == NULL );
147
    
148
    vectorOperators[i].push_back(op);
149
    vectorFactors[i].push_back(factor);
150
  }
151
152
153
154
155
156
157
  
  template <class Traits>
  void ProblemStatSeq<Traits>::addVectorOperator(OperatorType& op, 
                                                 int i,
                                                 double* factor, 
                                                 double* estFactor)
  {
158
    addVectorOperator(make_shared<OperatorType>(op), i, factor, estFactor);
159
  }
160

161
  
162
  template <class Traits>
163
  void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, shared_ptr<OperatorType> > ops)
164
165
166
167
  {
    for (auto op : ops)
      addVectorOperator(op.second, op.first);
  }
168
169
  
  /// Adds a Dirichlet boundary condition
170
  template <class Traits>
171
    template <class Predicate, class Values>
172
173
174
  void ProblemStatSeq<Traits>::addDirichletBC(Predicate const& predicate, 
                                              int row, int col, 
                                              Values const& values)
175
176
177
178
179
180
181
  {
    static_assert(Dune::Functions::Concept::isCallable<Predicate, WorldVector>(), 
      "Function passed to addDirichletBC for predicate does not model the Callable<WorldVector> concept");
    
    static_assert(Dune::Functions::Concept::isCallable<Values, WorldVector>(), 
      "Function passed to addDirichletBC for values does not model the Callable<WorldVector> concept");

182
183
    AMDIS_TEST_EXIT(row >= 0 && row < nComponents, "row number out of range: " << row);
    AMDIS_TEST_EXIT(col >= 0 && col < nComponents, "col number out of range: " << col);
184
    
185
//     boundaryConditionSet = true;
186

187
188
    using BcType = DirichletBC<WorldVector>;      
    dirichletBc[{row, col}].emplace_back( new BcType{predicate, values} );
189
190
191
  }
  
  
192
193
194
195
  template <class Traits>
  void ProblemStatSeq<Traits>::solve(AdaptInfo& adaptInfo, 
                                     bool createMatrixData, 
                                     bool storeMatrixData)
196
197
  {
    Timer t;
198
199
200
201
202
    
    SolverInfo solverInfo(name + "->solver");
    solverInfo.setCreateMatrixData(createMatrixData);
    solverInfo.setStoreMatrixData(storeMatrixData);
    
203
    solution->compress();
204
205
206
    linearSolver->solve(systemMatrix->getMatrix(), 
                        solution->getVector(), rhs->getVector(),
                        solverInfo);
207
    
208
    if (solverInfo.getInfo() > 0) {
209
210
      AMDIS_MSG("solution of discrete system needed " << t.elapsed() << " seconds");
    
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
      if (solverInfo.getAbsResidual() >= 0.0) {
        if (solverInfo.getRelResidual() >= 0.0)
          AMDIS_MSG("Residual norm: ||b-Ax|| = " << solverInfo.getAbsResidual()
                        << ", ||b-Ax||/||b|| = " << solverInfo.getRelResidual());
        else
          AMDIS_MSG("Residual norm: ||b-Ax|| = " << solverInfo.getAbsResidual());
      }
    }
    
    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() << " ";
      AMDIS_ERROR_EXIT("Tolerance " << tol_str.str() << " could not be reached!");
    }
230
231
232
  }
  
  
233
234
235
  template <class Traits>
  void ProblemStatSeq<Traits>::buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag,
                                                 bool asmMatrix, bool asmVector)
236
237
  {
    Timer t;
238
    
239
    For<0, nComponents>::loop([this](auto const _r) {
240
241
242
	this->getNodeFactory(_r).initializeIndices();
    });
    
243

244
245
    size_t nnz = 0;
    For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector](auto const _r) 
246
    {
247
      AMDIS_MSG(this->getFeSpace(_r).size() << " DOFs for FeSpace[" << _r << "]");
248
      
249
      if (asmVector) {
250
251
	rhs->compress(_r);
        rhs->getVector(_r) = 0.0;
252
253
      }

254
      For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector, _r](auto const _c) 
255
256
      {
        // The DOFMatrix which should be assembled
257
        auto& dofmatrix    = (*systemMatrix)(_r, _c);
258
259
        auto& solution_vec = (*solution)[_c];
        auto& rhs_vec      = (*rhs)[_r];
260
	  
261
        int r = 0, c = 0;
262
263
264
	if (asmMatrix) {

	  // init boundary condition
265
266
267
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
268
269
              for (auto bc : bc_list.second)
                bc->init(c == int(_c), dofmatrix, solution_vec, rhs_vec);
270
271
272
	    }
	  }
	  
273
	  // assemble the DOFMatrix block and the corresponding rhs vector, of r==c
274
	  this->assemble({int(_r), int(_c)}, &dofmatrix, (_r == _c && asmVector ? &rhs_vec : NULL));
275
276

	  // finish boundary condition
277
278
279
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
280
281
              for (auto bc : bc_list.second)
                bc->finish(c == int(_c), dofmatrix, solution_vec, rhs_vec);
282
283
	    }
	  }
284
	  
285
	  nnz += dofmatrix.getNnz();
286
	}
287
288
      });
    });
289
    
290
291
    AMDIS_MSG("fillin of assembled matrix: " << nnz);
    
292
293
294
295
    AMDIS_MSG("buildAfterCoarsen needed " << t.elapsed() << " seconds");
  }
  

296
297
  template <class Traits>
  void ProblemStatSeq<Traits>::writeFiles(AdaptInfo& adaptInfo, bool force)
298
  {
299
    Timer t;    
300
    filewriter->write(adaptInfo.getTime(), *solution);
301
    AMDIS_MSG("writeFiles needed " << t.elapsed() << " seconds");
302
303
304
  }


305
306
307
308
309
  template <class Traits>
    template <class Matrix, class Vector>
  void ProblemStatSeq<Traits>::assemble(std::pair<int, int> row_col,
                                        Matrix* dofmatrix,
                                        Vector* rhs)
310
  {
311
312
313
    auto const& rowFeSpace = dofmatrix->getRowFeSpace();
    auto const& colFeSpace = dofmatrix->getColFeSpace();
    
314
315
316
317
318
319
320
    dofmatrix->init();
    auto matrixOp = matrixOperators[row_col];
    auto matrixOpFac = matrixFactors[row_col];
    auto vectorOp = vectorOperators[row_col.first];
    auto vectorOpFac = vectorFactors[row_col.first];
    
    for (auto op : matrixOp)
321
      op->init(rowFeSpace, colFeSpace);
322
    for (auto op : vectorOp)
323
324
325
326
327
328
329
330
331
332
333
      op->init(rowFeSpace, colFeSpace);
            
    auto rowLocalView = rowFeSpace.localView();
    auto rowIndexSet = rowFeSpace.localIndexSet();
    
    auto colLocalView = colFeSpace.localView();
    auto colIndexSet = colFeSpace.localIndexSet();
    
    for (const auto& element : elements(*meshView)) {
      rowLocalView.bind(element);
      rowIndexSet.bind(rowLocalView);
334
      
335
336
      colLocalView.bind(element);
      colIndexSet.bind(colLocalView);
337
      
338
339
      if (dofmatrix) {
        ElementMatrix elementMatrix;
340
341
        bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix, 
                                    matrixOp, matrixOpFac);
342
343
344
345
346
347
        if (add)
          addElementMatrix(*dofmatrix, rowIndexSet, colIndexSet, elementMatrix);
      }
      
      if (rhs) {
        ElementVector elementVector;
348
349
        bool add = getElementVector(rowLocalView, elementVector, 
                                    vectorOp, vectorOpFac);
350
351
        if (add)
          addElementVector(*rhs, rowIndexSet, elementVector);
352
      }
353
    }
354
    dofmatrix->finish();
355
356
357
  }
  
  
358
  template <class Traits>
359
    template <class RowView, class ColView>
360
361
362
  bool ProblemStatSeq<Traits>::getElementMatrix(RowView const& rowLocalView,
                                                ColView const& colLocalView,
                                                ElementMatrix& elementMatrix,
363
                                                std::list<shared_ptr<OperatorType>>& operators,
364
                                                std::list<double*> const& factors)
365
  {
366
367
    const auto nRows = rowLocalView.tree().finiteElement().size();
    const auto nCols = colLocalView.tree().finiteElement().size();
368

369
370
371
372
    AMDIS_TEST_EXIT_DBG(operators.size() == factors.size(),
                        "Nr. of operators and factors must be consistent!"); 

    // fills the entire matrix with zeroes
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
373
    elementMatrix.setSize(nRows, nCols);
374
    elementMatrix = 0.0;
375

376
377
378
379
380
381
382
383
    bool add = false;
    
    auto op_it = operators.begin();
    auto fac_it = factors.begin();
    for (; op_it != operators.end(); ++op_it, ++fac_it) {
      bool add_op = (*op_it)->getElementMatrix(rowLocalView, colLocalView, elementMatrix, *fac_it);    
      add = add || add_op;
    }
384
    
385
    return add;
386
387
388
  }
  
  
389
  template <class Traits>
390
    template <class RowView>
391
392
  bool ProblemStatSeq<Traits>::getElementVector(RowView const& rowLocalView,
                                                ElementVector& elementVector,
393
                                                std::list<shared_ptr<OperatorType>>& operators,
394
                                                std::list<double*> const& factors)
395
  {
396
    const auto nRows = rowLocalView.tree().finiteElement().size();
397

398
399
400
    AMDIS_TEST_EXIT_DBG(operators.size() == factors.size(),
                        "Nr. of operators and factors must be consistent!"); 

401
    // Set all entries to zero
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
402
    elementVector.resize(nRows);
403
404
    elementVector = 0.0;
    
405
406
407
408
409
410
411
412
413
414
    bool add = false;
    
    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;
    }
    
    return add;
415
  }
416
  
417
  
418
  template <class Traits>
419
    template <class Matrix, class RowIndexSet, class ColIndexSet>
420
421
422
423
  void ProblemStatSeq<Traits>::addElementMatrix(Matrix& dofmatrix,
                                                RowIndexSet const& rowIndexSet,
                                                ColIndexSet const& colIndexSet,
                                                ElementMatrix const& elementMatrix)
424
425
  {
    for (size_t i = 0; i < elementMatrix.N(); ++i) {
426
427
428
429
430
431
432
      // The global index of the i−th vertex of the element
      const auto row = rowIndexSet.index(i);
      for (size_t j = 0; j < elementMatrix.M(); ++j) {
        // The global index of the j−th vertex of the element
        const auto col = colIndexSet.index(j);
        dofmatrix(row, col) += elementMatrix[i][j];
      }
433
434
435
    }
  }
  
436
  
437
  template <class Traits>
438
    template <class Vector, class IndexSet>
439
440
441
  void ProblemStatSeq<Traits>::addElementVector(Vector& dofvector,
                                                IndexSet const& indexSet,
                                                ElementVector const& elementVector)
442
443
  {
    for (size_t i = 0; i < elementVector.size(); ++i) {
444
445
446
      // The global index of the i-th vertex
      const auto row = indexSet.index(i);
      dofvector[row] += elementVector[i];
447
448
    }
  }
449
450

} // end namespace AMDiS