My Project
NonlinearSolverEbos.hpp
1 /*
2  Copyright 2015 SINTEF ICT, Applied Mathematics.
3  Copyright 2015 Statoil ASA.
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef OPM_NON_LINEAR_SOLVER_EBOS_HPP
22 #define OPM_NON_LINEAR_SOLVER_EBOS_HPP
23 
24 #include <opm/simulators/timestepping/SimulatorReport.hpp>
25 #include <opm/common/ErrorMacros.hpp>
26 #include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
27 
28 #include <opm/models/utils/parametersystem.hh>
29 #include <opm/models/utils/propertysystem.hh>
30 #include <opm/models/utils/basicproperties.hh>
31 #include <opm/models/nonlinear/newtonmethodproperties.hh>
32 #include <opm/common/Exceptions.hpp>
33 
34 #include <dune/common/fmatrix.hh>
35 #include <dune/istl/bcrsmatrix.hh>
36 #include <memory>
37 
38 namespace Opm::Properties {
39 
40 namespace TTag {
42 }
43 
44 template<class TypeTag, class MyTypeTag>
46  using type = UndefinedProperty;
47 };
48 
49 // we are reusing NewtonMaxIterations from opm-models
50 // See opm/models/nonlinear/newtonmethodproperties.hh
51 
52 template<class TypeTag, class MyTypeTag>
54  using type = UndefinedProperty;
55 };
56 template<class TypeTag, class MyTypeTag>
58  using type = UndefinedProperty;
59 };
60 
61 template<class TypeTag>
62 struct NewtonMaxRelax<TypeTag, TTag::FlowNonLinearSolver> {
63  using type = GetPropType<TypeTag, Scalar>;
64  static constexpr type value = 0.5;
65 };
66 template<class TypeTag>
67 struct NewtonMaxIterations<TypeTag, TTag::FlowNonLinearSolver> {
68  static constexpr int value = 20;
69 };
70 template<class TypeTag>
71 struct NewtonMinIterations<TypeTag, TTag::FlowNonLinearSolver> {
72  static constexpr int value = 1;
73 };
74 template<class TypeTag>
75 struct NewtonRelaxationType<TypeTag, TTag::FlowNonLinearSolver> {
76  static constexpr auto value = "dampen";
77 };
78 
79 } // namespace Opm::Properties
80 
81 namespace Opm {
82 
83 class WellState;
84 
87  template <class TypeTag, class PhysicalModel>
89  {
90  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
91 
92  public:
93  // Available relaxation scheme types.
94  enum RelaxType {
95  Dampen,
96  SOR
97  };
98 
99  // Solver parameters controlling nonlinear process.
101  {
102  RelaxType relaxType_;
103  double relaxMax_;
104  double relaxIncrement_;
105  double relaxRelTol_;
106  int maxIter_; // max nonlinear iterations
107  int minIter_; // min nonlinear iterations
108 
110  {
111  // set default values
112  reset();
113 
114  // overload with given parameters
115  relaxMax_ = EWOMS_GET_PARAM(TypeTag, Scalar, NewtonMaxRelax);
116  maxIter_ = EWOMS_GET_PARAM(TypeTag, int, NewtonMaxIterations);
117  minIter_ = EWOMS_GET_PARAM(TypeTag, int, NewtonMinIterations);
118 
119  const auto& relaxationTypeString = EWOMS_GET_PARAM(TypeTag, std::string, NewtonRelaxationType);
120  if (relaxationTypeString == "dampen") {
121  relaxType_ = Dampen;
122  } else if (relaxationTypeString == "sor") {
123  relaxType_ = SOR;
124  } else {
125  OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxationTypeString);
126  }
127  }
128 
129  static void registerParameters()
130  {
131  EWOMS_REGISTER_PARAM(TypeTag, Scalar, NewtonMaxRelax, "The maximum relaxation factor of a Newton iteration");
132  EWOMS_REGISTER_PARAM(TypeTag, int, NewtonMaxIterations, "The maximum number of Newton iterations per time step");
133  EWOMS_REGISTER_PARAM(TypeTag, int, NewtonMinIterations, "The minimum number of Newton iterations per time step");
134  EWOMS_REGISTER_PARAM(TypeTag, std::string, NewtonRelaxationType, "The type of relaxation used by Newton method");
135  }
136 
137  void reset()
138  {
139  // default values for the solver parameters
140  relaxType_ = Dampen;
141  relaxMax_ = 0.5;
142  relaxIncrement_ = 0.1;
143  relaxRelTol_ = 0.2;
144  maxIter_ = 10;
145  minIter_ = 1;
146  }
147 
148  };
149 
150  // Forwarding types from PhysicalModel.
151  //typedef typename PhysicalModel::WellState WellState;
152 
153  // --------- Public methods ---------
154 
163  std::unique_ptr<PhysicalModel> model)
164  : param_(param)
165  , model_(std::move(model))
166  , linearizations_(0)
167  , nonlinearIterations_(0)
168  , linearIterations_(0)
169  , wellIterations_(0)
170  , nonlinearIterationsLast_(0)
171  , linearIterationsLast_(0)
172  , wellIterationsLast_(0)
173  {
174  if (!model_) {
175  OPM_THROW(std::logic_error, "Must provide a non-null model argument for NonlinearSolver.");
176  }
177  }
178 
179 
181  {
182  SimulatorReportSingle report;
183  report.global_time = timer.simulationTimeElapsed();
184  report.timestep_length = timer.currentStepLength();
185 
186  // Do model-specific once-per-step calculations.
187  report += model_->prepareStep(timer);
188 
189  int iteration = 0;
190 
191  // Let the model do one nonlinear iteration.
192 
193  // Set up for main solver loop.
194  bool converged = false;
195 
196  // ---------- Main nonlinear solver loop ----------
197  do {
198  try {
199  // Do the nonlinear step. If we are in a converged state, the
200  // model will usually do an early return without an expensive
201  // solve, unless the minIter() count has not been reached yet.
202  auto iterReport = model_->nonlinearIteration(iteration, timer, *this);
203  iterReport.global_time = timer.simulationTimeElapsed();
204  report += iterReport;
205  report.converged = iterReport.converged;
206 
207  converged = report.converged;
208  iteration += 1;
209  }
210  catch (...) {
211  // if an iteration fails during a time step, all previous iterations
212  // count as a failure as well
213  failureReport_ = report;
214  failureReport_ += model_->failureReport();
215  throw;
216  }
217  }
218  while ( (!converged && (iteration <= maxIter())) || (iteration <= minIter()));
219 
220  if (!converged) {
221  failureReport_ = report;
222 
223  std::string msg = "Solver convergence failure - Failed to complete a time step within " + std::to_string(maxIter()) + " iterations.";
224  OPM_THROW_NOLOG(TooManyIterations, msg);
225  }
226 
227  // Do model-specific post-step actions.
228  report += model_->afterStep(timer);
229  report.converged = true;
230  return report;
231  }
232 
235  { return failureReport_; }
236 
238  int linearizations() const
239  { return linearizations_; }
240 
243  { return nonlinearIterations_; }
244 
246  int linearIterations() const
247  { return linearIterations_; }
248 
250  int wellIterations() const
251  { return wellIterations_; }
252 
255  { return nonlinearIterationsLast_; }
256 
259  { return linearIterationsLast_; }
260 
263  { return wellIterationsLast_; }
264 
265  std::vector<std::vector<double> >
266  computeFluidInPlace(const std::vector<int>& fipnum) const
267  { return model_->computeFluidInPlace(fipnum); }
268 
270  const PhysicalModel& model() const
271  { return *model_; }
272 
274  PhysicalModel& model()
275  { return *model_; }
276 
278  void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
279  const int it, bool& oscillate, bool& stagnate) const
280  {
281  // The detection of oscillation in two primary variable results in the report of the detection
282  // of oscillation for the solver.
283  // Only the saturations are used for oscillation detection for the black oil model.
284  // Stagnate is not used for any treatment here.
285 
286  if ( it < 2 ) {
287  oscillate = false;
288  stagnate = false;
289  return;
290  }
291 
292  stagnate = true;
293  int oscillatePhase = 0;
294  const std::vector<double>& F0 = residualHistory[it];
295  const std::vector<double>& F1 = residualHistory[it - 1];
296  const std::vector<double>& F2 = residualHistory[it - 2];
297  for (int p= 0; p < model_->numPhases(); ++p){
298  const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
299  const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);
300 
301  oscillatePhase += (d1 < relaxRelTol()) && (relaxRelTol() < d2);
302 
303  // Process is 'stagnate' unless at least one phase
304  // exhibits significant residual change.
305  stagnate = (stagnate && !(std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3));
306  }
307 
308  oscillate = (oscillatePhase > 1);
309  }
310 
311 
314  template <class BVector>
315  void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
316  {
317  // The dxOld is updated with dx.
318  // If omega is equal to 1., no relaxtion will be appiled.
319 
320  BVector tempDxOld = dxOld;
321  dxOld = dx;
322 
323  switch (relaxType()) {
324  case Dampen: {
325  if (omega == 1.) {
326  return;
327  }
328  auto i = dx.size();
329  for (i = 0; i < dx.size(); ++i) {
330  dx[i] *= omega;
331  }
332  return;
333  }
334  case SOR: {
335  if (omega == 1.) {
336  return;
337  }
338  auto i = dx.size();
339  for (i = 0; i < dx.size(); ++i) {
340  dx[i] *= omega;
341  tempDxOld[i] *= (1.-omega);
342  dx[i] += tempDxOld[i];
343  }
344  return;
345  }
346  default:
347  OPM_THROW(std::runtime_error, "Can only handle Dampen and SOR relaxation type.");
348  }
349 
350  return;
351  }
352 
354  double relaxMax() const
355  { return param_.relaxMax_; }
356 
358  double relaxIncrement() const
359  { return param_.relaxIncrement_; }
360 
362  enum RelaxType relaxType() const
363  { return param_.relaxType_; }
364 
366  double relaxRelTol() const
367  { return param_.relaxRelTol_; }
368 
370  int maxIter() const
371  { return param_.maxIter_; }
372 
374  int minIter() const
375  { return param_.minIter_; }
376 
378  void setParameters(const SolverParameters& param)
379  { param_ = param; }
380 
381  private:
382  // --------- Data members ---------
383  SimulatorReportSingle failureReport_;
384  SolverParameters param_;
385  std::unique_ptr<PhysicalModel> model_;
386  int linearizations_;
387  int nonlinearIterations_;
388  int linearIterations_;
389  int wellIterations_;
390  int nonlinearIterationsLast_;
391  int linearIterationsLast_;
392  int wellIterationsLast_;
393  };
394 } // namespace Opm
395 
396 #endif // OPM_NON_LINEAR_SOLVER_EBOS_HPP
A nonlinear solver class suitable for general fully-implicit models, as well as pressure,...
Definition: NonlinearSolverEbos.hpp:89
double relaxMax() const
The greatest relaxation factor (i.e. smallest factor) allowed.
Definition: NonlinearSolverEbos.hpp:354
int maxIter() const
The maximum number of nonlinear iterations allowed.
Definition: NonlinearSolverEbos.hpp:370
void detectOscillations(const std::vector< std::vector< double >> &residualHistory, const int it, bool &oscillate, bool &stagnate) const
Detect oscillation or stagnation in a given residual history.
Definition: NonlinearSolverEbos.hpp:278
PhysicalModel & model()
Mutable reference to physical model.
Definition: NonlinearSolverEbos.hpp:274
int wellIterationsLastStep() const
Number of well iterations used in all calls to step().
Definition: NonlinearSolverEbos.hpp:262
const SimulatorReportSingle & failureReport() const
return the statistics if the step() method failed
Definition: NonlinearSolverEbos.hpp:234
double relaxIncrement() const
The step-change size for the relaxation factor.
Definition: NonlinearSolverEbos.hpp:358
int linearIterationsLastStep() const
Number of linear solver iterations used in the last call to step().
Definition: NonlinearSolverEbos.hpp:258
int linearIterations() const
Number of linear solver iterations used in all calls to step().
Definition: NonlinearSolverEbos.hpp:246
enum RelaxType relaxType() const
The relaxation type (Dampen or SOR).
Definition: NonlinearSolverEbos.hpp:362
int linearizations() const
Number of linearizations used in all calls to step().
Definition: NonlinearSolverEbos.hpp:238
void stabilizeNonlinearUpdate(BVector &dx, BVector &dxOld, const double omega) const
Apply a stabilization to dx, depending on dxOld and relaxation parameters.
Definition: NonlinearSolverEbos.hpp:315
int minIter() const
The minimum number of nonlinear iterations allowed.
Definition: NonlinearSolverEbos.hpp:374
int nonlinearIterationsLastStep() const
Number of nonlinear solver iterations used in the last call to step().
Definition: NonlinearSolverEbos.hpp:254
double relaxRelTol() const
The relaxation relative tolerance.
Definition: NonlinearSolverEbos.hpp:366
int wellIterations() const
Number of well iterations used in all calls to step().
Definition: NonlinearSolverEbos.hpp:250
const PhysicalModel & model() const
Reference to physical model.
Definition: NonlinearSolverEbos.hpp:270
NonlinearSolverEbos(const SolverParameters &param, std::unique_ptr< PhysicalModel > model)
Construct solver for a given model.
Definition: NonlinearSolverEbos.hpp:162
void setParameters(const SolverParameters &param)
Set parameters to override those given at construction time.
Definition: NonlinearSolverEbos.hpp:378
int nonlinearIterations() const
Number of full nonlinear solver iterations used in all calls to step().
Definition: NonlinearSolverEbos.hpp:242
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:34
virtual double currentStepLength() const =0
Current step length.
virtual double simulationTimeElapsed() const =0
Time elapsed since the start of the simulation until the beginning of the current time step [s].
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellState.hpp:56
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
Definition: NonlinearSolverEbos.hpp:101
Definition: NonlinearSolverEbos.hpp:45
Definition: NonlinearSolverEbos.hpp:53
Definition: NonlinearSolverEbos.hpp:57
Definition: NonlinearSolverEbos.hpp:41
A struct for returning timing data from a simulator to its caller.
Definition: SimulatorReport.hpp:31