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/common/Exceptions.hpp>
32
33#include <dune/common/fmatrix.hh>
34#include <dune/istl/bcrsmatrix.hh>
35#include <memory>
36
37namespace Opm::Properties {
38
39namespace TTag {
41}
42
43template<class TypeTag, class MyTypeTag>
45 using type = UndefinedProperty;
46};
47template<class TypeTag, class MyTypeTag>
49 using type = UndefinedProperty;
50};
51template<class TypeTag, class MyTypeTag>
53 using type = UndefinedProperty;
54};
55template<class TypeTag, class MyTypeTag>
57 using type = UndefinedProperty;
58};
59
60template<class TypeTag>
61struct NewtonMaxRelax<TypeTag, TTag::FlowNonLinearSolver> {
62 using type = GetPropType<TypeTag, Scalar>;
63 static constexpr type value = 0.5;
64};
65template<class TypeTag>
66struct FlowNewtonMaxIterations<TypeTag, TTag::FlowNonLinearSolver> {
67 static constexpr int value = 20;
68};
69template<class TypeTag>
70struct NewtonMinIterations<TypeTag, TTag::FlowNonLinearSolver> {
71 static constexpr int value = 1;
72};
73template<class TypeTag>
74struct NewtonRelaxationType<TypeTag, TTag::FlowNonLinearSolver> {
75 static constexpr auto value = "dampen";
76};
77
78} // namespace Opm::Properties
79
80namespace Opm {
81
82class WellState;
83
86 template <class TypeTag, class PhysicalModel>
88 {
89 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
90
91 public:
92 // Available relaxation scheme types.
93 enum RelaxType {
94 Dampen,
95 SOR
96 };
97
98 // Solver parameters controlling nonlinear process.
100 {
101 RelaxType relaxType_;
102 double relaxMax_;
103 double relaxIncrement_;
104 double relaxRelTol_;
105 int maxIter_; // max nonlinear iterations
106 int minIter_; // min nonlinear iterations
107
109 {
110 // set default values
111 reset();
112
113 // overload with given parameters
114 relaxMax_ = EWOMS_GET_PARAM(TypeTag, Scalar, NewtonMaxRelax);
115 maxIter_ = EWOMS_GET_PARAM(TypeTag, int, FlowNewtonMaxIterations);
116 minIter_ = EWOMS_GET_PARAM(TypeTag, int, NewtonMinIterations);
117
118 const auto& relaxationTypeString = EWOMS_GET_PARAM(TypeTag, std::string, NewtonRelaxationType);
119 if (relaxationTypeString == "dampen") {
120 relaxType_ = Dampen;
121 } else if (relaxationTypeString == "sor") {
122 relaxType_ = SOR;
123 } else {
124 OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxationTypeString);
125 }
126 }
127
128 static void registerParameters()
129 {
130 EWOMS_REGISTER_PARAM(TypeTag, Scalar, NewtonMaxRelax, "The maximum relaxation factor of a Newton iteration used by flow");
131 EWOMS_REGISTER_PARAM(TypeTag, int, FlowNewtonMaxIterations, "The maximum number of Newton iterations per time step used by flow");
132 EWOMS_REGISTER_PARAM(TypeTag, int, NewtonMinIterations, "The minimum number of Newton iterations per time step used by flow");
133 EWOMS_REGISTER_PARAM(TypeTag, std::string, NewtonRelaxationType, "The type of relaxation used by flow's Newton method");
134 }
135
136 void reset()
137 {
138 // default values for the solver parameters
139 relaxType_ = Dampen;
140 relaxMax_ = 0.5;
141 relaxIncrement_ = 0.1;
142 relaxRelTol_ = 0.2;
143 maxIter_ = 10;
144 minIter_ = 1;
145 }
146
147 };
148
149 // Forwarding types from PhysicalModel.
150 //typedef typename PhysicalModel::WellState WellState;
151
152 // --------- Public methods ---------
153
162 std::unique_ptr<PhysicalModel> model)
163 : param_(param)
164 , model_(std::move(model))
165 , linearizations_(0)
166 , nonlinearIterations_(0)
167 , linearIterations_(0)
168 , wellIterations_(0)
169 , nonlinearIterationsLast_(0)
170 , linearIterationsLast_(0)
171 , wellIterationsLast_(0)
172 {
173 if (!model_) {
174 OPM_THROW(std::logic_error, "Must provide a non-null model argument for NonlinearSolver.");
175 }
176 }
177
178
180 {
182 report.global_time = timer.simulationTimeElapsed();
183 report.timestep_length = timer.currentStepLength();
184
185 // Do model-specific once-per-step calculations.
186 report += model_->prepareStep(timer);
187
188 int iteration = 0;
189
190 // Let the model do one nonlinear iteration.
191
192 // Set up for main solver loop.
193 bool converged = false;
194
195 // ---------- Main nonlinear solver loop ----------
196 do {
197 try {
198 // Do the nonlinear step. If we are in a converged state, the
199 // model will usually do an early return without an expensive
200 // solve, unless the minIter() count has not been reached yet.
201 auto iterReport = model_->nonlinearIteration(iteration, timer, *this);
202 iterReport.global_time = timer.simulationTimeElapsed();
203 report += iterReport;
204 report.converged = iterReport.converged;
205
206 converged = report.converged;
207 iteration += 1;
208 }
209 catch (...) {
210 // if an iteration fails during a time step, all previous iterations
211 // count as a failure as well
212 failureReport_ = report;
213 failureReport_ += model_->failureReport();
214 throw;
215 }
216 }
217 while ( (!converged && (iteration <= maxIter())) || (iteration <= minIter()));
218
219 if (!converged) {
220 failureReport_ = report;
221
222 std::string msg = "Solver convergence failure - Failed to complete a time step within " + std::to_string(maxIter()) + " iterations.";
223 OPM_THROW_NOLOG(TooManyIterations, msg);
224 }
225
226 // Do model-specific post-step actions.
227 report += model_->afterStep(timer);
228 report.converged = true;
229 return report;
230 }
231
234 { return failureReport_; }
235
237 int linearizations() const
238 { return linearizations_; }
239
242 { return nonlinearIterations_; }
243
246 { return linearIterations_; }
247
249 int wellIterations() const
250 { return wellIterations_; }
251
254 { return nonlinearIterationsLast_; }
255
258 { return linearIterationsLast_; }
259
262 { return wellIterationsLast_; }
263
264 std::vector<std::vector<double> >
265 computeFluidInPlace(const std::vector<int>& fipnum) const
266 { return model_->computeFluidInPlace(fipnum); }
267
269 const PhysicalModel& model() const
270 { return *model_; }
271
273 PhysicalModel& model()
274 { return *model_; }
275
277 void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
278 const int it, bool& oscillate, bool& stagnate) const
279 {
280 // The detection of oscillation in two primary variable results in the report of the detection
281 // of oscillation for the solver.
282 // Only the saturations are used for oscillation detection for the black oil model.
283 // Stagnate is not used for any treatment here.
284
285 if ( it < 2 ) {
286 oscillate = false;
287 stagnate = false;
288 return;
289 }
290
291 stagnate = true;
292 int oscillatePhase = 0;
293 const std::vector<double>& F0 = residualHistory[it];
294 const std::vector<double>& F1 = residualHistory[it - 1];
295 const std::vector<double>& F2 = residualHistory[it - 2];
296 for (int p= 0; p < model_->numPhases(); ++p){
297 const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
298 const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);
299
300 oscillatePhase += (d1 < relaxRelTol()) && (relaxRelTol() < d2);
301
302 // Process is 'stagnate' unless at least one phase
303 // exhibits significant residual change.
304 stagnate = (stagnate && !(std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3));
305 }
306
307 oscillate = (oscillatePhase > 1);
308 }
309
310
313 template <class BVector>
314 void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
315 {
316 // The dxOld is updated with dx.
317 // If omega is equal to 1., no relaxtion will be appiled.
318
319 BVector tempDxOld = dxOld;
320 dxOld = dx;
321
322 switch (relaxType()) {
323 case Dampen: {
324 if (omega == 1.) {
325 return;
326 }
327 auto i = dx.size();
328 for (i = 0; i < dx.size(); ++i) {
329 dx[i] *= omega;
330 }
331 return;
332 }
333 case SOR: {
334 if (omega == 1.) {
335 return;
336 }
337 auto i = dx.size();
338 for (i = 0; i < dx.size(); ++i) {
339 dx[i] *= omega;
340 tempDxOld[i] *= (1.-omega);
341 dx[i] += tempDxOld[i];
342 }
343 return;
344 }
345 default:
346 OPM_THROW(std::runtime_error, "Can only handle Dampen and SOR relaxation type.");
347 }
348
349 return;
350 }
351
353 double relaxMax() const
354 { return param_.relaxMax_; }
355
357 double relaxIncrement() const
358 { return param_.relaxIncrement_; }
359
361 enum RelaxType relaxType() const
362 { return param_.relaxType_; }
363
365 double relaxRelTol() const
366 { return param_.relaxRelTol_; }
367
369 int maxIter() const
370 { return param_.maxIter_; }
371
373 int minIter() const
374 { return param_.minIter_; }
375
378 { param_ = param; }
379
380 private:
381 // --------- Data members ---------
382 SimulatorReportSingle failureReport_;
383 SolverParameters param_;
384 std::unique_ptr<PhysicalModel> model_;
385 int linearizations_;
386 int nonlinearIterations_;
387 int linearIterations_;
388 int wellIterations_;
389 int nonlinearIterationsLast_;
390 int linearIterationsLast_;
391 int wellIterationsLast_;
392 };
393} // namespace Opm
394
395#endif // OPM_NON_LINEAR_SOLVER_EBOS_HPP
A nonlinear solver class suitable for general fully-implicit models, as well as pressure,...
Definition: NonlinearSolverEbos.hpp:88
double relaxMax() const
The greatest relaxation factor (i.e. smallest factor) allowed.
Definition: NonlinearSolverEbos.hpp:353
int maxIter() const
The maximum number of nonlinear iterations allowed.
Definition: NonlinearSolverEbos.hpp:369
const SimulatorReportSingle & failureReport() const
return the statistics if the step() method failed
Definition: NonlinearSolverEbos.hpp:233
PhysicalModel & model()
Mutable reference to physical model.
Definition: NonlinearSolverEbos.hpp:273
const PhysicalModel & model() const
Reference to physical model.
Definition: NonlinearSolverEbos.hpp:269
int wellIterationsLastStep() const
Number of well iterations used in all calls to step().
Definition: NonlinearSolverEbos.hpp:261
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:277
double relaxIncrement() const
The step-change size for the relaxation factor.
Definition: NonlinearSolverEbos.hpp:357
int linearIterationsLastStep() const
Number of linear solver iterations used in the last call to step().
Definition: NonlinearSolverEbos.hpp:257
int linearIterations() const
Number of linear solver iterations used in all calls to step().
Definition: NonlinearSolverEbos.hpp:245
enum RelaxType relaxType() const
The relaxation type (Dampen or SOR).
Definition: NonlinearSolverEbos.hpp:361
int linearizations() const
Number of linearizations used in all calls to step().
Definition: NonlinearSolverEbos.hpp:237
void stabilizeNonlinearUpdate(BVector &dx, BVector &dxOld, const double omega) const
Apply a stabilization to dx, depending on dxOld and relaxation parameters.
Definition: NonlinearSolverEbos.hpp:314
int minIter() const
The minimum number of nonlinear iterations allowed.
Definition: NonlinearSolverEbos.hpp:373
int nonlinearIterationsLastStep() const
Number of nonlinear solver iterations used in the last call to step().
Definition: NonlinearSolverEbos.hpp:253
double relaxRelTol() const
The relaxation relative tolerance.
Definition: NonlinearSolverEbos.hpp:365
int wellIterations() const
Number of well iterations used in all calls to step().
Definition: NonlinearSolverEbos.hpp:249
NonlinearSolverEbos(const SolverParameters &param, std::unique_ptr< PhysicalModel > model)
Construct solver for a given model.
Definition: NonlinearSolverEbos.hpp:161
void setParameters(const SolverParameters &param)
Set parameters to override those given at construction time.
Definition: NonlinearSolverEbos.hpp:377
int nonlinearIterations() const
Number of full nonlinear solver iterations used in all calls to step().
Definition: NonlinearSolverEbos.hpp:241
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:100
Definition: NonlinearSolverEbos.hpp:48
Definition: NonlinearSolverEbos.hpp:44
Definition: NonlinearSolverEbos.hpp:52
Definition: NonlinearSolverEbos.hpp:56
Definition: NonlinearSolverEbos.hpp:40
A struct for returning timing data from a simulator to its caller.
Definition: SimulatorReport.hpp:31