28#ifndef EWOMS_FV_BASE_LINEARIZER_HH
29#define EWOMS_FV_BASE_LINEARIZER_HH
34#include <opm/common/Exceptions.hpp>
35#include <opm/common/TimingMacros.hpp>
36#include <opm/grid/utility/SparseTable.hpp>
43#include <dune/common/version.hh>
44#include <dune/common/fvector.hh>
45#include <dune/common/fmatrix.hh>
57template<
class TypeTag>
58class EcfvDiscretization;
69template<
class TypeTag>
96 using Element =
typename GridView::template Codim<0>::Entity;
97 using ElementIterator =
typename GridView::template Codim<0>::Iterator;
99 using Vector = GlobalEqVector;
101 using IstlMatrix =
typename SparseMatrixAdapter::IstlMatrix;
106 using MatrixBlock =
typename SparseMatrixAdapter::MatrixBlock;
107 using VectorBlock = Dune::FieldVector<Scalar, numEq>;
124 auto it = elementCtx_.begin();
125 const auto&
endIt = elementCtx_.end();
126 for (; it !=
endIt; ++it)
147 simulatorPtr_ = &simulator;
149 auto it = elementCtx_.begin();
150 const auto&
endIt = elementCtx_.end();
151 for (; it !=
endIt; ++it){
154 elementCtx_.resize(0);
155 fullDomain_ = std::make_unique<FullDomain>(simulator.
gridView());
200 template <
class SubDomainType>
208 initFirstIteration_();
211 if (
static_cast<std::size_t
>(
domain.view.size(0)) == model_().numTotalDof()) {
223 catch (
const std::exception&
e)
225 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
226 <<
" caught an exception while linearizing:" <<
e.what()
227 <<
"\n" << std::flush;
232 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
233 <<
" caught an exception while linearizing"
234 <<
"\n" << std::flush;
240 throw NumericalProblem(
"A process did not succeed in linearizing the system");
244 { jacobian_->finalize(); }
256 auto& model = model_();
257 const auto& comm = simulator_().
gridView().comm();
261 model.auxiliaryModule(
auxModIdx)->linearize(*jacobian_, residual_);
263 catch (
const std::exception&
e) {
266 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
267 <<
" caught an exception while linearizing:" <<
e.what()
268 <<
"\n" << std::flush;
282 {
return *jacobian_; }
285 {
return *jacobian_; }
291 {
return residual_; }
294 {
return residual_; }
296 void setLinearizationType(LinearizationType linearizationType){
297 linearizationType_ = linearizationType;
300 const LinearizationType& getLinearizationType()
const{
301 return linearizationType_;
304 void updateDiscretizationParameters()
309 void updateBoundaryConditionData()
314 void updateFlowsInfo() {
324 {
return constraintsMap_; }
340 {
return floresInfo_;}
342 template <
class SubDomainType>
346 initFirstIteration_();
362 ElementContext&
elemCtx = *elementCtx_[threadId];
370 residual_[
globI] = 0.0;
371 jacobian_->clearRow(
globI, 0.0);
378 Simulator& simulator_()
379 {
return *simulatorPtr_; }
380 const Simulator& simulator_()
const
381 {
return *simulatorPtr_; }
384 {
return simulator_().
problem(); }
385 const Problem& problem_()
const
386 {
return simulator_().
problem(); }
389 {
return simulator_().
model(); }
390 const Model& model_()
const
391 {
return simulator_().
model(); }
393 const GridView& gridView_()
const
394 {
return problem_().gridView(); }
396 const ElementMapper& elementMapper_()
const
397 {
return model_().elementMapper(); }
399 const DofMapper& dofMapper_()
const
400 {
return model_().dofMapper(); }
402 void initFirstIteration_()
408 residual_.resize(model_().numTotalDof());
414 elementCtx_[threadId] =
new ElementContext(simulator_());
420 const auto& model = model_();
421 Stencil stencil(gridView_(), model_().dofMapper());
425 sparsityPattern_.clear();
426 sparsityPattern_.resize(model.numTotalDof());
429 stencil.update(
elem);
443 size_t numAuxMod = model.numAuxiliaryModules();
445 model.auxiliaryModule(
auxModIdx)->addNeighbors(sparsityPattern_);
448 jacobian_.reset(
new SparseMatrixAdapter(simulator_()));
451 jacobian_->reserve(sparsityPattern_);
464 void updateConstraintsMap_()
466 if (!enableConstraints_())
470 constraintsMap_.clear();
484 ElementContext&
elemCtx = *elementCtx_[threadId];
493 Constraints constraints;
494 elemCtx.problem().constraints(constraints,
498 if (constraints.isActive()) {
500 constraintsMap_[
globI] = constraints;
509 template <
class SubDomainType>
521 if (model_().newtonMethod().numIterations() == 0)
522 updateConstraintsMap_();
524 applyConstraintsToSolution_();
552 ||
nextElem.partitionType() == Dune::InteriorEntity)
563 linearizeElement_(
elem);
589 applyConstraintsToLinearization_();
594 template <
class ElementType>
599 ElementContext *
elementCtx = elementCtx_[threadId];
600 auto& localLinearizer = model_().localLinearizer(threadId);
607 globalMatrixMutex_.lock();
609 size_t numPrimaryDof =
elementCtx->numPrimaryDof(0);
625 globalMatrixMutex_.unlock();
630 void applyConstraintsToSolution_()
632 if (!enableConstraints_())
636 auto&
sol = model_().solution(0);
637 auto&
oldSol = model_().solution(1);
639 auto it = constraintsMap_.begin();
640 const auto&
endIt = constraintsMap_.end();
641 for (; it !=
endIt; ++it) {
642 sol[it->first] = it->second;
643 oldSol[it->first] = it->second;
649 void applyConstraintsToLinearization_()
651 if (!enableConstraints_())
654 auto it = constraintsMap_.begin();
655 const auto&
endIt = constraintsMap_.end();
656 for (; it !=
endIt; ++it) {
668 static bool enableConstraints_()
671 Simulator *simulatorPtr_;
672 std::vector<ElementContext*> elementCtx_;
676 std::map<unsigned, Constraints> constraintsMap_;
689 std::unique_ptr<SparseMatrixAdapter> jacobian_;
692 GlobalEqVector residual_;
694 LinearizationType linearizationType_;
696 std::mutex globalMatrixMutex_;
698 std::vector<std::set<unsigned int>> sparsityPattern_;
702 explicit FullDomain(
const GridView&
v) : view (
v) {}
704 std::vector<bool> interior;
709 std::unique_ptr<FullDomain> fullDomain_;
Base class for specifying auxiliary equations.
The common code for the linearizers of non-linear systems of equations.
Definition fvbaselinearizer.hh:71
const auto & getFlowsInfo() const
Return constant reference to the flowsInfo.
Definition fvbaselinearizer.hh:331
const GlobalEqVector & residual() const
Return constant reference to global residual vector.
Definition fvbaselinearizer.hh:290
static void registerParameters()
Register all run-time parameters for the Jacobian linearizer.
Definition fvbaselinearizer.hh:133
void linearize()
Linearize the full system of non-linear equations.
Definition fvbaselinearizer.hh:179
void init(Simulator &simulator)
Initialize the linearizer.
Definition fvbaselinearizer.hh:145
void linearizeDomain()
Linearize the part of the non-linear system of equations that is associated with the spatial domain.
Definition fvbaselinearizer.hh:195
const SparseMatrixAdapter & jacobian() const
Return constant reference to global Jacobian matrix backend.
Definition fvbaselinearizer.hh:281
const auto & getFloresInfo() const
Return constant reference to the floresInfo.
Definition fvbaselinearizer.hh:339
void linearizeAuxiliaryEquations()
Linearize the part of the non-linear system of equations that is associated with the spatial domain.
Definition fvbaselinearizer.hh:250
const std::map< unsigned, Constraints > & constraintsMap() const
Returns the map of constraint degrees of freedom.
Definition fvbaselinearizer.hh:323
void eraseMatrix()
Causes the Jacobian matrix to be recreated from scratch before the next iteration.
Definition fvbaselinearizer.hh:165
Definition matrixblock.hh:227
Manages the initializing and running of time dependent problems.
Definition simulator.hh:92
Problem & problem()
Return the object which specifies the pysical setup of the simulation.
Definition simulator.hh:291
const GridView & gridView() const
Return the grid view for which the simulation is done.
Definition simulator.hh:272
Model & model()
Return the physical model used in the simulation.
Definition simulator.hh:278
Simplifies multi-threaded capabilities.
Definition threadmanager.hpp:36
static unsigned maxThreads()
Return the maximum number of threads of the current process.
Definition threadmanager.hpp:66
static unsigned threadId()
Return the index of the current OpenMP thread.
Definition threadmanager.cpp:84
Declare the properties used by the infrastructure code of the finite volume discretizations.
Provides data handles for parallel communication which operate on DOFs.
The common code for the linearizers of non-linear systems of equations.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
ElementType
The types of reference elements available.
Definition vcfvstencil.hh:52
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
Provides an STL-iterator like interface to iterate over the enties of a GridView in OpenMP threaded a...
Simplifies multi-threaded capabilities.