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 <dune/istl/owneroverlapcopy.hh>
26 
27 #include <ebos/eclbasevanguard.hh>
28 
29 #include <opm/common/ErrorMacros.hpp>
30 
31 #include <opm/models/discretization/common/fvbaseproperties.hh>
32 #include <opm/models/common/multiphasebaseproperties.hh>
33 #include <opm/models/utils/parametersystem.hh>
34 #include <opm/models/utils/propertysystem.hh>
35 #include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
36 #include <opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp>
37 #include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
38 #include <opm/simulators/linalg/matrixblock.hh>
39 #include <opm/simulators/linalg/istlsparsematrixadapter.hh>
40 #include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
41 #include <opm/simulators/linalg/WellOperators.hpp>
42 #include <opm/simulators/linalg/WriteSystemMatrixHelper.hpp>
43 #include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
44 #include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
45 #include <opm/simulators/linalg/setupPropertyTree.hpp>
46 
47 #include <any>
48 #include <cstddef>
49 #include <functional>
50 #include <memory>
51 #include <set>
52 #include <sstream>
53 #include <string>
54 #include <tuple>
55 #include <vector>
56 
57 namespace Opm::Properties {
58 
59 namespace TTag {
61  using InheritsFrom = std::tuple<FlowIstlSolverParams>;
62 };
63 }
64 
65 template <class TypeTag, class MyTypeTag>
66 struct EclWellModel;
67 
70 template<class TypeTag>
71 struct SparseMatrixAdapter<TypeTag, TTag::FlowIstlSolver>
72 {
73 private:
74  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
75  enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
77 
78 public:
79  using type = typename Linear::IstlSparseMatrixAdapter<Block>;
80 };
81 
82 } // namespace Opm::Properties
83 
84 namespace Opm
85 {
86 
87 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
88 template<class Matrix, class Vector, int block_size> class BdaBridge;
89 class WellContributions;
90 #endif
91 
92 namespace detail
93 {
94 
95 template<class Matrix, class Vector, class Comm>
97 {
98  using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
99  using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
101 
102  void create(const Matrix& matrix,
103  bool parallel,
104  const PropertyTree& prm,
105  size_t pressureIndex,
106  std::function<Vector()> trueFunc,
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  size_t interiorCellNum_ = 0;
114 };
115 
116 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
117 template<class Matrix, class Vector>
118 struct BdaSolverInfo
119 {
120  using WellContribFunc = std::function<void(WellContributions&)>;
122 
123  BdaSolverInfo(const std::string& accelerator_mode,
124  const std::string& fpga_bitstream,
125  const int linear_solver_verbosity,
126  const int maxit,
127  const double tolerance,
128  const int platformID,
129  const int deviceID,
130  const std::string& opencl_ilu_reorder,
131  const std::string& linsolver);
132 
133  ~BdaSolverInfo();
134 
135  template<class Grid>
136  void prepare(const Grid& grid,
137  const Dune::CartesianIndexMapper<Grid>& cartMapper,
138  const std::vector<Well>& wellsForConn,
139  const std::vector<int>& cellPartition,
140  const size_t nonzeroes,
141  const bool useWellConn);
142 
143  bool apply(Vector& rhs,
144  const bool useWellConn,
145  WellContribFunc getContribs,
146  const int rank,
147  Matrix& matrix,
148  Vector& x,
149  Dune::InverseOperatorResult& result);
150 
151  int numJacobiBlocks_ = 0;
152 
153 private:
156  template<class Grid>
157  void blockJacobiAdjacency(const Grid& grid,
158  const std::vector<int>& cell_part,
159  size_t nonzeroes);
160 
161  void copyMatToBlockJac(const Matrix& mat, Matrix& blockJac);
162 
163  std::unique_ptr<Bridge> bridge_;
164  std::string accelerator_mode_;
165  std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
166  std::vector<std::set<int>> wellConnectionsGraph_;
167 };
168 #endif
169 
170 #ifdef HAVE_MPI
172 void copyParValues(std::any& parallelInformation, size_t size,
173  Dune::OwnerOverlapCopyCommunication<int,int>& comm);
174 #endif
175 
178 template<class Matrix>
179 void makeOverlapRowsInvalid(Matrix& matrix,
180  const std::vector<int>& overlapRows);
181 
184 template<class Matrix, class Grid>
185 std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
186  const std::vector<int>& cell_part,
187  size_t nonzeroes,
188  const std::vector<std::set<int>>& wellConnectionsGraph);
189 }
190 
195  template <class TypeTag>
197  {
198  protected:
199  using GridView = GetPropType<TypeTag, Properties::GridView>;
200  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
201  using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
202  using Vector = GetPropType<TypeTag, Properties::GlobalEqVector>;
203  using Indices = GetPropType<TypeTag, Properties::Indices>;
204  using WellModel = GetPropType<TypeTag, Properties::EclWellModel>;
205  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
206  using Matrix = typename SparseMatrixAdapter::IstlMatrix;
207  using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
208  using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
209  using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
210  using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
213  using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
214  constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
215 
216 #if HAVE_MPI
217  using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
218 #else
219  using CommunicationType = Dune::CollectiveCommunication<int>;
220 #endif
221 
222  public:
223  using AssembledLinearOperatorType = Dune::AssembledLinearOperator< Matrix, Vector, Vector >;
224 
225  static void registerParameters()
226  {
227  FlowLinearSolverParameters::registerParameters<TypeTag>();
228  }
229 
233  explicit ISTLSolverEbos(const Simulator& simulator)
234  : simulator_(simulator),
235  iterations_( 0 ),
236  calls_( 0 ),
237  converged_(false),
238  matrix_()
239  {
240  const bool on_io_rank = (simulator.gridView().comm().rank() == 0);
241 #if HAVE_MPI
242  comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) );
243 #endif
244  parameters_.template init<TypeTag>();
245  prm_ = setupPropertyTree(parameters_,
246  EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter),
247  EWOMS_PARAM_IS_SET(TypeTag, int, CprMaxEllIter));
248 
249 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
250  {
251  std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
252  if ((simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode != "none")) {
253  if (on_io_rank) {
254  OpmLog::warning("Cannot use GPU or FPGA with MPI, GPU/FPGA are disabled");
255  }
256  accelerator_mode = "none";
257  }
258  const int platformID = EWOMS_GET_PARAM(TypeTag, int, OpenclPlatformId);
259  const int deviceID = EWOMS_GET_PARAM(TypeTag, int, BdaDeviceId);
260  const int maxit = EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIter);
261  const double tolerance = EWOMS_GET_PARAM(TypeTag, double, LinearSolverReduction);
262  const std::string opencl_ilu_reorder = EWOMS_GET_PARAM(TypeTag, std::string, OpenclIluReorder);
263  const int linear_solver_verbosity = parameters_.linear_solver_verbosity_;
264  std::string fpga_bitstream = EWOMS_GET_PARAM(TypeTag, std::string, FpgaBitstream);
265  std::string linsolver = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolver);
266  bdaBridge = std::make_unique<detail::BdaSolverInfo<Matrix,Vector>>(accelerator_mode,
267  fpga_bitstream,
268  linear_solver_verbosity,
269  maxit,
270  tolerance,
271  platformID,
272  deviceID,
273  opencl_ilu_reorder,
274  linsolver);
275  }
276 #else
277  if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) != "none") {
278  OPM_THROW(std::logic_error,"Cannot use accelerated solver since CUDA, OpenCL and amgcl were not found by cmake and FPGA was not enabled");
279  }
280 #endif
281  extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
282 
283  // For some reason simulator_.model().elementMapper() is not initialized at this stage
284  //const auto& elemMapper = simulator_.model().elementMapper(); //does not work.
285  // Set it up manually
286  ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
287  detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
288  useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
289 #if HAVE_FPGA
290  // check usage of MatrixAddWellContributions: for FPGA they must be included
291  if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) == "fpga" && !useWellConn_) {
292  OPM_THROW(std::logic_error,"fpgaSolver needs --matrix-add-well-contributions=true");
293  }
294 #endif
295  const bool ownersFirst = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst);
296  if (!ownersFirst) {
297  const std::string msg = "The linear solver no longer supports --owner-cells-first=false.";
298  if (on_io_rank) {
299  OpmLog::error(msg);
300  }
301  OPM_THROW_NOLOG(std::runtime_error, msg);
302  }
303 
304  flexibleSolver_.interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true);
305 
306  // Print parameters to PRT/DBG logs.
307  if (on_io_rank) {
308  std::ostringstream os;
309  os << "Property tree for linear solver:\n";
310  prm_.write_json(os, true);
311  OpmLog::note(os.str());
312  }
313  }
314 
315  // nothing to clean here
316  void eraseMatrix()
317  {
318  }
319 
320  void prepare(const SparseMatrixAdapter& M, Vector& b)
321  {
322  static bool firstcall = true;
323 #if HAVE_MPI
324  if (firstcall) {
325  const size_t size = M.istlMatrix().N();
326  detail::copyParValues(parallelInformation_, size, *comm_);
327  }
328 #endif
329 
330  // update matrix entries for solvers.
331  if (firstcall) {
332  // ebos will not change the matrix object. Hence simply store a pointer
333  // to the original one with a deleter that does nothing.
334  // Outch! We need to be able to scale the linear system! Hence const_cast
335  matrix_ = const_cast<Matrix*>(&M.istlMatrix());
336 
337  useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
338  // setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver)
339 #if HAVE_OPENCL
340  bdaBridge->numJacobiBlocks_ = EWOMS_GET_PARAM(TypeTag, int, NumJacobiBlocks);
341  bdaBridge->prepare(simulator_.vanguard().grid(),
342  simulator_.vanguard().cartesianIndexMapper(),
343  simulator_.vanguard().schedule().getWellsatEnd(),
344  simulator_.vanguard().cellPartition(),
345  getMatrix().nonzeroes(), useWellConn_);
346 #endif
347  } else {
348  // Pointers should not change
349  if ( &(M.istlMatrix()) != matrix_ ) {
350  OPM_THROW(std::logic_error, "Matrix objects are expected to be reused when reassembling!"
351  <<" old pointer was " << matrix_ << ", new one is " << (&M.istlMatrix()) );
352  }
353  }
354  rhs_ = &b;
355 
356  if (isParallel() && prm_.get<std::string>("preconditioner.type") != "ParOverILU0") {
357  detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
358  }
359  prepareFlexibleSolver();
360  firstcall = false;
361  }
362 
363 
364  void setResidual(Vector& /* b */)
365  {
366  // rhs_ = &b; // Must be handled in prepare() instead.
367  }
368 
369  void getResidual(Vector& b) const
370  {
371  b = *rhs_;
372  }
373 
374  void setMatrix(const SparseMatrixAdapter& /* M */)
375  {
376  // matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead.
377  }
378 
379  bool solve(Vector& x)
380  {
381  calls_ += 1;
382  // Write linear system if asked for.
383  const int verbosity = prm_.get<int>("verbosity", 0);
384  const bool write_matrix = verbosity > 10;
385  if (write_matrix) {
386  Helper::writeSystem(simulator_, //simulator is only used to get names
387  getMatrix(),
388  *rhs_,
389  comm_.get());
390  }
391 
392  // Solve system.
393  Dune::InverseOperatorResult result;
394 
395  // Use GPU if: available, chosen by user, and successful.
396  // Use FPGA if: support compiled, chosen by user, and successful.
397 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
398  std::function<void(WellContributions&)> getContribs =
399  [this](WellContributions& w)
400  {
401  this->simulator_.problem().wellModel().getWellContributions(w);
402  };
403  if (!bdaBridge->apply(*rhs_, useWellConn_, getContribs,
404  simulator_.gridView().comm().rank(),
405  const_cast<Matrix&>(this->getMatrix()),
406  x, result))
407 #endif
408  {
409  assert(flexibleSolver_.solver_);
410  flexibleSolver_.solver_->apply(x, *rhs_, result);
411  }
412 
413  // Check convergence, iterations etc.
414  checkConvergence(result);
415 
416  return converged_;
417  }
418 
419 
425 
427  int iterations () const { return iterations_; }
428 
430  const std::any& parallelInformation() const { return parallelInformation_; }
431 
432  protected:
433 #if HAVE_MPI
434  using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
435 #endif
436 
437  void checkConvergence( const Dune::InverseOperatorResult& result ) const
438  {
439  // store number of iterations
440  iterations_ = result.iterations;
441  converged_ = result.converged;
442 
443  // Check for failure of linear solver.
444  if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
445  const std::string msg("Convergence failure for linear solver.");
446  OPM_THROW_NOLOG(NumericalIssue, msg);
447  }
448  }
449  protected:
450 
451  bool isParallel() const {
452 #if HAVE_MPI
453  return comm_->communicator().size() > 1;
454 #else
455  return false;
456 #endif
457  }
458 
459  void prepareFlexibleSolver()
460  {
461 
462  if (shouldCreateSolver()) {
463  std::function<Vector()> trueFunc =
464  [this]
465  {
466  return this->getTrueImpesWeights(pressureIndex);
467  };
468 
469  if (!useWellConn_) {
470  auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
471  flexibleSolver_.wellOperator_ = std::move(wellOp);
472  }
473 
474  flexibleSolver_.create(getMatrix(),
475  isParallel(),
476  prm_,
477  pressureIndex,
478  trueFunc,
479  *comm_);
480  }
481  else
482  {
483  flexibleSolver_.pre_->update();
484  }
485  }
486 
487 
490  bool shouldCreateSolver() const
491  {
492  // Decide if we should recreate the solver or just do
493  // a minimal preconditioner update.
494  if (!flexibleSolver_.solver_) {
495  return true;
496  }
497  if (this->parameters_.cpr_reuse_setup_ == 0) {
498  // Always recreate solver.
499  return true;
500  }
501  if (this->parameters_.cpr_reuse_setup_ == 1) {
502  // Recreate solver on the first iteration of every timestep.
503  const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
504  return newton_iteration == 0;
505  }
506  if (this->parameters_.cpr_reuse_setup_ == 2) {
507  // Recreate solver if the last solve used more than 10 iterations.
508  return this->iterations() > 10;
509  }
510  if (this->parameters_.cpr_reuse_setup_ == 3) {
511  // Recreate solver if the last solve used more than 10 iterations.
512  return false;
513  }
514  if (this->parameters_.cpr_reuse_setup_ == 4) {
515  // Recreate solver every 'step' solve calls.
516  const int step = this->parameters_.cpr_reuse_interval_;
517  const bool create = ((calls_ % step) == 0);
518  return create;
519  }
520 
521  // If here, we have an invalid parameter.
522  const bool on_io_rank = (simulator_.gridView().comm().rank() == 0);
523  std::string msg = "Invalid value: " + std::to_string(this->parameters_.cpr_reuse_setup_)
524  + " for --cpr-reuse-setup parameter, run with --help to see allowed values.";
525  if (on_io_rank) {
526  OpmLog::error(msg);
527  }
528  throw std::runtime_error(msg);
529 
530  // Never reached.
531  return false;
532  }
533 
534 
535  // Weights to make approximate pressure equations.
536  // Calculated from the storage terms (only) of the
537  // conservation equations, ignoring all other terms.
538  Vector getTrueImpesWeights(int pressureVarIndex) const
539  {
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 calls_;
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 #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
570  std::unique_ptr<detail::BdaSolverInfo<Matrix, Vector>> bdaBridge;
571 #endif
572 
573  detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType> flexibleSolver_;
574  std::vector<int> overlapRows_;
575  std::vector<int> interiorRows_;
576 
577  bool useWellConn_;
578 
579  FlowLinearSolverParameters parameters_;
580  PropertyTree prm_;
581 
582  std::shared_ptr< CommunicationType > comm_;
583  }; // end ISTLSolver
584 
585 } // namespace Opm
586 #endif
Definition: findOverlapRowsAndColumns.hpp:29
Interface class adding the update() method to the preconditioner interface.
Definition: PreconditionerWithUpdate.hpp:32
BdaBridge acts as interface between opm-simulators with the BdaSolvers.
Definition: BdaBridge.hpp:40
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition: ISTLSolverEbos.hpp:197
const std::any & parallelInformation() const
Definition: ISTLSolverEbos.hpp:430
int iterations() const
Solve the system of linear equations Ax = b, with A being the combined derivative matrix of the resid...
Definition: ISTLSolverEbos.hpp:427
ISTLSolverEbos(const Simulator &simulator)
Construct a system solver.
Definition: ISTLSolverEbos.hpp:233
bool shouldCreateSolver() const
Return true if we should (re)create the whole solver, instead of just calling update() on the precond...
Definition: ISTLSolverEbos.hpp:490
Definition: MatrixMarketSpecializations.hpp:18
Definition: PropertyTree.hpp:37
This class serves to eliminate the need to include the WellContributions into the matrix (with –matri...
Definition: WellContributions.hpp:53
Definition: WellOperators.hpp:60
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
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:40
Definition: ISTLSolverEbos.hpp:66
Definition: ISTLSolverEbos.hpp:60
Definition: ISTLSolverEbos.hpp:97