My Project
ISTLSolverEbos.hpp
1/*
2 Copyright 2016 IRIS AS
3 Copyright 2019, 2020 Equinor ASA
4 Copyright 2020 SINTEF Digital, Mathematics and Cybernetics
5
6 This file is part of the Open Porous Media project (OPM).
7
8 OPM is free software: you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12
13 OPM is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with OPM. If not, see <http://www.gnu.org/licenses/>.
20*/
21
22#ifndef OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED
23#define OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED
24
25#include <dune/istl/owneroverlapcopy.hh>
26
27#include <ebos/eclbasevanguard.hh>
28
29#include <opm/common/ErrorMacros.hpp>
30
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>
46
47#include <any>
48#include <cstddef>
49#include <functional>
50#include <memory>
51#include <set>
52#include <sstream>
53#include <string>
54#include <tuple>
55#include <vector>
56
57namespace Opm::Properties {
58
59namespace TTag {
61 using InheritsFrom = std::tuple<FlowIstlSolverParams>;
62};
63}
64
65template <class TypeTag, class MyTypeTag>
67
70template<class TypeTag>
71struct SparseMatrixAdapter<TypeTag, TTag::FlowIstlSolver>
72{
73private:
74 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
75 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
77
78public:
79 using type = typename Linear::IstlSparseMatrixAdapter<Block>;
80};
81
82} // namespace Opm::Properties
83
84namespace Opm
85{
86
87#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
88template<class Matrix, class Vector, int block_size> class BdaBridge;
90#endif
91
92namespace detail
93{
94
95template<class Matrix, class Vector, class Comm>
97{
98 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
99 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
101
102 void create(const Matrix& matrix,
103 bool parallel,
104 const PropertyTree& prm,
105 size_t pressureIndex,
106 std::function<Vector()> trueFunc,
107 Comm& comm);
108
109 std::unique_ptr<AbstractSolverType> solver_;
110 std::unique_ptr<AbstractOperatorType> op_;
111 std::unique_ptr<LinearOperatorExtra<Vector,Vector>> wellOperator_;
112 AbstractPreconditionerType* pre_ = nullptr;
113 size_t interiorCellNum_ = 0;
114};
115
116#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
117template<class Matrix, class Vector>
118struct BdaSolverInfo
119{
120 using WellContribFunc = std::function<void(WellContributions&)>;
122
123 BdaSolverInfo(const std::string& accelerator_mode,
124 const std::string& fpga_bitstream,
125 const int linear_solver_verbosity,
126 const int maxit,
127 const double tolerance,
128 const int platformID,
129 const int deviceID,
130 const std::string& opencl_ilu_reorder,
131 const std::string& linsolver);
132
133 ~BdaSolverInfo();
134
135 template<class Grid>
136 void prepare(const Grid& grid,
137 const Dune::CartesianIndexMapper<Grid>& cartMapper,
138 const std::vector<Well>& wellsForConn,
139 const std::vector<int>& cellPartition,
140 const size_t nonzeroes,
141 const bool useWellConn);
142
143 bool apply(Vector& rhs,
144 const bool useWellConn,
145 WellContribFunc getContribs,
146 const int rank,
147 Matrix& matrix,
148 Vector& x,
149 Dune::InverseOperatorResult& result);
150
151 int numJacobiBlocks_ = 0;
152
153private:
156 template<class Grid>
157 void blockJacobiAdjacency(const Grid& grid,
158 const std::vector<int>& cell_part,
159 size_t nonzeroes);
160
161 void copyMatToBlockJac(const Matrix& mat, Matrix& blockJac);
162
163 std::unique_ptr<Bridge> bridge_;
164 std::string accelerator_mode_;
165 std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
166 std::vector<std::set<int>> wellConnectionsGraph_;
167};
168#endif
169
170#ifdef HAVE_MPI
172void copyParValues(std::any& parallelInformation, size_t size,
173 Dune::OwnerOverlapCopyCommunication<int,int>& comm);
174#endif
175
178template<class Matrix>
179void makeOverlapRowsInvalid(Matrix& matrix,
180 const std::vector<int>& overlapRows);
181
184template<class Matrix, class Grid>
185std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
186 const std::vector<int>& cell_part,
187 size_t nonzeroes,
188 const std::vector<std::set<int>>& wellConnectionsGraph);
189}
190
195 template <class TypeTag>
197 {
198 protected:
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;
215
216#if HAVE_MPI
217 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
218#else
219 using CommunicationType = Dune::CollectiveCommunication<int>;
220#endif
221
222 public:
223 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
224
225 static void registerParameters()
226 {
227 FlowLinearSolverParameters::registerParameters<TypeTag>();
228 }
229
233 explicit ISTLSolverEbos(const Simulator& simulator)
234 : simulator_(simulator),
235 iterations_( 0 ),
236 calls_( 0 ),
237 converged_(false),
238 matrix_()
239 {
240 const bool on_io_rank = (simulator.gridView().comm().rank() == 0);
241#if HAVE_MPI
242 comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) );
243#endif
244 parameters_.template init<TypeTag>();
245 prm_ = setupPropertyTree(parameters_,
246 EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter),
247 EWOMS_PARAM_IS_SET(TypeTag, int, CprMaxEllIter));
248
249#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
250 {
251 std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
252 if ((simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode != "none")) {
253 if (on_io_rank) {
254 OpmLog::warning("Cannot use GPU or FPGA with MPI, GPU/FPGA are disabled");
255 }
256 accelerator_mode = "none";
257 }
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,
267 fpga_bitstream,
268 linear_solver_verbosity,
269 maxit,
270 tolerance,
271 platformID,
272 deviceID,
273 opencl_ilu_reorder,
274 linsolver);
275 }
276#else
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");
279 }
280#endif
281 extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
282
283 // For some reason simulator_.model().elementMapper() is not initialized at this stage
284 //const auto& elemMapper = simulator_.model().elementMapper(); //does not work.
285 // Set it up manually
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);
289#if HAVE_FPGA
290 // check usage of MatrixAddWellContributions: for FPGA they must be included
291 if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) == "fpga" && !useWellConn_) {
292 OPM_THROW(std::logic_error,"fpgaSolver needs --matrix-add-well-contributions=true");
293 }
294#endif
295 const bool ownersFirst = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst);
296 if (!ownersFirst) {
297 const std::string msg = "The linear solver no longer supports --owner-cells-first=false.";
298 if (on_io_rank) {
299 OpmLog::error(msg);
300 }
301 OPM_THROW_NOLOG(std::runtime_error, msg);
302 }
303
304 flexibleSolver_.interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true);
305
306 // Print parameters to PRT/DBG logs.
307 if (on_io_rank) {
308 std::ostringstream os;
309 os << "Property tree for linear solver:\n";
310 prm_.write_json(os, true);
311 OpmLog::note(os.str());
312 }
313 }
314
315 // nothing to clean here
316 void eraseMatrix()
317 {
318 }
319
320 void prepare(const SparseMatrixAdapter& M, Vector& b)
321 {
322 static bool firstcall = true;
323#if HAVE_MPI
324 if (firstcall) {
325 const size_t size = M.istlMatrix().N();
326 detail::copyParValues(parallelInformation_, size, *comm_);
327 }
328#endif
329
330 // update matrix entries for solvers.
331 if (firstcall) {
332 // ebos will not change the matrix object. Hence simply store a pointer
333 // to the original one with a deleter that does nothing.
334 // Outch! We need to be able to scale the linear system! Hence const_cast
335 matrix_ = const_cast<Matrix*>(&M.istlMatrix());
336
337 useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
338 // setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver)
339#if HAVE_OPENCL
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_);
346#endif
347 } else {
348 // Pointers should not change
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()) );
352 }
353 }
354 rhs_ = &b;
355
356 if (isParallel() && prm_.get<std::string>("preconditioner.type") != "ParOverILU0") {
357 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
358 }
359 prepareFlexibleSolver();
360 firstcall = false;
361 }
362
363
364 void setResidual(Vector& /* b */)
365 {
366 // rhs_ = &b; // Must be handled in prepare() instead.
367 }
368
369 void getResidual(Vector& b) const
370 {
371 b = *rhs_;
372 }
373
374 void setMatrix(const SparseMatrixAdapter& /* M */)
375 {
376 // matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead.
377 }
378
379 bool solve(Vector& x)
380 {
381 calls_ += 1;
382 // Write linear system if asked for.
383 const int verbosity = prm_.get<int>("verbosity", 0);
384 const bool write_matrix = verbosity > 10;
385 if (write_matrix) {
386 Helper::writeSystem(simulator_, //simulator is only used to get names
387 getMatrix(),
388 *rhs_,
389 comm_.get());
390 }
391
392 // Solve system.
393 Dune::InverseOperatorResult result;
394
395 // Use GPU if: available, chosen by user, and successful.
396 // Use FPGA if: support compiled, chosen by user, and successful.
397#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
398 std::function<void(WellContributions&)> getContribs =
399 [this](WellContributions& w)
400 {
401 this->simulator_.problem().wellModel().getWellContributions(w);
402 };
403 if (!bdaBridge->apply(*rhs_, useWellConn_, getContribs,
404 simulator_.gridView().comm().rank(),
405 const_cast<Matrix&>(this->getMatrix()),
406 x, result))
407#endif
408 {
409 assert(flexibleSolver_.solver_);
410 flexibleSolver_.solver_->apply(x, *rhs_, result);
411 }
412
413 // Check convergence, iterations etc.
414 checkConvergence(result);
415
416 return converged_;
417 }
418
419
425
427 int iterations () const { return iterations_; }
428
430 const std::any& parallelInformation() const { return parallelInformation_; }
431
432 protected:
433#if HAVE_MPI
434 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
435#endif
436
437 void checkConvergence( const Dune::InverseOperatorResult& result ) const
438 {
439 // store number of iterations
440 iterations_ = result.iterations;
441 converged_ = result.converged;
442
443 // Check for failure of linear solver.
444 if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
445 const std::string msg("Convergence failure for linear solver.");
446 OPM_THROW_NOLOG(NumericalIssue, msg);
447 }
448 }
449 protected:
450
451 bool isParallel() const {
452#if HAVE_MPI
453 return comm_->communicator().size() > 1;
454#else
455 return false;
456#endif
457 }
458
459 void prepareFlexibleSolver()
460 {
461
462 if (shouldCreateSolver()) {
463 std::function<Vector()> trueFunc =
464 [this]
465 {
466 return this->getTrueImpesWeights(pressureIndex);
467 };
468
469 if (!useWellConn_) {
470 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
471 flexibleSolver_.wellOperator_ = std::move(wellOp);
472 }
473
474 flexibleSolver_.create(getMatrix(),
475 isParallel(),
476 prm_,
477 pressureIndex,
478 trueFunc,
479 *comm_);
480 }
481 else
482 {
483 flexibleSolver_.pre_->update();
484 }
485 }
486
487
491 {
492 // Decide if we should recreate the solver or just do
493 // a minimal preconditioner update.
494 if (!flexibleSolver_.solver_) {
495 return true;
496 }
497 if (this->parameters_.cpr_reuse_setup_ == 0) {
498 // Always recreate solver.
499 return true;
500 }
501 if (this->parameters_.cpr_reuse_setup_ == 1) {
502 // Recreate solver on the first iteration of every timestep.
503 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
504 return newton_iteration == 0;
505 }
506 if (this->parameters_.cpr_reuse_setup_ == 2) {
507 // Recreate solver if the last solve used more than 10 iterations.
508 return this->iterations() > 10;
509 }
510 if (this->parameters_.cpr_reuse_setup_ == 3) {
511 // Recreate solver if the last solve used more than 10 iterations.
512 return false;
513 }
514 if (this->parameters_.cpr_reuse_setup_ == 4) {
515 // Recreate solver every 'step' solve calls.
516 const int step = this->parameters_.cpr_reuse_interval_;
517 const bool create = ((calls_ % step) == 0);
518 return create;
519 }
520
521 // If here, we have an invalid parameter.
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.";
525 if (on_io_rank) {
526 OpmLog::error(msg);
527 }
528 throw std::runtime_error(msg);
529
530 // Never reached.
531 return false;
532 }
533
534
535 // Weights to make approximate pressure equations.
536 // Calculated from the storage terms (only) of the
537 // conservation equations, ignoring all other terms.
538 Vector getTrueImpesWeights(int pressureVarIndex) const
539 {
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());
546 return weights;
547 }
548
549 Matrix& getMatrix()
550 {
551 return *matrix_;
552 }
553
554 const Matrix& getMatrix() const
555 {
556 return *matrix_;
557 }
558
559 const Simulator& simulator_;
560 mutable int iterations_;
561 mutable int calls_;
562 mutable bool converged_;
563 std::any parallelInformation_;
564
565 // non-const to be able to scale the linear system
566 Matrix* matrix_;
567 Vector *rhs_;
568
569#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
570 std::unique_ptr<detail::BdaSolverInfo<Matrix, Vector>> bdaBridge;
571#endif
572
573 detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType> flexibleSolver_;
574 std::vector<int> overlapRows_;
575 std::vector<int> interiorRows_;
576
577 bool useWellConn_;
578
579 FlowLinearSolverParameters parameters_;
580 PropertyTree prm_;
581
582 std::shared_ptr< CommunicationType > comm_;
583 }; // end ISTLSolver
584
585} // namespace Opm
586#endif
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