ProblemStat.inc.hpp 13.7 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
226
	auto const& rowFeSpace = this->getFeSpace(_r);
        auto const& colFeSpace = this->getFeSpace(_c);
227
228
	
        // The DOFMatrix which should be assembled
229
        auto& dofmatrix    = (*systemMatrix)(_r, _c);
230
231
        auto& solution_vec = (*solution)[_c];
        auto& rhs_vec      = (*rhs)[_r];
232
	  
233
        int r = 0, c = 0;
234
235
236
	if (asmMatrix) {

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

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

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


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

341
342
343
344
    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
345
    elementMatrix.setSize(nRows, nCols);
346
    elementMatrix = 0.0;
347

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

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

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

} // end namespace AMDiS