28 #ifndef EWOMS_FV_BASE_LINEARIZER_HH 29 #define EWOMS_FV_BASE_LINEARIZER_HH 38 #include <opm/common/ErrorMacros.hpp> 39 #include <opm/common/Exceptions.hpp> 41 #include <dune/common/version.hh> 42 #include <dune/common/fvector.hh> 43 #include <dune/common/fmatrix.hh> 45 #include <type_traits> 52 template<
class TypeTag>
53 class EcfvDiscretization;
64 template<
class TypeTag>
69 typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
74 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
76 typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
77 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
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;
83 typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints;
87 typedef typename GET_PROP_TYPE(TypeTag, GridCommHandleFactory) GridCommHandleFactory;
89 typedef Opm::MathToolbox<Evaluation> Toolbox;
91 typedef typename GridView::template Codim<0>::Entity Element;
92 typedef typename GridView::template Codim<0>::Iterator ElementIterator;
94 typedef GlobalEqVector Vector;
95 typedef JacobianMatrix Matrix;
98 enum { historySize =
GET_PROP_VALUE(TypeTag, TimeDiscHistorySize) };
100 typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
101 typedef Dune::FieldVector<Scalar, numEq> VectorBlock;
103 static const bool linearizeNonLocalElements =
GET_PROP_VALUE(TypeTag, LinearizeNonLocalElements);
120 auto it = elementCtx_.begin();
121 const auto& endIt = elementCtx_.end();
122 for (; it != endIt; ++it)
143 simulatorPtr_ = &simulator;
176 initFirstIteration_();
183 catch (
const std::exception& e)
185 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
186 <<
" caught an exception while linearizing:" << e.what()
187 <<
"\n" << std::flush;
190 #if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,5) 191 catch (
const Dune::Exception& e)
193 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
194 <<
" caught an exception while linearizing:" << e.what()
195 <<
"\n" << std::flush;
201 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
202 <<
" caught an exception while linearizing" 203 <<
"\n" << std::flush;
206 succeeded = gridView_().comm().min(succeeded);
209 OPM_THROW(Opm::NumericalProblem,
210 "A process did not succeed in linearizing the system");
227 {
return residual_; }
230 {
return residual_; }
238 {
return constraintsMap_; }
242 {
return *simulatorPtr_; }
243 const Simulator& simulator_()
const 244 {
return *simulatorPtr_; }
247 {
return simulator_().
problem(); }
248 const Problem& problem_()
const 249 {
return simulator_().
problem(); }
252 {
return simulator_().
model(); }
253 const Model& model_()
const 254 {
return simulator_().
model(); }
256 const GridView& gridView_()
const 257 {
return problem_().gridView(); }
259 const ElementMapper& elementMapper_()
const 260 {
return model_().elementMapper(); }
262 const DofMapper& dofMapper_()
const 263 {
return model_().dofMapper(); }
265 void initFirstIteration_()
272 residual_.resize(model_().numTotalDof());
278 elementCtx_[threadId] =
new ElementContext(simulator_());
284 size_t numAllDof = model_().numTotalDof();
287 matrix_ =
new Matrix(numAllDof, numAllDof, Matrix::random);
289 Stencil stencil(gridView_(), model_().dofMapper() );
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);
301 for (
unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
302 unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
304 for (
unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
305 unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
306 neighbors[myIdx].insert(neighborIdx);
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);
319 for (
unsigned dofIdx = 0; dofIdx < numAllDof; ++ dofIdx)
320 matrix_->setrowsize(dofIdx, neighbors[dofIdx].size());
321 matrix_->endrowsizes();
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);
332 matrix_->endindices();
344 void updateConstraintsMap_()
346 if (!enableConstraints_())
351 constraintsMap_.clear();
354 ThreadedEntityIterator<GridView, 0> threadedElemIt(gridView_());
359 ElementIterator elemIt = threadedElemIt.beginParallel();
360 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
363 const Element& elem = *elemIt;
364 ElementContext& elemCtx = *elementCtx_[threadId];
365 elemCtx.updateStencil(elem);
369 for (
unsigned primaryDofIdx = 0;
370 primaryDofIdx < elemCtx.numPrimaryDof(0);
373 Constraints constraints;
374 elemCtx.problem().constraints(constraints,
378 if (constraints.isActive()) {
379 unsigned globI = elemCtx.globalSpaceIndex(primaryDofIdx, 0);
380 constraintsMap_[globI] = constraints;
396 if (model_().newtonMethod().numIterations() == 0)
397 updateConstraintsMap_();
399 applyConstraintsToSolution_();
404 ThreadedEntityIterator<GridView, 0> threadedElemIt(gridView_());
409 ElementIterator elemIt = threadedElemIt.beginParallel();
410 ElementIterator nextElemIt = elemIt;
411 for (; !threadedElemIt.isFinished(elemIt); elemIt = nextElemIt) {
414 nextElemIt = threadedElemIt.increment();
415 if (!threadedElemIt.isFinished(nextElemIt)) {
416 const auto& nextElem = *nextElemIt;
417 if (linearizeNonLocalElements
418 || nextElem.partitionType() == Dune::InteriorEntity)
420 model_().prefetch(nextElem);
421 problem_().prefetch(nextElem);
425 const Element& elem = *elemIt;
426 if (!linearizeNonLocalElements && elem.partitionType() != Dune::InteriorEntity)
429 linearizeElement_(elem);
433 applyConstraintsToLinearization_();
435 linearizeAuxiliaryEquations_();
439 void linearizeElement_(
const Element& elem)
443 ElementContext *elementCtx = elementCtx_[threadId];
444 auto& localLinearizer = model_().localLinearizer(threadId);
447 localLinearizer.linearize(*elementCtx, elem);
451 globalMatrixMutex_.lock();
453 size_t numPrimaryDof = elementCtx->numPrimaryDof(0);
454 for (
unsigned primaryDofIdx = 0; primaryDofIdx < numPrimaryDof; ++ primaryDofIdx) {
455 unsigned globI = elementCtx->globalSpaceIndex(primaryDofIdx, 0);
458 residual_[globI] += localLinearizer.residual(primaryDofIdx);
461 for (
unsigned dofIdx = 0; dofIdx < elementCtx->numDof(0); ++ dofIdx) {
462 unsigned globJ = elementCtx->globalSpaceIndex(dofIdx, 0);
464 (*matrix_)[globJ][globI] += localLinearizer.jacobian(dofIdx, primaryDofIdx);
469 globalMatrixMutex_.unlock();
472 void linearizeAuxiliaryEquations_()
474 auto& model = model_();
475 for (
unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx)
476 model.auxiliaryModule(auxModIdx)->linearize(*matrix_, residual_);
481 void applyConstraintsToSolution_()
483 if (!enableConstraints_())
487 auto& sol = model_().solution(0);
488 auto& oldSol = model_().solution(1);
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;
500 void applyConstraintsToLinearization_()
502 if (!enableConstraints_())
505 MatrixBlock idBlock = 0.0;
506 for (
unsigned i = 0; i < numEq; ++i)
509 auto it = constraintsMap_.begin();
510 const auto& endIt = constraintsMap_.end();
511 for (; it != endIt; ++it) {
512 unsigned constraintDofIdx = it->first;
515 auto colIt = (*matrix_)[constraintDofIdx].begin();
516 const auto& colEndIt = (*matrix_)[constraintDofIdx].end();
517 for (; colIt != colEndIt; ++colIt)
521 (*matrix_)[constraintDofIdx][constraintDofIdx] = idBlock;
524 residual_[constraintDofIdx] = 0.0;
528 static bool enableConstraints_()
531 Simulator *simulatorPtr_;
532 std::vector<ElementContext*> elementCtx_;
536 std::map<unsigned, Constraints> constraintsMap_;
541 GlobalEqVector residual_;
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
Definition: baseauxiliarymodule.hh:37
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:189
Base class for specifying auxiliary equations.
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
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
const std::map< unsigned, Constraints > & constraintsMap() const
Returns the map of constraint degrees of freedom.
Definition: fvbaselinearizer.hh:237
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 GridView & gridView() const
Return the grid view for which the simulation is done.
Definition: simulator.hh:183
static unsigned maxThreads()
Return the maximum number of threads of the current process.
Definition: threadmanager.hh:113
const GlobalEqVector & residual() const
Return constant reference to global residual vector.
Definition: fvbaselinearizer.hh:226
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 Matrix & matrix() const
Return constant reference to global Jacobian matrix.
Definition: fvbaselinearizer.hh:217