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/CriticalError.hpp>
29#include <opm/common/ErrorMacros.hpp>
30#include <opm/common/Exceptions.hpp>
31#include <opm/common/TimingMacros.hpp>
32
37#include <opm/simulators/flow/BlackoilModelParameters.hpp>
40#include <opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp>
41#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
42#include <opm/simulators/linalg/matrixblock.hh>
44#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
45#include <opm/simulators/linalg/WellOperators.hpp>
46#include <opm/simulators/linalg/WriteSystemMatrixHelper.hpp>
47#include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
48#include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
49#include <opm/simulators/linalg/setupPropertyTree.hpp>
50
51#include <any>
52#include <cstddef>
53#include <functional>
54#include <memory>
55#include <set>
56#include <sstream>
57#include <string>
58#include <tuple>
59#include <vector>
60
61namespace Opm::Properties {
62
63namespace TTag {
65 using InheritsFrom = std::tuple<FlowIstlSolverParams>;
66};
67}
68
69template <class TypeTag, class MyTypeTag>
70struct WellModel;
71
74template<class TypeTag>
75struct SparseMatrixAdapter<TypeTag, TTag::FlowIstlSolver>
76{
77private:
81
82public:
83 using type = typename Linear::IstlSparseMatrixAdapter<Block>;
84};
85
86} // namespace Opm::Properties
87
88namespace Opm
89{
90
91
92namespace detail
93{
94
95template<class Matrix, class Vector, class Comm>
96struct FlexibleSolverInfo
97{
98 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
99 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
100 using AbstractPreconditionerType = Dune::PreconditionerWithUpdate<Vector, Vector>;
101
102 void create(const Matrix& matrix,
103 bool parallel,
104 const PropertyTree& prm,
105 std::size_t pressureIndex,
106 std::function<Vector()> weightCalculator,
107 const bool forceSerial,
108 Comm& comm);
109
110 std::unique_ptr<AbstractSolverType> solver_;
111 std::unique_ptr<AbstractOperatorType> op_;
112 std::unique_ptr<LinearOperatorExtra<Vector,Vector>> wellOperator_;
113 AbstractPreconditionerType* pre_ = nullptr;
114 std::size_t interiorCellNum_ = 0;
115};
116
117
118#ifdef HAVE_MPI
120void copyParValues(std::any& parallelInformation, std::size_t size,
121 Dune::OwnerOverlapCopyCommunication<int,int>& comm);
122#endif
123
126template<class Matrix>
127void makeOverlapRowsInvalid(Matrix& matrix,
128 const std::vector<int>& overlapRows);
129
132template<class Matrix, class Grid>
133std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
134 const std::vector<int>& cell_part,
135 std::size_t nonzeroes,
136 const std::vector<std::set<int>>& wellConnectionsGraph);
137}
138
143 template <class TypeTag>
145 {
146 protected:
154 using Matrix = typename SparseMatrixAdapter::IstlMatrix;
157 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
158 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
162 constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
163
164 enum { enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>() };
165 constexpr static bool isIncompatibleWithCprw = enablePolymerMolarWeight;
166
167#if HAVE_MPI
168 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
169#else
170 using CommunicationType = Dune::Communication<int>;
171#endif
172
173 public:
174 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
175
176 static void registerParameters()
177 {
178 FlowLinearSolverParameters::registerParameters();
179 }
180
188 ISTLSolver(const Simulator& simulator,
190 bool forceSerial = false)
191 : simulator_(simulator),
192 iterations_( 0 ),
193 converged_(false),
194 matrix_(nullptr),
195 parameters_{parameters},
196 forceSerial_(forceSerial)
197 {
198 initialize();
199 }
200
203 explicit ISTLSolver(const Simulator& simulator)
204 : simulator_(simulator),
205 iterations_( 0 ),
206 solveCount_(0),
207 converged_(false),
208 matrix_(nullptr)
209 {
210 parameters_.resize(1);
211 parameters_[0].init(simulator_.vanguard().eclState().getSimulationConfig().useCPR());
212 initialize();
213 }
214
215 void initialize()
216 {
218
219 if (isIncompatibleWithCprw) {
220 // Polymer injectivity is incompatible with the CPRW linear solver.
221 if (parameters_[0].linsolver_ == "cprw" || parameters_[0].linsolver_ == "hybrid") {
222 OPM_THROW(std::runtime_error,
223 "The polymer injectivity model is incompatible with the CPRW linear solver.\n"
224 "Choose a different option, for example --linear-solver=ilu0");
225 }
226 }
227
228 if (parameters_[0].linsolver_ == "hybrid") {
229 // Experimental hybrid configuration.
230 // When chosen, will set up two solvers, one with CPRW
231 // and the other with ILU0 preconditioner. More general
232 // options may be added later.
233 prm_.clear();
234 parameters_.clear();
235 {
236 FlowLinearSolverParameters para;
237 para.init(false);
238 para.linsolver_ = "cprw";
239 parameters_.push_back(para);
240 prm_.push_back(setupPropertyTree(parameters_[0],
241 Parameters::IsSet<Parameters::LinearSolverMaxIter>(),
242 Parameters::IsSet<Parameters::LinearSolverReduction>()));
243 }
244 {
245 FlowLinearSolverParameters para;
246 para.init(false);
247 para.linsolver_ = "ilu0";
248 parameters_.push_back(para);
249 prm_.push_back(setupPropertyTree(parameters_[1],
250 Parameters::IsSet<Parameters::LinearSolverMaxIter>(),
251 Parameters::IsSet<Parameters::LinearSolverReduction>()));
252 }
253 // ------------
254 } else {
255 assert(parameters_.size() == 1);
256 assert(prm_.empty());
257
258 // Do a normal linear solver setup.
259 if (parameters_[0].is_nldd_local_solver_) {
260 prm_.push_back(setupPropertyTree(parameters_[0],
261 Parameters::IsSet<Parameters::NlddLocalLinearSolverMaxIter>(),
262 Parameters::IsSet<Parameters::NlddLocalLinearSolverReduction>()));
263 }
264 else {
265 prm_.push_back(setupPropertyTree(parameters_[0],
266 Parameters::IsSet<Parameters::LinearSolverMaxIter>(),
267 Parameters::IsSet<Parameters::LinearSolverReduction>()));
268 }
269 }
270 flexibleSolver_.resize(prm_.size());
271
272 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
273#if HAVE_MPI
274 comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) );
275#endif
276 extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
277
278 // For some reason simulator_.model().elementMapper() is not initialized at this stage
279 //const auto& elemMapper = simulator_.model().elementMapper(); //does not work.
280 // Set it up manually
281 ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
282 detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
283 useWellConn_ = Parameters::Get<Parameters::MatrixAddWellContributions>();
284 const bool ownersFirst = Parameters::Get<Parameters::OwnerCellsFirst>();
285 if (!ownersFirst) {
286 const std::string msg = "The linear solver no longer supports --owner-cells-first=false.";
287 if (on_io_rank) {
288 OpmLog::error(msg);
289 }
290 OPM_THROW_NOLOG(std::runtime_error, msg);
291 }
292
293 const int interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true);
294 for (auto& f : flexibleSolver_) {
295 f.interiorCellNum_ = interiorCellNum_;
296 }
297
298#if HAVE_MPI
299 if (isParallel()) {
300 const std::size_t size = simulator_.vanguard().grid().leafGridView().size(0);
301 detail::copyParValues(parallelInformation_, size, *comm_);
302 }
303#endif
304
305 // Print parameters to PRT/DBG logs.
306 if (on_io_rank && parameters_[activeSolverNum_].linear_solver_print_json_definition_) {
307 std::ostringstream os;
308 os << "\nProperty tree for linear solvers:\n";
309 for (std::size_t i = 0; i<prm_.size(); i++) {
310 prm_[i].write_json(os, true);
311 }
312 OpmLog::note(os.str());
313 }
314 }
315
316 // nothing to clean here
317 void eraseMatrix()
318 {
319 }
320
321 void setActiveSolver(const int num)
322 {
323 if (num > static_cast<int>(prm_.size()) - 1) {
324 OPM_THROW(std::logic_error, "Solver number " + std::to_string(num) + " not available.");
325 }
326 activeSolverNum_ = num;
327 if (simulator_.gridView().comm().rank() == 0) {
328 OpmLog::debug("Active solver = " + std::to_string(activeSolverNum_)
329 + " (" + parameters_[activeSolverNum_].linsolver_ + ")");
330 }
331 }
332
333 int numAvailableSolvers()
334 {
335 return flexibleSolver_.size();
336 }
337
338 void initPrepare(const Matrix& M, Vector& b)
339 {
340 const bool firstcall = (matrix_ == nullptr);
341
342 // update matrix entries for solvers.
343 if (firstcall) {
344 // model will not change the matrix object. Hence simply store a pointer
345 // to the original one with a deleter that does nothing.
346 // Outch! We need to be able to scale the linear system! Hence const_cast
347 matrix_ = const_cast<Matrix*>(&M);
348
349 useWellConn_ = Parameters::Get<Parameters::MatrixAddWellContributions>();
350 // setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver)
351 } else {
352 // Pointers should not change
353 if ( &M != matrix_ ) {
354 OPM_THROW(std::logic_error,
355 "Matrix objects are expected to be reused when reassembling!");
356 }
357 }
358 rhs_ = &b;
359
360 // TODO: check all solvers, not just one.
361 if (isParallel() && prm_[activeSolverNum_].template get<std::string>("preconditioner.type") != "ParOverILU0") {
362 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
363 }
364 }
365
366 void prepare(const SparseMatrixAdapter& M, Vector& b)
367 {
368 prepare(M.istlMatrix(), b);
369 }
370
371 void prepare(const Matrix& M, Vector& b)
372 {
374 try {
375 initPrepare(M,b);
376
377 prepareFlexibleSolver();
378 } OPM_CATCH_AND_RETHROW_AS_CRITICAL_ERROR("This is likely due to a faulty linear solver JSON specification. Check for errors related to missing nodes.");
379 }
380
381
382 void setResidual(Vector& /* b */)
383 {
384 // rhs_ = &b; // Must be handled in prepare() instead.
385 }
386
387 void getResidual(Vector& b) const
388 {
389 b = *rhs_;
390 }
391
392 void setMatrix(const SparseMatrixAdapter& /* M */)
393 {
394 // matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead.
395 }
396
397 int getSolveCount() const {
398 return solveCount_;
399 }
400
401 void resetSolveCount() {
402 solveCount_ = 0;
403 }
404
405 bool solve(Vector& x)
406 {
408 ++solveCount_;
409 // Write linear system if asked for.
410 const int verbosity = prm_[activeSolverNum_].get("verbosity", 0);
411 const bool write_matrix = verbosity > 10;
412 if (write_matrix) {
413 Helper::writeSystem(simulator_, //simulator is only used to get names
414 getMatrix(),
415 *rhs_,
416 comm_.get());
417 }
418
419 // Solve system.
420 Dune::InverseOperatorResult result;
421 {
423 assert(flexibleSolver_[activeSolverNum_].solver_);
424 flexibleSolver_[activeSolverNum_].solver_->apply(x, *rhs_, result);
425 }
426
427 // Check convergence, iterations etc.
428 checkConvergence(result);
429
430 return converged_;
431 }
432
433
439
441 int iterations () const { return iterations_; }
442
444 const std::any& parallelInformation() const { return parallelInformation_; }
445
446 const CommunicationType* comm() const { return comm_.get(); }
447
448 void setDomainIndex(const int index)
449 {
450 domainIndex_ = index;
451 }
452
453 bool isNlddLocalSolver() const
454 {
455 return parameters_[activeSolverNum_].is_nldd_local_solver_;
456 }
457
458 protected:
459#if HAVE_MPI
460 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
461#endif
462
463 void checkConvergence( const Dune::InverseOperatorResult& result ) const
464 {
465 // store number of iterations
466 iterations_ = result.iterations;
467 converged_ = result.converged;
468 if(!converged_){
469 if(result.reduction < parameters_[activeSolverNum_].relaxed_linear_solver_reduction_){
470 std::stringstream ss;
471 ss<< "Full linear solver tolerance not achieved. The reduction is:" << result.reduction
472 << " after " << result.iterations << " iterations ";
473 OpmLog::warning(ss.str());
474 converged_ = true;
475 }
476 }
477 // Check for failure of linear solver.
478 if (!parameters_[activeSolverNum_].ignoreConvergenceFailure_ && !converged_) {
479 const std::string msg("Convergence failure for linear solver.");
481 }
482 }
483 protected:
484
485 bool isParallel() const {
486#if HAVE_MPI
487 return !forceSerial_ && comm_->communicator().size() > 1;
488#else
489 return false;
490#endif
491 }
492
493 void prepareFlexibleSolver()
494 {
496 if (shouldCreateSolver()) {
497 if (!useWellConn_) {
498 if (isNlddLocalSolver()) {
499 auto wellOp = std::make_unique<DomainWellModelAsLinearOperator<WellModel, Vector, Vector>>(simulator_.problem().wellModel());
500 wellOp->setDomainIndex(domainIndex_);
501 flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp);
502 }
503 else {
504 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
505 flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp);
506 }
507 }
508 std::function<Vector()> weightCalculator = this->getWeightsCalculator(prm_[activeSolverNum_], getMatrix(), pressureIndex);
510 flexibleSolver_[activeSolverNum_].create(getMatrix(),
511 isParallel(),
512 prm_[activeSolverNum_],
513 pressureIndex,
515 forceSerial_,
516 *comm_);
517 }
518 else
519 {
521 flexibleSolver_[activeSolverNum_].pre_->update();
522 }
523 }
524
525
529 {
530 // Decide if we should recreate the solver or just do
531 // a minimal preconditioner update.
532 if (flexibleSolver_.empty()) {
533 return true;
534 }
535 if (!flexibleSolver_[activeSolverNum_].solver_) {
536 return true;
537 }
538
539 if (flexibleSolver_[activeSolverNum_].pre_->hasPerfectUpdate()) {
540 return false;
541 }
542
543 // For AMG based preconditioners, the hierarchy depends on the matrix values
544 // so it is recreated at certain intervals
545 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 0) {
546 // Always recreate solver.
547 return true;
548 }
549 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 1) {
550 // Recreate solver on the first iteration of every timestep.
551 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
552 return newton_iteration == 0;
553 }
554 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 2) {
555 // Recreate solver if the last solve used more than 10 iterations.
556 return this->iterations() > 10;
557 }
558 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 3) {
559 // Never recreate the solver
560 return false;
561 }
562 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 4) {
563 // Recreate solver every 'step' solve calls.
564 const int step = this->parameters_[activeSolverNum_].cpr_reuse_interval_;
565 const bool create = ((solveCount_ % step) == 0);
566 return create;
567 }
568 // If here, we have an invalid parameter.
569 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
570 std::string msg = "Invalid value: " + std::to_string(this->parameters_[activeSolverNum_].cpr_reuse_setup_)
571 + " for --cpr-reuse-setup parameter, run with --help to see allowed values.";
572 if (on_io_rank) {
573 OpmLog::error(msg);
574 }
575 throw std::runtime_error(msg);
576
577 return false;
578 }
579
580
581 // Weights to make approximate pressure equations.
582 // Calculated from the storage terms (only) of the
583 // conservation equations, ignoring all other terms.
584 std::function<Vector()> getWeightsCalculator(const PropertyTree& prm,
585 const Matrix& matrix,
586 std::size_t pressIndex) const
587 {
588 std::function<Vector()> weightsCalculator;
589
590 using namespace std::string_literals;
591
592 auto preconditionerType = prm.get("preconditioner.type"s, "cpr"s);
593 if (preconditionerType == "cpr" || preconditionerType == "cprt"
594 || preconditionerType == "cprw" || preconditionerType == "cprwt") {
595 const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
596 const auto weightsType = prm.get("preconditioner.weight_type"s, "quasiimpes"s);
597 if (weightsType == "quasiimpes") {
598 // weights will be created as default in the solver
599 // assignment p = pressureIndex prevent compiler warning about
600 // capturing variable with non-automatic storage duration
601 weightsCalculator = [matrix, transpose, pressIndex]() {
602 return Amg::getQuasiImpesWeights<Matrix, Vector>(matrix,
604 transpose);
605 };
606 } else if ( weightsType == "trueimpes" ) {
608 [this, pressIndex]
609 {
610 Vector weights(rhs_->size());
611 ElementContext elemCtx(simulator_);
612 Amg::getTrueImpesWeights(pressIndex, weights,
613 simulator_.vanguard().gridView(),
614 elemCtx, simulator_.model(),
616 return weights;
617 };
618 } else if (weightsType == "trueimpesanalytic" ) {
620 [this, pressIndex]
621 {
622 Vector weights(rhs_->size());
623 ElementContext elemCtx(simulator_);
624 Amg::getTrueImpesWeightsAnalytic(pressIndex, weights,
625 simulator_.vanguard().gridView(),
626 elemCtx, simulator_.model(),
628 return weights;
629 };
630 } else {
631 OPM_THROW(std::invalid_argument,
632 "Weights type " + weightsType +
633 "not implemented for cpr."
634 " Please use quasiimpes, trueimpes or trueimpesanalytic.");
635 }
636 }
637 return weightsCalculator;
638 }
639
640
641 Matrix& getMatrix()
642 {
643 return *matrix_;
644 }
645
646 const Matrix& getMatrix() const
647 {
648 return *matrix_;
649 }
650
651 const Simulator& simulator_;
652 mutable int iterations_;
653 mutable int solveCount_;
654 mutable bool converged_;
655 std::any parallelInformation_;
656
657 // non-const to be able to scale the linear system
658 Matrix* matrix_;
659 Vector *rhs_;
660
661 int activeSolverNum_ = 0;
662 std::vector<detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType>> flexibleSolver_;
663 std::vector<int> overlapRows_;
664 std::vector<int> interiorRows_;
665
666 int domainIndex_ = -1;
667
668 bool useWellConn_;
669
670 std::vector<FlowLinearSolverParameters> parameters_;
671 bool forceSerial_ = false;
672 std::vector<PropertyTree> prm_;
673
674 std::shared_ptr< CommunicationType > comm_;
675 }; // end ISTLSolver
676
677} // namespace Opm
678
679#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:145
ISTLSolver(const Simulator &simulator, const FlowLinearSolverParameters &parameters, bool forceSerial=false)
Construct a system solver.
Definition ISTLSolver.hpp:188
int iterations() const
Solve the system of linear equations Ax = b, with A being the combined derivative matrix of the resid...
Definition ISTLSolver.hpp:441
bool shouldCreateSolver() const
Return true if we should (re)create the whole solver, instead of just calling update() on the precond...
Definition ISTLSolver.hpp:528
const std::any & parallelInformation() const
Definition ISTLSolver.hpp:444
ISTLSolver(const Simulator &simulator)
Construct a system solver.
Definition ISTLSolver.hpp:203
A sparse matrix interface backend for BCRSMatrix from dune-istl.
Definition istlsparsematrixadapter.hh:43
Definition matrixblock.hh:227
Hierarchical collection of key/value pairs.
Definition PropertyTree.hpp:39
T get(const std::string &key) const
Retrieve property value given hierarchical property key.
Definition PropertyTree.cpp:59
static unsigned threadId()
Return the index of the current OpenMP thread.
Definition threadmanager.cpp:84
Definition WellOperators.hpp:70
Declare the properties used by the infrastructure code of the finite volume discretizations.
A sparse matrix interface backend for BCRSMatrix from dune-istl.
Defines the common properties required by the porous medium multi-phase models.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
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 file provides the infrastructure to retrieve run-time parameters.
The Opm property system, traits with inheritance.
This class carries all parameters for the NewtonIterationBlackoilInterleaved class.
Definition FlowLinearSolverParameters.hpp:95
The class that allows to manipulate sparse matrices.
Definition linalgproperties.hh:50
Definition ISTLSolver.hpp:64
Definition FlowBaseProblemProperties.hpp:87