20 #ifndef OPM_ISTLSOLVEREBOSFLEXIBLE_HEADER_INCLUDED
21 #define OPM_ISTLSOLVEREBOSFLEXIBLE_HEADER_INCLUDED
23 #include <opm/simulators/linalg/matrixblock.hh>
24 #include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
25 #include <opm/simulators/linalg/FlexibleSolver.hpp>
26 #include <opm/simulators/linalg/setupPropertyTree.hpp>
27 #include <opm/simulators/linalg/WriteSystemMatrixHelper.hpp>
29 #include <opm/common/ErrorMacros.hpp>
34 namespace Opm::Properties {
38 using InheritsFrom = std::tuple<FlowIstlSolverParams>;
57 template <
class TypeTag>
60 using GridView = GetPropType<TypeTag, Properties::GridView>;
61 using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
62 using VectorType = GetPropType<TypeTag, Properties::GlobalEqVector>;
63 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
64 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
65 using MatrixType =
typename SparseMatrixAdapter::IstlMatrix;
66 using WellModel = GetPropType<TypeTag, Properties::EclWellModel>;
68 using Communication = Dune::OwnerOverlapCopyCommunication<int, int>;
70 using Communication = int;
72 using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
77 using Vector = GetPropType<TypeTag, Properties::GlobalEqVector>;
78 using Indices = GetPropType<TypeTag, Properties::Indices>;
79 typedef typename SparseMatrixAdapter::IstlMatrix Matrix;
80 typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType;
81 typedef typename Vector::block_type BlockVector;
82 using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
83 using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
84 typedef typename GridView::template Codim<0>::Entity Element;
85 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
86 constexpr
static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
89 static void registerParameters()
91 FlowLinearSolverParameters::registerParameters<TypeTag>();
95 : simulator_(simulator)
96 , ownersFirst_(EWOMS_GET_PARAM(TypeTag,
bool, OwnerCellsFirst))
97 , matrixAddWellContributions_(EWOMS_GET_PARAM(TypeTag,
bool, MatrixAddWellContributions))
98 , interiorCellNum_(detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), ownersFirst_))
100 parameters_.template init<TypeTag>();
102 EWOMS_PARAM_IS_SET(TypeTag,
int, LinearSolverMaxIter),
103 EWOMS_PARAM_IS_SET(TypeTag,
int, CprMaxEllIter));
109 using ElementMapper =
110 Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
111 ElementMapper elemMapper(simulator_.gridView(), Dune::mcmgElementLayout());
112 detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
114 if (parallelInformation_.type() ==
typeid(ParallelISTLInformation)) {
116 const ParallelISTLInformation* parinfo = std::any_cast<ParallelISTLInformation>(¶llelInformation_);
118 comm_.reset(
new Communication(parinfo->communicator()));
122 if (simulator.gridView().comm().rank() == 0) {
123 std::ostringstream os;
124 os <<
"Property tree for linear solver:\n";
125 prm_.write_json(os,
true);
126 OpmLog::note(os.str());
134 void prepare(SparseMatrixAdapter& mat, VectorType& b)
137 static bool firstcall =
true;
138 if (firstcall && parallelInformation_.type() ==
typeid(ParallelISTLInformation)) {
140 const ParallelISTLInformation* parinfo = std::any_cast<ParallelISTLInformation>(¶llelInformation_);
142 const size_t size = mat.istlMatrix().N();
143 parinfo->copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), size, 1);
146 if (isParallel() && matrixAddWellContributions_) {
151 matrix_ = &mat.istlMatrix();
152 std::function<VectorType()> weightsCalculator = getWeightsCalculator(mat.istlMatrix(), b);
154 if (shouldCreateSolver()) {
157 if (matrixAddWellContributions_) {
158 using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Communication>;
159 auto op = std::make_unique<ParOperatorType>(mat.istlMatrix(), *comm_);
160 auto sol = std::make_unique<SolverType>(*op, *comm_, prm_, weightsCalculator,
162 solver_ = std::move(sol);
163 linear_operator_ = std::move(op);
166 OPM_THROW(std::runtime_error,
"In parallel, the flexible solver requires "
167 "--owner-cells-first=true when --matrix-add-well-contributions=false is used.");
170 auto well_op = std::make_unique<WellModelOpType>(simulator_.problem().wellModel());
171 auto op = std::make_unique<ParOperatorType>(mat.istlMatrix(), *well_op, interiorCellNum_);
172 auto sol = std::make_unique<SolverType>(*op, *comm_, prm_, weightsCalculator,
174 solver_ = std::move(sol);
175 linear_operator_ = std::move(op);
176 well_operator_ = std::move(well_op);
180 if (matrixAddWellContributions_) {
181 using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
182 auto op = std::make_unique<SeqOperatorType>(mat.istlMatrix());
183 auto sol = std::make_unique<SolverType>(*op, prm_, weightsCalculator,
185 solver_ = std::move(sol);
186 linear_operator_ = std::move(op);
189 auto well_op = std::make_unique<WellModelOpType>(simulator_.problem().wellModel());
190 auto op = std::make_unique<SeqOperatorType>(mat.istlMatrix(), *well_op);
191 auto sol = std::make_unique<SolverType>(*op, prm_, weightsCalculator,
193 solver_ = std::move(sol);
194 linear_operator_ = std::move(op);
195 well_operator_ = std::move(well_op);
200 solver_->preconditioner().update();
205 bool solve(VectorType& x)
207 solver_->apply(x, rhs_, res_);
209 return res_.converged;
212 bool isParallel()
const
215 return parallelInformation_.type() ==
typeid(ParallelISTLInformation);
221 int iterations()
const
223 return res_.iterations;
226 void setResidual(VectorType& )
231 void setMatrix(
const SparseMatrixAdapter& )
238 bool shouldCreateSolver()
const
245 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
246 bool recreate_solver =
false;
247 if (this->parameters_.cpr_reuse_setup_ == 0) {
249 recreate_solver =
true;
250 }
else if (this->parameters_.cpr_reuse_setup_ == 1) {
252 if (newton_iteration == 0) {
253 recreate_solver =
true;
255 }
else if (this->parameters_.cpr_reuse_setup_ == 2) {
257 if (this->iterations() > 10) {
258 recreate_solver =
true;
261 assert(this->parameters_.cpr_reuse_setup_ == 3);
262 assert(recreate_solver ==
false);
265 return recreate_solver;
268 std::function<VectorType()> getWeightsCalculator(
const MatrixType& mat,
const VectorType& b)
const
270 std::function<VectorType()> weightsCalculator;
272 using namespace std::string_literals;
274 auto preconditionerType = prm_.get(
"preconditioner.type",
"cpr"s);
275 if (preconditionerType ==
"cpr" || preconditionerType ==
"cprt") {
276 const bool transpose = preconditionerType ==
"cprt";
277 const auto weightsType = prm_.get(
"preconditioner.weight_type",
"quasiimpes"s);
278 if (weightsType ==
"quasiimpes") {
280 weightsCalculator = [&mat, transpose, p = pressureIndex]() {
281 return Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(mat, p, transpose);
283 }
else if (weightsType ==
"trueimpes") {
284 weightsCalculator = [
this, &b, p = pressureIndex]() {
285 return this->getTrueImpesWeights(b, p);
288 OPM_THROW(std::invalid_argument,
289 "Weights type " << weightsType <<
"not implemented for cpr."
290 <<
" Please use quasiimpes or trueimpes.");
293 return weightsCalculator;
301 const int numEq = MatrixType::block_type::rows;
302 typename MatrixType::block_type diag_block(0.0);
303 for (
int eq = 0; eq < numEq; ++eq)
304 diag_block[eq][eq] = 1.0;
307 for (
auto row = overlapRows_.begin(); row != overlapRows_.end(); row++ )
314 matrix[lcell][lcell] = diag_block;
318 VectorType getTrueImpesWeights(
const VectorType& b,
const int pressureVarIndex)
const
320 VectorType weights(b.size());
321 ElementContext elemCtx(simulator_);
322 Opm::Amg::getTrueImpesWeights(pressureVarIndex, weights, simulator_.vanguard().gridView(),
323 elemCtx, simulator_.model(),
324 ThreadManager::threadId());
330 const int verbosity = prm_.get<
int>(
"verbosity");
331 const bool write_matrix = verbosity > 10;
333 Opm::Helper::writeSystem(this->simulator_,
341 const Simulator& simulator_;
343 std::unique_ptr<WellModelOpType> well_operator_;
344 std::unique_ptr<AbstractOperatorType> linear_operator_;
345 std::unique_ptr<SolverType> solver_;
346 FlowLinearSolverParameters parameters_;
349 Dune::InverseOperatorResult res_;
350 std::any parallelInformation_;
352 bool matrixAddWellContributions_;
353 int interiorCellNum_;
354 std::unique_ptr<Communication> comm_;
355 std::vector<int> overlapRows_;
356 std::vector<int> interiorRows_;
A solver class that encapsulates all needed objects for a linear solver (operator,...
Definition: FlexibleSolver.hpp:40
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition: ISTLSolverEbosFlexible.hpp:59
void makeOverlapRowsInvalid(MatrixType &matrix) const
Zero out off-diagonal blocks on rows corresponding to overlap cells Diagonal blocks on ovelap rows ar...
Definition: ISTLSolverEbosFlexible.hpp:298
Linear operator wrapper for well model.
Definition: WellOperators.hpp:50
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition: WellOperators.hpp:173
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition: WellOperators.hpp:100
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:26
void extractParallelGridInformationToISTL(std::any &anyComm, const UnstructuredGrid &grid)
Extracts the information about the data decomposition from the grid for dune-istl.
Definition: ParallelIstlInformation.hpp:676
PropertyTree setupPropertyTree(FlowLinearSolverParameters p, bool LinearSolverMaxIterSet, bool CprMaxEllIterSet)
Set up a property tree intended for FlexibleSolver by either reading the tree from a JSON file or cre...
Definition: setupPropertyTree.cpp:37
Definition: ISTLSolverEbosFlexible.hpp:37