22#ifndef OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED
23#define OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED
25#include <dune/istl/owneroverlapcopy.hh>
27#include <ebos/eclbasevanguard.hh>
29#include <opm/common/ErrorMacros.hpp>
30#include <opm/common/Exceptions.hpp>
31#include <opm/common/TimingMacros.hpp>
33#include <opm/models/discretization/common/fvbaseproperties.hh>
34#include <opm/models/common/multiphasebaseproperties.hh>
35#include <opm/models/utils/parametersystem.hh>
36#include <opm/models/utils/propertysystem.hh>
37#include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
38#include <opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp>
39#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
40#include <opm/simulators/linalg/matrixblock.hh>
41#include <opm/simulators/linalg/istlsparsematrixadapter.hh>
42#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
43#include <opm/simulators/linalg/WellOperators.hpp>
44#include <opm/simulators/linalg/WriteSystemMatrixHelper.hpp>
45#include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
46#include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
47#include <opm/simulators/linalg/setupPropertyTree.hpp>
59namespace Opm::Properties {
63 using InheritsFrom = std::tuple<FlowIstlSolverParams>;
67template <
class TypeTag,
class MyTypeTag>
72template<
class TypeTag>
73struct SparseMatrixAdapter<TypeTag, TTag::FlowIstlSolver>
76 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
77 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
81 using type =
typename Linear::IstlSparseMatrixAdapter<Block>;
90template<
class Matrix,
class Vector,
int block_size>
class BdaBridge;
97template<
class Matrix,
class Vector,
class Comm>
100 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
101 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
104 void create(
const Matrix& matrix,
107 size_t pressureIndex,
108 std::function<Vector()> trueFunc,
111 std::unique_ptr<AbstractSolverType> solver_;
112 std::unique_ptr<AbstractOperatorType> op_;
113 std::unique_ptr<LinearOperatorExtra<Vector,Vector>> wellOperator_;
115 size_t interiorCellNum_ = 0;
118#if COMPILE_BDA_BRIDGE
119template<
class Matrix,
class Vector>
125 BdaSolverInfo(
const std::string& accelerator_mode,
126 const int linear_solver_verbosity,
128 const double tolerance,
129 const int platformID,
131 const bool opencl_ilu_parallel,
132 const std::string& linsolver);
137 void prepare(
const Grid& grid,
139 const std::vector<Well>& wellsForConn,
140 const std::vector<int>& cellPartition,
141 const size_t nonzeroes,
142 const bool useWellConn);
144 bool apply(Vector& rhs,
145 const bool useWellConn,
146 WellContribFunc getContribs,
150 Dune::InverseOperatorResult& result);
154 int numJacobiBlocks_ = 0;
160 void blockJacobiAdjacency(
const Grid& grid,
161 const std::vector<int>& cell_part,
164 void copyMatToBlockJac(
const Matrix& mat, Matrix& blockJac);
166 std::unique_ptr<Bridge> bridge_;
167 std::string accelerator_mode_;
168 std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
169 std::vector<std::set<int>> wellConnectionsGraph_;
175void copyParValues(std::any& parallelInformation,
size_t size,
176 Dune::OwnerOverlapCopyCommunication<int,int>& comm);
181template<
class Matrix>
182void makeOverlapRowsInvalid(Matrix& matrix,
183 const std::vector<int>& overlapRows);
187template<
class Matrix,
class Gr
id>
188std::unique_ptr<Matrix> blockJacobiAdjacency(
const Grid& grid,
189 const std::vector<int>& cell_part,
191 const std::vector<std::set<int>>& wellConnectionsGraph);
198 template <
class TypeTag>
202 using GridView = GetPropType<TypeTag, Properties::GridView>;
203 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
204 using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
205 using Vector = GetPropType<TypeTag, Properties::GlobalEqVector>;
206 using Indices = GetPropType<TypeTag, Properties::Indices>;
207 using WellModel = GetPropType<TypeTag, Properties::EclWellModel>;
208 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
209 using Matrix =
typename SparseMatrixAdapter::IstlMatrix;
210 using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
211 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
212 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
213 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
216 using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
217 constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
220 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
222 using CommunicationType = Dune::CollectiveCommunication<int>;
226 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
228 static void registerParameters()
230 FlowLinearSolverParameters::registerParameters<TypeTag>();
237 : simulator_(simulator),
243 OPM_TIMEBLOCK(IstlSolverEbos);
244 const bool on_io_rank = (simulator.gridView().comm().rank() == 0);
246 comm_.reset(
new CommunicationType( simulator_.vanguard().grid().comm() ) );
248 parameters_.template init<TypeTag>(simulator_.vanguard().eclState().getSimulationConfig().useCPR());
250 EWOMS_PARAM_IS_SET(TypeTag,
int, LinearSolverMaxIter),
251 EWOMS_PARAM_IS_SET(TypeTag,
double, LinearSolverReduction));
253#if COMPILE_BDA_BRIDGE
255 std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
256 if ((simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode !=
"none")) {
258 OpmLog::warning(
"Cannot use GPU with MPI, GPU are disabled");
260 accelerator_mode =
"none";
262 const int platformID = EWOMS_GET_PARAM(TypeTag,
int, OpenclPlatformId);
263 const int deviceID = EWOMS_GET_PARAM(TypeTag,
int, BdaDeviceId);
264 const int maxit = EWOMS_GET_PARAM(TypeTag,
int, LinearSolverMaxIter);
265 const double tolerance = EWOMS_GET_PARAM(TypeTag,
double, LinearSolverReduction);
266 const bool opencl_ilu_parallel = EWOMS_GET_PARAM(TypeTag,
bool, OpenclIluParallel);
267 const int linear_solver_verbosity = parameters_.linear_solver_verbosity_;
268 std::string linsolver = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolver);
269 bdaBridge = std::make_unique<detail::BdaSolverInfo<Matrix,Vector>>(accelerator_mode,
270 linear_solver_verbosity,
279 if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) !=
"none") {
280 OPM_THROW(std::logic_error,
"Cannot use accelerated solver since CUDA, OpenCL and amgcl were not found by cmake");
283 extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
288 ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
289 detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
290 useWellConn_ = EWOMS_GET_PARAM(TypeTag,
bool, MatrixAddWellContributions);
291 const bool ownersFirst = EWOMS_GET_PARAM(TypeTag,
bool, OwnerCellsFirst);
293 const std::string msg =
"The linear solver no longer supports --owner-cells-first=false.";
297 OPM_THROW_NOLOG(std::runtime_error, msg);
300 flexibleSolver_.interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(),
true);
304 std::ostringstream os;
305 os <<
"Property tree for linear solver:\n";
306 prm_.write_json(os,
true);
307 OpmLog::note(os.str());
316 void prepare(
const SparseMatrixAdapter& M, Vector& b)
318 OPM_TIMEBLOCK(istlSolverEbosPrepare);
319 static bool firstcall =
true;
322 const size_t size = M.istlMatrix().N();
323 detail::copyParValues(parallelInformation_, size, *comm_);
332 matrix_ =
const_cast<Matrix*
>(&M.istlMatrix());
334 useWellConn_ = EWOMS_GET_PARAM(TypeTag,
bool, MatrixAddWellContributions);
337 bdaBridge->numJacobiBlocks_ = EWOMS_GET_PARAM(TypeTag,
int, NumJacobiBlocks);
338 bdaBridge->prepare(simulator_.vanguard().grid(),
339 simulator_.vanguard().cartesianIndexMapper(),
340 simulator_.vanguard().schedule().getWellsatEnd(),
341 simulator_.vanguard().cellPartition(),
342 getMatrix().nonzeroes(), useWellConn_);
346 if ( &(M.istlMatrix()) != matrix_ ) {
347 OPM_THROW(std::logic_error,
348 "Matrix objects are expected to be reused when reassembling!");
353 if (isParallel() && prm_.get<std::string>(
"preconditioner.type") !=
"ParOverILU0") {
354 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
356#if COMPILE_BDA_BRIDGE
357 if(!bdaBridge->gpuActive()){
358 prepareFlexibleSolver();
361 prepareFlexibleSolver();
368 void setResidual(Vector& )
373 void getResidual(Vector& b)
const
378 void setMatrix(
const SparseMatrixAdapter& )
383 bool solve(Vector& x)
385 OPM_TIMEBLOCK(istlSolverEbosSolve);
388 const int verbosity = prm_.get<
int>(
"verbosity", 0);
389 const bool write_matrix = verbosity > 10;
391 Helper::writeSystem(simulator_,
398 Dune::InverseOperatorResult result;
400#if COMPILE_BDA_BRIDGE
401 std::function<void(WellContributions&)> getContribs =
402 [
this](WellContributions& w)
404 this->simulator_.problem().wellModel().getWellContributions(w);
406 if (!bdaBridge->apply(*rhs_, useWellConn_, getContribs,
407 simulator_.gridView().comm().rank(),
408 const_cast<Matrix&
>(this->getMatrix()),
412 OPM_TIMEBLOCK(flexibleSolverApply);
413#if COMPILE_BDA_BRIDGE
414 if(bdaBridge->gpuActive()){
415 prepareFlexibleSolver();
418 assert(flexibleSolver_.solver_);
419 flexibleSolver_.solver_->apply(x, *rhs_, result);
423 checkConvergence(result);
443 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
446 void checkConvergence(
const Dune::InverseOperatorResult& result )
const
449 iterations_ = result.iterations;
450 converged_ = result.converged;
452 if(result.reduction < parameters_.relaxed_linear_solver_reduction_){
453 std::stringstream ss;
454 ss<<
"Full linear solver tolerance not achieved. The reduction is:" << result.reduction
455 <<
" after " << result.iterations <<
" iterations ";
456 OpmLog::warning(ss.str());
461 if (!parameters_.ignoreConvergenceFailure_ && !converged_) {
462 const std::string msg(
"Convergence failure for linear solver.");
463 OPM_THROW_NOLOG(NumericalProblem, msg);
468 bool isParallel()
const {
470 return comm_->communicator().size() > 1;
476 void prepareFlexibleSolver()
478 OPM_TIMEBLOCK(flexibleSolverPrepare);
480 std::function<Vector()> trueFunc =
483 return this->getTrueImpesWeights(pressureIndex);
487 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
488 flexibleSolver_.wellOperator_ = std::move(wellOp);
490 OPM_TIMEBLOCK(flexibleSolverCreate);
491 flexibleSolver_.create(getMatrix(),
500 OPM_TIMEBLOCK(flexibleSolverUpdate);
501 flexibleSolver_.pre_->update();
512 if (!flexibleSolver_.solver_) {
515 if (this->parameters_.cpr_reuse_setup_ == 0) {
519 if (this->parameters_.cpr_reuse_setup_ == 1) {
521 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
522 return newton_iteration == 0;
524 if (this->parameters_.cpr_reuse_setup_ == 2) {
528 if (this->parameters_.cpr_reuse_setup_ == 3) {
532 if (this->parameters_.cpr_reuse_setup_ == 4) {
534 const int step = this->parameters_.cpr_reuse_interval_;
535 const bool create = ((calls_ % step) == 0);
540 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
541 std::string msg =
"Invalid value: " + std::to_string(this->parameters_.cpr_reuse_setup_)
542 +
" for --cpr-reuse-setup parameter, run with --help to see allowed values.";
546 throw std::runtime_error(msg);
556 Vector getTrueImpesWeights(
int pressureVarIndex)
const
558 OPM_TIMEBLOCK(getTrueImpesWeights);
559 Vector weights(rhs_->size());
560 ElementContext elemCtx(simulator_);
561 Amg::getTrueImpesWeights(pressureVarIndex, weights,
562 simulator_.vanguard().gridView(),
563 elemCtx, simulator_.model(),
564 ThreadManager::threadId());
573 const Matrix& getMatrix()
const
578 const Simulator& simulator_;
579 mutable int iterations_;
581 mutable bool converged_;
582 std::any parallelInformation_;
588#if COMPILE_BDA_BRIDGE
589 std::unique_ptr<detail::BdaSolverInfo<Matrix, Vector>> bdaBridge;
592 detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType> flexibleSolver_;
593 std::vector<int> overlapRows_;
594 std::vector<int> interiorRows_;
598 FlowLinearSolverParameters parameters_;
601 std::shared_ptr< CommunicationType > comm_;
Definition: findOverlapRowsAndColumns.hpp:29
Interface class adding the update() method to the preconditioner interface.
Definition: PreconditionerWithUpdate.hpp:32
BdaBridge acts as interface between opm-simulators with the BdaSolvers.
Definition: BdaBridge.hpp:37
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition: ISTLSolverEbos.hpp:200
int iterations() const
Solve the system of linear equations Ax = b, with A being the combined derivative matrix of the resid...
Definition: ISTLSolverEbos.hpp:436
const std::any & parallelInformation() const
Definition: ISTLSolverEbos.hpp:439
ISTLSolverEbos(const Simulator &simulator)
Construct a system solver.
Definition: ISTLSolverEbos.hpp:236
bool shouldCreateSolver() const
Return true if we should (re)create the whole solver, instead of just calling update() on the precond...
Definition: ISTLSolverEbos.hpp:508
Definition: MatrixMarketSpecializations.hpp:17
Definition: PropertyTree.hpp:37
This class serves to eliminate the need to include the WellContributions into the matrix (with –matri...
Definition: WellContributions.hpp:52
Definition: WellOperators.hpp:60
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
PropertyTree setupPropertyTree(FlowLinearSolverParameters p, bool linearSolverMaxIterSet, bool linearSolverReductionSet)
Set up a property tree intended for FlexibleSolver by either reading the tree from a JSON file or cre...
Definition: setupPropertyTree.cpp:40
Definition: ISTLSolverEbos.hpp:68
Definition: ISTLSolverEbos.hpp:62
Definition: ISTLSolverEbos.hpp:99