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>
31#include <opm/models/discretization/common/fvbaseproperties.hh>
32#include <opm/models/common/multiphasebaseproperties.hh>
33#include <opm/models/utils/parametersystem.hh>
34#include <opm/models/utils/propertysystem.hh>
35#include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
36#include <opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp>
37#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
38#include <opm/simulators/linalg/matrixblock.hh>
39#include <opm/simulators/linalg/istlsparsematrixadapter.hh>
40#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
41#include <opm/simulators/linalg/WellOperators.hpp>
42#include <opm/simulators/linalg/WriteSystemMatrixHelper.hpp>
43#include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
44#include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
45#include <opm/simulators/linalg/setupPropertyTree.hpp>
57namespace Opm::Properties {
61 using InheritsFrom = std::tuple<FlowIstlSolverParams>;
65template <
class TypeTag,
class MyTypeTag>
70template<
class TypeTag>
71struct SparseMatrixAdapter<TypeTag, TTag::FlowIstlSolver>
74 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
75 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
79 using type =
typename Linear::IstlSparseMatrixAdapter<Block>;
87#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
88template<
class Matrix,
class Vector,
int block_size>
class BdaBridge;
95template<
class Matrix,
class Vector,
class Comm>
98 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
99 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
102 void create(
const Matrix& matrix,
105 size_t pressureIndex,
106 std::function<Vector()> trueFunc,
109 std::unique_ptr<AbstractSolverType> solver_;
110 std::unique_ptr<AbstractOperatorType> op_;
111 std::unique_ptr<LinearOperatorExtra<Vector,Vector>> wellOperator_;
113 size_t interiorCellNum_ = 0;
116#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
117template<
class Matrix,
class Vector>
123 BdaSolverInfo(
const std::string& accelerator_mode,
124 const std::string& fpga_bitstream,
125 const int linear_solver_verbosity,
127 const double tolerance,
128 const int platformID,
130 const std::string& opencl_ilu_reorder,
131 const std::string& linsolver);
136 void prepare(
const Grid& grid,
138 const std::vector<Well>& wellsForConn,
139 const std::vector<int>& cellPartition,
140 const size_t nonzeroes,
141 const bool useWellConn);
143 bool apply(Vector& rhs,
144 const bool useWellConn,
145 WellContribFunc getContribs,
149 Dune::InverseOperatorResult& result);
151 int numJacobiBlocks_ = 0;
157 void blockJacobiAdjacency(
const Grid& grid,
158 const std::vector<int>& cell_part,
161 void copyMatToBlockJac(
const Matrix& mat, Matrix& blockJac);
163 std::unique_ptr<Bridge> bridge_;
164 std::string accelerator_mode_;
165 std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
166 std::vector<std::set<int>> wellConnectionsGraph_;
172void copyParValues(std::any& parallelInformation,
size_t size,
173 Dune::OwnerOverlapCopyCommunication<int,int>& comm);
178template<
class Matrix>
179void makeOverlapRowsInvalid(Matrix& matrix,
180 const std::vector<int>& overlapRows);
184template<
class Matrix,
class Gr
id>
185std::unique_ptr<Matrix> blockJacobiAdjacency(
const Grid& grid,
186 const std::vector<int>& cell_part,
188 const std::vector<std::set<int>>& wellConnectionsGraph);
195 template <
class TypeTag>
199 using GridView = GetPropType<TypeTag, Properties::GridView>;
200 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
201 using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
202 using Vector = GetPropType<TypeTag, Properties::GlobalEqVector>;
203 using Indices = GetPropType<TypeTag, Properties::Indices>;
204 using WellModel = GetPropType<TypeTag, Properties::EclWellModel>;
205 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
206 using Matrix =
typename SparseMatrixAdapter::IstlMatrix;
207 using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
208 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
209 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
210 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
213 using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
214 constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
217 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
219 using CommunicationType = Dune::CollectiveCommunication<int>;
223 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
225 static void registerParameters()
227 FlowLinearSolverParameters::registerParameters<TypeTag>();
234 : simulator_(simulator),
240 const bool on_io_rank = (simulator.gridView().comm().rank() == 0);
242 comm_.reset(
new CommunicationType( simulator_.vanguard().grid().comm() ) );
244 parameters_.template init<TypeTag>();
246 EWOMS_PARAM_IS_SET(TypeTag,
int, LinearSolverMaxIter),
247 EWOMS_PARAM_IS_SET(TypeTag,
int, CprMaxEllIter));
249#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
251 std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
252 if ((simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode !=
"none")) {
254 OpmLog::warning(
"Cannot use GPU or FPGA with MPI, GPU/FPGA are disabled");
256 accelerator_mode =
"none";
258 const int platformID = EWOMS_GET_PARAM(TypeTag,
int, OpenclPlatformId);
259 const int deviceID = EWOMS_GET_PARAM(TypeTag,
int, BdaDeviceId);
260 const int maxit = EWOMS_GET_PARAM(TypeTag,
int, LinearSolverMaxIter);
261 const double tolerance = EWOMS_GET_PARAM(TypeTag,
double, LinearSolverReduction);
262 const std::string opencl_ilu_reorder = EWOMS_GET_PARAM(TypeTag, std::string, OpenclIluReorder);
263 const int linear_solver_verbosity = parameters_.linear_solver_verbosity_;
264 std::string fpga_bitstream = EWOMS_GET_PARAM(TypeTag, std::string, FpgaBitstream);
265 std::string linsolver = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolver);
266 bdaBridge = std::make_unique<detail::BdaSolverInfo<Matrix,Vector>>(accelerator_mode,
268 linear_solver_verbosity,
277 if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) !=
"none") {
278 OPM_THROW(std::logic_error,
"Cannot use accelerated solver since CUDA, OpenCL and amgcl were not found by cmake and FPGA was not enabled");
281 extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
286 ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
287 detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
288 useWellConn_ = EWOMS_GET_PARAM(TypeTag,
bool, MatrixAddWellContributions);
291 if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) ==
"fpga" && !useWellConn_) {
292 OPM_THROW(std::logic_error,
"fpgaSolver needs --matrix-add-well-contributions=true");
295 const bool ownersFirst = EWOMS_GET_PARAM(TypeTag,
bool, OwnerCellsFirst);
297 const std::string msg =
"The linear solver no longer supports --owner-cells-first=false.";
301 OPM_THROW_NOLOG(std::runtime_error, msg);
304 flexibleSolver_.interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(),
true);
308 std::ostringstream os;
309 os <<
"Property tree for linear solver:\n";
310 prm_.write_json(os,
true);
311 OpmLog::note(os.str());
320 void prepare(
const SparseMatrixAdapter& M, Vector& b)
322 static bool firstcall =
true;
325 const size_t size = M.istlMatrix().N();
326 detail::copyParValues(parallelInformation_, size, *comm_);
335 matrix_ =
const_cast<Matrix*
>(&M.istlMatrix());
337 useWellConn_ = EWOMS_GET_PARAM(TypeTag,
bool, MatrixAddWellContributions);
340 bdaBridge->numJacobiBlocks_ = EWOMS_GET_PARAM(TypeTag,
int, NumJacobiBlocks);
341 bdaBridge->prepare(simulator_.vanguard().grid(),
342 simulator_.vanguard().cartesianIndexMapper(),
343 simulator_.vanguard().schedule().getWellsatEnd(),
344 simulator_.vanguard().cellPartition(),
345 getMatrix().nonzeroes(), useWellConn_);
349 if ( &(M.istlMatrix()) != matrix_ ) {
350 OPM_THROW(std::logic_error,
"Matrix objects are expected to be reused when reassembling!"
351 <<
" old pointer was " << matrix_ <<
", new one is " << (&M.istlMatrix()) );
356 if (isParallel() && prm_.get<std::string>(
"preconditioner.type") !=
"ParOverILU0") {
357 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
359 prepareFlexibleSolver();
364 void setResidual(Vector& )
369 void getResidual(Vector& b)
const
374 void setMatrix(
const SparseMatrixAdapter& )
379 bool solve(Vector& x)
383 const int verbosity = prm_.get<
int>(
"verbosity", 0);
384 const bool write_matrix = verbosity > 10;
386 Helper::writeSystem(simulator_,
393 Dune::InverseOperatorResult result;
397#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
398 std::function<void(WellContributions&)> getContribs =
399 [
this](WellContributions& w)
401 this->simulator_.problem().wellModel().getWellContributions(w);
403 if (!bdaBridge->apply(*rhs_, useWellConn_, getContribs,
404 simulator_.gridView().comm().rank(),
405 const_cast<Matrix&
>(this->getMatrix()),
409 assert(flexibleSolver_.solver_);
410 flexibleSolver_.solver_->apply(x, *rhs_, result);
414 checkConvergence(result);
434 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
437 void checkConvergence(
const Dune::InverseOperatorResult& result )
const
440 iterations_ = result.iterations;
441 converged_ = result.converged;
444 if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
445 const std::string msg(
"Convergence failure for linear solver.");
446 OPM_THROW_NOLOG(NumericalIssue, msg);
451 bool isParallel()
const {
453 return comm_->communicator().size() > 1;
459 void prepareFlexibleSolver()
463 std::function<Vector()> trueFunc =
466 return this->getTrueImpesWeights(pressureIndex);
470 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
471 flexibleSolver_.wellOperator_ = std::move(wellOp);
474 flexibleSolver_.create(getMatrix(),
483 flexibleSolver_.pre_->update();
494 if (!flexibleSolver_.solver_) {
497 if (this->parameters_.cpr_reuse_setup_ == 0) {
501 if (this->parameters_.cpr_reuse_setup_ == 1) {
503 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
504 return newton_iteration == 0;
506 if (this->parameters_.cpr_reuse_setup_ == 2) {
510 if (this->parameters_.cpr_reuse_setup_ == 3) {
514 if (this->parameters_.cpr_reuse_setup_ == 4) {
516 const int step = this->parameters_.cpr_reuse_interval_;
517 const bool create = ((calls_ % step) == 0);
522 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
523 std::string msg =
"Invalid value: " + std::to_string(this->parameters_.cpr_reuse_setup_)
524 +
" for --cpr-reuse-setup parameter, run with --help to see allowed values.";
528 throw std::runtime_error(msg);
538 Vector getTrueImpesWeights(
int pressureVarIndex)
const
540 Vector weights(rhs_->size());
541 ElementContext elemCtx(simulator_);
542 Amg::getTrueImpesWeights(pressureVarIndex, weights,
543 simulator_.vanguard().gridView(),
544 elemCtx, simulator_.model(),
545 ThreadManager::threadId());
554 const Matrix& getMatrix()
const
559 const Simulator& simulator_;
560 mutable int iterations_;
562 mutable bool converged_;
563 std::any parallelInformation_;
569#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
570 std::unique_ptr<detail::BdaSolverInfo<Matrix, Vector>> bdaBridge;
573 detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType> flexibleSolver_;
574 std::vector<int> overlapRows_;
575 std::vector<int> interiorRows_;
579 FlowLinearSolverParameters parameters_;
582 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:40
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition: ISTLSolverEbos.hpp:197
int iterations() const
Solve the system of linear equations Ax = b, with A being the combined derivative matrix of the resid...
Definition: ISTLSolverEbos.hpp:427
const std::any & parallelInformation() const
Definition: ISTLSolverEbos.hpp:430
ISTLSolverEbos(const Simulator &simulator)
Construct a system solver.
Definition: ISTLSolverEbos.hpp:233
bool shouldCreateSolver() const
Return true if we should (re)create the whole solver, instead of just calling update() on the precond...
Definition: ISTLSolverEbos.hpp:490
Definition: MatrixMarketSpecializations.hpp:18
Definition: PropertyTree.hpp:37
This class serves to eliminate the need to include the WellContributions into the matrix (with –matri...
Definition: WellContributions.hpp:53
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 CprMaxEllIterSet)
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:66
Definition: ISTLSolverEbos.hpp:60
Definition: ISTLSolverEbos.hpp:97