My Project
Loading...
Searching...
No Matches
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#include <opm/common/Exceptions.hpp>
31#include <opm/common/TimingMacros.hpp>
32
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>
48
49#include <any>
50#include <cstddef>
51#include <functional>
52#include <memory>
53#include <set>
54#include <sstream>
55#include <string>
56#include <tuple>
57#include <vector>
58
59namespace Opm::Properties {
60
61namespace TTag {
63 using InheritsFrom = std::tuple<FlowIstlSolverParams>;
64};
65}
66
67template <class TypeTag, class MyTypeTag>
69
72template<class TypeTag>
73struct SparseMatrixAdapter<TypeTag, TTag::FlowIstlSolver>
74{
75private:
79
80public:
81 using type = typename Linear::IstlSparseMatrixAdapter<Block>;
82};
83
84} // namespace Opm::Properties
85
86namespace Opm
87{
88
89#if COMPILE_BDA_BRIDGE
90template<class Matrix, class Vector, int block_size> class BdaBridge;
92#endif
93
94namespace detail
95{
96
97template<class Matrix, class Vector, class Comm>
98struct FlexibleSolverInfo
99{
100 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
101 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
102 using AbstractPreconditionerType = Dune::PreconditionerWithUpdate<Vector, Vector>;
103
104 void create(const Matrix& matrix,
105 bool parallel,
106 const PropertyTree& prm,
107 size_t pressureIndex,
108 std::function<Vector()> trueFunc,
109 Comm& comm);
110
111 std::unique_ptr<AbstractSolverType> solver_;
112 std::unique_ptr<AbstractOperatorType> op_;
113 std::unique_ptr<LinearOperatorExtra<Vector,Vector>> wellOperator_;
114 AbstractPreconditionerType* pre_ = nullptr;
115 size_t interiorCellNum_ = 0;
116};
117
118#if COMPILE_BDA_BRIDGE
119template<class Matrix, class Vector>
120struct BdaSolverInfo
121{
122 using WellContribFunc = std::function<void(WellContributions&)>;
123 using Bridge = BdaBridge<Matrix,Vector,Matrix::block_type::rows>;
124
125 BdaSolverInfo(const std::string& accelerator_mode,
126 const int linear_solver_verbosity,
127 const int maxit,
128 const double tolerance,
129 const int platformID,
130 const int deviceID,
131 const bool opencl_ilu_parallel,
132 const std::string& linsolver);
133
134 ~BdaSolverInfo();
135
136 template<class Grid>
137 void prepare(const Grid& grid,
138 const Dune::CartesianIndexMapper<Grid>& cartMapper,
139 const std::vector<Well>& wellsForConn,
140 const std::vector<int>& cellPartition,
141 const size_t nonzeroes,
142 const bool useWellConn);
143
144 bool apply(Vector& rhs,
145 const bool useWellConn,
146 WellContribFunc getContribs,
147 const int rank,
148 Matrix& matrix,
149 Vector& x,
150 Dune::InverseOperatorResult& result);
151
152 bool gpuActive();
153
154 int numJacobiBlocks_ = 0;
155
156private:
159 template<class Grid>
160 void blockJacobiAdjacency(const Grid& grid,
161 const std::vector<int>& cell_part,
162 size_t nonzeroes);
163
164 void copyMatToBlockJac(const Matrix& mat, Matrix& blockJac);
165
166 std::unique_ptr<Bridge> bridge_;
167 std::string accelerator_mode_;
168 std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
169 std::vector<std::set<int>> wellConnectionsGraph_;
170};
171#endif
172
173#ifdef HAVE_MPI
175void copyParValues(std::any& parallelInformation, size_t size,
176 Dune::OwnerOverlapCopyCommunication<int,int>& comm);
177#endif
178
181template<class Matrix>
182void makeOverlapRowsInvalid(Matrix& matrix,
183 const std::vector<int>& overlapRows);
184
187template<class Matrix, class Grid>
188std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
189 const std::vector<int>& cell_part,
190 size_t nonzeroes,
191 const std::vector<std::set<int>>& wellConnectionsGraph);
192}
193
198 template <class TypeTag>
200 {
201 protected:
209 using Matrix = typename SparseMatrixAdapter::IstlMatrix;
212 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
213 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
217 constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
218
219#if HAVE_MPI
220 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
221#else
222 using CommunicationType = Dune::CollectiveCommunication<int>;
223#endif
224
225 public:
226 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
227
228 static void registerParameters()
229 {
230 FlowLinearSolverParameters::registerParameters<TypeTag>();
231 }
232
236 explicit ISTLSolverEbos(const Simulator& simulator)
237 : simulator_(simulator),
238 iterations_( 0 ),
239 calls_( 0 ),
240 converged_(false),
241 matrix_()
242 {
244 const bool on_io_rank = (simulator.gridView().comm().rank() == 0);
245#if HAVE_MPI
246 comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) );
247#endif
248 parameters_.template init<TypeTag>(simulator_.vanguard().eclState().getSimulationConfig().useCPR());
249 prm_ = setupPropertyTree(parameters_,
250 EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter),
251 EWOMS_PARAM_IS_SET(TypeTag, double, LinearSolverReduction));
252
253#if COMPILE_BDA_BRIDGE
254 {
255 std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
256 if ((simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode != "none")) {
257 if (on_io_rank) {
258 OpmLog::warning("Cannot use GPU with MPI, GPU are disabled");
259 }
260 accelerator_mode = "none";
261 }
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,
271 maxit,
272 tolerance,
273 platformID,
274 deviceID,
275 opencl_ilu_parallel,
276 linsolver);
277 }
278#else
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");
281 }
282#endif
283 extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
284
285 // For some reason simulator_.model().elementMapper() is not initialized at this stage
286 //const auto& elemMapper = simulator_.model().elementMapper(); //does not work.
287 // Set it up manually
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);
292 if (!ownersFirst) {
293 const std::string msg = "The linear solver no longer supports --owner-cells-first=false.";
294 if (on_io_rank) {
295 OpmLog::error(msg);
296 }
297 OPM_THROW_NOLOG(std::runtime_error, msg);
298 }
299
300 flexibleSolver_.interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true);
301
302 // Print parameters to PRT/DBG logs.
303 if (on_io_rank) {
304 std::ostringstream os;
305 os << "Property tree for linear solver:\n";
306 prm_.write_json(os, true);
307 OpmLog::note(os.str());
308 }
309 }
310
311 // nothing to clean here
312 void eraseMatrix()
313 {
314 }
315
316 void prepare(const SparseMatrixAdapter& M, Vector& b)
317 {
318 OPM_TIMEBLOCK(istlSolverEbosPrepare);
319 static bool firstcall = true;
320#if HAVE_MPI
321 if (firstcall) {
322 const size_t size = M.istlMatrix().N();
323 detail::copyParValues(parallelInformation_, size, *comm_);
324 }
325#endif
326
327 // update matrix entries for solvers.
328 if (firstcall) {
329 // ebos will not change the matrix object. Hence simply store a pointer
330 // to the original one with a deleter that does nothing.
331 // Outch! We need to be able to scale the linear system! Hence const_cast
332 matrix_ = const_cast<Matrix*>(&M.istlMatrix());
333
334 useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
335 // setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver)
336#if HAVE_OPENCL
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_);
343#endif
344 } else {
345 // Pointers should not change
346 if ( &(M.istlMatrix()) != matrix_ ) {
347 OPM_THROW(std::logic_error,
348 "Matrix objects are expected to be reused when reassembling!");
349 }
350 }
351 rhs_ = &b;
352
353 if (isParallel() && prm_.get<std::string>("preconditioner.type") != "ParOverILU0") {
354 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
355 }
356#if COMPILE_BDA_BRIDGE
357 if(!bdaBridge->gpuActive()){
358 prepareFlexibleSolver();
359 }
360#else
361 prepareFlexibleSolver();
362#endif
363
364 firstcall = false;
365 }
366
367
368 void setResidual(Vector& /* b */)
369 {
370 // rhs_ = &b; // Must be handled in prepare() instead.
371 }
372
373 void getResidual(Vector& b) const
374 {
375 b = *rhs_;
376 }
377
378 void setMatrix(const SparseMatrixAdapter& /* M */)
379 {
380 // matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead.
381 }
382
383 bool solve(Vector& x)
384 {
385 OPM_TIMEBLOCK(istlSolverEbosSolve);
386 calls_ += 1;
387 // Write linear system if asked for.
388 const int verbosity = prm_.get<int>("verbosity", 0);
389 const bool write_matrix = verbosity > 10;
390 if (write_matrix) {
391 Helper::writeSystem(simulator_, //simulator is only used to get names
392 getMatrix(),
393 *rhs_,
394 comm_.get());
395 }
396
397 // Solve system.
398 Dune::InverseOperatorResult result;
399
400#if COMPILE_BDA_BRIDGE
401 std::function<void(WellContributions&)> getContribs =
402 [this](WellContributions& w)
403 {
404 this->simulator_.problem().wellModel().getWellContributions(w);
405 };
406 if (!bdaBridge->apply(*rhs_, useWellConn_, getContribs,
407 simulator_.gridView().comm().rank(),
408 const_cast<Matrix&>(this->getMatrix()),
409 x, result))
410#endif
411 {
412 OPM_TIMEBLOCK(flexibleSolverApply);
413#if COMPILE_BDA_BRIDGE
414 if(bdaBridge->gpuActive()){
415 prepareFlexibleSolver();
416 }
417#endif
418 assert(flexibleSolver_.solver_);
419 flexibleSolver_.solver_->apply(x, *rhs_, result);
420 }
421
422 // Check convergence, iterations etc.
423 checkConvergence(result);
424
425 return converged_;
426 }
427
428
434
436 int iterations () const { return iterations_; }
437
439 const std::any& parallelInformation() const { return parallelInformation_; }
440
441 protected:
442#if HAVE_MPI
443 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
444#endif
445
446 void checkConvergence( const Dune::InverseOperatorResult& result ) const
447 {
448 // store number of iterations
449 iterations_ = result.iterations;
450 converged_ = result.converged;
451 if(!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());
457 converged_ = true;
458 }
459 }
460 // Check for failure of linear solver.
461 if (!parameters_.ignoreConvergenceFailure_ && !converged_) {
462 const std::string msg("Convergence failure for linear solver.");
463 OPM_THROW_NOLOG(NumericalProblem, msg);
464 }
465 }
466 protected:
467
468 bool isParallel() const {
469#if HAVE_MPI
470 return comm_->communicator().size() > 1;
471#else
472 return false;
473#endif
474 }
475
476 void prepareFlexibleSolver()
477 {
478 OPM_TIMEBLOCK(flexibleSolverPrepare);
479 if (shouldCreateSolver()) {
480 std::function<Vector()> trueFunc =
481 [this]
482 {
483 return this->getTrueImpesWeights(pressureIndex);
484 };
485
486 if (!useWellConn_) {
487 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
488 flexibleSolver_.wellOperator_ = std::move(wellOp);
489 }
490 OPM_TIMEBLOCK(flexibleSolverCreate);
491 flexibleSolver_.create(getMatrix(),
492 isParallel(),
493 prm_,
494 pressureIndex,
495 trueFunc,
496 *comm_);
497 }
498 else
499 {
500 OPM_TIMEBLOCK(flexibleSolverUpdate);
501 flexibleSolver_.pre_->update();
502 }
503 }
504
505
509 {
510 // Decide if we should recreate the solver or just do
511 // a minimal preconditioner update.
512 if (!flexibleSolver_.solver_) {
513 return true;
514 }
515 if (this->parameters_.cpr_reuse_setup_ == 0) {
516 // Always recreate solver.
517 return true;
518 }
519 if (this->parameters_.cpr_reuse_setup_ == 1) {
520 // Recreate solver on the first iteration of every timestep.
521 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
522 return newton_iteration == 0;
523 }
524 if (this->parameters_.cpr_reuse_setup_ == 2) {
525 // Recreate solver if the last solve used more than 10 iterations.
526 return this->iterations() > 10;
527 }
528 if (this->parameters_.cpr_reuse_setup_ == 3) {
529 // Recreate solver if the last solve used more than 10 iterations.
530 return false;
531 }
532 if (this->parameters_.cpr_reuse_setup_ == 4) {
533 // Recreate solver every 'step' solve calls.
534 const int step = this->parameters_.cpr_reuse_interval_;
535 const bool create = ((calls_ % step) == 0);
536 return create;
537 }
538
539 // If here, we have an invalid parameter.
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.";
543 if (on_io_rank) {
544 OpmLog::error(msg);
545 }
546 throw std::runtime_error(msg);
547
548 // Never reached.
549 return false;
550 }
551
552
553 // Weights to make approximate pressure equations.
554 // Calculated from the storage terms (only) of the
555 // conservation equations, ignoring all other terms.
556 Vector getTrueImpesWeights(int pressureVarIndex) const
557 {
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());
565 return weights;
566 }
567
568 Matrix& getMatrix()
569 {
570 return *matrix_;
571 }
572
573 const Matrix& getMatrix() const
574 {
575 return *matrix_;
576 }
577
578 const Simulator& simulator_;
579 mutable int iterations_;
580 mutable int calls_;
581 mutable bool converged_;
582 std::any parallelInformation_;
583
584 // non-const to be able to scale the linear system
585 Matrix* matrix_;
586 Vector *rhs_;
587
588#if COMPILE_BDA_BRIDGE
589 std::unique_ptr<detail::BdaSolverInfo<Matrix, Vector>> bdaBridge;
590#endif
591
592 detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType> flexibleSolver_;
593 std::vector<int> overlapRows_;
594 std::vector<int> interiorRows_;
595
596 bool useWellConn_;
597
598 FlowLinearSolverParameters parameters_;
599 PropertyTree prm_;
600
601 std::shared_ptr< CommunicationType > comm_;
602 }; // end ISTLSolver
603
604} // namespace Opm
605#endif
Definition findOverlapRowsAndColumns.hpp:29
Interface class adding the update() method to the preconditioner interface.
Definition PreconditionerWithUpdate.hpp:32
Definition AquiferInterface.hpp:35
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 PropertyTree.hpp:37
This class serves to eliminate the need to include the WellContributions into the matrix (with –matri...
Definition WellContributions.hpp:52
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