ProblemStat.inc.hpp 13.6 KB
Newer Older
1
2
#pragma once

3
4
#include "AdaptInfo.hpp"
#include "Loops.hpp"
5
6
7
8
#include "Timer.hpp"

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
110
111
112
113
114
115

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


  template <class Traits>
  void ProblemStatSeq<Traits>::addMatrixOperator(OperatorType& op, 
                                                 int i, int j,
                                                 double* factor, 
                                                 double* estFactor)
  {
    // TODO: currently the factors are ignored
116
    assert( estFactor == NULL );
117
    
118
119
    matrixOperators[{i,j}].emplace_back(new OperatorType{op});
    matrixFactors[{i,j}].push_back(factor);
120
121
122
123
124
125
126
127
128
  }

  
  template <class Traits>
  void ProblemStatSeq<Traits>::addVectorOperator(OperatorType& op, 
                                                 int i,
                                                 double* factor, 
                                                 double* estFactor)
  {
129
130
    // TODO: currently the est-factors are ignored
    assert( estFactor == NULL );
131
    
132
133
    vectorOperators[i].emplace_back(new OperatorType{op});
    vectorFactors[i].push_back(factor);
134
135
  }

136
137
138
  
  
  /// Adds a Dirichlet boundary condition
139
  template <class Traits>
140
    template <class Predicate, class Values>
141
142
143
  void ProblemStatSeq<Traits>::addDirichletBC(Predicate const& predicate, 
                                              int row, int col, 
                                              Values const& values)
144
145
146
147
148
149
150
  {
    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");

151
152
    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);
153
    
154
//     boundaryConditionSet = true;
155

156
157
    using BcType = DirichletBC<WorldVector>;      
    dirichletBc[{row, col}].emplace_back( new BcType{predicate, values} );
158
159
160
  }
  
  
161
162
163
164
  template <class Traits>
  void ProblemStatSeq<Traits>::solve(AdaptInfo& adaptInfo, 
                                     bool createMatrixData, 
                                     bool storeMatrixData)
165
166
  {
    Timer t;
167
168
169
170
171
    
    SolverInfo solverInfo(name + "->solver");
    solverInfo.setCreateMatrixData(createMatrixData);
    solverInfo.setStoreMatrixData(storeMatrixData);
    
172
    solution->compress();
173
174
175
    linearSolver->solve(systemMatrix->getMatrix(), 
                        solution->getVector(), rhs->getVector(),
                        solverInfo);
176
    
177
    AMDIS_MSG("solution of discrete system needed " << t.elapsed() << " seconds");
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
    
    if (solverInfo.getInfo() > 0) {
      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!");
    }
199
200
201
  }
  
  
202
203
204
  template <class Traits>
  void ProblemStatSeq<Traits>::buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag,
                                                 bool asmMatrix, bool asmVector)
205
206
  {
    Timer t;
207
    
208
    For<0, nComponents>::loop([this](auto const _r) {
209
210
211
	this->getNodeFactory(_r).initializeIndices();
    });
    
212

213
214
    size_t nnz = 0;
    For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector](auto const _r) 
215
    {
216
      AMDIS_MSG(this->getFeSpace(_r).size() << " DOFs for FeSpace[" << _r << "]");
217
      
218
      if (asmVector) {
219
220
	rhs->compress(_r);
        rhs->getVector(_r) = 0.0;
221
222
      }

223
      For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector, _r](auto const _c) 
224
225
      {
        // The DOFMatrix which should be assembled
226
        auto& dofmatrix    = (*systemMatrix)(_r, _c);
227
228
        auto& solution_vec = (*solution)[_c];
        auto& rhs_vec      = (*rhs)[_r];
229
	  
230
        int r = 0, c = 0;
231
232
233
	if (asmMatrix) {

	  // init boundary condition
234
235
236
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
237
238
              for (auto bc : bc_list.second)
                bc->init(c == int(_c), dofmatrix, solution_vec, rhs_vec);
239
240
241
	    }
	  }
	  
242
	  // assemble the DOFMatrix block and the corresponding rhs vector, of r==c
243
	  this->assemble({int(_r), int(_c)}, &dofmatrix, (_r == _c && asmVector ? &rhs_vec : NULL));
244
245

	  // finish boundary condition
246
247
248
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
249
250
              for (auto bc : bc_list.second)
                bc->finish(c == int(_c), dofmatrix, solution_vec, rhs_vec);
251
252
	    }
	  }
253
	  
254
	  nnz += dofmatrix.getNnz();
255
	}
256
257
      });
    });
258
    
259
260
    AMDIS_MSG("fillin of assembled matrix: " << nnz);
    
261
262
263
264
    AMDIS_MSG("buildAfterCoarsen needed " << t.elapsed() << " seconds");
  }
  

265
266
  template <class Traits>
  void ProblemStatSeq<Traits>::writeFiles(AdaptInfo& adaptInfo, bool force)
267
  {
268
    Timer t;    
269
    filewriter->write(adaptInfo.getTime(), *solution);
270
    AMDIS_MSG("writeFiles needed " << t.elapsed() << " seconds");
271
272
273
  }


274
275
276
277
278
  template <class Traits>
    template <class Matrix, class Vector>
  void ProblemStatSeq<Traits>::assemble(std::pair<int, int> row_col,
                                        Matrix* dofmatrix,
                                        Vector* rhs)
279
  {
280
281
282
    auto const& rowFeSpace = dofmatrix->getRowFeSpace();
    auto const& colFeSpace = dofmatrix->getColFeSpace();
    
283
284
285
286
287
288
289
    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)
290
      op->init(rowFeSpace, colFeSpace);
291
    for (auto op : vectorOp)
292
293
294
295
296
297
298
299
300
301
302
      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);
303
      
304
305
      colLocalView.bind(element);
      colIndexSet.bind(colLocalView);
306
      
307
308
      if (dofmatrix) {
        ElementMatrix elementMatrix;
309
310
        bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix, 
                                    matrixOp, matrixOpFac);
311
312
313
314
315
316
        if (add)
          addElementMatrix(*dofmatrix, rowIndexSet, colIndexSet, elementMatrix);
      }
      
      if (rhs) {
        ElementVector elementVector;
317
318
        bool add = getElementVector(rowLocalView, elementVector, 
                                    vectorOp, vectorOpFac);
319
320
        if (add)
          addElementVector(*rhs, rowIndexSet, elementVector);
321
      }
322
    }
323
    dofmatrix->finish();
324
325
326
  }
  
  
327
  template <class Traits>
328
    template <class RowView, class ColView>
329
330
331
  bool ProblemStatSeq<Traits>::getElementMatrix(RowView const& rowLocalView,
                                                ColView const& colLocalView,
                                                ElementMatrix& elementMatrix,
332
333
                                                std::list<std::shared_ptr<OperatorType>>& operators,
                                                std::list<double*> const& factors)
334
  {
335
336
    const auto nRows = rowLocalView.tree().finiteElement().size();
    const auto nCols = colLocalView.tree().finiteElement().size();
337

338
339
340
341
    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
342
    elementMatrix.setSize(nRows, nCols);
343
    elementMatrix = 0.0;
344

345
346
347
348
349
350
351
352
    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;
    }
353
    
354
    return add;
355
356
357
  }
  
  
358
  template <class Traits>
359
    template <class RowView>
360
361
  bool ProblemStatSeq<Traits>::getElementVector(RowView const& rowLocalView,
                                                ElementVector& elementVector,
362
363
                                                std::list<std::shared_ptr<OperatorType>>& operators,
                                                std::list<double*> const& factors)
364
  {
365
    const auto nRows = rowLocalView.tree().finiteElement().size();
366

367
368
369
    AMDIS_TEST_EXIT_DBG(operators.size() == factors.size(),
                        "Nr. of operators and factors must be consistent!"); 

370
    // Set all entries to zero
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
371
    elementVector.resize(nRows);
372
373
    elementVector = 0.0;
    
374
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)->getElementVector(rowLocalView, elementVector, *fac_it);
      add = add || add_op;
    }
    
    return add;
384
  }
385
  
386
  
387
  template <class Traits>
388
    template <class Matrix, class RowIndexSet, class ColIndexSet>
389
390
391
392
  void ProblemStatSeq<Traits>::addElementMatrix(Matrix& dofmatrix,
                                                RowIndexSet const& rowIndexSet,
                                                ColIndexSet const& colIndexSet,
                                                ElementMatrix const& elementMatrix)
393
394
  {
    for (size_t i = 0; i < elementMatrix.N(); ++i) {
395
396
397
398
399
400
401
      // 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];
      }
402
403
404
    }
  }
  
405
  
406
  template <class Traits>
407
    template <class Vector, class IndexSet>
408
409
410
  void ProblemStatSeq<Traits>::addElementVector(Vector& dofvector,
                                                IndexSet const& indexSet,
                                                ElementVector const& elementVector)
411
412
  {
    for (size_t i = 0; i < elementVector.size(); ++i) {
413
414
415
      // The global index of the i-th vertex
      const auto row = indexSet.index(i);
      dofvector[row] += elementVector[i];
416
417
    }
  }
418
419

} // end namespace AMDiS