ProblemStat.inc.hpp 11.2 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
72
	
	solution->compress();
	rhs->compress();
73
74
75
    }
  
  
76
77
78
79
80
    template <class Param>
    void ProblemStat<Param>::addMatrixOperator(OperatorType& op, 
					       int i, int j,
					       double* factor, 
					       double* estFactor)
81
82
83
84
85
86
87
88
    {
	// TODO: currently the factors are ignored
	assert( factor == NULL && estFactor == NULL );
	
	matrixOperators[std::make_pair(i,j)].push_back(std::make_shared<OperatorType>(op));
    }

    
89
90
91
92
93
    template <class Param>
    void ProblemStat<Param>::addVectorOperator(OperatorType& op, 
					       int i,
					       double* factor, 
					       double* estFactor)
94
95
96
97
98
99
100
101
102
103
    {
	// 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
104
  template <class Param>
105
    template <class Predicate, class Values>
106
107
108
  void ProblemStat<Param>::addDirichletBC(Predicate const& predicate, 
					  int row, int col, 
					  Values const& values)
109
110
111
112
113
114
115
  {
    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");

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

    using BcType = DirichletBC<WorldVector>;
    std::shared_ptr<BcType> dirichlet = std::make_shared<BcType>(predicate, values);
      
    dirichletBc[std::make_pair(row, col)].push_back( dirichlet );
  }
  
  
128
129
130
131
  template <class Param>
  void ProblemStat<Param>::solve(AdaptInfo& adaptInfo, 
				 bool createMatrixData, 
				 bool storeMatrixData)
132
133
134
135
136
137
138
139
140
  {
//     if (!solver) {
//       WARNING("no solver\n");
//       return;
//     }

    Timer t;
    
    // NOTE: just to show the priciple, with initial solution 0
141
142
143
144
145
//     For<0, nComponents>::loop([this](auto const _i) {
// 	solution->getVector(_i) = 0.0;
//     });
    
    using VectorType = std::tuple_element_t<0, typename SystemVectorType::BaseVectors>;
146
    constexpr index_<0> _0 = index_<0>();
147
148
    
    // Technicality : turn the matrix into a linear operator
149
    Dune::MatrixAdapter<DOFMatrix, VectorType, VectorType> op(systemMatrix[_0][_0]);
150
151

    // Sequential incomplete LU decomposition as the preconditioner
152
    Dune::SeqILU0<DOFMatrix, VectorType, VectorType> ilu0(systemMatrix[_0][_0], 1.0);
153
154

    // Preconditioned conjugate−gradient solver
155
    Dune::CGSolver<VectorType> cg(op, ilu0, 1.e-4, 500, 2);
156
157

    // storing some statistics
158
    Dune::InverseOperatorResult statistics;
159
160

    // solve the linear system
161
    solution->compress();
162
    cg.apply(solution->getVector(_0), rhs->getVector(_0), statistics);
163
    
164
    AMDIS_MSG("solution of discrete system needed " << t.elapsed() << " seconds");
165

166
167
168
169
    adaptInfo.setSolverIterations(statistics.iterations);
//     adaptInfo.setMaxSolverIterations(solver->getMaxIterations());
//     adaptInfo.setSolverTolerance(solver->getTolerance());
    adaptInfo.setSolverResidual(statistics.reduction);
170
171
172
  }
  
  
173
174
175
  template <class Param>
  void ProblemStat<Param>::buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag,
					     bool asmMatrix, bool asmVector)
176
177
  {
    Timer t;
178
179
180
181
182
183
    
    For<0, nComponents>::loop([this](auto const _r) 
    {
	this->getNodeFactory(_r).initializeIndices();
    });
    
184
185
186
187

    For<0, nComponents>::loop([this, asmMatrix, asmVector](auto const _r) 
    {
      if (asmVector) {
188
189
	rhs->compress(_r);
        rhs->getVector(_r) = 0.0;
190
191
      }

192
      For<0, nComponents>::loop([this, asmMatrix, asmVector, _r](auto const _c) 
193
      {
194
195
	auto const& rowFeSpace = this->getFeSpace(_r);
        auto const& colFeSpace = this->getFeSpace(_c);
196
	
197
198
	int r = int(_r), c = int(_c);
	auto row_col = std::make_pair(r, c);
199
200
	
        // The DOFMatrix which should be assembled
201
        DOFMatrix& matrix = systemMatrix[_r][_c];
202
203
204
	  
	if (asmMatrix) {
	  // TODO: calc occupationPattern only once for each feSpace
205
206
	  Dune::MatrixIndexSet occupationPattern;
	  this->getOccupationPattern(rowFeSpace, colFeSpace, occupationPattern);
207
208
209
210
211
212
	  occupationPattern.exportIdx(matrix);

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

	  // init boundary condition
213
214
215
216
217
218
219
220
	  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));
221
222
223
	    }
	  }
	  
224
225
	  this->assemble(row_col, rowFeSpace, colFeSpace, &matrix,
		  (_r == _c && asmVector ? &rhs->getVector(_r) : NULL));
226
227

	  // finish boundary condition
228
229
230
231
232
233
234
235
	  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));
236
237
238
	    }
	  }
	}
239
240
      });
    });
241
    
242
243
244
245
246
247
248
    AMDIS_MSG("buildAfterCoarsen needed " << t.elapsed() << " seconds");
  }
  

  template <class Param>
  void ProblemStat<Param>::writeFiles(AdaptInfo& adaptInfo, double time)
  {
249
250
    Timer t;    
    filewriter->write(time, *solution);
251
    AMDIS_MSG("writeFiles needed " << t.elapsed() << " seconds");
252
253
254
255
  }
  
  
  // Get the occupation pattern
256
  template <class Param>
257
    template <class RowFeSpace, class ColFeSpace>
258
259
260
  void ProblemStat<Param>::getOccupationPattern(const RowFeSpace& rowFeSpace, 
						const ColFeSpace& colFeSpace, 
						Dune::MatrixIndexSet& nb)
261
262
263
264
  {
      nb.resize(rowFeSpace.size(), colFeSpace.size());

      // A loop over all elements of the grid
265
      auto rowLocalView = rowFeSpace.localView();
266
267
      auto rowIndexSet = rowFeSpace.localIndexSet();
      
268
      auto colLocalView = colFeSpace.localView();
269
270
271
      auto colIndexSet = colFeSpace.localIndexSet();

      for (const auto& element : elements(*meshView)) {
272
273
	  rowLocalView.bind(element);
	  rowIndexSet.bind(rowLocalView); // TODO: can this be moved out of the loop?
274
	  
275
276
	  colLocalView.bind(element);
	  colIndexSet.bind(colLocalView); // TODO: can this be moved out of the loop?
277
278
279
280
281
282
283
284
285
286
287
288
289
290
	  
	  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);
	      }
	  }
      }
  }


291
292
293
294
295
296
297
  template <class Param>
    template <class RowFeSpace, class ColFeSpace, class Vector>
  void ProblemStat<Param>::assemble(std::pair<int, int> row_col,
				    const RowFeSpace& rowFeSpace, 
				    const ColFeSpace& colFeSpace, 
				    DOFMatrix* matrix,
				    Vector* rhs)
298
  {
299
300
301
302
303
      for (auto op : matrixOperators[row_col])
	op->init(rowFeSpace, colFeSpace);
      for (auto op : vectorOperators[row_col.first])
	op->init(rowFeSpace, colFeSpace);
	      
304
      auto rowLocalView = rowFeSpace.localView();
305
306
      auto rowIndexSet = rowFeSpace.localIndexSet();
      
307
      auto colLocalView = colFeSpace.localView();
308
309
      auto colIndexSet = colFeSpace.localIndexSet();
      
310
      for (const auto& element : elements(*meshView)) {
311
312
	  rowLocalView.bind(element);
	  rowIndexSet.bind(rowLocalView);
313
	  
314
315
	  colLocalView.bind(element);
	  colIndexSet.bind(colLocalView);
316
	  
317
	  if (matrix) {
318
319
	      ElementMatrix elementMatrix;

320
321
	      bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix, matrixOperators[row_col]);
	      
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
	      if (add) {
		  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];
		      }
		  }
	      }
	  }
	  
	  if (rhs) {
	      ElementVector elementVector;
	      
338
	      bool add = getElementVector(rowLocalView, elementVector, vectorOperators[row_col.first]);
339
340
341
342
343
344
345
346
347
348
349
350
351
	      
	      if (add) {
		  for (size_t i = 0; i < elementVector.size(); ++i) {
		      // The global index of the i-th vertex
		      const auto row = rowIndexSet.index(i);
		      (*rhs)[row] += elementVector[i];
		  }
	      }
	  }
      }
  }
  
  
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
390
391
    
    return !operators.empty();
  }

} // end namespace AMDiS