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#include <dune/istl/solver.hh>
27
28#include <ebos/eclbasevanguard.hh>
29
30#include <opm/common/ErrorMacros.hpp>
31#include <opm/common/Exceptions.hpp>
32#include <opm/common/TimingMacros.hpp>
33
34#include <opm/models/discretization/common/fvbaseproperties.hh>
35#include <opm/models/common/multiphasebaseproperties.hh>
36#include <opm/models/utils/parametersystem.hh>
37#include <opm/models/utils/propertysystem.hh>
38#include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
39#include <opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp>
40#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
41#include <opm/simulators/linalg/matrixblock.hh>
42#include <opm/simulators/linalg/istlsparsematrixadapter.hh>
43#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
44#include <opm/simulators/linalg/WellOperators.hpp>
45#include <opm/simulators/linalg/WriteSystemMatrixHelper.hpp>
46#include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
47#include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
48#include <opm/simulators/linalg/setupPropertyTree.hpp>
49
50#include <any>
51#include <cstddef>
52#include <functional>
53#include <memory>
54#include <set>
55#include <sstream>
56#include <string>
57#include <tuple>
58#include <vector>
59
60namespace Opm::Properties {
61
62namespace TTag {
64 using InheritsFrom = std::tuple<FlowIstlSolverParams>;
65};
66}
67
68template <class TypeTag, class MyTypeTag>
70
73template<class TypeTag>
74struct SparseMatrixAdapter<TypeTag, TTag::FlowIstlSolver>
75{
76private:
80
81public:
82 using type = typename Linear::IstlSparseMatrixAdapter<Block>;
83};
84
85} // namespace Opm::Properties
86
87namespace Opm
88{
89
90
91namespace detail
92{
93
94template<class Matrix, class Vector, class Comm>
95struct FlexibleSolverInfo
96{
97 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
98 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
99 using AbstractPreconditionerType = Dune::PreconditionerWithUpdate<Vector, Vector>;
100
101 void create(const Matrix& matrix,
102 bool parallel,
103 const PropertyTree& prm,
104 std::size_t pressureIndex,
105 std::function<Vector()> trueFunc,
106 const bool forceSerial,
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 std::size_t interiorCellNum_ = 0;
114};
115
116
117#ifdef HAVE_MPI
119void copyParValues(std::any& parallelInformation, std::size_t size,
120 Dune::OwnerOverlapCopyCommunication<int,int>& comm);
121#endif
122
125template<class Matrix>
126void makeOverlapRowsInvalid(Matrix& matrix,
127 const std::vector<int>& overlapRows);
128
131template<class Matrix, class Grid>
132std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
133 const std::vector<int>& cell_part,
134 std::size_t nonzeroes,
135 const std::vector<std::set<int>>& wellConnectionsGraph);
136}
137
142 template <class TypeTag>
144 {
145 protected:
153 using Matrix = typename SparseMatrixAdapter::IstlMatrix;
156 using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
157 using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
161 constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
162
163#if HAVE_MPI
164 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
165#else
166 using CommunicationType = Dune::CollectiveCommunication<int>;
167#endif
168
169 public:
170 using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
171
172 static void registerParameters()
173 {
174 FlowLinearSolverParameters::registerParameters<TypeTag>();
175 }
176
185 : simulator_(simulator),
186 iterations_( 0 ),
187 converged_(false),
188 matrix_(nullptr),
189 parameters_{parameters},
190 forceSerial_(forceSerial)
191 {
192 initialize();
193 }
194
197 explicit ISTLSolverEbos(const Simulator& simulator)
198 : simulator_(simulator),
199 iterations_( 0 ),
200 solveCount_(0),
201 converged_(false),
202 matrix_(nullptr)
203 {
204 parameters_.resize(1);
205 parameters_[0].template init<TypeTag>(simulator_.vanguard().eclState().getSimulationConfig().useCPR());
206 initialize();
207 }
208
209 void initialize()
210 {
212
213 if (parameters_[0].linsolver_ == "hybrid") {
214 // Experimental hybrid configuration.
215 // When chosen, will set up two solvers, one with CPRW
216 // and the other with ILU0 preconditioner. More general
217 // options may be added later.
218 prm_.clear();
219 parameters_.clear();
220 {
222 para.init<TypeTag>(false);
223 para.linsolver_ = "cprw";
224 parameters_.push_back(para);
225 prm_.push_back(setupPropertyTree(parameters_[0],
226 EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter),
227 EWOMS_PARAM_IS_SET(TypeTag, double, LinearSolverReduction)));
228 }
229 {
230 FlowLinearSolverParameters para;
231 para.init<TypeTag>(false);
232 para.linsolver_ = "ilu0";
233 parameters_.push_back(para);
234 prm_.push_back(setupPropertyTree(parameters_[1],
235 EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter),
236 EWOMS_PARAM_IS_SET(TypeTag, double, LinearSolverReduction)));
237 }
238 // ------------
239 } else {
240 // Do a normal linear solver setup.
241 assert(parameters_.size() == 1);
242 assert(prm_.empty());
243 prm_.push_back(setupPropertyTree(parameters_[0],
244 EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter),
245 EWOMS_PARAM_IS_SET(TypeTag, double, LinearSolverReduction)));
246 }
247 flexibleSolver_.resize(prm_.size());
248
249 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
250#if HAVE_MPI
251 comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) );
252#endif
253 extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
254
255 // For some reason simulator_.model().elementMapper() is not initialized at this stage
256 //const auto& elemMapper = simulator_.model().elementMapper(); //does not work.
257 // Set it up manually
258 ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
259 detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
260 useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
261 const bool ownersFirst = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst);
262 if (!ownersFirst) {
263 const std::string msg = "The linear solver no longer supports --owner-cells-first=false.";
264 if (on_io_rank) {
265 OpmLog::error(msg);
266 }
267 OPM_THROW_NOLOG(std::runtime_error, msg);
268 }
269
270 const int interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true);
271 for (auto& f : flexibleSolver_) {
272 f.interiorCellNum_ = interiorCellNum_;
273 }
274
275 // Print parameters to PRT/DBG logs.
276 if (on_io_rank && parameters_[activeSolverNum_].linear_solver_print_json_definition_) {
277 std::ostringstream os;
278 os << "Property tree for linear solvers:\n";
279 for (std::size_t i = 0; i<prm_.size(); i++) {
280 prm_[i].write_json(os, true);
281 std::cerr<< "debug: ["<<i<<"] : " << os.str() <<std::endl;
282 }
283 OpmLog::note(os.str());
284 }
285 }
286
287 // nothing to clean here
288 void eraseMatrix()
289 {
290 }
291
292 void setActiveSolver(const int num)
293 {
294 if (num > static_cast<int>(prm_.size()) - 1) {
295 OPM_THROW(std::logic_error, "Solver number " + std::to_string(num) + " not available.");
296 }
297 activeSolverNum_ = num;
298 auto cc = Dune::MPIHelper::getCollectiveCommunication();
299 if (cc.rank() == 0) {
300 OpmLog::debug("Active solver = " + std::to_string(activeSolverNum_)
301 + " (" + parameters_[activeSolverNum_].linsolver_ + ")");
302 }
303 }
304
305 int numAvailableSolvers()
306 {
307 return flexibleSolver_.size();
308 }
309
310 void prepare(const SparseMatrixAdapter& M, Vector& b)
311 {
312 prepare(M.istlMatrix(), b);
313 }
314
315 void prepare(const Matrix& M, Vector& b)
316 {
317 OPM_TIMEBLOCK(istlSolverEbosPrepare);
318 const bool firstcall = (matrix_ == nullptr);
319#if HAVE_MPI
320 if (firstcall && isParallel()) {
321 const std::size_t size = M.N();
322 detail::copyParValues(parallelInformation_, size, *comm_);
323 }
324#endif
325
326 // update matrix entries for solvers.
327 if (firstcall) {
328 // ebos will not change the matrix object. Hence simply store a pointer
329 // to the original one with a deleter that does nothing.
330 // Outch! We need to be able to scale the linear system! Hence const_cast
331 matrix_ = const_cast<Matrix*>(&M);
332
333 useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
334 // setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver)
335 } else {
336 // Pointers should not change
337 if ( &M != matrix_ ) {
338 OPM_THROW(std::logic_error,
339 "Matrix objects are expected to be reused when reassembling!");
340 }
341 }
342 rhs_ = &b;
343
344 // TODO: check all solvers, not just one.
345 if (isParallel() && prm_[activeSolverNum_].template get<std::string>("preconditioner.type") != "ParOverILU0") {
346 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
347 }
348 prepareFlexibleSolver();
349 }
350
351
352 void setResidual(Vector& /* b */)
353 {
354 // rhs_ = &b; // Must be handled in prepare() instead.
355 }
356
357 void getResidual(Vector& b) const
358 {
359 b = *rhs_;
360 }
361
362 void setMatrix(const SparseMatrixAdapter& /* M */)
363 {
364 // matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead.
365 }
366
367 int getSolveCount() const {
368 return solveCount_;
369 }
370
371 void resetSolveCount() {
372 solveCount_ = 0;
373 }
374
375 bool solve(Vector& x)
376 {
377 OPM_TIMEBLOCK(istlSolverEbosSolve);
378 ++solveCount_;
379 // Write linear system if asked for.
380 const int verbosity = prm_[activeSolverNum_].get("verbosity", 0);
381 const bool write_matrix = verbosity > 10;
382 if (write_matrix) {
383 Helper::writeSystem(simulator_, //simulator is only used to get names
384 getMatrix(),
385 *rhs_,
386 comm_.get());
387 }
388
389 // Solve system.
390 Dune::InverseOperatorResult result;
391 {
392 OPM_TIMEBLOCK(flexibleSolverApply);
393 assert(flexibleSolver_[activeSolverNum_].solver_);
394 flexibleSolver_[activeSolverNum_].solver_->apply(x, *rhs_, result);
395 }
396
397 // Check convergence, iterations etc.
398 checkConvergence(result);
399
400 return converged_;
401 }
402
403
409
411 int iterations () const { return iterations_; }
412
414 const std::any& parallelInformation() const { return parallelInformation_; }
415
416 const CommunicationType* comm() const { return comm_.get(); }
417
418 protected:
419#if HAVE_MPI
420 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
421#endif
422
423 void checkConvergence( const Dune::InverseOperatorResult& result ) const
424 {
425 // store number of iterations
426 iterations_ = result.iterations;
427 converged_ = result.converged;
428 if(!converged_){
429 if(result.reduction < parameters_[activeSolverNum_].relaxed_linear_solver_reduction_){
430 std::stringstream ss;
431 ss<< "Full linear solver tolerance not achieved. The reduction is:" << result.reduction
432 << " after " << result.iterations << " iterations ";
433 OpmLog::warning(ss.str());
434 converged_ = true;
435 }
436 }
437 // Check for failure of linear solver.
438 if (!parameters_[activeSolverNum_].ignoreConvergenceFailure_ && !converged_) {
439 const std::string msg("Convergence failure for linear solver.");
440 OPM_THROW_NOLOG(NumericalProblem, msg);
441 }
442 }
443 protected:
444
445 bool isParallel() const {
446#if HAVE_MPI
447 return !forceSerial_ && comm_->communicator().size() > 1;
448#else
449 return false;
450#endif
451 }
452
453 void prepareFlexibleSolver()
454 {
455 OPM_TIMEBLOCK(flexibleSolverPrepare);
456 if (shouldCreateSolver()) {
457 std::function<Vector()> trueFunc =
458 [this]
459 {
460 return this->getTrueImpesWeights(pressureIndex);
461 };
462
463 if (!useWellConn_) {
464 auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
465 flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp);
466 }
467 OPM_TIMEBLOCK(flexibleSolverCreate);
468 flexibleSolver_[activeSolverNum_].create(getMatrix(),
469 isParallel(),
470 prm_[activeSolverNum_],
471 pressureIndex,
472 trueFunc,
473 forceSerial_,
474 *comm_);
475 }
476 else
477 {
478 OPM_TIMEBLOCK(flexibleSolverUpdate);
479 flexibleSolver_[activeSolverNum_].pre_->update();
480 }
481 }
482
483
487 {
488 // Decide if we should recreate the solver or just do
489 // a minimal preconditioner update.
490 if (flexibleSolver_.empty()) {
491 return true;
492 }
493 if (!flexibleSolver_[activeSolverNum_].solver_) {
494 return true;
495 }
496 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 0) {
497 // Always recreate solver.
498 return true;
499 }
500 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 1) {
501 // Recreate solver on the first iteration of every timestep.
502 const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
503 return newton_iteration == 0;
504 }
505 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 2) {
506 // Recreate solver if the last solve used more than 10 iterations.
507 return this->iterations() > 10;
508 }
509 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 3) {
510 // Recreate solver if the last solve used more than 10 iterations.
511 return false;
512 }
513 if (this->parameters_[activeSolverNum_].cpr_reuse_setup_ == 4) {
514 // Recreate solver every 'step' solve calls.
515 const int step = this->parameters_[activeSolverNum_].cpr_reuse_interval_;
516 const bool create = ((solveCount_ % step) == 0);
517 return create;
518 }
519
520 // If here, we have an invalid parameter.
521 const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
522 std::string msg = "Invalid value: " + std::to_string(this->parameters_[activeSolverNum_].cpr_reuse_setup_)
523 + " for --cpr-reuse-setup parameter, run with --help to see allowed values.";
524 if (on_io_rank) {
525 OpmLog::error(msg);
526 }
527 throw std::runtime_error(msg);
528
529 // Never reached.
530 return false;
531 }
532
533
534 // Weights to make approximate pressure equations.
535 // Calculated from the storage terms (only) of the
536 // conservation equations, ignoring all other terms.
537 Vector getTrueImpesWeights(int pressureVarIndex) const
538 {
539 OPM_TIMEBLOCK(getTrueImpesWeights);
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 solveCount_;
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 int activeSolverNum_ = 0;
570 std::vector<detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType>> flexibleSolver_;
571 std::vector<int> overlapRows_;
572 std::vector<int> interiorRows_;
573
574 bool useWellConn_;
575
576 std::vector<FlowLinearSolverParameters> parameters_;
577 bool forceSerial_ = false;
578 std::vector<PropertyTree> prm_;
579
580 std::shared_ptr< CommunicationType > comm_;
581 }; // end ISTLSolver
582
583} // namespace Opm
584#endif
Interface class adding the update() method to the preconditioner interface.
Definition PreconditionerWithUpdate.hpp:32
Definition AquiferInterface.hpp:35
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition ISTLSolverEbos.hpp:144
ISTLSolverEbos(const Simulator &simulator, const FlowLinearSolverParameters &parameters, bool forceSerial=false)
Construct a system solver.
Definition ISTLSolverEbos.hpp:184
int iterations() const
Solve the system of linear equations Ax = b, with A being the combined derivative matrix of the resid...
Definition ISTLSolverEbos.hpp:411
const std::any & parallelInformation() const
Definition ISTLSolverEbos.hpp:414
ISTLSolverEbos(const Simulator &simulator)
Construct a system solver.
Definition ISTLSolverEbos.hpp:197
bool shouldCreateSolver() const
Return true if we should (re)create the whole solver, instead of just calling update() on the precond...
Definition ISTLSolverEbos.hpp:486
Definition PropertyTree.hpp:37
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 ISTLSolverEbos.hpp:69
Definition ISTLSolverEbos.hpp:63