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>
38 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
39 #include <opm/simulators/linalg/bda/BdaBridge.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);
152 if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) !=
"none") {
153 OPM_THROW(std::logic_error,
"Cannot use accelerated solver since CUDA, OpenCL and amgcl were not found by cmake and FPGA was not enabled");
161 ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
162 detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
164 useWellConn_ = EWOMS_GET_PARAM(TypeTag,
bool, MatrixAddWellContributions);
167 if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) ==
"fpga" && !useWellConn_) {
168 OPM_THROW(std::logic_error,
"fpgaSolver needs --matrix-add-well-contributions=true");
171 const bool ownersFirst = EWOMS_GET_PARAM(TypeTag,
bool, OwnerCellsFirst);
173 const std::string msg =
"The linear solver no longer supports --owner-cells-first=false.";
177 OPM_THROW_NOLOG(std::runtime_error, msg);
180 interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(),
true);
184 std::ostringstream os;
185 os <<
"Property tree for linear solver:\n";
186 prm_.write_json(os,
true);
187 OpmLog::note(os.str());
195 void prepare(
const SparseMatrixAdapter& M, Vector& b)
197 static bool firstcall =
true;
199 if (firstcall && parallelInformation_.type() ==
typeid(ParallelISTLInformation)) {
201 const ParallelISTLInformation* parinfo = std::any_cast<ParallelISTLInformation>(¶llelInformation_);
203 const size_t size = M.istlMatrix().N();
204 parinfo->copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), size, 1);
213 matrix_ =
const_cast<Matrix*
>(&M.istlMatrix());
216 if ( &(M.istlMatrix()) != matrix_ ) {
217 OPM_THROW(std::logic_error,
"Matrix objects are expected to be reused when reassembling!"
218 <<
" old pointer was " << matrix_ <<
", new one is " << (&M.istlMatrix()) );
223 if (isParallel() && prm_.get<std::string>(
"preconditioner.type") !=
"ParOverILU0") {
226 prepareFlexibleSolver();
231 void setResidual(Vector& ) {
235 void getResidual(Vector& b)
const {
239 void setMatrix(
const SparseMatrixAdapter& ) {
243 bool solve(Vector& x) {
245 const int verbosity = prm_.get<
int>(
"verbosity", 0);
246 const bool write_matrix = verbosity > 10;
248 Helper::writeSystem(simulator_,
255 Dune::InverseOperatorResult result;
256 bool accelerator_was_used =
false;
260 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
261 bool use_gpu = bdaBridge->getUseGpu();
262 bool use_fpga = bdaBridge->getUseFpga();
263 if (use_gpu || use_fpga) {
264 const std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
265 WellContributions wellContribs(accelerator_mode, useWellConn_);
266 bdaBridge->initWellContributions(wellContribs);
269 #if HAVE_CUDA || HAVE_OPENCL
271 simulator_.problem().wellModel().getWellContributions(wellContribs);
276 bdaBridge->solve_system(
const_cast<Matrix*
>(&getMatrix()), *rhs_, wellContribs, result);
277 if (result.converged) {
279 bdaBridge->get_result(x);
280 accelerator_was_used =
true;
286 if (simulator_.gridView().comm().rank() == 0) {
287 OpmLog::warning(bdaBridge->getAccleratorName() +
" did not converge, now trying Dune to solve current linear system...");
294 if (!accelerator_was_used) {
295 assert(flexibleSolver_);
296 flexibleSolver_->apply(x, *rhs_, result);
300 checkConvergence(result);
323 Matrix::block_type::rows,
324 Matrix::block_type::cols> >,
325 Vector, Vector> SeqPreconditioner;
328 typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
334 void checkConvergence(
const Dune::InverseOperatorResult& result )
const
337 iterations_ = result.iterations;
338 converged_ = result.converged;
341 if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
342 const std::string msg(
"Convergence failure for linear solver.");
343 OPM_THROW_NOLOG(NumericalIssue, msg);
348 bool isParallel()
const {
350 return comm_->communicator().size() > 1;
356 void prepareFlexibleSolver()
365 using ParOperatorType = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
366 linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
367 flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
370 using ParOperatorType = WellModelGhostLastMatrixAdapter<Matrix, Vector, Vector, true>;
371 wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
372 linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *wellOperator_, interiorCellNum_);
373 flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
379 using SeqOperatorType = Dune::MatrixAdapter<Matrix, Vector, Vector>;
380 linearOperatorForFlexibleSolver_ = std::make_unique<SeqOperatorType>(getMatrix());
381 flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator,
384 using SeqOperatorType = WellModelMatrixAdapter<Matrix, Vector, Vector, false>;
385 wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
386 linearOperatorForFlexibleSolver_ = std::make_unique<SeqOperatorType>(getMatrix(), *wellOperator_);
387 flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator,
394 flexibleSolver_->preconditioner().update();
405 if (!flexibleSolver_) {
408 if (this->parameters_.cpr_reuse_setup_ == 0) {
412 if (this->parameters_.cpr_reuse_setup_ == 1) {
414 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
415 return newton_iteration == 0;
417 if (this->parameters_.cpr_reuse_setup_ == 2) {
423 assert(this->parameters_.cpr_reuse_setup_ == 3);
432 std::function<Vector()> weightsCalculator;
434 using namespace std::string_literals;
436 auto preconditionerType = prm_.get(
"preconditioner.type"s,
"cpr"s);
437 if (preconditionerType ==
"cpr" || preconditionerType ==
"cprt") {
438 const bool transpose = preconditionerType ==
"cprt";
439 const auto weightsType = prm_.get(
"preconditioner.weight_type"s,
"quasiimpes"s);
440 if (weightsType ==
"quasiimpes") {
444 weightsCalculator = [
this, transpose, p = pressureIndex]() {
445 return Amg::getQuasiImpesWeights<Matrix, Vector>(this->getMatrix(), p, transpose);
447 }
else if (weightsType ==
"trueimpes") {
450 weightsCalculator = [
this, p = pressureIndex]() {
451 return this->getTrueImpesWeights(p);
454 OPM_THROW(std::invalid_argument,
455 "Weights type " << weightsType <<
"not implemented for cpr."
456 <<
" Please use quasiimpes or trueimpes.");
459 return weightsCalculator;
466 Vector getTrueImpesWeights(
int pressureVarIndex)
const
468 Vector weights(rhs_->size());
469 ElementContext elemCtx(simulator_);
470 Amg::getTrueImpesWeights(pressureVarIndex, weights, simulator_.vanguard().gridView(),
471 elemCtx, simulator_.model(),
472 ThreadManager::threadId());
482 const int numEq = Matrix::block_type::rows;
483 typename Matrix::block_type diag_block(0.0);
484 for (
int eq = 0; eq < numEq; ++eq)
485 diag_block[eq][eq] = 1.0;
488 for (
auto row = overlapRows_.begin(); row != overlapRows_.end(); row++ )
495 matrix[lcell][lcell] = diag_block;
505 const Matrix& getMatrix()
const
510 const Simulator& simulator_;
511 mutable int iterations_;
512 mutable bool converged_;
513 std::any parallelInformation_;
519 std::unique_ptr<FlexibleSolverType> flexibleSolver_;
520 std::unique_ptr<AbstractOperatorType> linearOperatorForFlexibleSolver_;
521 std::unique_ptr<WellModelAsLinearOperator<WellModel, Vector, Vector>> wellOperator_;
522 std::vector<int> overlapRows_;
523 std::vector<int> interiorRows_;
524 std::vector<std::set<int>> wellConnectionsGraph_;
527 size_t interiorCellNum_;
529 FlowLinearSolverParameters parameters_;
531 bool scale_variables_;
533 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:46
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:430
const std::any & parallelInformation() const
Definition: ISTLSolverEbos.hpp:316
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:479
int iterations() const
Solve the system of linear equations Ax = b, with A being the combined derivative matrix of the resid...
Definition: ISTLSolverEbos.hpp:313
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:401
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: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: ISTLSolverEbos.hpp:51
Definition: ISTLSolverEbos.hpp:45