22 #ifndef OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED
23 #define OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED
25 #include <opm/models/utils/parametersystem.hh>
26 #include <opm/models/utils/propertysystem.hh>
27 #include <opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp>
28 #include <opm/simulators/linalg/FlexibleSolver.hpp>
29 #include <opm/simulators/linalg/MatrixBlock.hpp>
30 #include <opm/simulators/linalg/ParallelIstlInformation.hpp>
31 #include <opm/simulators/linalg/WellOperators.hpp>
32 #include <opm/simulators/linalg/WriteSystemMatrixHelper.hpp>
33 #include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
34 #include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
35 #include <opm/simulators/linalg/setupPropertyTree.hpp>
37 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
38 #include <opm/simulators/linalg/bda/BdaBridge.hpp>
39 #include <opm/simulators/linalg/bda/WellContributions.hpp>
42 namespace Opm::Properties {
46 using InheritsFrom = std::tuple<FlowIstlSolverParams>;
50 template <
class TypeTag,
class MyTypeTag>
55 template<
class TypeTag>
56 struct SparseMatrixAdapter<TypeTag, TTag::FlowIstlSolver>
59 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
60 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
64 typedef typename Linear::IstlSparseMatrixAdapter<Block> type;
76 template <
class TypeTag>
80 using GridView = GetPropType<TypeTag, Properties::GridView>;
81 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
82 using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
83 using Vector = GetPropType<TypeTag, Properties::GlobalEqVector>;
84 using Indices = GetPropType<TypeTag, Properties::Indices>;
85 using WellModel = GetPropType<TypeTag, Properties::EclWellModel>;
86 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
87 using Matrix =
typename SparseMatrixAdapter::IstlMatrix;
88 using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
89 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
91 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
93 using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
94 constexpr
static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
96 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
97 static const unsigned int block_size = Matrix::block_type::rows;
98 std::unique_ptr<BdaBridge<Matrix, Vector, block_size>> bdaBridge;
102 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
104 using CommunicationType = Dune::CollectiveCommunication< int >;
108 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
110 static void registerParameters()
112 FlowLinearSolverParameters::registerParameters<TypeTag>();
119 : simulator_(simulator),
124 const bool on_io_rank = (simulator.gridView().comm().rank() == 0);
126 comm_.reset(
new CommunicationType( simulator_.vanguard().grid().comm() ) );
128 parameters_.template init<TypeTag>();
130 EWOMS_PARAM_IS_SET(TypeTag,
int, LinearSolverMaxIter),
131 EWOMS_PARAM_IS_SET(TypeTag,
int, CprMaxEllIter));
133 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
135 std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
136 if ((simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode !=
"none")) {
138 OpmLog::warning(
"Cannot use GPU or FPGA with MPI, GPU/FPGA are disabled");
140 accelerator_mode =
"none";
142 const int platformID = EWOMS_GET_PARAM(TypeTag,
int, OpenclPlatformId);
143 const int deviceID = EWOMS_GET_PARAM(TypeTag,
int, BdaDeviceId);
144 const int maxit = EWOMS_GET_PARAM(TypeTag,
int, LinearSolverMaxIter);
145 const double tolerance = EWOMS_GET_PARAM(TypeTag,
double, LinearSolverReduction);
146 const std::string opencl_ilu_reorder = EWOMS_GET_PARAM(TypeTag, std::string, OpenclIluReorder);
147 const int linear_solver_verbosity = parameters_.linear_solver_verbosity_;
148 std::string fpga_bitstream = EWOMS_GET_PARAM(TypeTag, std::string, FpgaBitstream);
149 std::string linsolver = EWOMS_GET_PARAM(TypeTag, std::string, Linsolver);
150 bdaBridge.reset(
new BdaBridge<Matrix, Vector, block_size>(accelerator_mode, fpga_bitstream, linear_solver_verbosity, maxit, tolerance, platformID, deviceID, opencl_ilu_reorder, linsolver));
153 if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) !=
"none") {
154 OPM_THROW(std::logic_error,
"Cannot use accelerated solver since CUDA, OpenCL and amgcl were not found by cmake and FPGA was not enabled");
162 ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
163 detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
165 useWellConn_ = EWOMS_GET_PARAM(TypeTag,
bool, MatrixAddWellContributions);
168 if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) ==
"fpga" && !useWellConn_) {
169 OPM_THROW(std::logic_error,
"fpgaSolver needs --matrix-add-well-contributions=true");
172 const bool ownersFirst = EWOMS_GET_PARAM(TypeTag,
bool, OwnerCellsFirst);
174 const std::string msg =
"The linear solver no longer supports --owner-cells-first=false.";
178 OPM_THROW_NOLOG(std::runtime_error, msg);
181 interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(),
true);
185 std::ostringstream os;
186 os <<
"Property tree for linear solver:\n";
187 prm_.write_json(os,
true);
188 OpmLog::note(os.str());
196 void prepare(
const SparseMatrixAdapter& M, Vector& b)
198 static bool firstcall =
true;
200 if (firstcall && parallelInformation_.type() ==
typeid(ParallelISTLInformation)) {
202 const ParallelISTLInformation* parinfo = std::any_cast<ParallelISTLInformation>(¶llelInformation_);
204 const size_t size = M.istlMatrix().N();
205 parinfo->copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), size, 1);
214 matrix_ =
const_cast<Matrix*
>(&M.istlMatrix());
217 if ( &(M.istlMatrix()) != matrix_ ) {
218 OPM_THROW(std::logic_error,
"Matrix objects are expected to be reused when reassembling!"
219 <<
" old pointer was " << matrix_ <<
", new one is " << (&M.istlMatrix()) );
224 if (isParallel() && prm_.get<std::string>(
"preconditioner.type") !=
"ParOverILU0") {
227 prepareFlexibleSolver();
232 void setResidual(Vector& ) {
236 void getResidual(Vector& b)
const {
240 void setMatrix(
const SparseMatrixAdapter& ) {
244 bool solve(Vector& x) {
246 const int verbosity = prm_.get<
int>(
"verbosity", 0);
247 const bool write_matrix = verbosity > 10;
249 Helper::writeSystem(simulator_,
256 Dune::InverseOperatorResult result;
257 bool accelerator_was_used =
false;
261 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
262 bool use_gpu = bdaBridge->getUseGpu();
263 bool use_fpga = bdaBridge->getUseFpga();
264 if (use_gpu || use_fpga) {
265 const std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
266 auto wellContribs = WellContributions::create(accelerator_mode, useWellConn_);
267 bdaBridge->initWellContributions(*wellContribs);
270 #if HAVE_CUDA || HAVE_OPENCL
272 simulator_.problem().wellModel().getWellContributions(*wellContribs);
277 bdaBridge->solve_system(
const_cast<Matrix*
>(&getMatrix()), *rhs_, *wellContribs, result);
278 if (result.converged) {
280 bdaBridge->get_result(x);
281 accelerator_was_used =
true;
287 if (simulator_.gridView().comm().rank() == 0) {
288 OpmLog::warning(bdaBridge->getAccleratorName() +
" did not converge, now trying Dune to solve current linear system...");
295 if (!accelerator_was_used) {
296 assert(flexibleSolver_);
297 flexibleSolver_->apply(x, *rhs_, result);
301 checkConvergence(result);
324 Matrix::block_type::rows,
325 Matrix::block_type::cols> >,
326 Vector, Vector> SeqPreconditioner;
329 typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
335 void checkConvergence(
const Dune::InverseOperatorResult& result )
const
338 iterations_ = result.iterations;
339 converged_ = result.converged;
342 if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
343 const std::string msg(
"Convergence failure for linear solver.");
344 OPM_THROW_NOLOG(NumericalIssue, msg);
349 bool isParallel()
const {
351 return comm_->communicator().size() > 1;
357 void prepareFlexibleSolver()
366 using ParOperatorType = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
367 linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
368 flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
371 using ParOperatorType = WellModelGhostLastMatrixAdapter<Matrix, Vector, Vector, true>;
372 wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
373 linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *wellOperator_, interiorCellNum_);
374 flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
380 using SeqOperatorType = Dune::MatrixAdapter<Matrix, Vector, Vector>;
381 linearOperatorForFlexibleSolver_ = std::make_unique<SeqOperatorType>(getMatrix());
382 flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator,
385 using SeqOperatorType = WellModelMatrixAdapter<Matrix, Vector, Vector, false>;
386 wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
387 linearOperatorForFlexibleSolver_ = std::make_unique<SeqOperatorType>(getMatrix(), *wellOperator_);
388 flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator,
395 flexibleSolver_->preconditioner().update();
406 if (!flexibleSolver_) {
409 if (this->parameters_.cpr_reuse_setup_ == 0) {
413 if (this->parameters_.cpr_reuse_setup_ == 1) {
415 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
416 return newton_iteration == 0;
418 if (this->parameters_.cpr_reuse_setup_ == 2) {
424 assert(this->parameters_.cpr_reuse_setup_ == 3);
433 std::function<Vector()> weightsCalculator;
435 using namespace std::string_literals;
437 auto preconditionerType = prm_.get(
"preconditioner.type"s,
"cpr"s);
438 if (preconditionerType ==
"cpr" || preconditionerType ==
"cprt") {
439 const bool transpose = preconditionerType ==
"cprt";
440 const auto weightsType = prm_.get(
"preconditioner.weight_type"s,
"quasiimpes"s);
441 if (weightsType ==
"quasiimpes") {
445 weightsCalculator = [
this, transpose, p = pressureIndex]() {
446 return Amg::getQuasiImpesWeights<Matrix, Vector>(this->getMatrix(), p, transpose);
448 }
else if (weightsType ==
"trueimpes") {
451 weightsCalculator = [
this, p = pressureIndex]() {
452 return this->getTrueImpesWeights(p);
455 OPM_THROW(std::invalid_argument,
456 "Weights type " << weightsType <<
"not implemented for cpr."
457 <<
" Please use quasiimpes or trueimpes.");
460 return weightsCalculator;
467 Vector getTrueImpesWeights(
int pressureVarIndex)
const
469 Vector weights(rhs_->size());
470 ElementContext elemCtx(simulator_);
471 Amg::getTrueImpesWeights(pressureVarIndex, weights, simulator_.vanguard().gridView(),
472 elemCtx, simulator_.model(),
473 ThreadManager::threadId());
483 const int numEq = Matrix::block_type::rows;
484 typename Matrix::block_type diag_block(0.0);
485 for (
int eq = 0; eq < numEq; ++eq)
486 diag_block[eq][eq] = 1.0;
489 for (
auto row = overlapRows_.begin(); row != overlapRows_.end(); row++ )
496 matrix[lcell][lcell] = diag_block;
506 const Matrix& getMatrix()
const
511 const Simulator& simulator_;
512 mutable int iterations_;
513 mutable bool converged_;
514 std::any parallelInformation_;
520 std::unique_ptr<FlexibleSolverType> flexibleSolver_;
521 std::unique_ptr<AbstractOperatorType> linearOperatorForFlexibleSolver_;
522 std::unique_ptr<WellModelAsLinearOperator<WellModel, Vector, Vector>> wellOperator_;
523 std::vector<int> overlapRows_;
524 std::vector<int> interiorRows_;
525 std::vector<std::set<int>> wellConnectionsGraph_;
528 size_t interiorCellNum_;
530 FlowLinearSolverParameters parameters_;
532 bool scale_variables_;
534 std::shared_ptr< CommunicationType > comm_;
A solver class that encapsulates all needed objects for a linear solver (operator,...
Definition: FlexibleSolver.hpp:40
Definition: MatrixBlock.hpp:370
BdaBridge acts as interface between opm-simulators with the BdaSolvers.
Definition: BdaBridge.hpp:39
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition: ISTLSolverEbos.hpp:78
std::function< Vector()> getWeightsCalculator() const
Return an appropriate weight function if a cpr preconditioner is asked for.
Definition: ISTLSolverEbos.hpp:431
const std::any & parallelInformation() const
Definition: ISTLSolverEbos.hpp:317
void makeOverlapRowsInvalid(Matrix &matrix) const
Zero out off-diagonal blocks on rows corresponding to overlap cells Diagonal blocks on ovelap rows ar...
Definition: ISTLSolverEbos.hpp:480
int iterations() const
Solve the system of linear equations Ax = b, with A being the combined derivative matrix of the resid...
Definition: ISTLSolverEbos.hpp:314
ISTLSolverEbos(const Simulator &simulator)
Construct a system solver.
Definition: ISTLSolverEbos.hpp:118
bool shouldCreateSolver() const
Return true if we should (re)create the whole solver, instead of just calling update() on the precond...
Definition: ISTLSolverEbos.hpp:402
Definition: MatrixMarketSpecializations.hpp:17
A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as.
Definition: ParallelOverlappingILU0.hpp:611
Linear operator wrapper for well model.
Definition: WellOperators.hpp:50
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
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:36
Definition: ISTLSolverEbos.hpp:51
Definition: ISTLSolverEbos.hpp:45