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 <opm/models/utils/parametersystem.hh>
26 #include <opm/models/utils/propertysystem.hh>
27 #include <opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp>
28 #include <opm/simulators/linalg/FlexibleSolver.hpp>
29 #include <opm/simulators/linalg/MatrixBlock.hpp>
30 #include <opm/simulators/linalg/ParallelIstlInformation.hpp>
31 #include <opm/simulators/linalg/WellOperators.hpp>
32 #include <opm/simulators/linalg/WriteSystemMatrixHelper.hpp>
33 #include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
34 #include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
35 #include <opm/simulators/linalg/setupPropertyTree.hpp>
36 
37 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
38 #include <opm/simulators/linalg/bda/BdaBridge.hpp>
39 #include <opm/simulators/linalg/bda/WellContributions.hpp>
40 #endif
41 
42 namespace Opm::Properties {
43 
44 namespace TTag {
46  using InheritsFrom = std::tuple<FlowIstlSolverParams>;
47 };
48 }
49 
50 template <class TypeTag, class MyTypeTag>
51 struct EclWellModel;
52 
55 template<class TypeTag>
56 struct SparseMatrixAdapter<TypeTag, TTag::FlowIstlSolver>
57 {
58 private:
59  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
60  enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
62 
63 public:
64  typedef typename Linear::IstlSparseMatrixAdapter<Block> type;
65 };
66 
67 } // namespace Opm::Properties
68 
69 namespace Opm
70 {
71 
76  template <class TypeTag>
78  {
79  protected:
80  using GridView = GetPropType<TypeTag, Properties::GridView>;
81  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
82  using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
83  using Vector = GetPropType<TypeTag, Properties::GlobalEqVector>;
84  using Indices = GetPropType<TypeTag, Properties::Indices>;
85  using WellModel = GetPropType<TypeTag, Properties::EclWellModel>;
86  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
87  using Matrix = typename SparseMatrixAdapter::IstlMatrix;
88  using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
89  using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
91  using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
93  using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
94  constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
95 
96 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
97  static const unsigned int block_size = Matrix::block_type::rows;
98  std::unique_ptr<BdaBridge<Matrix, Vector, block_size>> bdaBridge;
99 #endif
100 
101 #if HAVE_MPI
102  using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
103 #else
104  using CommunicationType = Dune::CollectiveCommunication< int >;
105 #endif
106 
107  public:
108  using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
109 
110  static void registerParameters()
111  {
112  FlowLinearSolverParameters::registerParameters<TypeTag>();
113  }
114 
118  explicit ISTLSolverEbos(const Simulator& simulator)
119  : simulator_(simulator),
120  iterations_( 0 ),
121  converged_(false),
122  matrix_()
123  {
124  const bool on_io_rank = (simulator.gridView().comm().rank() == 0);
125 #if HAVE_MPI
126  comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) );
127 #endif
128  parameters_.template init<TypeTag>();
129  prm_ = setupPropertyTree(parameters_,
130  EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter),
131  EWOMS_PARAM_IS_SET(TypeTag, int, CprMaxEllIter));
132 
133 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
134  {
135  std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
136  if ((simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode != "none")) {
137  if (on_io_rank) {
138  OpmLog::warning("Cannot use GPU or FPGA with MPI, GPU/FPGA are disabled");
139  }
140  accelerator_mode = "none";
141  }
142  const int platformID = EWOMS_GET_PARAM(TypeTag, int, OpenclPlatformId);
143  const int deviceID = EWOMS_GET_PARAM(TypeTag, int, BdaDeviceId);
144  const int maxit = EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIter);
145  const double tolerance = EWOMS_GET_PARAM(TypeTag, double, LinearSolverReduction);
146  const std::string opencl_ilu_reorder = EWOMS_GET_PARAM(TypeTag, std::string, OpenclIluReorder);
147  const int linear_solver_verbosity = parameters_.linear_solver_verbosity_;
148  std::string fpga_bitstream = EWOMS_GET_PARAM(TypeTag, std::string, FpgaBitstream);
149  std::string linsolver = EWOMS_GET_PARAM(TypeTag, std::string, Linsolver);
150  bdaBridge.reset(new BdaBridge<Matrix, Vector, block_size>(accelerator_mode, fpga_bitstream, linear_solver_verbosity, maxit, tolerance, platformID, deviceID, opencl_ilu_reorder, linsolver));
151  }
152 #else
153  if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) != "none") {
154  OPM_THROW(std::logic_error,"Cannot use accelerated solver since CUDA, OpenCL and amgcl were not found by cmake and FPGA was not enabled");
155  }
156 #endif
157  extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
158 
159  // For some reason simulator_.model().elementMapper() is not initialized at this stage
160  // Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work.
161  // Set it up manually
162  ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
163  detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
164 
165  useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
166 #if HAVE_FPGA
167  // check usage of MatrixAddWellContributions: for FPGA they must be included
168  if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) == "fpga" && !useWellConn_) {
169  OPM_THROW(std::logic_error,"fpgaSolver needs --matrix-add-well-contributions=true");
170  }
171 #endif
172  const bool ownersFirst = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst);
173  if (!ownersFirst) {
174  const std::string msg = "The linear solver no longer supports --owner-cells-first=false.";
175  if (on_io_rank) {
176  OpmLog::error(msg);
177  }
178  OPM_THROW_NOLOG(std::runtime_error, msg);
179  }
180 
181  interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true);
182 
183  // Print parameters to PRT/DBG logs.
184  if (on_io_rank) {
185  std::ostringstream os;
186  os << "Property tree for linear solver:\n";
187  prm_.write_json(os, true);
188  OpmLog::note(os.str());
189  }
190  }
191 
192  // nothing to clean here
193  void eraseMatrix() {
194  }
195 
196  void prepare(const SparseMatrixAdapter& M, Vector& b)
197  {
198  static bool firstcall = true;
199 #if HAVE_MPI
200  if (firstcall && parallelInformation_.type() == typeid(ParallelISTLInformation)) {
201  // Parallel case.
202  const ParallelISTLInformation* parinfo = std::any_cast<ParallelISTLInformation>(&parallelInformation_);
203  assert(parinfo);
204  const size_t size = M.istlMatrix().N();
205  parinfo->copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), size, 1);
206  }
207 #endif
208 
209  // update matrix entries for solvers.
210  if (firstcall) {
211  // ebos will not change the matrix object. Hence simply store a pointer
212  // to the original one with a deleter that does nothing.
213  // Outch! We need to be able to scale the linear system! Hence const_cast
214  matrix_ = const_cast<Matrix*>(&M.istlMatrix());
215  } else {
216  // Pointers should not change
217  if ( &(M.istlMatrix()) != matrix_ ) {
218  OPM_THROW(std::logic_error, "Matrix objects are expected to be reused when reassembling!"
219  <<" old pointer was " << matrix_ << ", new one is " << (&M.istlMatrix()) );
220  }
221  }
222  rhs_ = &b;
223 
224  if (isParallel() && prm_.get<std::string>("preconditioner.type") != "ParOverILU0") {
225  makeOverlapRowsInvalid(getMatrix());
226  }
227  prepareFlexibleSolver();
228  firstcall = false;
229  }
230 
231 
232  void setResidual(Vector& /* b */) {
233  // rhs_ = &b; // Must be handled in prepare() instead.
234  }
235 
236  void getResidual(Vector& b) const {
237  b = *rhs_;
238  }
239 
240  void setMatrix(const SparseMatrixAdapter& /* M */) {
241  // matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead.
242  }
243 
244  bool solve(Vector& x) {
245  // Write linear system if asked for.
246  const int verbosity = prm_.get<int>("verbosity", 0);
247  const bool write_matrix = verbosity > 10;
248  if (write_matrix) {
249  Helper::writeSystem(simulator_, //simulator is only used to get names
250  getMatrix(),
251  *rhs_,
252  comm_.get());
253  }
254 
255  // Solve system.
256  Dune::InverseOperatorResult result;
257  bool accelerator_was_used = false;
258 
259  // Use GPU if: available, chosen by user, and successful.
260  // Use FPGA if: support compiled, chosen by user, and successful.
261 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
262  bool use_gpu = bdaBridge->getUseGpu();
263  bool use_fpga = bdaBridge->getUseFpga();
264  if (use_gpu || use_fpga) {
265  const std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
266  auto wellContribs = WellContributions::create(accelerator_mode, useWellConn_);
267  bdaBridge->initWellContributions(*wellContribs);
268 
269  // the WellContributions can only be applied separately with CUDA or OpenCL, not with an FPGA or amgcl
270 #if HAVE_CUDA || HAVE_OPENCL
271  if (!useWellConn_) {
272  simulator_.problem().wellModel().getWellContributions(*wellContribs);
273  }
274 #endif
275 
276  // Const_cast needed since the CUDA stuff overwrites values for better matrix condition..
277  bdaBridge->solve_system(const_cast<Matrix*>(&getMatrix()), *rhs_, *wellContribs, result);
278  if (result.converged) {
279  // get result vector x from non-Dune backend, iff solve was successful
280  bdaBridge->get_result(x);
281  accelerator_was_used = true;
282  } else {
283  // warn about CPU fallback
284  // BdaBridge might have disabled its BdaSolver for this simulation due to some error
285  // in that case the BdaBridge is disabled and flexibleSolver is always used
286  // or maybe the BdaSolver did not converge in time, then it will be used next linear solve
287  if (simulator_.gridView().comm().rank() == 0) {
288  OpmLog::warning(bdaBridge->getAccleratorName() + " did not converge, now trying Dune to solve current linear system...");
289  }
290  }
291  }
292 #endif
293 
294  // Otherwise, use flexible istl solver.
295  if (!accelerator_was_used) {
296  assert(flexibleSolver_);
297  flexibleSolver_->apply(x, *rhs_, result);
298  }
299 
300  // Check convergence, iterations etc.
301  checkConvergence(result);
302 
303  return converged_;
304  }
305 
306 
312 
314  int iterations () const { return iterations_; }
315 
317  const std::any& parallelInformation() const { return parallelInformation_; }
318 
319  protected:
320  // 3x3 matrix block inversion was unstable at least 2.3 until and including
321  // 2.5.0. There may still be some issue with the 4x4 matrix block inversion
322  // we therefore still use the block inversion in OPM
323  typedef ParallelOverlappingILU0<Dune::BCRSMatrix<Dune::MatrixBlock<typename Matrix::field_type,
324  Matrix::block_type::rows,
325  Matrix::block_type::cols> >,
326  Vector, Vector> SeqPreconditioner;
327 
328 #if HAVE_MPI
329  typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
330  // 3x3 matrix block inversion was unstable from at least 2.3 until and
331  // including 2.5.0
332  typedef ParallelOverlappingILU0<Matrix,Vector,Vector,Comm> ParPreconditioner;
333 #endif
334 
335  void checkConvergence( const Dune::InverseOperatorResult& result ) const
336  {
337  // store number of iterations
338  iterations_ = result.iterations;
339  converged_ = result.converged;
340 
341  // Check for failure of linear solver.
342  if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
343  const std::string msg("Convergence failure for linear solver.");
344  OPM_THROW_NOLOG(NumericalIssue, msg);
345  }
346  }
347  protected:
348 
349  bool isParallel() const {
350 #if HAVE_MPI
351  return comm_->communicator().size() > 1;
352 #else
353  return false;
354 #endif
355  }
356 
357  void prepareFlexibleSolver()
358  {
359 
360  std::function<Vector()> weightsCalculator = getWeightsCalculator();
361 
362  if (shouldCreateSolver()) {
363  if (isParallel()) {
364 #if HAVE_MPI
365  if (useWellConn_) {
366  using ParOperatorType = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
367  linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
368  flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
369  pressureIndex);
370  } else {
371  using ParOperatorType = WellModelGhostLastMatrixAdapter<Matrix, Vector, Vector, true>;
372  wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
373  linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *wellOperator_, interiorCellNum_);
374  flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
375  pressureIndex);
376  }
377 #endif
378  } else {
379  if (useWellConn_) {
380  using SeqOperatorType = Dune::MatrixAdapter<Matrix, Vector, Vector>;
381  linearOperatorForFlexibleSolver_ = std::make_unique<SeqOperatorType>(getMatrix());
382  flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator,
383  pressureIndex);
384  } else {
385  using SeqOperatorType = WellModelMatrixAdapter<Matrix, Vector, Vector, false>;
386  wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
387  linearOperatorForFlexibleSolver_ = std::make_unique<SeqOperatorType>(getMatrix(), *wellOperator_);
388  flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator,
389  pressureIndex);
390  }
391  }
392  }
393  else
394  {
395  flexibleSolver_->preconditioner().update();
396  }
397  }
398 
399 
402  bool shouldCreateSolver() const
403  {
404  // Decide if we should recreate the solver or just do
405  // a minimal preconditioner update.
406  if (!flexibleSolver_) {
407  return true;
408  }
409  if (this->parameters_.cpr_reuse_setup_ == 0) {
410  // Always recreate solver.
411  return true;
412  }
413  if (this->parameters_.cpr_reuse_setup_ == 1) {
414  // Recreate solver on the first iteration of every timestep.
415  const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
416  return newton_iteration == 0;
417  }
418  if (this->parameters_.cpr_reuse_setup_ == 2) {
419  // Recreate solver if the last solve used more than 10 iterations.
420  return this->iterations() > 10;
421  }
422 
423  // Otherwise, do not recreate solver.
424  assert(this->parameters_.cpr_reuse_setup_ == 3);
425 
426  return false;
427  }
428 
429 
431  std::function<Vector()> getWeightsCalculator() const
432  {
433  std::function<Vector()> weightsCalculator;
434 
435  using namespace std::string_literals;
436 
437  auto preconditionerType = prm_.get("preconditioner.type"s, "cpr"s);
438  if (preconditionerType == "cpr" || preconditionerType == "cprt") {
439  const bool transpose = preconditionerType == "cprt";
440  const auto weightsType = prm_.get("preconditioner.weight_type"s, "quasiimpes"s);
441  if (weightsType == "quasiimpes") {
442  // weights will be created as default in the solver
443  // assignment p = pressureIndex prevent compiler warning about
444  // capturing variable with non-automatic storage duration
445  weightsCalculator = [this, transpose, p = pressureIndex]() {
446  return Amg::getQuasiImpesWeights<Matrix, Vector>(this->getMatrix(), p, transpose);
447  };
448  } else if (weightsType == "trueimpes") {
449  // assignment p = pressureIndex prevent compiler warning about
450  // capturing variable with non-automatic storage duration
451  weightsCalculator = [this, p = pressureIndex]() {
452  return this->getTrueImpesWeights(p);
453  };
454  } else {
455  OPM_THROW(std::invalid_argument,
456  "Weights type " << weightsType << "not implemented for cpr."
457  << " Please use quasiimpes or trueimpes.");
458  }
459  }
460  return weightsCalculator;
461  }
462 
463 
464  // Weights to make approximate pressure equations.
465  // Calculated from the storage terms (only) of the
466  // conservation equations, ignoring all other terms.
467  Vector getTrueImpesWeights(int pressureVarIndex) const
468  {
469  Vector weights(rhs_->size());
470  ElementContext elemCtx(simulator_);
471  Amg::getTrueImpesWeights(pressureVarIndex, weights, simulator_.vanguard().gridView(),
472  elemCtx, simulator_.model(),
473  ThreadManager::threadId());
474  return weights;
475  }
476 
477 
480  void makeOverlapRowsInvalid(Matrix& matrix) const
481  {
482  //value to set on diagonal
483  const int numEq = Matrix::block_type::rows;
484  typename Matrix::block_type diag_block(0.0);
485  for (int eq = 0; eq < numEq; ++eq)
486  diag_block[eq][eq] = 1.0;
487 
488  //loop over precalculated overlap rows and columns
489  for (auto row = overlapRows_.begin(); row != overlapRows_.end(); row++ )
490  {
491  int lcell = *row;
492  // Zero out row.
493  matrix[lcell] = 0.0;
494 
495  //diagonal block set to diag(1.0).
496  matrix[lcell][lcell] = diag_block;
497  }
498  }
499 
500 
501  Matrix& getMatrix()
502  {
503  return *matrix_;
504  }
505 
506  const Matrix& getMatrix() const
507  {
508  return *matrix_;
509  }
510 
511  const Simulator& simulator_;
512  mutable int iterations_;
513  mutable bool converged_;
514  std::any parallelInformation_;
515 
516  // non-const to be able to scale the linear system
517  Matrix* matrix_;
518  Vector *rhs_;
519 
520  std::unique_ptr<FlexibleSolverType> flexibleSolver_;
521  std::unique_ptr<AbstractOperatorType> linearOperatorForFlexibleSolver_;
522  std::unique_ptr<WellModelAsLinearOperator<WellModel, Vector, Vector>> wellOperator_;
523  std::vector<int> overlapRows_;
524  std::vector<int> interiorRows_;
525  std::vector<std::set<int>> wellConnectionsGraph_;
526 
527  bool useWellConn_;
528  size_t interiorCellNum_;
529 
530  FlowLinearSolverParameters parameters_;
531  PropertyTree prm_;
532  bool scale_variables_;
533 
534  std::shared_ptr< CommunicationType > comm_;
535  }; // end ISTLSolver
536 
537 } // namespace Opm
538 #endif
A solver class that encapsulates all needed objects for a linear solver (operator,...
Definition: FlexibleSolver.hpp:40
Definition: MatrixBlock.hpp:370
BdaBridge acts as interface between opm-simulators with the BdaSolvers.
Definition: BdaBridge.hpp:39
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition: ISTLSolverEbos.hpp:78
std::function< Vector()> getWeightsCalculator() const
Return an appropriate weight function if a cpr preconditioner is asked for.
Definition: ISTLSolverEbos.hpp:431
const std::any & parallelInformation() const
Definition: ISTLSolverEbos.hpp:317
void makeOverlapRowsInvalid(Matrix &matrix) const
Zero out off-diagonal blocks on rows corresponding to overlap cells Diagonal blocks on ovelap rows ar...
Definition: ISTLSolverEbos.hpp:480
int iterations() const
Solve the system of linear equations Ax = b, with A being the combined derivative matrix of the resid...
Definition: ISTLSolverEbos.hpp:314
ISTLSolverEbos(const Simulator &simulator)
Construct a system solver.
Definition: ISTLSolverEbos.hpp:118
bool shouldCreateSolver() const
Return true if we should (re)create the whole solver, instead of just calling update() on the precond...
Definition: ISTLSolverEbos.hpp:402
Definition: MatrixMarketSpecializations.hpp:17
A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as.
Definition: ParallelOverlappingILU0.hpp:611
Linear operator wrapper for well model.
Definition: WellOperators.hpp:50
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
void extractParallelGridInformationToISTL(std::any &anyComm, const UnstructuredGrid &grid)
Extracts the information about the data decomposition from the grid for dune-istl.
Definition: ParallelIstlInformation.hpp:676
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:36
Definition: ISTLSolverEbos.hpp:51
Definition: ISTLSolverEbos.hpp:45