27 #ifndef EWOMS_FV_BASE_LINEARIZER_HH
28 #define EWOMS_FV_BASE_LINEARIZER_HH
37 #include <dune/common/fvector.hh>
38 #include <dune/common/fmatrix.hh>
40 #include <type_traits>
47 template<
class TypeTag>
48 class EcfvDiscretization;
59 template<
class TypeTag>
64 typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
69 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
71 typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
72 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
74 typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
75 typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
76 typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
81 typedef typename GET_PROP_TYPE(TypeTag, GridCommHandleFactory) GridCommHandleFactory;
83 typedef Opm::MathToolbox<Evaluation> Toolbox;
85 typedef typename GridView::template Codim<0>::Entity Element;
86 typedef typename GridView::template Codim<0>::Iterator ElementIterator;
88 typedef GlobalEqVector Vector;
89 typedef JacobianMatrix Matrix;
92 enum { historySize =
GET_PROP_VALUE(TypeTag, TimeDiscHistorySize) };
94 typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
95 typedef Dune::FieldVector<Scalar, numEq> VectorBlock;
97 static const bool linearizeNonLocalElements =
GET_PROP_VALUE(TypeTag, LinearizeNonLocalElements);
136 relinearizationAccuracy_ = 0.0;
142 auto it = elementCtx_.begin();
143 const auto &endIt = elementCtx_.end();
144 for (; it != endIt; ++it)
154 "Re-use of the linearized system of equations at the first "
155 "iteration of the next time step");
157 "relinearize only those degrees of freedom that have changed "
158 "'sufficiently' between two Newton iterations");
170 void init(Simulator& simulator)
172 simulatorPtr_ = &simulator;
201 initFirstIteration_();
206 bool linearizationReused = reuseLinearization_;
210 int curNumRelin = numTotalElems_ - numGreenElems_;
211 Scalar curRelAcc = relinearizationAccuracy_;
217 succeeded = gridView_().comm().min(succeeded);
219 catch (Opm::NumericalProblem &e)
221 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
222 <<
" caught an exception while linearizing:" << e.what()
223 <<
"\n" << std::flush;
225 succeeded = gridView_().comm().min(succeeded);
229 OPM_THROW(Opm::NumericalProblem,
230 "A process did not succeed in linearizing the system");
233 if (!linearizationReused && enablePartialRelinearization_()) {
234 model_().newtonMethod().endIterMsg()
235 <<
", relinearized " << curNumRelin <<
" of " << numTotalElems_
236 <<
" elements (" << 100*Scalar(curNumRelin)/numTotalElems_ <<
"%)"
237 <<
" and achieved a accuracy of " << curRelAcc;
250 if (enableLinearizationRecycling_())
251 reuseLinearization_ = yesno;
261 reuseLinearization_ =
false;
264 relinearizationAccuracy_ = 0.0;
265 relinearizationTolerance_ = 0.0;
267 if (enablePartialRelinearization_()) {
268 std::fill(dofError_.begin(), dofError_.end(), 1e100);
269 std::fill(dofColor_.begin(), dofColor_.end(),
Red);
270 std::fill(elementColor_.begin(), elementColor_.end(),
Red);
274 auto &model = model_();
275 int numGridDof = model.numGridDof();
276 for (
int timeIdx = 0; timeIdx < historySize; ++timeIdx) {
277 for (
int dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
278 model.setIntensiveQuantitiesCacheEntryValidity(dofIdx, timeIdx,
false);
294 {
return relinearizationAccuracy_; }
300 {
return maxDofError_; }
314 if (!enablePartialRelinearization_())
317 unsigned numGridDof = model_().numGridDof();
322 for (
unsigned globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
323 if (!model_().isLocalDof(globalDofIdx)) {
325 dofError_[globalDofIdx] = 0;
332 const auto &d = uDelta[globalDofIdx];
334 for (
unsigned pvIdx = 0; pvIdx < d.size(); ++pvIdx) {
335 Scalar tmp = std::abs(d[pvIdx]*model_().primaryVarWeight(globalDofIdx, pvIdx));
336 distRel = std::max(distRel, tmp);
338 Valgrind::CheckDefined(distRel);
339 dofError_[globalDofIdx] += distRel;
340 maxDofError_ = std::max(maxDofError_, dofError_[globalDofIdx]);
352 if (enablePartialRelinearization_())
353 dofError_[dofIdx] = 1e100;
354 this->model_().setIntensiveQuantitiesCacheEntryValidity(dofIdx, 0,
false);
362 if (!enablePartialRelinearization_())
364 return dofColor_[dofIdx];
372 if (!enablePartialRelinearization_())
374 return elementColor_[elemIdx];
381 {
return relinearizationTolerance_; }
387 { relinearizationTolerance_ = tolerance; }
393 {
return dofError_[dofIdx]; }
405 {
return residual_; }
413 void computeColors_()
415 if (!enablePartialRelinearization_())
420 for (
unsigned dofIdx = 0; dofIdx < dofColor_.size(); ++dofIdx) {
421 if (dofError_[dofIdx] > relinearizationTolerance_) {
424 dofColor_[dofIdx] =
Red;
425 dofError_[dofIdx] = 0.0;
427 else if (!model_().isLocalDof(dofIdx)) {
429 dofColor_[dofIdx] =
Red;
430 dofError_[dofIdx] = 0.0;
433 dofColor_[dofIdx] =
Green;
437 int numGridDof = this->model_().numGridDof();
438 for (
int dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
439 if (dofColor_[dofIdx] !=
Red)
442 typedef typename JacobianMatrix::ColIterator ColIterator;
443 ColIterator colIt = (*matrix_)[dofIdx].begin();
444 const ColIterator &colEndIt = (*matrix_)[dofIdx].end();
445 for (; colIt != colEndIt; ++colIt) {
446 int dof2Idx = colIt.index();
447 if (dof2Idx >= numGridDof)
450 if (dofColor_[dof2Idx] !=
Red)
451 dofColor_[dof2Idx] =
Yellow;
459 if (std::is_same<Discretization, EcfvDiscretization<TypeTag> >::value)
460 elementColor_ = dofColor_;
462 const auto& gridView = this->model_().gridView();
463 const auto& elemMapper = this->model_().elementMapper();
464 Stencil stencil(gridView);
465 ElementIterator elemIt = gridView.template begin<0>();
466 const ElementIterator& elemEndIt = gridView.template end<0>();
467 for (; elemIt != elemEndIt; ++elemIt) {
468 const Element& elem = *elemIt;
469 stencil.updateTopology(elem);
471 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
472 int elemIdx = elemMapper.index(elem);
474 int elemIdx = elemMapper.map(elem);
476 int numElemDof = stencil.numDof();
480 elementColor_[elemIdx] =
Green;
481 for (
int elemDofIdx = 0; elemDofIdx < numElemDof; ++elemDofIdx) {
482 int dofIdx = stencil.globalSpaceIndex(elemDofIdx);
483 if (dofColor_[dofIdx] !=
Green) {
484 elementColor_[elemIdx] =
Red;
492 relinearizationAccuracy_ = 0;
494 for (
unsigned dofIdx = 0; dofIdx < dofColor_.size(); ++dofIdx) {
495 if (dofColor_[dofIdx] !=
Red) {
496 relinearizationAccuracy_ =
497 std::max(relinearizationAccuracy_, dofError_[dofIdx]);
499 if (dofColor_[dofIdx] ==
Green)
505 relinearizationAccuracy_ = gridView_().comm().max(relinearizationAccuracy_);
506 numGreenElems_ = gridView_().comm().sum(numGreenElems_);
509 static bool enableLinearizationRecycling_()
510 {
return EWOMS_GET_PARAM(TypeTag,
bool, EnableLinearizationRecycling); }
511 static bool enablePartialRelinearization_()
512 {
return EWOMS_GET_PARAM(TypeTag,
bool, EnablePartialRelinearization); }
514 Simulator &simulator_()
515 {
return *simulatorPtr_; }
516 const Simulator &simulator_()
const
517 {
return *simulatorPtr_; }
520 {
return simulator_().
problem(); }
521 const Problem &problem_()
const
522 {
return simulator_().
problem(); }
525 {
return simulator_().
model(); }
526 const Model &model_()
const
527 {
return simulator_().
model(); }
529 const GridView &gridView_()
const
530 {
return problem_().gridView(); }
532 const ElementMapper &elementMapper_()
const
533 {
return model_().elementMapper(); }
535 const DofMapper &dofMapper_()
const
536 {
return model_().dofMapper(); }
538 void initFirstIteration_()
542 int numGridDof = model_().numGridDof();
543 int numAllDof = model_().numTotalDof();
544 int numElems = gridView_().size(0);
545 numTotalElems_ = gridView_().comm().sum(numElems);
553 residual_.resize(numAllDof);
559 elementCtx_[threadId] =
new ElementContext(simulator_());
563 if (enableLinearizationRecycling_()) {
564 storageJacobian_.resize(numGridDof);
565 storageTerm_.resize(numGridDof);
569 reuseLinearization_ =
false;
570 if (enablePartialRelinearization_()) {
571 dofColor_.resize(numGridDof);
572 dofError_.resize(numGridDof);
573 elementColor_.resize(numElems);
581 int numAllDof = model_().numTotalDof();
584 matrix_ =
new Matrix(numAllDof, numAllDof, Matrix::random);
586 Stencil stencil(gridView_());
590 typedef std::set<int> NeighborSet;
591 std::vector<NeighborSet> neighbors(numAllDof);
592 ElementIterator elemIt = gridView_().template begin<0>();
593 const ElementIterator elemEndIt = gridView_().template end<0>();
594 for (; elemIt != elemEndIt; ++elemIt) {
595 const Element &elem = *elemIt;
596 stencil.update(elem);
598 for (
int primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
599 int myIdx = stencil.globalSpaceIndex(primaryDofIdx);
601 for (
int dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
602 int neighborIdx = stencil.globalSpaceIndex(dofIdx);
603 neighbors[myIdx].insert(neighborIdx);
610 const auto& model = model_();
611 int numAuxMod = model.numAuxiliaryModules();
612 for (
int auxModIdx = 0; auxModIdx < numAuxMod; ++auxModIdx)
613 model.auxiliaryModule(auxModIdx)->addNeighbors(neighbors);
616 for (
int dofIdx = 0; dofIdx < numAllDof; ++ dofIdx)
617 matrix_->setrowsize(dofIdx, neighbors[dofIdx].size());
618 matrix_->endrowsizes();
623 for (
int dofIdx = 0; dofIdx < numAllDof; ++ dofIdx) {
624 typename NeighborSet::iterator nIt = neighbors[dofIdx].begin();
625 typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end();
626 for (; nIt != nEndIt; ++nIt)
627 matrix_->addindex(dofIdx, *nIt);
629 matrix_->endindices();
637 size_t numGridDof = model_().numGridDof();
638 size_t numTotalDof = model_().numTotalDof();
640 if (!enablePartialRelinearization_()) {
646 if (enableLinearizationRecycling_()) {
647 for (
unsigned dofIdx=0; dofIdx < numGridDof; ++ dofIdx) {
648 storageTerm_[dofIdx] = 0.0;
649 storageJacobian_[dofIdx] = 0.0;
658 if (enableLinearizationRecycling_())
659 for (
unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx)
660 storageTerm_[dofIdx] = 0.0;
663 for (
unsigned dofIdx = numGridDof; dofIdx < numTotalDof; ++dofIdx) {
665 typedef typename JacobianMatrix::ColIterator ColIterator;
666 ColIterator colIt = (*matrix_)[dofIdx].begin();
667 const ColIterator &colEndIt = (*matrix_)[dofIdx].end();
668 for (; colIt != colEndIt; ++colIt) {
674 for (
unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
684 typedef typename JacobianMatrix::ColIterator ColIterator;
685 ColIterator colIt = (*matrix_)[dofIdx].beforeEnd();
686 const ColIterator colBeginIt = (*matrix_)[dofIdx].beforeBegin();
687 for (; colIt != colBeginIt; -- colIt) {
688 if (colIt.index() < numGridDof)
699 if (enableLinearizationRecycling_())
700 storageJacobian_[dofIdx] = 0.0;
704 typedef typename JacobianMatrix::ColIterator ColIterator;
705 ColIterator colIt = (*matrix_)[dofIdx].begin();
706 const ColIterator &colEndIt = (*matrix_)[dofIdx].end();
707 for (; colIt != colEndIt; ++colIt) {
708 if (colIt.index() >= numGridDof || dofColor_[colIt.index()] !=
Green)
722 Scalar curDt = problem_().simulator().timeStepSize();
723 if (reuseLinearization_) {
724 int numGridDof = model_().numGridDof();
725 for (
int dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
730 residual_[dofIdx] -= storageTerm_[dofIdx];
733 MatrixBlock &J_ii = (*matrix_)[dofIdx][dofIdx];
735 J_ii -= storageJacobian_[dofIdx];
736 storageJacobian_[dofIdx] *= oldDt_/curDt;
737 J_ii += storageJacobian_[dofIdx];
740 reuseLinearization_ =
false;
743 linearizeAuxiliaryEquations_();
745 problem_().newtonMethod().endIterMsg()
746 <<
", linear system of equations rescaled using previous time step";
755 ThreadedEntityIterator<GridView, 0> threadedElemIt(gridView_());
760 ElementIterator elemIt = gridView_().template begin<0>();
761 for (threadedElemIt.beginParallel(elemIt);
762 !threadedElemIt.isFinished(elemIt);
763 threadedElemIt.increment(elemIt))
765 const Element &elem = *elemIt;
767 if (!linearizeNonLocalElements && elem.partitionType() != Dune::InteriorEntity)
770 linearizeElement_(elem);
774 linearizeAuxiliaryEquations_();
778 void linearizeElement_(
const Element &elem)
780 if (enablePartialRelinearization_()) {
781 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
782 int globalElemIdx = model_().elementMapper().index(elem);
784 int globalElemIdx = model_().elementMapper().map(elem);
787 linearizeGreenElement_(elem);
794 ElementContext *elementCtx = elementCtx_[threadId];
795 auto &localLinearizer = model_().localLinearizer(threadId);
798 elementCtx->updateAll(elem);
799 localLinearizer.linearize(*elementCtx);
803 int numPrimaryDof = elementCtx->numPrimaryDof(0);
804 for (
int primaryDofIdx = 0; primaryDofIdx < numPrimaryDof; ++ primaryDofIdx) {
805 int globI = elementCtx->globalSpaceIndex(primaryDofIdx, 0);
814 residual_[globI] += localLinearizer.residual(primaryDofIdx);
816 if (enableLinearizationRecycling_()) {
817 storageTerm_[globI] += localLinearizer.residualStorage(primaryDofIdx);
818 storageJacobian_[globI] += localLinearizer.jacobianStorage(primaryDofIdx);
822 for (
int dofIdx = 0; dofIdx < elementCtx->numDof(0); ++ dofIdx) {
823 int globJ = elementCtx->globalSpaceIndex(dofIdx, 0);
830 (*matrix_)[globJ][globI] += localLinearizer.jacobian(dofIdx, primaryDofIdx);
838 void linearizeGreenElement_(
const Element &elem)
841 ElementContext *elementCtx = elementCtx_[threadId];
843 elementCtx->updateAll(elem);
844 auto& localResidual = model_().localResidual(threadId);
845 localResidual.eval(*elementCtx);
848 for (
int dofIdx=0; dofIdx < elementCtx->numPrimaryDof(0); ++ dofIdx) {
849 int globI = elementCtx->globalSpaceIndex(dofIdx, 0);
852 for (
int eqIdx = 0; eqIdx < numEq; ++ eqIdx)
853 residual_[globI][eqIdx] += Toolbox::value(localResidual.residual(dofIdx)[eqIdx]);
854 if (enableLinearizationRecycling_())
855 for (
int eqIdx = 0; eqIdx < numEq; ++ eqIdx)
856 storageTerm_[globI][eqIdx] += Toolbox::value(localResidual.storageTerm(dofIdx)[eqIdx]);
861 void linearizeAuxiliaryEquations_()
863 auto& model = model_();
864 for (
unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx)
865 model.auxiliaryModule(auxModIdx)->linearize(*matrix_, residual_);
868 Simulator *simulatorPtr_;
869 std::vector<ElementContext*> elementCtx_;
874 GlobalEqVector residual_;
877 bool reuseLinearization_;
879 std::vector<MatrixBlock> storageJacobian_;
880 std::vector<VectorBlock> storageTerm_;
885 std::vector<EntityColor> dofColor_;
886 std::vector<Scalar> dofError_;
887 std::vector<EntityColor> elementColor_;
892 Scalar relinearizationTolerance_;
893 Scalar relinearizationAccuracy_;
FvBaseLinearizer()
Definition: fvbaselinearizer.hh:127
void relinearizeAll()
If partial relinearization is enabled, this method causes all elements to be relinearized in the next...
Definition: fvbaselinearizer.hh:258
Scalar relinearizationTolerance() const
Returns the maximum error for which a degree of freedom is not relinearized.
Definition: fvbaselinearizer.hh:380
Problem & problem()
Return the object which specifies the pysical setup of the simulation.
Definition: simulator.hh:189
EntityColor dofColor(int dofIdx) const
Returns the "relinearization color" of a degree of freedom.
Definition: fvbaselinearizer.hh:360
Declare the properties used by the infrastructure code of the finite volume discretizations.
void linearize()
Linearize the global non-linear system of equations.
Definition: fvbaselinearizer.hh:195
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:468
The common code for the linearizers of non-linear systems of equations.
Definition: fvbaselinearizer.hh:60
void init(Simulator &simulator)
Initialize the linearizer.
Definition: fvbaselinearizer.hh:170
const GlobalEqVector & residual() const
Return constant reference to global residual vector.
Definition: fvbaselinearizer.hh:404
void setLinearizationReusable(bool yesno=true)
If linearization recycling is enabled, this method specifies whether the next call to linearize() jus...
Definition: fvbaselinearizer.hh:248
Definition: fvbaselinearizer.hh:124
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:485
This class implements an exception-safe scoped lock-keeper.
Definition: locks.hh:63
const Matrix & matrix() const
Return constant reference to global Jacobian matrix.
Definition: fvbaselinearizer.hh:398
Simplifies multi-threaded capabilities.
Definition: threadmanager.hh:48
void updateRelinearizationErrors(const GlobalEqVector &uDelta, const GlobalEqVector &resid)
Update the distance where the non-linear system was originally insistently linearized and the point w...
Definition: fvbaselinearizer.hh:312
void recreateMatrix()
Causes the Jacobian matrix to be recreated in the next iteration.
Definition: fvbaselinearizer.hh:180
Definition: fvbaselinearizer.hh:119
void setRelinearizationTolerance(Scalar tolerance)
Sets the maximum error for which a degree of freedom is not relinearized.
Definition: fvbaselinearizer.hh:386
void markDofRed(int dofIdx)
Ensure that a given degree of freedom is relinarized in the next iteration.
Definition: fvbaselinearizer.hh:350
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:73
Scalar dofError(int dofIdx) const
Returns the error for a given degree of freedom after the last iteration.
Definition: fvbaselinearizer.hh:392
Definition: baseauxiliarymodule.hh:35
Scalar relinearizationAccuracy() const
Returns the largest error of a "green" degree of freedom for the most recent call of the linearize() ...
Definition: fvbaselinearizer.hh:293
EntityColor elementColor(int elemIdx) const
Returns the "relinearization color" of an element.
Definition: fvbaselinearizer.hh:370
Scalar maxDofError() const
The maximum deflection seen for any DOF after an Newton update.
Definition: fvbaselinearizer.hh:299
Provides data handles for parallel communication which operate on DOFs.
static int maxThreads()
Return the maximum number of threads of the current process.
Definition: threadmanager.hh:117
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:64
static int threadId()
Return the index of the current OpenMP thread.
Definition: threadmanager.hh:123
EntityColor
The colors of elements and degrees of freedom required for partial relinearization.
Definition: fvbaselinearizer.hh:108
static void registerParameters()
Register all run-time parameters for the Jacobian linearizer.
Definition: fvbaselinearizer.hh:151
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:176
~FvBaseLinearizer()
Definition: fvbaselinearizer.hh:139
Implements a shallow wrapper around the "raw" locks provided by OpenMP.
Definition: locks.hh:35
GridView & gridView()
Return the grid view for which the simulation is done.
Definition: simulator.hh:164
Definition: fvbaselinearizer.hh:113
Base class for specifying auxiliary equations.
Simplifies multi-threaded capabilities.
Provides an STL-iterator like interface to iterate over the enties of a GridView in OpenMP threaded a...
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:95