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


4
5
6
#include <dune/istl/solvers.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
7
#include <dune/istl/matrixmarket.hh>
8
9
10

#include "AdaptInfo.hpp"
#include "Loops.hpp"
11
12
13
14
#include "Timer.hpp"

namespace AMDiS
{
15
16
17
18
    template <class Param>
    void ProblemStat<Param>::initialize(Flag initFlag,
					ProblemStat* adoptProblem,
					Flag adoptFlag) 
19
20
21
    {
	// === create meshes ===
	if (mesh) {
22
	    AMDIS_WARNING("mesh already created");
23
24
25
26
	}
	else {
	    if (initFlag.isSet(CREATE_MESH) ||
		(!adoptFlag.isSet(INIT_MESH) &&
27
28
		(initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) 
	    {
29
30
31
32
33
34
		createMesh();
	    }

	    if (adoptProblem &&
		(adoptFlag.isSet(INIT_MESH) ||
		adoptFlag.isSet(INIT_SYSTEM) ||
35
36
		adoptFlag.isSet(INIT_FE_SPACE))) 
	    {
37
38
39
40
41
42
43
		mesh = adoptProblem->getMesh();
	    }
	    
	    componentMeshes.resize(nComponents, mesh.get());
	}

	if (!mesh) {
44
	    AMDIS_WARNING("no mesh created");
45
46
	}
	
47
48
49
	int globalRefinements = 0;
	Parameters::get(meshName + "->global refinements", globalRefinements);
	mesh->globalRefine(globalRefinements);
50
51
52
	
	// === create fespace ===
	if (feSpaces) {
53
	    AMDIS_WARNING("feSpaces already created");
54
55
56
	}
	else {
	    if (initFlag.isSet(INIT_FE_SPACE) ||
57
58
		(initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) 
	    {
59
60
61
62
		createFeSpaces();
	    }

	    if (adoptProblem &&
63
64
		(adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) 
	    {
65
66
67
68
		feSpaces = adoptProblem->getFeSpaces();
	    }
	}
	
69
	filewriter = std::make_shared<FileWriter<Param>>(name + "->output", *meshView, componentNames);
70
71
	
	solution->compress();
72
73
74
    }
  
  
75
76
77
78
79
    template <class Param>
    void ProblemStat<Param>::addMatrixOperator(OperatorType& op, 
					       int i, int j,
					       double* factor, 
					       double* estFactor)
80
81
82
83
84
85
86
87
    {
	// TODO: currently the factors are ignored
	assert( factor == NULL && estFactor == NULL );
	
	matrixOperators[std::make_pair(i,j)].push_back(std::make_shared<OperatorType>(op));
    }

    
88
89
90
91
92
    template <class Param>
    void ProblemStat<Param>::addVectorOperator(OperatorType& op, 
					       int i,
					       double* factor, 
					       double* estFactor)
93
94
95
96
97
98
99
100
101
102
    {
	// TODO: currently the factors are ignored
	assert( factor == NULL && estFactor == NULL );
	
	vectorOperators[i].push_back(std::make_shared<OperatorType>(op));
    }

  
  
  /// Adds a Dirichlet boundary condition
103
  template <class Param>
104
    template <class Predicate, class Values>
105
106
107
  void ProblemStat<Param>::addDirichletBC(Predicate const& predicate, 
					  int row, int col, 
					  Values const& values)
108
109
110
111
112
113
114
  {
    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");

115
116
    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);
117
    
118
//     boundaryConditionSet = true;
119
120
121
122
123
124
125
126

    using BcType = DirichletBC<WorldVector>;
    std::shared_ptr<BcType> dirichlet = std::make_shared<BcType>(predicate, values);
      
    dirichletBc[std::make_pair(row, col)].push_back( dirichlet );
  }
  
  
127
128
129
130
  template <class Param>
  void ProblemStat<Param>::solve(AdaptInfo& adaptInfo, 
				 bool createMatrixData, 
				 bool storeMatrixData)
131
  {
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
      using Vector = typename SystemVectorType::MultiVector;
      Dune::Richardson<Vector, Vector> precon(1.0);
      
      solveAdvanced(adaptInfo, precon);
  }
  
  template <class Param>
  void ProblemStat<Param>::solveAdvanced(AdaptInfo& adaptInfo, 
					 BasePreconditioner& precon)
  {
    
    size_t maxIter = 500;
    Parameters::get(name + "->solver->max iteration", maxIter);
    
    double relTol = 1.e-6;
//     Parameters::get(name + "->solver->relative tolerance", relTol);
    Parameters::get(name + "->solver->reduction", relTol);
    
    int info = 2;
    Parameters::get(name + "->solver->info", info);
    
    int restart = 30;
    Parameters::get(name + "->solver->restart", restart);
155
156

    Timer t;
157
158
    using Matrix = typename SystemMatrixType::MultiMatrix;
    using Vector = typename SystemVectorType::MultiVector;
159
    
160
161
    // Technicality : turn the matrix into a linear operator
    Dune::MatrixAdapter<Matrix, Vector, Vector> op(systemMatrix->getMatrix());
162
    
163
    PreconAdapter<BasePreconditioner> precon_adapter(precon);
164
165
    
    // Preconditioned conjugate−gradient solver
166
    Dune::RestartedGMResSolver<Vector> solver(op, precon_adapter, relTol, restart, maxIter, info);
167
168

    // storing some statistics
169
    Dune::InverseOperatorResult statistics;
170
171

    // solve the linear system
172
    solution->compress();
173
    solver.apply(solution->getVector(), rhs->getVector(), statistics);
174
    
175
    AMDIS_MSG("solution of discrete system needed " << t.elapsed() << " seconds");
176

177
178
179
180
    adaptInfo.setSolverIterations(statistics.iterations);
//     adaptInfo.setMaxSolverIterations(solver->getMaxIterations());
//     adaptInfo.setSolverTolerance(solver->getTolerance());
    adaptInfo.setSolverResidual(statistics.reduction);
181
182
183
  }
  
  
184
185
186
  template <class Param>
  void ProblemStat<Param>::buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag,
					     bool asmMatrix, bool asmVector)
187
188
  {
    Timer t;
189
190
191
192
193
194
    
    For<0, nComponents>::loop([this](auto const _r) 
    {
	this->getNodeFactory(_r).initializeIndices();
    });
    
195

196
197
    size_t nnz = 0;
    For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector](auto const _r) 
198
    {
199
200
      AMDIS_MSG(this->getFeSpace(_r).size() << " DOFs for " << "FeSpace[" << _r << "]");
      
201
      if (asmVector) {
202
203
	rhs->compress(_r);
        rhs->getVector(_r) = 0.0;
204
205
      }

206
      For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector, _r](auto const _c) 
207
      {
208
209
	auto const& rowFeSpace = this->getFeSpace(_r);
        auto const& colFeSpace = this->getFeSpace(_c);
210
	
211
212
	int r = int(_r), c = int(_c);
	auto row_col = std::make_pair(r, c);
213
214
	
        // The DOFMatrix which should be assembled
215
        auto& matrix = systemMatrix->getMatrix(_r, _c);
216
217
218
	  
	if (asmMatrix) {
	  // TODO: calc occupationPattern only once for each feSpace
219
220
	  Dune::MatrixIndexSet occupationPattern;
	  this->getOccupationPattern(rowFeSpace, colFeSpace, occupationPattern);
221
222
223
224
225
226
	  occupationPattern.exportIdx(matrix);

	  // TODO: only set to zero if not keep the matrix
	  matrix = 0.0;

	  // init boundary condition
227
228
229
230
231
232
233
234
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
		for (auto bc : bc_list.second)
		    bc->init(c == int(_c), 
			     rowFeSpace, colFeSpace, &matrix, 
			     &rhs->getVector(_r), 
			     &solution->getVector(_c));
235
236
237
	    }
	  }
	  
238
239
	  this->assemble(row_col, rowFeSpace, colFeSpace, &matrix,
		  (_r == _c && asmVector ? &rhs->getVector(_r) : NULL));
240
241

	  // finish boundary condition
242
243
244
245
246
247
248
249
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
		for (auto bc : bc_list.second)
		    bc->finish(c == int(_c), 
			       rowFeSpace, colFeSpace, &matrix, 
			       &rhs->getVector(_r), 
			       &solution->getVector(_c));
250
251
	    }
	  }
252
253
	  
	  nnz += matrix.nonzeroes();
254
	}
255
256
      });
    });
257
    
258
259
    AMDIS_MSG("fillin of assembled matrix: " << nnz);
    
260
261
262
263
264
265
266
    AMDIS_MSG("buildAfterCoarsen needed " << t.elapsed() << " seconds");
  }
  

  template <class Param>
  void ProblemStat<Param>::writeFiles(AdaptInfo& adaptInfo, double time)
  {
267
268
    Timer t;    
    filewriter->write(time, *solution);
269
    AMDIS_MSG("writeFiles needed " << t.elapsed() << " seconds");
270
271
272
273
  }
  
  
  // Get the occupation pattern
274
  template <class Param>
275
    template <class RowFeSpace, class ColFeSpace>
276
277
278
  void ProblemStat<Param>::getOccupationPattern(const RowFeSpace& rowFeSpace, 
						const ColFeSpace& colFeSpace, 
						Dune::MatrixIndexSet& nb)
279
280
281
282
  {
      nb.resize(rowFeSpace.size(), colFeSpace.size());

      // A loop over all elements of the grid
283
      auto rowLocalView = rowFeSpace.localView();
284
285
      auto rowIndexSet = rowFeSpace.localIndexSet();
      
286
      auto colLocalView = colFeSpace.localView();
287
288
289
      auto colIndexSet = colFeSpace.localIndexSet();

      for (const auto& element : elements(*meshView)) {
290
291
	  rowLocalView.bind(element);
	  rowIndexSet.bind(rowLocalView); // TODO: can this be moved out of the loop?
292
	  
293
294
	  colLocalView.bind(element);
	  colIndexSet.bind(colLocalView); // TODO: can this be moved out of the loop?
295
296
297
298
299
300
301
302
303
304
305
306
307
308
	  
	  for (size_t i = 0; i < rowIndexSet.size(); ++i) {
	      // The global index of the i-th vertex of the element
	      auto row = rowIndexSet.index(i);
	      for (size_t j = 0; j < colIndexSet.size(); ++j) {
		  // The global index of the j-th vertex of the element
		  auto col = colIndexSet.index(j);
		  nb.add(row, col);
	      }
	  }
      }
  }


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

Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
362
    elementMatrix.setSize(nRows, nCols);
363
364
365
    elementMatrix = 0.0; // fills the entire matrix with zeroes

    for (auto op : operators)
366
      op->getElementMatrix(rowLocalView, colLocalView, elementMatrix);
367
368
369
370
371
372
    
    return !operators.empty();
  }
  
  
  
373
  template <class Param>
374
    template <class RowView>
375
  bool ProblemStat<Param>::getElementVector(RowView const& rowLocalView,
376
377
					    ElementVector& elementVector,
					    std::list<std::shared_ptr<OperatorType>>& operators)
378
  {
379
    const auto nRows = rowLocalView.tree().finiteElement().size();
380
381

    // Set all entries to zero
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
382
    elementVector.resize(nRows);
383
384
385
    elementVector = 0.0;

    for (auto op : operators)
386
      op->getElementVector(rowLocalView, elementVector);
387
388
389
    
    return !operators.empty();
  }
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
  
  template <class Param>
    template <class Matrix, class RowIndexSet, class ColIndexSet>
  void ProblemStat<Param>::addElementMatrix(Matrix& matrix,
					    RowIndexSet const& rowIndexSet,
					    ColIndexSet const& colIndexSet,
					    ElementMatrix const& elementMatrix)
  {
    for (size_t i = 0; i < elementMatrix.N(); ++i) {
	// 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);
	    matrix[row][col] += elementMatrix[i][j];
	}
    }
  }
  
  template <class Param>
    template <class Vector, class IndexSet>
  void ProblemStat<Param>::addElementVector(Vector& vector,
					    IndexSet const& indexSet,
					    ElementVector const& elementVector)
  {
    for (size_t i = 0; i < elementVector.size(); ++i) {
	// The global index of the i-th vertex
	const auto row = indexSet.index(i);
	vector[row] += elementVector[i];
    }
  }
421
422

} // end namespace AMDiS