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 
38 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
39 #include <opm/simulators/linalg/bda/BdaBridge.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  bdaBridge.reset(new BdaBridge<Matrix, Vector, block_size>(accelerator_mode, fpga_bitstream, linear_solver_verbosity, maxit, tolerance, platformID, deviceID, opencl_ilu_reorder));
150  }
151 #else
152  if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) != "none") {
153  OPM_THROW(std::logic_error,"Cannot use accelerated solver since CUDA, OpenCL and amgcl were not found by cmake and FPGA was not enabled");
154  }
155 #endif
156  extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
157 
158  // For some reason simulator_.model().elementMapper() is not initialized at this stage
159  // Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work.
160  // Set it up manually
161  ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
162  detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
163 
164  useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
165 #if HAVE_FPGA
166  // check usage of MatrixAddWellContributions: for FPGA they must be included
167  if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) == "fpga" && !useWellConn_) {
168  OPM_THROW(std::logic_error,"fpgaSolver needs --matrix-add-well-contributions=true");
169  }
170 #endif
171  const bool ownersFirst = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst);
172  if (!ownersFirst) {
173  const std::string msg = "The linear solver no longer supports --owner-cells-first=false.";
174  if (on_io_rank) {
175  OpmLog::error(msg);
176  }
177  OPM_THROW_NOLOG(std::runtime_error, msg);
178  }
179 
180  interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true);
181 
182  // Print parameters to PRT/DBG logs.
183  if (on_io_rank) {
184  std::ostringstream os;
185  os << "Property tree for linear solver:\n";
186  prm_.write_json(os, true);
187  OpmLog::note(os.str());
188  }
189  }
190 
191  // nothing to clean here
192  void eraseMatrix() {
193  }
194 
195  void prepare(const SparseMatrixAdapter& M, Vector& b)
196  {
197  static bool firstcall = true;
198 #if HAVE_MPI
199  if (firstcall && parallelInformation_.type() == typeid(ParallelISTLInformation)) {
200  // Parallel case.
201  const ParallelISTLInformation* parinfo = std::any_cast<ParallelISTLInformation>(&parallelInformation_);
202  assert(parinfo);
203  const size_t size = M.istlMatrix().N();
204  parinfo->copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), size, 1);
205  }
206 #endif
207 
208  // update matrix entries for solvers.
209  if (firstcall) {
210  // ebos will not change the matrix object. Hence simply store a pointer
211  // to the original one with a deleter that does nothing.
212  // Outch! We need to be able to scale the linear system! Hence const_cast
213  matrix_ = const_cast<Matrix*>(&M.istlMatrix());
214  } else {
215  // Pointers should not change
216  if ( &(M.istlMatrix()) != matrix_ ) {
217  OPM_THROW(std::logic_error, "Matrix objects are expected to be reused when reassembling!"
218  <<" old pointer was " << matrix_ << ", new one is " << (&M.istlMatrix()) );
219  }
220  }
221  rhs_ = &b;
222 
223  if (isParallel() && prm_.get<std::string>("preconditioner.type") != "ParOverILU0") {
224  makeOverlapRowsInvalid(getMatrix());
225  }
226  prepareFlexibleSolver();
227  firstcall = false;
228  }
229 
230 
231  void setResidual(Vector& /* b */) {
232  // rhs_ = &b; // Must be handled in prepare() instead.
233  }
234 
235  void getResidual(Vector& b) const {
236  b = *rhs_;
237  }
238 
239  void setMatrix(const SparseMatrixAdapter& /* M */) {
240  // matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead.
241  }
242 
243  bool solve(Vector& x) {
244  // Write linear system if asked for.
245  const int verbosity = prm_.get<int>("verbosity", 0);
246  const bool write_matrix = verbosity > 10;
247  if (write_matrix) {
248  Helper::writeSystem(simulator_, //simulator is only used to get names
249  getMatrix(),
250  *rhs_,
251  comm_.get());
252  }
253 
254  // Solve system.
255  Dune::InverseOperatorResult result;
256  bool accelerator_was_used = false;
257 
258  // Use GPU if: available, chosen by user, and successful.
259  // Use FPGA if: support compiled, chosen by user, and successful.
260 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
261  bool use_gpu = bdaBridge->getUseGpu();
262  bool use_fpga = bdaBridge->getUseFpga();
263  if (use_gpu || use_fpga) {
264  const std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
265  WellContributions wellContribs(accelerator_mode, useWellConn_);
266  bdaBridge->initWellContributions(wellContribs);
267 
268  // the WellContributions can only be applied separately with CUDA or OpenCL, not with an FPGA or amgcl
269 #if HAVE_CUDA || HAVE_OPENCL
270  if (!useWellConn_) {
271  simulator_.problem().wellModel().getWellContributions(wellContribs);
272  }
273 #endif
274 
275  // Const_cast needed since the CUDA stuff overwrites values for better matrix condition..
276  bdaBridge->solve_system(const_cast<Matrix*>(&getMatrix()), *rhs_, wellContribs, result);
277  if (result.converged) {
278  // get result vector x from non-Dune backend, iff solve was successful
279  bdaBridge->get_result(x);
280  accelerator_was_used = true;
281  } else {
282  // warn about CPU fallback
283  // BdaBridge might have disabled its BdaSolver for this simulation due to some error
284  // in that case the BdaBridge is disabled and flexibleSolver is always used
285  // or maybe the BdaSolver did not converge in time, then it will be used next linear solve
286  if (simulator_.gridView().comm().rank() == 0) {
287  OpmLog::warning(bdaBridge->getAccleratorName() + " did not converge, now trying Dune to solve current linear system...");
288  }
289  }
290  }
291 #endif
292 
293  // Otherwise, use flexible istl solver.
294  if (!accelerator_was_used) {
295  assert(flexibleSolver_);
296  flexibleSolver_->apply(x, *rhs_, result);
297  }
298 
299  // Check convergence, iterations etc.
300  checkConvergence(result);
301 
302  return converged_;
303  }
304 
305 
311 
313  int iterations () const { return iterations_; }
314 
316  const std::any& parallelInformation() const { return parallelInformation_; }
317 
318  protected:
319  // 3x3 matrix block inversion was unstable at least 2.3 until and including
320  // 2.5.0. There may still be some issue with the 4x4 matrix block inversion
321  // we therefore still use the block inversion in OPM
322  typedef ParallelOverlappingILU0<Dune::BCRSMatrix<Dune::MatrixBlock<typename Matrix::field_type,
323  Matrix::block_type::rows,
324  Matrix::block_type::cols> >,
325  Vector, Vector> SeqPreconditioner;
326 
327 #if HAVE_MPI
328  typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
329  // 3x3 matrix block inversion was unstable from at least 2.3 until and
330  // including 2.5.0
331  typedef ParallelOverlappingILU0<Matrix,Vector,Vector,Comm> ParPreconditioner;
332 #endif
333 
334  void checkConvergence( const Dune::InverseOperatorResult& result ) const
335  {
336  // store number of iterations
337  iterations_ = result.iterations;
338  converged_ = result.converged;
339 
340  // Check for failure of linear solver.
341  if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
342  const std::string msg("Convergence failure for linear solver.");
343  OPM_THROW_NOLOG(NumericalIssue, msg);
344  }
345  }
346  protected:
347 
348  bool isParallel() const {
349 #if HAVE_MPI
350  return comm_->communicator().size() > 1;
351 #else
352  return false;
353 #endif
354  }
355 
356  void prepareFlexibleSolver()
357  {
358 
359  std::function<Vector()> weightsCalculator = getWeightsCalculator();
360 
361  if (shouldCreateSolver()) {
362  if (isParallel()) {
363 #if HAVE_MPI
364  if (useWellConn_) {
365  using ParOperatorType = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
366  linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
367  flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
368  pressureIndex);
369  } else {
370  using ParOperatorType = WellModelGhostLastMatrixAdapter<Matrix, Vector, Vector, true>;
371  wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
372  linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *wellOperator_, interiorCellNum_);
373  flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
374  pressureIndex);
375  }
376 #endif
377  } else {
378  if (useWellConn_) {
379  using SeqOperatorType = Dune::MatrixAdapter<Matrix, Vector, Vector>;
380  linearOperatorForFlexibleSolver_ = std::make_unique<SeqOperatorType>(getMatrix());
381  flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator,
382  pressureIndex);
383  } else {
384  using SeqOperatorType = WellModelMatrixAdapter<Matrix, Vector, Vector, false>;
385  wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
386  linearOperatorForFlexibleSolver_ = std::make_unique<SeqOperatorType>(getMatrix(), *wellOperator_);
387  flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator,
388  pressureIndex);
389  }
390  }
391  }
392  else
393  {
394  flexibleSolver_->preconditioner().update();
395  }
396  }
397 
398 
401  bool shouldCreateSolver() const
402  {
403  // Decide if we should recreate the solver or just do
404  // a minimal preconditioner update.
405  if (!flexibleSolver_) {
406  return true;
407  }
408  if (this->parameters_.cpr_reuse_setup_ == 0) {
409  // Always recreate solver.
410  return true;
411  }
412  if (this->parameters_.cpr_reuse_setup_ == 1) {
413  // Recreate solver on the first iteration of every timestep.
414  const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
415  return newton_iteration == 0;
416  }
417  if (this->parameters_.cpr_reuse_setup_ == 2) {
418  // Recreate solver if the last solve used more than 10 iterations.
419  return this->iterations() > 10;
420  }
421 
422  // Otherwise, do not recreate solver.
423  assert(this->parameters_.cpr_reuse_setup_ == 3);
424 
425  return false;
426  }
427 
428 
430  std::function<Vector()> getWeightsCalculator() const
431  {
432  std::function<Vector()> weightsCalculator;
433 
434  using namespace std::string_literals;
435 
436  auto preconditionerType = prm_.get("preconditioner.type"s, "cpr"s);
437  if (preconditionerType == "cpr" || preconditionerType == "cprt") {
438  const bool transpose = preconditionerType == "cprt";
439  const auto weightsType = prm_.get("preconditioner.weight_type"s, "quasiimpes"s);
440  if (weightsType == "quasiimpes") {
441  // weights will be created as default in the solver
442  // assignment p = pressureIndex prevent compiler warning about
443  // capturing variable with non-automatic storage duration
444  weightsCalculator = [this, transpose, p = pressureIndex]() {
445  return Amg::getQuasiImpesWeights<Matrix, Vector>(this->getMatrix(), p, transpose);
446  };
447  } else if (weightsType == "trueimpes") {
448  // assignment p = pressureIndex prevent compiler warning about
449  // capturing variable with non-automatic storage duration
450  weightsCalculator = [this, p = pressureIndex]() {
451  return this->getTrueImpesWeights(p);
452  };
453  } else {
454  OPM_THROW(std::invalid_argument,
455  "Weights type " << weightsType << "not implemented for cpr."
456  << " Please use quasiimpes or trueimpes.");
457  }
458  }
459  return weightsCalculator;
460  }
461 
462 
463  // Weights to make approximate pressure equations.
464  // Calculated from the storage terms (only) of the
465  // conservation equations, ignoring all other terms.
466  Vector getTrueImpesWeights(int pressureVarIndex) const
467  {
468  Vector weights(rhs_->size());
469  ElementContext elemCtx(simulator_);
470  Amg::getTrueImpesWeights(pressureVarIndex, weights, simulator_.vanguard().gridView(),
471  elemCtx, simulator_.model(),
472  ThreadManager::threadId());
473  return weights;
474  }
475 
476 
479  void makeOverlapRowsInvalid(Matrix& matrix) const
480  {
481  //value to set on diagonal
482  const int numEq = Matrix::block_type::rows;
483  typename Matrix::block_type diag_block(0.0);
484  for (int eq = 0; eq < numEq; ++eq)
485  diag_block[eq][eq] = 1.0;
486 
487  //loop over precalculated overlap rows and columns
488  for (auto row = overlapRows_.begin(); row != overlapRows_.end(); row++ )
489  {
490  int lcell = *row;
491  // Zero out row.
492  matrix[lcell] = 0.0;
493 
494  //diagonal block set to diag(1.0).
495  matrix[lcell][lcell] = diag_block;
496  }
497  }
498 
499 
500  Matrix& getMatrix()
501  {
502  return *matrix_;
503  }
504 
505  const Matrix& getMatrix() const
506  {
507  return *matrix_;
508  }
509 
510  const Simulator& simulator_;
511  mutable int iterations_;
512  mutable bool converged_;
513  std::any parallelInformation_;
514 
515  // non-const to be able to scale the linear system
516  Matrix* matrix_;
517  Vector *rhs_;
518 
519  std::unique_ptr<FlexibleSolverType> flexibleSolver_;
520  std::unique_ptr<AbstractOperatorType> linearOperatorForFlexibleSolver_;
521  std::unique_ptr<WellModelAsLinearOperator<WellModel, Vector, Vector>> wellOperator_;
522  std::vector<int> overlapRows_;
523  std::vector<int> interiorRows_;
524  std::vector<std::set<int>> wellConnectionsGraph_;
525 
526  bool useWellConn_;
527  size_t interiorCellNum_;
528 
529  FlowLinearSolverParameters parameters_;
530  PropertyTree prm_;
531  bool scale_variables_;
532 
533  std::shared_ptr< CommunicationType > comm_;
534  }; // end ISTLSolver
535 
536 } // namespace Opm
537 #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:46
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:430
const std::any & parallelInformation() const
Definition: ISTLSolverEbos.hpp:316
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:479
int iterations() const
Solve the system of linear equations Ax = b, with A being the combined derivative matrix of the resid...
Definition: ISTLSolverEbos.hpp:313
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:401
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:26
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:37
Definition: ISTLSolverEbos.hpp:51
Definition: ISTLSolverEbos.hpp:45