All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
fvbaselinearizer.hh
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
28 #ifndef EWOMS_FV_BASE_LINEARIZER_HH
29 #define EWOMS_FV_BASE_LINEARIZER_HH
30 
31 #include "fvbaseproperties.hh"
32 
37 
38 #include <opm/common/ErrorMacros.hpp>
39 #include <opm/common/Exceptions.hpp>
40 
41 #include <dune/common/version.hh>
42 #include <dune/common/fvector.hh>
43 #include <dune/common/fmatrix.hh>
44 
45 #include <type_traits>
46 #include <iostream>
47 #include <vector>
48 #include <set>
49 
50 namespace Ewoms {
51 // forward declarations
52 template<class TypeTag>
53 class EcfvDiscretization;
54 
64 template<class TypeTag>
66 {
68  typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
69  typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
70  typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
71  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
72  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
73  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
74  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
75  typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
76  typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
77  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
78 
79  typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
80  typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
81  typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
82  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
83  typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints;
84  typedef typename GET_PROP_TYPE(TypeTag, Stencil) Stencil;
85  typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager;
86 
87  typedef typename GET_PROP_TYPE(TypeTag, GridCommHandleFactory) GridCommHandleFactory;
88 
89  typedef Opm::MathToolbox<Evaluation> Toolbox;
90 
91  typedef typename GridView::template Codim<0>::Entity Element;
92  typedef typename GridView::template Codim<0>::Iterator ElementIterator;
93 
94  typedef GlobalEqVector Vector;
95  typedef JacobianMatrix Matrix;
96 
97  enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
98  enum { historySize = GET_PROP_VALUE(TypeTag, TimeDiscHistorySize) };
99 
100  typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
101  typedef Dune::FieldVector<Scalar, numEq> VectorBlock;
102 
103  static const bool linearizeNonLocalElements = GET_PROP_VALUE(TypeTag, LinearizeNonLocalElements);
104 
105  // copying the linearizer is not a good idea
108 
109 public:
111  {
112  simulatorPtr_ = 0;
113 
114  matrix_ = 0;
115  }
116 
118  {
119  delete matrix_;
120  auto it = elementCtx_.begin();
121  const auto& endIt = elementCtx_.end();
122  for (; it != endIt; ++it)
123  delete *it;
124  }
125 
129  static void registerParameters()
130  { }
131 
141  void init(Simulator& simulator)
142  {
143  simulatorPtr_ = &simulator;
144  delete matrix_; // <- note that this even works for nullpointers!
145  matrix_ = 0;
146  }
147 
155  void eraseMatrix()
156  {
157  delete matrix_; // <- note that this even works for nullpointers!
158  matrix_ = 0;
159  }
160 
170  void linearize()
171  {
172  // we defer the initialization of the Jacobian matrix until here because the
173  // auxiliary modules usually assume the problem, model and grid to be fully
174  // initialized...
175  if (!matrix_)
176  initFirstIteration_();
177 
178  int succeeded;
179  try {
180  linearize_();
181  succeeded = 1;
182  }
183  catch (const std::exception& e)
184  {
185  std::cout << "rank " << simulator_().gridView().comm().rank()
186  << " caught an exception while linearizing:" << e.what()
187  << "\n" << std::flush;
188  succeeded = 0;
189  }
190 #if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,5)
191  catch (const Dune::Exception& e)
192  {
193  std::cout << "rank " << simulator_().gridView().comm().rank()
194  << " caught an exception while linearizing:" << e.what()
195  << "\n" << std::flush;
196  succeeded = 0;
197  }
198 #endif
199  catch (...)
200  {
201  std::cout << "rank " << simulator_().gridView().comm().rank()
202  << " caught an exception while linearizing"
203  << "\n" << std::flush;
204  succeeded = 0;
205  }
206  succeeded = gridView_().comm().min(succeeded);
207 
208  if (!succeeded) {
209  OPM_THROW(Opm::NumericalProblem,
210  "A process did not succeed in linearizing the system");
211  }
212  }
213 
217  const Matrix& matrix() const
218  { return *matrix_; }
219 
220  Matrix& matrix()
221  { return *matrix_; }
222 
226  const GlobalEqVector& residual() const
227  { return residual_; }
228 
229  GlobalEqVector& residual()
230  { return residual_; }
231 
237  const std::map<unsigned, Constraints>& constraintsMap() const
238  { return constraintsMap_; }
239 
240 private:
241  Simulator& simulator_()
242  { return *simulatorPtr_; }
243  const Simulator& simulator_() const
244  { return *simulatorPtr_; }
245 
246  Problem& problem_()
247  { return simulator_().problem(); }
248  const Problem& problem_() const
249  { return simulator_().problem(); }
250 
251  Model& model_()
252  { return simulator_().model(); }
253  const Model& model_() const
254  { return simulator_().model(); }
255 
256  const GridView& gridView_() const
257  { return problem_().gridView(); }
258 
259  const ElementMapper& elementMapper_() const
260  { return model_().elementMapper(); }
261 
262  const DofMapper& dofMapper_() const
263  { return model_().dofMapper(); }
264 
265  void initFirstIteration_()
266  {
267  // initialize the BCRS matrix for the Jacobian of the residual function
268  createMatrix_();
269 
270  // initialize the Jacobian matrix and the vector for the residual function
271  (*matrix_) = 0.0;
272  residual_.resize(model_().numTotalDof());
273  residual_ = 0;
274 
275  // create the per-thread context objects
276  elementCtx_.resize(ThreadManager::maxThreads());
277  for (unsigned threadId = 0; threadId != ThreadManager::maxThreads(); ++ threadId)
278  elementCtx_[threadId] = new ElementContext(simulator_());
279  }
280 
281  // Construct the BCRS matrix for the Jacobian of the residual function
282  void createMatrix_()
283  {
284  size_t numAllDof = model_().numTotalDof();
285 
286  // allocate raw matrix
287  matrix_ = new Matrix(numAllDof, numAllDof, Matrix::random);
288 
289  Stencil stencil(gridView_(), model_().dofMapper() );
290 
291  // for the main model, find out the global indices of the neighboring degrees of
292  // freedom of each primary degree of freedom
293  typedef std::set<unsigned> NeighborSet;
294  std::vector<NeighborSet> neighbors(numAllDof);
295  ElementIterator elemIt = gridView_().template begin<0>();
296  const ElementIterator elemEndIt = gridView_().template end<0>();
297  for (; elemIt != elemEndIt; ++elemIt) {
298  const Element& elem = *elemIt;
299  stencil.update(elem);
300 
301  for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
302  unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
303 
304  for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
305  unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
306  neighbors[myIdx].insert(neighborIdx);
307  }
308  }
309  }
310 
311  // add the additional neighbors and degrees of freedom caused by the auxiliary
312  // equations
313  const auto& model = model_();
314  size_t numAuxMod = model.numAuxiliaryModules();
315  for (unsigned auxModIdx = 0; auxModIdx < numAuxMod; ++auxModIdx)
316  model.auxiliaryModule(auxModIdx)->addNeighbors(neighbors);
317 
318  // allocate space for the rows of the matrix
319  for (unsigned dofIdx = 0; dofIdx < numAllDof; ++ dofIdx)
320  matrix_->setrowsize(dofIdx, neighbors[dofIdx].size());
321  matrix_->endrowsizes();
322 
323  // fill the rows with indices. each degree of freedom talks to
324  // all of its neighbors. (it also talks to itself since
325  // degrees of freedom are sometimes quite egocentric.)
326  for (unsigned dofIdx = 0; dofIdx < numAllDof; ++ dofIdx) {
327  typename NeighborSet::iterator nIt = neighbors[dofIdx].begin();
328  typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end();
329  for (; nIt != nEndIt; ++nIt)
330  matrix_->addindex(dofIdx, *nIt);
331  }
332  matrix_->endindices();
333  }
334 
335  // reset the global linear system of equations.
336  void resetSystem_()
337  {
338  residual_ = 0.0;
339  (*matrix_) = 0.0;
340  }
341 
342  // query the problem for all constraint degrees of freedom. note that this method is
343  // quite involved and is thus relatively slow.
344  void updateConstraintsMap_()
345  {
346  if (!enableConstraints_())
347  // constraints are not explictly enabled, so we don't need to consider them!
348  return;
349 
350  unsigned threadId = ThreadManager::threadId();
351  constraintsMap_.clear();
352 
353  // loop over all elements...
354  ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_());
355 #ifdef _OPENMP
356 #pragma omp parallel
357 #endif
358  {
359  ElementIterator elemIt = threadedElemIt.beginParallel();
360  for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
361  // create an element context (the solution-based quantities are not
362  // available here!)
363  const Element& elem = *elemIt;
364  ElementContext& elemCtx = *elementCtx_[threadId];
365  elemCtx.updateStencil(elem);
366 
367  // check if the problem wants to constrain any degree of the current
368  // element's freedom. if yes, add the constraint to the map.
369  for (unsigned primaryDofIdx = 0;
370  primaryDofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0);
371  ++ primaryDofIdx)
372  {
373  Constraints constraints;
374  elemCtx.problem().constraints(constraints,
375  elemCtx,
376  primaryDofIdx,
377  /*timeIdx=*/0);
378  if (constraints.isActive()) {
379  unsigned globI = elemCtx.globalSpaceIndex(primaryDofIdx, /*timeIdx=*/0);
380  constraintsMap_[globI] = constraints;
381  continue;
382  }
383  }
384  }
385  }
386  }
387 
388  // linearize the whole system
389  void linearize_()
390  {
391  resetSystem_();
392 
393  // before the first iteration of each time step, we need to update the
394  // constraints. (i.e., we assume that constraints can be time dependent, but they
395  // can't depend on the solution.)
396  if (model_().newtonMethod().numIterations() == 0)
397  updateConstraintsMap_();
398 
399  applyConstraintsToSolution_();
400 
401  *matrix_ = 0.0;
402 
403  // relinearize the elements...
404  ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_());
405 #ifdef _OPENMP
406 #pragma omp parallel
407 #endif
408  {
409  ElementIterator elemIt = threadedElemIt.beginParallel();
410  ElementIterator nextElemIt = elemIt;
411  for (; !threadedElemIt.isFinished(elemIt); elemIt = nextElemIt) {
412  // give the model and the problem a chance to prefetch the data required
413  // to linearize the next element, but only if we need to consider it
414  nextElemIt = threadedElemIt.increment();
415  if (!threadedElemIt.isFinished(nextElemIt)) {
416  const auto& nextElem = *nextElemIt;
417  if (linearizeNonLocalElements
418  || nextElem.partitionType() == Dune::InteriorEntity)
419  {
420  model_().prefetch(nextElem);
421  problem_().prefetch(nextElem);
422  }
423  }
424 
425  const Element& elem = *elemIt;
426  if (!linearizeNonLocalElements && elem.partitionType() != Dune::InteriorEntity)
427  continue;
428 
429  linearizeElement_(elem);
430  }
431  }
432 
433  applyConstraintsToLinearization_();
434 
435  linearizeAuxiliaryEquations_();
436  }
437 
438  // linearize an element in the interior of the process' grid partition
439  void linearizeElement_(const Element& elem)
440  {
441  unsigned threadId = ThreadManager::threadId();
442 
443  ElementContext *elementCtx = elementCtx_[threadId];
444  auto& localLinearizer = model_().localLinearizer(threadId);
445 
446  // the actual work of linearization is done by the local linearizer class
447  localLinearizer.linearize(*elementCtx, elem);
448 
449  // update the right hand side and the Jacobian matrix
450  if (GET_PROP_VALUE(TypeTag, UseLinearizationLock))
451  globalMatrixMutex_.lock();
452 
453  size_t numPrimaryDof = elementCtx->numPrimaryDof(/*timeIdx=*/0);
454  for (unsigned primaryDofIdx = 0; primaryDofIdx < numPrimaryDof; ++ primaryDofIdx) {
455  unsigned globI = elementCtx->globalSpaceIndex(/*spaceIdx=*/primaryDofIdx, /*timeIdx=*/0);
456 
457  // update the right hand side
458  residual_[globI] += localLinearizer.residual(primaryDofIdx);
459 
460  // update the global Jacobian matrix
461  for (unsigned dofIdx = 0; dofIdx < elementCtx->numDof(/*timeIdx=*/0); ++ dofIdx) {
462  unsigned globJ = elementCtx->globalSpaceIndex(/*spaceIdx=*/dofIdx, /*timeIdx=*/0);
463 
464  (*matrix_)[globJ][globI] += localLinearizer.jacobian(dofIdx, primaryDofIdx);
465  }
466  }
467 
468  if (GET_PROP_VALUE(TypeTag, UseLinearizationLock))
469  globalMatrixMutex_.unlock();
470  }
471 
472  void linearizeAuxiliaryEquations_()
473  {
474  auto& model = model_();
475  for (unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx)
476  model.auxiliaryModule(auxModIdx)->linearize(*matrix_, residual_);
477  }
478 
479  // apply the constraints to the solution. (i.e., the solution of constraint degrees
480  // of freedom is set to the value of the constraint.)
481  void applyConstraintsToSolution_()
482  {
483  if (!enableConstraints_())
484  return;
485 
486  // TODO: assuming a history size of 2 only works for Euler time discretizations!
487  auto& sol = model_().solution(/*timeIdx=*/0);
488  auto& oldSol = model_().solution(/*timeIdx=*/1);
489 
490  auto it = constraintsMap_.begin();
491  const auto& endIt = constraintsMap_.end();
492  for (; it != endIt; ++it) {
493  sol[it->first] = it->second;
494  oldSol[it->first] = it->second;
495  }
496  }
497 
498  // apply the constraints to the linearization. (i.e., for constrain degrees of
499  // freedom the Jacobian matrix maps to identity and the residual is zero)
500  void applyConstraintsToLinearization_()
501  {
502  if (!enableConstraints_())
503  return;
504 
505  MatrixBlock idBlock = 0.0;
506  for (unsigned i = 0; i < numEq; ++i)
507  idBlock[i][i] = 1.0;
508 
509  auto it = constraintsMap_.begin();
510  const auto& endIt = constraintsMap_.end();
511  for (; it != endIt; ++it) {
512  unsigned constraintDofIdx = it->first;
513 
514  // reset the column of the Jacobian matrix
515  auto colIt = (*matrix_)[constraintDofIdx].begin();
516  const auto& colEndIt = (*matrix_)[constraintDofIdx].end();
517  for (; colIt != colEndIt; ++colIt)
518  *colIt = 0.0;
519 
520  // put an identity matrix on the main diagonal of the Jacobian
521  (*matrix_)[constraintDofIdx][constraintDofIdx] = idBlock;
522 
523  // make the right-hand side of constraint DOFs zero
524  residual_[constraintDofIdx] = 0.0;
525  }
526  }
527 
528  static bool enableConstraints_()
529  { return GET_PROP_VALUE(TypeTag, EnableConstraints); }
530 
531  Simulator *simulatorPtr_;
532  std::vector<ElementContext*> elementCtx_;
533 
534  // The constraint equations (only non-empty if the
535  // EnableConstraints property is true)
536  std::map<unsigned, Constraints> constraintsMap_;
537 
538  // the jacobian matrix
539  Matrix *matrix_;
540  // the right-hand side
541  GlobalEqVector residual_;
542 
543 
544  OmpMutex globalMatrixMutex_;
545 };
546 
547 } // namespace Ewoms
548 
549 #endif
The common code for the linearizers of non-linear systems of equations.
Definition: fvbaselinearizer.hh:65
Simplifies multi-threaded capabilities.
Definition: threadmanager.hh:52
static void registerParameters()
Register all run-time parameters for the Jacobian linearizer.
Definition: fvbaselinearizer.hh:129
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:189
Base class for specifying auxiliary equations.
Definition: locks.hh:54
const std::map< unsigned, Constraints > & constraintsMap() const
Returns the map of constraint degrees of freedom.
Definition: fvbaselinearizer.hh:237
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
const Matrix & matrix() const
Return constant reference to global Jacobian matrix.
Definition: fvbaselinearizer.hh:217
Simplifies multi-threaded capabilities.
Problem & problem()
Return the object which specifies the pysical setup of the simulation.
Definition: simulator.hh:202
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
static unsigned threadId()
Return the index of the current OpenMP thread.
Definition: threadmanager.hh:119
void eraseMatrix()
Causes the Jacobian matrix to be recreated from scratch before the next iteration.
Definition: fvbaselinearizer.hh:155
Declare the properties used by the infrastructure code of the finite volume discretizations.
Provides data handles for parallel communication which operate on DOFs.
const GlobalEqVector & residual() const
Return constant reference to global residual vector.
Definition: fvbaselinearizer.hh:226
static unsigned maxThreads()
Return the maximum number of threads of the current process.
Definition: threadmanager.hh:113
void init(Simulator &simulator)
Initialize the linearizer.
Definition: fvbaselinearizer.hh:141
Provides an STL-iterator like interface to iterate over the enties of a GridView in OpenMP threaded a...
void linearize()
Linearize the global non-linear system of equations.
Definition: fvbaselinearizer.hh:170
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
const GridView & gridView() const
Return the grid view for which the simulation is done.
Definition: simulator.hh:183