ProblemStat.inc.hpp 16.2 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
    matrixChanging[{i,j}] = true;
121
122
  }

123
124
125
126
127
128
129

  template <class Traits>
  void ProblemStatSeq<Traits>::addMatrixOperator(OperatorType& op, 
                                                 int i, int j,
                                                 double* factor, 
                                                 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
154
155
156
                                                 int i,
                                                 double* factor, 
                                                 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
166
167
168
169
170
  
  template <class Traits>
  void ProblemStatSeq<Traits>::addVectorOperator(OperatorType& op, 
                                                 int i,
                                                 double* factor, 
                                                 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
190
  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;
    }
  }
  
191
  /// Adds a Dirichlet boundary condition
192
  template <class Traits>
193
    template <class Predicate, class Values>
194
195
196
  void ProblemStatSeq<Traits>::addDirichletBC(Predicate const& predicate, 
                                              int row, int col, 
                                              Values const& values)
197
  {
198
    static_assert( Concepts::Callable<Predicate, WorldVector>, 
199
200
      "Function passed to addDirichletBC for predicate does not model the Callable<WorldVector> concept");
    
201
    static_assert( Concepts::Callable<Values, WorldVector>, 
202
203
      "Function passed to addDirichletBC for values does not model the Callable<WorldVector> concept");

204
205
    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);
206
    
207
//     boundaryConditionSet = true;
208

209
210
    using BcType = DirichletBC<WorldVector>;      
    dirichletBc[{row, col}].emplace_back( new BcType{predicate, values} );
211
212
213
  }
  
  
214
215
216
217
  template <class Traits>
  void ProblemStatSeq<Traits>::solve(AdaptInfo& adaptInfo, 
                                     bool createMatrixData, 
                                     bool storeMatrixData)
218
219
  {
    Timer t;
220
221
222
223
224
    
    SolverInfo solverInfo(name + "->solver");
    solverInfo.setCreateMatrixData(createMatrixData);
    solverInfo.setStoreMatrixData(storeMatrixData);
    
225
    solution->compress();
226
227
228
    linearSolver->solve(systemMatrix->getMatrix(), 
                        solution->getVector(), rhs->getVector(),
                        solverInfo);
229
    
230
    if (solverInfo.getInfo() > 0) {
231
232
      AMDIS_MSG("solution of discrete system needed " << t.elapsed() << " seconds");
    
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
      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!");
    }
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
263
264
	this->getNodeFactory(_r).initializeIndices();
    });
    
265

266
    size_t nnz = 0;
267
    For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector_](auto const _r) 
268
    {
269
      AMDIS_MSG(this->getFeSpace(_r).size() << " DOFs for FeSpace[" << _r << "]");
270
      
271
272
      bool asmVector = asmVector_ && (!vectorAssembled[int(_r)] || vectorChanging[int(_r)]);
      
273
      if (asmVector) {
274
275
	rhs->compress(_r);
        rhs->getVector(_r) = 0.0;
276
277
      }

278
      For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector, _r](auto const _c) 
279
280
      {
        // The DOFMatrix which should be assembled
281
        auto& dofmatrix    = (*systemMatrix)(_r, _c);
282
283
        auto& solution_vec = (*solution)[_c];
        auto& rhs_vec      = (*rhs)[_r];
284
	  
285
286
287
        auto row_col = std::make_pair(int(_r), int(_c));
        bool asmMatrix = asmMatrix_ && (!matrixAssembled[row_col] || matrixChanging[row_col]);
        
288
        int r = 0, c = 0;
289
290
	if (asmMatrix) {
	  // init boundary condition
291
292
293
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
294
295
              for (auto bc : bc_list.second)
                bc->init(c == int(_c), dofmatrix, solution_vec, rhs_vec);
296
297
	    }
	  }
298
        }
299
	  
300
301
        // assemble the DOFMatrix block and the corresponding rhs vector, of r==c
        this->assemble(row_col, dofmatrix, asmMatrix, rhs_vec, (_r == _c && asmVector));
302

303
        if (asmMatrix) {
304
	  // finish boundary condition
305
306
307
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
308
309
              for (auto bc : bc_list.second)
                bc->finish(c == int(_c), dofmatrix, solution_vec, rhs_vec);
310
311
	    }
	  }
312
	  
313
	  nnz += dofmatrix.getNnz();
314
	}
315
316
      });
    });
317
    
318
319
    AMDIS_MSG("fillin of assembled matrix: " << nnz);
    
320
321
322
323
    AMDIS_MSG("buildAfterCoarsen needed " << t.elapsed() << " seconds");
  }
  

324
325
  template <class Traits>
  void ProblemStatSeq<Traits>::writeFiles(AdaptInfo& adaptInfo, bool force)
326
  {
327
    Timer t;    
328
    filewriter->write(adaptInfo.getTime(), *solution);
329
    AMDIS_MSG("writeFiles needed " << t.elapsed() << " seconds");
330
331
332
  }


333
334
335
  template <class Traits>
    template <class Matrix, class Vector>
  void ProblemStatSeq<Traits>::assemble(std::pair<int, int> row_col,
336
337
338
339
340
                                        Matrix& dofmatrix, bool asmMatrix,
                                        Vector& rhs,       bool asmVector)
  {    
    auto const& rowFeSpace = dofmatrix.getRowFeSpace();
    auto const& colFeSpace = dofmatrix.getColFeSpace();
341
    
342
343
344
345
346
    auto matrixOp = matrixOperators[row_col];
    auto matrixOpFac = matrixFactors[row_col];
    auto vectorOp = vectorOperators[row_col.first];
    auto vectorOpFac = vectorFactors[row_col.first];
    
347
348
349
350
351
352
353
354
355
356
    dofmatrix.init(asmMatrix);
    
    // return if nothing to to.
    if ((matrixOp.empty() || !asmMatrix) && 
        (vectorOp.empty()  || !asmVector))
    {
      dofmatrix.finish();
      return;
    }
    
357
    for (auto op : matrixOp)
358
      op->init(rowFeSpace, colFeSpace);
359
    for (auto op : vectorOp)
360
      op->init(rowFeSpace, colFeSpace);
361
    
362
363
364
365
366
367
    auto rowLocalView = rowFeSpace.localView();
    auto rowIndexSet = rowFeSpace.localIndexSet();
    
    auto colLocalView = colFeSpace.localView();
    auto colIndexSet = colFeSpace.localIndexSet();
    
368
    for (auto const& element : elements(*meshView)) {
369
370
      rowLocalView.bind(element);
      rowIndexSet.bind(rowLocalView);
371
      
372
373
      colLocalView.bind(element);
      colIndexSet.bind(colLocalView);
374
      
375
      if (asmMatrix) {
376
        ElementMatrix elementMatrix;
377
378
        bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix, 
                                    matrixOp, matrixOpFac);
379
        if (add)
380
          addElementMatrix(dofmatrix, rowIndexSet, colIndexSet, elementMatrix);
381
382
      }
      
383
      if (asmVector) {
384
        ElementVector elementVector;
385
386
        bool add = getElementVector(rowLocalView, elementVector, 
                                    vectorOp, vectorOpFac);
387
        if (add)
388
          addElementVector(rhs, rowIndexSet, elementVector);
389
      }
390
    }
391
392
393
394
395
396
397
398
    
    if (asmMatrix) {
      dofmatrix.finish();
      matrixAssembled[row_col] = true;
    }
    if (asmVector) {      
      vectorAssembled[row_col.first] = true;
    }
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
415
416
    AMDIS_TEST_EXIT_DBG(operators.size() == factors.size(),
                        "Nr. of operators and factors must be consistent!"); 

    // fills the entire matrix with zeroes
417
418
    elementMatrix.change_dim(nRows, nCols);
    set_to_zero(elementMatrix);
419

420
421
422
423
424
425
426
427
    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;
    }
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
444
    AMDIS_TEST_EXIT_DBG(operators.size() == factors.size(),
                        "Nr. of operators and factors must be consistent!"); 

445
    // Set all entries to zero
446
447
    elementVector.change_dim(nRows);
    set_to_zero(elementVector);
448
    
449
450
451
452
453
454
455
456
457
458
    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;
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