My Project
ISTLSolverEbosFlexible.hpp
1 /*
2  Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #ifndef OPM_ISTLSOLVEREBOSFLEXIBLE_HEADER_INCLUDED
21 #define OPM_ISTLSOLVEREBOSFLEXIBLE_HEADER_INCLUDED
22 
23 #include <opm/simulators/linalg/matrixblock.hh>
24 #include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
25 #include <opm/simulators/linalg/FlexibleSolver.hpp>
26 #include <opm/simulators/linalg/setupPropertyTree.hpp>
27 #include <opm/simulators/linalg/WriteSystemMatrixHelper.hpp>
28 
29 #include <opm/common/ErrorMacros.hpp>
30 
31 #include <memory>
32 #include <utility>
33 
34 namespace Opm::Properties {
35 
36 namespace TTag {
38  using InheritsFrom = std::tuple<FlowIstlSolverParams>;
39 };
40 }
41 
42 } // namespace Opm::Properties
43 
44 
45 namespace Opm
46 {
47 
48 //=====================================================================
49 // Implementation for ISTL-matrix based operator
50 //=====================================================================
57 template <class TypeTag>
59 {
60  using GridView = GetPropType<TypeTag, Properties::GridView>;
61  using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
62  using VectorType = GetPropType<TypeTag, Properties::GlobalEqVector>;
63  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
64  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
65  using MatrixType = typename SparseMatrixAdapter::IstlMatrix;
66  using WellModel = GetPropType<TypeTag, Properties::EclWellModel>;
67 #if HAVE_MPI
68  using Communication = Dune::OwnerOverlapCopyCommunication<int, int>;
69 #else
70  using Communication = int; // Dummy type.
71 #endif
72  using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
75 
76  // for quasiImpesWeights
77  using Vector = GetPropType<TypeTag, Properties::GlobalEqVector>;
78  using Indices = GetPropType<TypeTag, Properties::Indices>;
79  typedef typename SparseMatrixAdapter::IstlMatrix Matrix;
80  typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType;
81  typedef typename Vector::block_type BlockVector;
82  using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
83  using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
84  typedef typename GridView::template Codim<0>::Entity Element;
85  using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
86  constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
87 
88 public:
89  static void registerParameters()
90  {
91  FlowLinearSolverParameters::registerParameters<TypeTag>();
92  }
93 
94  explicit ISTLSolverEbosFlexible(const Simulator& simulator)
95  : simulator_(simulator)
96  , ownersFirst_(EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst))
97  , matrixAddWellContributions_(EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions))
98  , interiorCellNum_(detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), ownersFirst_))
99  {
100  parameters_.template init<TypeTag>();
101  prm_ = setupPropertyTree(parameters_,
102  EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter),
103  EWOMS_PARAM_IS_SET(TypeTag, int, CprMaxEllIter));
104 
105  extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
106  // For some reason simulator_.model().elementMapper() is not initialized at this stage
107  // Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work.
108  // Set it up manually
109  using ElementMapper =
110  Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
111  ElementMapper elemMapper(simulator_.gridView(), Dune::mcmgElementLayout());
112  detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
113 #if HAVE_MPI
114  if (parallelInformation_.type() == typeid(ParallelISTLInformation)) {
115  // Parallel case.
116  const ParallelISTLInformation* parinfo = std::any_cast<ParallelISTLInformation>(&parallelInformation_);
117  assert(parinfo);
118  comm_.reset(new Communication(parinfo->communicator()));
119  }
120 #endif
121  // Print parameters to PRT/DBG logs.
122  if (simulator.gridView().comm().rank() == 0) {
123  std::ostringstream os;
124  os << "Property tree for linear solver:\n";
125  prm_.write_json(os, true);
126  OpmLog::note(os.str());
127  }
128  }
129 
130  void eraseMatrix()
131  {
132  }
133 
134  void prepare(SparseMatrixAdapter& mat, VectorType& b)
135  {
136 #if HAVE_MPI
137  static bool firstcall = true;
138  if (firstcall && parallelInformation_.type() == typeid(ParallelISTLInformation)) {
139  // Parallel case.
140  const ParallelISTLInformation* parinfo = std::any_cast<ParallelISTLInformation>(&parallelInformation_);
141  assert(parinfo);
142  const size_t size = mat.istlMatrix().N();
143  parinfo->copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), size, 1);
144  firstcall = false;
145  }
146  if (isParallel() && matrixAddWellContributions_) {
147  makeOverlapRowsInvalid(mat.istlMatrix());
148  }
149 #endif
150 
151  matrix_ = &mat.istlMatrix(); // Store pointer for output if needed.
152  std::function<VectorType()> weightsCalculator = getWeightsCalculator(mat.istlMatrix(), b);
153 
154  if (shouldCreateSolver()) {
155  if (isParallel()) {
156 #if HAVE_MPI
157  if (matrixAddWellContributions_) {
158  using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Communication>;
159  auto op = std::make_unique<ParOperatorType>(mat.istlMatrix(), *comm_);
160  auto sol = std::make_unique<SolverType>(*op, *comm_, prm_, weightsCalculator,
161  pressureIndex);
162  solver_ = std::move(sol);
163  linear_operator_ = std::move(op);
164  } else {
165  if (!ownersFirst_) {
166  OPM_THROW(std::runtime_error, "In parallel, the flexible solver requires "
167  "--owner-cells-first=true when --matrix-add-well-contributions=false is used.");
168  }
170  auto well_op = std::make_unique<WellModelOpType>(simulator_.problem().wellModel());
171  auto op = std::make_unique<ParOperatorType>(mat.istlMatrix(), *well_op, interiorCellNum_);
172  auto sol = std::make_unique<SolverType>(*op, *comm_, prm_, weightsCalculator,
173  pressureIndex);
174  solver_ = std::move(sol);
175  linear_operator_ = std::move(op);
176  well_operator_ = std::move(well_op);
177  }
178 #endif // HAVE_MPI
179  } else {
180  if (matrixAddWellContributions_) {
181  using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
182  auto op = std::make_unique<SeqOperatorType>(mat.istlMatrix());
183  auto sol = std::make_unique<SolverType>(*op, prm_, weightsCalculator,
184  pressureIndex);
185  solver_ = std::move(sol);
186  linear_operator_ = std::move(op);
187  } else {
189  auto well_op = std::make_unique<WellModelOpType>(simulator_.problem().wellModel());
190  auto op = std::make_unique<SeqOperatorType>(mat.istlMatrix(), *well_op);
191  auto sol = std::make_unique<SolverType>(*op, prm_, weightsCalculator,
192  pressureIndex);
193  solver_ = std::move(sol);
194  linear_operator_ = std::move(op);
195  well_operator_ = std::move(well_op);
196  }
197  }
198  rhs_ = b;
199  } else {
200  solver_->preconditioner().update();
201  rhs_ = b;
202  }
203  }
204 
205  bool solve(VectorType& x)
206  {
207  solver_->apply(x, rhs_, res_);
208  this->writeMatrix();
209  return res_.converged;
210  }
211 
212  bool isParallel() const
213  {
214 #if HAVE_MPI
215  return parallelInformation_.type() == typeid(ParallelISTLInformation);
216 #else
217  return false;
218 #endif
219  }
220 
221  int iterations() const
222  {
223  return res_.iterations;
224  }
225 
226  void setResidual(VectorType& /* b */)
227  {
228  // rhs_ = &b; // Must be handled in prepare() instead.
229  }
230 
231  void setMatrix(const SparseMatrixAdapter& /* M */)
232  {
233  // matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead.
234  }
235 
236 protected:
237 
238  bool shouldCreateSolver() const
239  {
240  // Decide if we should recreate the solver or just do
241  // a minimal preconditioner update.
242  if (!solver_) {
243  return true;
244  }
245  const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
246  bool recreate_solver = false;
247  if (this->parameters_.cpr_reuse_setup_ == 0) {
248  // Always recreate solver.
249  recreate_solver = true;
250  } else if (this->parameters_.cpr_reuse_setup_ == 1) {
251  // Recreate solver on the first iteration of every timestep.
252  if (newton_iteration == 0) {
253  recreate_solver = true;
254  }
255  } else if (this->parameters_.cpr_reuse_setup_ == 2) {
256  // Recreate solver if the last solve used more than 10 iterations.
257  if (this->iterations() > 10) {
258  recreate_solver = true;
259  }
260  } else {
261  assert(this->parameters_.cpr_reuse_setup_ == 3);
262  assert(recreate_solver == false);
263  // Never recreate solver.
264  }
265  return recreate_solver;
266  }
267 
268  std::function<VectorType()> getWeightsCalculator(const MatrixType& mat, const VectorType& b) const
269  {
270  std::function<VectorType()> weightsCalculator;
271 
272  using namespace std::string_literals;
273 
274  auto preconditionerType = prm_.get("preconditioner.type", "cpr"s);
275  if (preconditionerType == "cpr" || preconditionerType == "cprt") {
276  const bool transpose = preconditionerType == "cprt";
277  const auto weightsType = prm_.get("preconditioner.weight_type", "quasiimpes"s);
278  if (weightsType == "quasiimpes") {
279  // weighs will be created as default in the solver
280  weightsCalculator = [&mat, transpose, p = pressureIndex]() {
281  return Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(mat, p, transpose);
282  };
283  } else if (weightsType == "trueimpes") {
284  weightsCalculator = [this, &b, p = pressureIndex]() {
285  return this->getTrueImpesWeights(b, p);
286  };
287  } else {
288  OPM_THROW(std::invalid_argument,
289  "Weights type " << weightsType << "not implemented for cpr."
290  << " Please use quasiimpes or trueimpes.");
291  }
292  }
293  return weightsCalculator;
294  }
295 
298  void makeOverlapRowsInvalid(MatrixType& matrix) const
299  {
300  //value to set on diagonal
301  const int numEq = MatrixType::block_type::rows;
302  typename MatrixType::block_type diag_block(0.0);
303  for (int eq = 0; eq < numEq; ++eq)
304  diag_block[eq][eq] = 1.0;
305 
306  //loop over precalculated overlap rows and columns
307  for (auto row = overlapRows_.begin(); row != overlapRows_.end(); row++ )
308  {
309  int lcell = *row;
310  // Zero out row.
311  matrix[lcell] = 0.0;
312 
313  //diagonal block set to diag(1.0).
314  matrix[lcell][lcell] = diag_block;
315  }
316  }
317 
318  VectorType getTrueImpesWeights(const VectorType& b, const int pressureVarIndex) const
319  {
320  VectorType weights(b.size());
321  ElementContext elemCtx(simulator_);
322  Opm::Amg::getTrueImpesWeights(pressureVarIndex, weights, simulator_.vanguard().gridView(),
323  elemCtx, simulator_.model(),
324  ThreadManager::threadId());
325  return weights;
326  }
327 
328  void writeMatrix()
329  {
330  const int verbosity = prm_.get<int>("verbosity");
331  const bool write_matrix = verbosity > 10;
332  if (write_matrix) {
333  Opm::Helper::writeSystem(this->simulator_, //simulator is only used to get names
334  *(this->matrix_),
335  this->rhs_,
336  comm_.get());
337  }
338  }
339 
340 
341  const Simulator& simulator_;
342  MatrixType* matrix_;
343  std::unique_ptr<WellModelOpType> well_operator_;
344  std::unique_ptr<AbstractOperatorType> linear_operator_;
345  std::unique_ptr<SolverType> solver_;
346  FlowLinearSolverParameters parameters_;
347  PropertyTree prm_;
348  VectorType rhs_;
349  Dune::InverseOperatorResult res_;
350  std::any parallelInformation_;
351  bool ownersFirst_;
352  bool matrixAddWellContributions_;
353  int interiorCellNum_;
354  std::unique_ptr<Communication> comm_;
355  std::vector<int> overlapRows_;
356  std::vector<int> interiorRows_;
357 }; // end ISTLSolverEbosFlexible
358 
359 } // namespace Opm
360 
361 #endif // OPM_ISTLSOLVEREBOSFLEXIBLE_HEADER_INCLUDED
A solver class that encapsulates all needed objects for a linear solver (operator,...
Definition: FlexibleSolver.hpp:40
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition: ISTLSolverEbosFlexible.hpp:59
void makeOverlapRowsInvalid(MatrixType &matrix) const
Zero out off-diagonal blocks on rows corresponding to overlap cells Diagonal blocks on ovelap rows ar...
Definition: ISTLSolverEbosFlexible.hpp:298
Linear operator wrapper for well model.
Definition: WellOperators.hpp:50
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition: WellOperators.hpp:173
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition: WellOperators.hpp:100
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: ISTLSolverEbosFlexible.hpp:37