ProblemStat.inc.hpp 12.5 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
21
22
23
24
  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) &&
          (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) 
      {
        createMesh();
      }
25

26
27
28
29
30
31
32
33
34
35
      if (adoptProblem &&
          (adoptFlag.isSet(INIT_MESH) ||
          adoptFlag.isSet(INIT_SYSTEM) ||
          adoptFlag.isSet(INIT_FE_SPACE))) 
      {
        mesh = adoptProblem->getMesh();
      }
      
      componentMeshes.resize(nComponents, mesh.get());
    }
36

37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
    if (!mesh)
      AMDIS_WARNING("no mesh created");
    
    
    int globalRefinements = 0;
    Parameters::get(meshName + "->global refinements", globalRefinements);
    mesh->globalRefine(globalRefinements);
    
    // === create fespace ===
    if (feSpaces) {
      AMDIS_WARNING("feSpaces already created");
    }
    else {
      if (initFlag.isSet(INIT_FE_SPACE) ||
          (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) 
      {
        createFeSpaces();
      }
55

56
57
58
59
60
      if (adoptProblem &&
          (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) 
      {
        feSpaces = adoptProblem->getFeSpaces();
      }
61
    }
62
63
64
65
66
67
68
69
70
71

    if (!feSpaces)
      AMDIS_WARNING("no feSpaces created\n");
    
    
    // === create system ===
    if (initFlag.isSet(INIT_SYSTEM))
      createMatricesAndVectors();

    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM))
72
    {
73
74
75
      solution = adoptProblem->getSolution();
      rhs = adoptProblem->getRhs();
      systemMatrix = adoptProblem->getSystemMatrix();
76
77
    }

78
79
    // === create solver ===
    if (linearSolver)
80
    {
81
82
83
84
85
86
87
88
89
90
91
92
      AMDIS_WARNING("solver already created\n");
    }
    else
    {
      if (initFlag.isSet(INIT_SOLVER))
        createSolver();

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

95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
    if (!linearSolver)
      AMDIS_WARNING("no solver created\n");

    // === 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
    assert( factor == NULL && estFactor == NULL );
    
    matrixOperators[std::make_pair(i,j)].push_back(std::make_shared<OperatorType>(op));
  }

  
  template <class Traits>
  void ProblemStatSeq<Traits>::addVectorOperator(OperatorType& op, 
                                                 int i,
                                                 double* factor, 
                                                 double* estFactor)
  {
    // TODO: currently the factors are ignored
    assert( factor == NULL && estFactor == NULL );
    
    vectorOperators[i].push_back(std::make_shared<OperatorType>(op));
  }

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

146
147
    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);
148
    
149
//     boundaryConditionSet = true;
150
151
152
153
154
155
156
157

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

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

222
      For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector, _r](auto const _c) 
223
      {
224
225
	auto const& rowFeSpace = this->getFeSpace(_r);
        auto const& colFeSpace = this->getFeSpace(_c);
226
	
227
228
	int r = int(_r), c = int(_c);
	auto row_col = std::make_pair(r, c);
229
230
	
        // The DOFMatrix which should be assembled
231
232
233
        auto& dofmatrix    = (*systemMatrix)(_r, _c);
        auto& solution_vec = solution->getDOFVector(_c);
        auto& rhs_vec      = rhs->getDOFVector(_r);
234
235
	  
	if (asmMatrix) {
236
	  dofmatrix.init();
237
238

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

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

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


280
281
282
283
284
  template <class Traits>
    template <class Matrix, class Vector>
  void ProblemStatSeq<Traits>::assemble(std::pair<int, int> row_col,
                                        Matrix* dofmatrix,
                                        Vector* rhs)
285
  {
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
    auto const& rowFeSpace = dofmatrix->getRowFeSpace();
    auto const& colFeSpace = dofmatrix->getColFeSpace();
    
    for (auto op : matrixOperators[row_col])
      op->init(rowFeSpace, colFeSpace);
    for (auto op : vectorOperators[row_col.first])
      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
309
310
311
312
313
314
315
316
317
318
      if (dofmatrix) {
        ElementMatrix elementMatrix;
        bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix, matrixOperators[row_col]);
        if (add)
          addElementMatrix(*dofmatrix, rowIndexSet, colIndexSet, elementMatrix);
      }
      
      if (rhs) {
        ElementVector elementVector;
        bool add = getElementVector(rowLocalView, elementVector, vectorOperators[row_col.first]);
        if (add)
          addElementVector(*rhs, rowIndexSet, elementVector);
319
      }
320
    }
321
322
323
  }
  
  
324
  template <class Traits>
325
    template <class RowView, class ColView>
326
327
328
329
  bool ProblemStatSeq<Traits>::getElementMatrix(RowView const& rowLocalView,
                                                ColView const& colLocalView,
                                                ElementMatrix& elementMatrix,
                                                std::list<std::shared_ptr<OperatorType>>& operators)
330
  {
331
332
    const auto nRows = rowLocalView.tree().finiteElement().size();
    const auto nCols = colLocalView.tree().finiteElement().size();
333

Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
334
    elementMatrix.setSize(nRows, nCols);
335
336
337
    elementMatrix = 0.0; // fills the entire matrix with zeroes

    for (auto op : operators)
338
      op->getElementMatrix(rowLocalView, colLocalView, elementMatrix);
339
340
341
342
343
344
    
    return !operators.empty();
  }
  
  
  
345
  template <class Traits>
346
    template <class RowView>
347
348
349
  bool ProblemStatSeq<Traits>::getElementVector(RowView const& rowLocalView,
                                                ElementVector& elementVector,
                                                std::list<std::shared_ptr<OperatorType>>& operators)
350
  {
351
    const auto nRows = rowLocalView.tree().finiteElement().size();
352
353

    // Set all entries to zero
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
354
    elementVector.resize(nRows);
355
356
357
    elementVector = 0.0;

    for (auto op : operators)
358
      op->getElementVector(rowLocalView, elementVector);
359
360
361
    
    return !operators.empty();
  }
362
  
363
  template <class Traits>
364
    template <class Matrix, class RowIndexSet, class ColIndexSet>
365
366
367
368
  void ProblemStatSeq<Traits>::addElementMatrix(Matrix& dofmatrix,
                                                RowIndexSet const& rowIndexSet,
                                                ColIndexSet const& colIndexSet,
                                                ElementMatrix const& elementMatrix)
369
370
  {
    for (size_t i = 0; i < elementMatrix.N(); ++i) {
371
372
373
374
375
376
377
      // 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];
      }
378
379
380
    }
  }
  
381
  template <class Traits>
382
    template <class Vector, class IndexSet>
383
384
385
  void ProblemStatSeq<Traits>::addElementVector(Vector& dofvector,
                                                IndexSet const& indexSet,
                                                ElementVector const& elementVector)
386
387
  {
    for (size_t i = 0; i < elementVector.size(); ++i) {
388
389
390
      // The global index of the i-th vertex
      const auto row = indexSet.index(i);
      dofvector[row] += elementVector[i];
391
392
    }
  }
393
394

} // end namespace AMDiS