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
38namespace Opm::Properties {
39
40namespace TTag {
42}
43
44template<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
52template<class TypeTag, class MyTypeTag>
54 using type = UndefinedProperty;
55};
56template<class TypeTag, class MyTypeTag>
58 using type = UndefinedProperty;
59};
60
61template<class TypeTag>
62struct NewtonMaxRelax<TypeTag, TTag::FlowNonLinearSolver> {
63 using type = GetPropType<TypeTag, Scalar>;
64 static constexpr type value = 0.5;
65};
66template<class TypeTag>
67struct NewtonMaxIterations<TypeTag, TTag::FlowNonLinearSolver> {
68 static constexpr int value = 20;
69};
70template<class TypeTag>
71struct NewtonMinIterations<TypeTag, TTag::FlowNonLinearSolver> {
72 static constexpr int value = 1;
73};
74template<class TypeTag>
75struct NewtonRelaxationType<TypeTag, TTag::FlowNonLinearSolver> {
76 static constexpr auto value = "dampen";
77};
78
79} // namespace Opm::Properties
80
81namespace Opm {
82
83class 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 {
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
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
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
const SimulatorReportSingle & failureReport() const
return the statistics if the step() method failed
Definition: NonlinearSolverEbos.hpp:234
PhysicalModel & model()
Mutable reference to physical model.
Definition: NonlinearSolverEbos.hpp:274
const PhysicalModel & model() const
Reference to physical model.
Definition: NonlinearSolverEbos.hpp:270
int wellIterationsLastStep() const
Number of well iterations used in all calls to step().
Definition: NonlinearSolverEbos.hpp:262
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
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
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