My Project
Loading...
Searching...
No Matches
ISTLSolver.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_HEADER_INCLUDED
23#define OPM_ISTLSOLVER_HEADER_INCLUDED
24
25#include <dune/istl/owneroverlapcopy.hh>
26#include <dune/istl/solver.hh>
27
28#include <opm/common/ErrorMacros.hpp>
29#include <opm/common/Exceptions.hpp>
30#include <opm/common/TimingMacros.hpp>
31
32#include <opm/models/discretization/common/fvbaseproperties.hh>
33#include <opm/models/common/multiphasebaseproperties.hh>
34#include <opm/models/utils/parametersystem.hh>
35#include <opm/models/utils/propertysystem.hh>
36#include <opm/simulators/flow/BlackoilModelParameters.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>
68struct WellModel;
69
72template<class TypeTag>
73struct SparseMatrixAdapter<TypeTag, TTag::FlowIstlSolver>
74{
75private:
76 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
77 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
79
80public:
81 using type = typename Linear::IstlSparseMatrixAdapter<Block>;
82};
83
84} // namespace Opm::Properties
85
86namespace Opm
87{
88
89
90namespace detail
91{
92
93template<class Matrix, class Vector, class Comm>
94struct FlexibleSolverInfo
95{
96 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
97 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
98 using AbstractPreconditionerType = Dune::PreconditionerWithUpdate<Vector, Vector>;
99
100 void create(const Matrix& matrix,
101 bool parallel,
102 const PropertyTree& prm,
103 std::size_t pressureIndex,
104 std::function<Vector()> weightCalculator,
105 const bool forceSerial,
106 Comm& comm);
107
108 std::unique_ptr<AbstractSolverType> solver_;
109 std::unique_ptr<AbstractOperatorType> op_;
110 std::unique_ptr<LinearOperatorExtra<Vector,Vector>> wellOperator_;
111 AbstractPreconditionerType* pre_ = nullptr;
112 std::size_t interiorCellNum_ = 0;
113};
114
115
116#ifdef HAVE_MPI
118void copyParValues(std::any& parallelInformation, std::size_t size,
119 Dune::OwnerOverlapCopyCommunication<int,int>& comm);
120#endif
121
124template<class Matrix>
125void makeOverlapRowsInvalid(Matrix& matrix,
126 const std::vector<int>& overlapRows);
127
130template<class Matrix, class Grid>
131std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
132 const std::vector<int>& cell_part,
133 std::size_t nonzeroes,
134 const std::vector<std::set<int>>& wellConnectionsGraph);
135}
136
141 template <class TypeTag>
143 {
144 protected:
145 using GridView = GetPropType<TypeTag, Properties::GridView>;
146 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
147 using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
148 using Vector = GetPropType<TypeTag, Properties::GlobalEqVector>;
149 using Indices = GetPropType<TypeTag, Properties::Indices>;
150 using WellModel = GetPropType<TypeTag, Properties::WellModel>;
151 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
152 using Matrix = typename SparseMatrixAdapter::IstlMatrix;
153 using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
154 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
155 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
156 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
159 using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
160 constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
161
162#if HAVE_MPI
163 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
164#else
165 using CommunicationType = Dune::CollectiveCommunication<int>;
166#endif
167
168 public:
169 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
170
171 static void registerParameters()
172 {
173 FlowLinearSolverParameters::registerParameters<TypeTag>();
174 }
175
183 ISTLSolver(const Simulator& simulator,
184 const FlowLinearSolverParameters& parameters,
185 bool forceSerial = false)
186 : simulator_(simulator),
187 iterations_( 0 ),
188 converged_(false),
189 matrix_(nullptr),
190 parameters_{parameters},
191 forceSerial_(forceSerial)
192 {
193 initialize();
194 }
195
198 explicit ISTLSolver(const Simulator& simulator)
199 : simulator_(simulator),
200 iterations_( 0 ),
201 solveCount_(0),
202 converged_(false),
203 matrix_(nullptr)
204 {
205 parameters_.resize(1);
206 parameters_[0].template init<TypeTag>(simulator_.vanguard().eclState().getSimulationConfig().useCPR());
207 initialize();
208 }
209
210 void initialize()
211 {
212 OPM_TIMEBLOCK(IstlSolver);
213
214 if (parameters_[0].linsolver_ == "hybrid") {
215 // Experimental hybrid configuration.
216 // When chosen, will set up two solvers, one with CPRW
217 // and the other with ILU0 preconditioner. More general
218 // options may be added later.
219 prm_.clear();
220 parameters_.clear();
221 {
223 para.init<TypeTag>(false);
224 para.linsolver_ = "cprw";
225 parameters_.push_back(para);
226 prm_.push_back(setupPropertyTree(parameters_[0],
227 Parameters::isSet<TypeTag,Properties::LinearSolverMaxIter>(),
228 Parameters::isSet<TypeTag,Properties::LinearSolverReduction>()));
229 }
230 {
231 FlowLinearSolverParameters para;
232 para.init<TypeTag>(false);
233 para.linsolver_ = "ilu0";
234 parameters_.push_back(para);
235 prm_.push_back(setupPropertyTree(parameters_[1],
236 Parameters::isSet<TypeTag,Properties::LinearSolverMaxIter>(),
237 Parameters::isSet<TypeTag,Properties::LinearSolverReduction>()));
238 }
239 // ------------
240 } else {
241 // Do a normal linear solver setup.
242 assert(parameters_.size() == 1);
243 assert(prm_.empty());
244 prm_.push_back(setupPropertyTree(parameters_[0],
245 Parameters::isSet<TypeTag,Properties::LinearSolverMaxIter>(),
246 Parameters::isSet<TypeTag,Properties::LinearSolverReduction>()));
247 }
248 flexibleSolver_.resize(prm_.size());
249
250 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
251#if HAVE_MPI
252 comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) );
253#endif
254 extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
255
256 // For some reason simulator_.model().elementMapper() is not initialized at this stage
257 //const auto& elemMapper = simulator_.model().elementMapper(); //does not work.
258 // Set it up manually
259 ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
260 detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
261 useWellConn_ = Parameters::get<TypeTag, Properties::MatrixAddWellContributions>();
262 const bool ownersFirst = Parameters::get<TypeTag, Properties::OwnerCellsFirst>();
263 if (!ownersFirst) {
264 const std::string msg = "The linear solver no longer supports --owner-cells-first=false.";
265 if (on_io_rank) {
266 OpmLog::error(msg);
267 }
268 OPM_THROW_NOLOG(std::runtime_error, msg);
269 }
270
271 const int interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true);
272 for (auto& f : flexibleSolver_) {
273 f.interiorCellNum_ = interiorCellNum_;
274 }
275
276#if HAVE_MPI
277 if (isParallel()) {
278 const std::size_t size = simulator_.vanguard().grid().leafGridView().size(0);
279 detail::copyParValues(parallelInformation_, size, *comm_);
280 }
281#endif
282
283 // Print parameters to PRT/DBG logs.
284 if (on_io_rank && parameters_[activeSolverNum_].linear_solver_print_json_definition_) {
285 std::ostringstream os;
286 os << "Property tree for linear solvers:\n";
287 for (std::size_t i = 0; i<prm_.size(); i++) {
288 prm_[i].write_json(os, true);
289 }
290 OpmLog::note(os.str());
291 }
292 }
293
294 // nothing to clean here
295 void eraseMatrix()
296 {
297 }
298
299 void setActiveSolver(const int num)
300 {
301 if (num > static_cast<int>(prm_.size()) - 1) {
302 OPM_THROW(std::logic_error, "Solver number " + std::to_string(num) + " not available.");
303 }
304 activeSolverNum_ = num;
305 if (simulator_.gridView().comm().rank() == 0) {
306 OpmLog::debug("Active solver = " + std::to_string(activeSolverNum_)
307 + " (" + parameters_[activeSolverNum_].linsolver_ + ")");
308 }
309 }
310
311 int numAvailableSolvers()
312 {
313 return flexibleSolver_.size();
314 }
315
316 void initPrepare(const Matrix& M, Vector& b)
317 {
318 const bool firstcall = (matrix_ == nullptr);
319
320 // update matrix entries for solvers.
321 if (firstcall) {
322 // model will not change the matrix object. Hence simply store a pointer
323 // to the original one with a deleter that does nothing.
324 // Outch! We need to be able to scale the linear system! Hence const_cast
325 matrix_ = const_cast<Matrix*>(&M);
326
327 useWellConn_ = Parameters::get<TypeTag, Properties::MatrixAddWellContributions>();
328 // setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver)
329 } else {
330 // Pointers should not change
331 if ( &M != matrix_ ) {
332 OPM_THROW(std::logic_error,
333 "Matrix objects are expected to be reused when reassembling!");
334 }
335 }
336 rhs_ = &b;
337
338 // TODO: check all solvers, not just one.
339 if (isParallel() && prm_[activeSolverNum_].template get<std::string>("preconditioner.type") != "ParOverILU0") {
340 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
341 }
342 }
343
344 void prepare(const SparseMatrixAdapter& M, Vector& b)
345 {
346 prepare(M.istlMatrix(), b);
347 }
348
349 void prepare(const Matrix& M, Vector& b)
350 {
351 OPM_TIMEBLOCK(istlSolverPrepare);
352
353 initPrepare(M,b);
354
355 prepareFlexibleSolver();
356 }
357
358
359 void setResidual(Vector& /* b */)
360 {
361 // rhs_ = &b; // Must be handled in prepare() instead.
362 }
363
364 void getResidual(Vector& b) const
365 {
366 b = *rhs_;
367 }
368
369 void setMatrix(const SparseMatrixAdapter& /* M */)
370 {
371 // matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead.
372 }
373
374 int getSolveCount() const {
375 return solveCount_;
376 }
377
378 void resetSolveCount() {
379 solveCount_ = 0;
380 }
381
382 bool solve(Vector& x)
383 {
384 OPM_TIMEBLOCK(istlSolverSolve);
385 ++solveCount_;
386 // Write linear system if asked for.
387 const int verbosity = prm_[activeSolverNum_].get("verbosity", 0);
388 const bool write_matrix = verbosity > 10;
389 if (write_matrix) {
390 Helper::writeSystem(simulator_, //simulator is only used to get names
391 getMatrix(),
392 *rhs_,
393 comm_.get());
394 }
395
396 // Solve system.
397 Dune::InverseOperatorResult result;
398 {
399 OPM_TIMEBLOCK(flexibleSolverApply);
400 assert(flexibleSolver_[activeSolverNum_].solver_);
401 flexibleSolver_[activeSolverNum_].solver_->apply(x, *rhs_, result);
402 }
403
404 // Check convergence, iterations etc.
405 checkConvergence(result);
406
407 return converged_;
408 }
409
410
416
418 int iterations () const { return iterations_; }
419
421 const std::any& parallelInformation() const { return parallelInformation_; }
422
423 const CommunicationType* comm() const { return comm_.get(); }
424
425 protected:
426#if HAVE_MPI
427 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
428#endif
429
430 void checkConvergence( const Dune::InverseOperatorResult& result ) const
431 {
432 // store number of iterations
433 iterations_ = result.iterations;
434 converged_ = result.converged;
435 if(!converged_){
436 if(result.reduction < parameters_[activeSolverNum_].relaxed_linear_solver_reduction_){
437 std::stringstream ss;
438 ss<< "Full linear solver tolerance not achieved. The reduction is:" << result.reduction
439 << " after " << result.iterations << " iterations ";
440 OpmLog::warning(ss.str());
441 converged_ = true;
442 }
443 }
444 // Check for failure of linear solver.
445 if (!parameters_[activeSolverNum_].ignoreConvergenceFailure_ && !converged_) {
446 const std::string msg("Convergence failure for linear solver.");
447 OPM_THROW_NOLOG(NumericalProblem, msg);
448 }
449 }
450 protected:
451
452 bool isParallel() const {
453#if HAVE_MPI
454 return !forceSerial_ && comm_->communicator().size() > 1;
455#else
456 return false;
457#endif
458 }
459
460 void prepareFlexibleSolver()
461 {
462 OPM_TIMEBLOCK(flexibleSolverPrepare);
463 if (shouldCreateSolver()) {
464 if (!useWellConn_) {
465 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
466 flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp);
467 }
468 std::function<Vector()> weightCalculator = this->getWeightsCalculator(prm_[activeSolverNum_], getMatrix(), pressureIndex);
469 OPM_TIMEBLOCK(flexibleSolverCreate);
470 flexibleSolver_[activeSolverNum_].create(getMatrix(),
471 isParallel(),
472 prm_[activeSolverNum_],
473 pressureIndex,
474 weightCalculator,
475 forceSerial_,
476 *comm_);
477 }
478 else
479 {
480 OPM_TIMEBLOCK(flexibleSolverUpdate);
481 flexibleSolver_[activeSolverNum_].pre_->update();
482 }
483 }
484
485
489 {
490 // Decide if we should recreate the solver or just do
491 // a minimal preconditioner update.
492 if (flexibleSolver_.empty()) {
493 return true;
494 }
495 if (!flexibleSolver_[activeSolverNum_].solver_) {
496 return true;
497 }
498 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 0) {
499 // Always recreate solver.
500 return true;
501 }
502 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 1) {
503 // Recreate solver on the first iteration of every timestep.
504 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
505 return newton_iteration == 0;
506 }
507 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 2) {
508 // Recreate solver if the last solve used more than 10 iterations.
509 return this->iterations() > 10;
510 }
511 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 3) {
512 // Recreate solver if the last solve used more than 10 iterations.
513 return false;
514 }
515 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 4) {
516 // Recreate solver every 'step' solve calls.
517 const int step = this->parameters_[activeSolverNum_].cpr_reuse_interval_;
518 const bool create = ((solveCount_ % step) == 0);
519 return create;
520 }
521
522 // If here, we have an invalid parameter.
523 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
524 std::string msg = "Invalid value: " + std::to_string(this->parameters_[activeSolverNum_].cpr_reuse_setup_)
525 + " for --cpr-reuse-setup parameter, run with --help to see allowed values.";
526 if (on_io_rank) {
527 OpmLog::error(msg);
528 }
529 throw std::runtime_error(msg);
530
531 // Never reached.
532 return false;
533 }
534
535
536 // Weights to make approximate pressure equations.
537 // Calculated from the storage terms (only) of the
538 // conservation equations, ignoring all other terms.
539 std::function<Vector()> getWeightsCalculator(const PropertyTree& prm,
540 const Matrix& matrix,
541 std::size_t pressIndex) const
542 {
543 std::function<Vector()> weightsCalculator;
544
545 using namespace std::string_literals;
546
547 auto preconditionerType = prm.get("preconditioner.type"s, "cpr"s);
548 if (preconditionerType == "cpr" || preconditionerType == "cprt"
549 || preconditionerType == "cprw" || preconditionerType == "cprwt") {
550 const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
551 const auto weightsType = prm.get("preconditioner.weight_type"s, "quasiimpes"s);
552 if (weightsType == "quasiimpes") {
553 // weights will be created as default in the solver
554 // assignment p = pressureIndex prevent compiler warning about
555 // capturing variable with non-automatic storage duration
556 weightsCalculator = [matrix, transpose, pressIndex]() {
557 return Amg::getQuasiImpesWeights<Matrix, Vector>(matrix,
558 pressIndex,
559 transpose);
560 };
561 } else if ( weightsType == "trueimpes" ) {
562 weightsCalculator =
563 [this, pressIndex]
564 {
565 Vector weights(rhs_->size());
566 ElementContext elemCtx(simulator_);
567 Amg::getTrueImpesWeights(pressIndex, weights,
568 simulator_.vanguard().gridView(),
569 elemCtx, simulator_.model(),
570 ThreadManager::threadId());
571 return weights;
572 };
573 } else if (weightsType == "trueimpesanalytic" ) {
574 weightsCalculator =
575 [this, pressIndex]
576 {
577 Vector weights(rhs_->size());
578 ElementContext elemCtx(simulator_);
579 Amg::getTrueImpesWeightsAnalytic(pressIndex, weights,
580 simulator_.vanguard().gridView(),
581 elemCtx, simulator_.model(),
582 ThreadManager::threadId());
583 return weights;
584 };
585 } else {
586 OPM_THROW(std::invalid_argument,
587 "Weights type " + weightsType +
588 "not implemented for cpr."
589 " Please use quasiimpes, trueimpes or trueimpesanalytic.");
590 }
591 }
592 return weightsCalculator;
593 }
594
595
596 Matrix& getMatrix()
597 {
598 return *matrix_;
599 }
600
601 const Matrix& getMatrix() const
602 {
603 return *matrix_;
604 }
605
606 const Simulator& simulator_;
607 mutable int iterations_;
608 mutable int solveCount_;
609 mutable bool converged_;
610 std::any parallelInformation_;
611
612 // non-const to be able to scale the linear system
613 Matrix* matrix_;
614 Vector *rhs_;
615
616 int activeSolverNum_ = 0;
617 std::vector<detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType>> flexibleSolver_;
618 std::vector<int> overlapRows_;
619 std::vector<int> interiorRows_;
620
621 bool useWellConn_;
622
623 std::vector<FlowLinearSolverParameters> parameters_;
624 bool forceSerial_ = false;
625 std::vector<PropertyTree> prm_;
626
627 std::shared_ptr< CommunicationType > comm_;
628 }; // end ISTLSolver
629
630} // namespace Opm
631
632#endif // OPM_ISTLSOLVER_HEADER_INCLUDED
Helper class for grid instantiation of ECL file-format using problems.
Interface class adding the update() method to the preconditioner interface.
Definition PreconditionerWithUpdate.hpp:32
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition ISTLSolver.hpp:143
ISTLSolver(const Simulator &simulator, const FlowLinearSolverParameters &parameters, bool forceSerial=false)
Construct a system solver.
Definition ISTLSolver.hpp:183
int iterations() const
Solve the system of linear equations Ax = b, with A being the combined derivative matrix of the resid...
Definition ISTLSolver.hpp:418
bool shouldCreateSolver() const
Return true if we should (re)create the whole solver, instead of just calling update() on the precond...
Definition ISTLSolver.hpp:488
const std::any & parallelInformation() const
Definition ISTLSolver.hpp:421
ISTLSolver(const Simulator &simulator)
Construct a system solver.
Definition ISTLSolver.hpp:198
Definition MatrixMarketSpecializations.hpp:17
Definition PropertyTree.hpp:37
Definition WellOperators.hpp:67
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
This class carries all parameters for the NewtonIterationBlackoilInterleaved class.
Definition FlowLinearSolverParameters.hpp:237
Definition ISTLSolver.hpp:62
Definition FlowProblemProperties.hpp:70