My Project
Loading...
Searching...
No Matches
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>
56template<class TypeTag, class MyTypeTag>
60
61template<class TypeTag>
62struct NewtonMaxRelax<TypeTag, TTag::FlowNonLinearSolver> {
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
85// Available relaxation scheme types.
86enum class NonlinearRelaxType {
87 Dampen,
88 SOR
89};
90
91namespace detail {
92
94void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
95 const int it, const int numPhases, const double relaxRelTol,
96 bool& oscillate, bool& stagnate);
97
100template <class BVector>
101void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
102 const double omega, NonlinearRelaxType relaxType);
103
104}
105
108 template <class TypeTag, class PhysicalModel>
110 {
112
113 public:
114 // Solver parameters controlling nonlinear process.
116 {
117 NonlinearRelaxType relaxType_;
118 double relaxMax_;
119 double relaxIncrement_;
120 double relaxRelTol_;
121 int maxIter_; // max nonlinear iterations
122 int minIter_; // min nonlinear iterations
123
125 {
126 // set default values
127 reset();
128
129 // overload with given parameters
130 relaxMax_ = EWOMS_GET_PARAM(TypeTag, Scalar, NewtonMaxRelax);
131 maxIter_ = EWOMS_GET_PARAM(TypeTag, int, NewtonMaxIterations);
132 minIter_ = EWOMS_GET_PARAM(TypeTag, int, NewtonMinIterations);
133
134 const auto& relaxationTypeString = EWOMS_GET_PARAM(TypeTag, std::string, NewtonRelaxationType);
135 if (relaxationTypeString == "dampen") {
136 relaxType_ = NonlinearRelaxType::Dampen;
137 } else if (relaxationTypeString == "sor") {
138 relaxType_ = NonlinearRelaxType::SOR;
139 } else {
140 OPM_THROW(std::runtime_error,
141 "Unknown Relaxtion Type " + relaxationTypeString);
142 }
143 }
144
145 static void registerParameters()
146 {
147 EWOMS_REGISTER_PARAM(TypeTag, Scalar, NewtonMaxRelax, "The maximum relaxation factor of a Newton iteration");
148 EWOMS_REGISTER_PARAM(TypeTag, int, NewtonMaxIterations, "The maximum number of Newton iterations per time step");
149 EWOMS_REGISTER_PARAM(TypeTag, int, NewtonMinIterations, "The minimum number of Newton iterations per time step");
150 EWOMS_REGISTER_PARAM(TypeTag, std::string, NewtonRelaxationType, "The type of relaxation used by Newton method");
151 }
152
153 void reset()
154 {
155 // default values for the solver parameters
156 relaxType_ = NonlinearRelaxType::Dampen;
157 relaxMax_ = 0.5;
158 relaxIncrement_ = 0.1;
159 relaxRelTol_ = 0.2;
160 maxIter_ = 10;
161 minIter_ = 1;
162 }
163
164 };
165
166 // Forwarding types from PhysicalModel.
167 //typedef typename PhysicalModel::WellState WellState;
168
169 // --------- Public methods ---------
170
179 std::unique_ptr<PhysicalModel> model)
180 : param_(param)
181 , model_(std::move(model))
182 , linearizations_(0)
183 , nonlinearIterations_(0)
184 , linearIterations_(0)
185 , wellIterations_(0)
186 , nonlinearIterationsLast_(0)
187 , linearIterationsLast_(0)
188 , wellIterationsLast_(0)
189 {
190 if (!model_) {
191 OPM_THROW(std::logic_error, "Must provide a non-null model argument for NonlinearSolver.");
192 }
193 }
194
195
197 {
199 report.global_time = timer.simulationTimeElapsed();
200 report.timestep_length = timer.currentStepLength();
201
202 // Do model-specific once-per-step calculations.
203 report += model_->prepareStep(timer);
204
205 int iteration = 0;
206
207 // Let the model do one nonlinear iteration.
208
209 // Set up for main solver loop.
210 bool converged = false;
211
212 // ---------- Main nonlinear solver loop ----------
213 do {
214 try {
215 // Do the nonlinear step. If we are in a converged state, the
216 // model will usually do an early return without an expensive
217 // solve, unless the minIter() count has not been reached yet.
218 auto iterReport = model_->nonlinearIteration(iteration, timer, *this);
219 iterReport.global_time = timer.simulationTimeElapsed();
220 report += iterReport;
221 report.converged = iterReport.converged;
222
223 converged = report.converged;
224 iteration += 1;
225 }
226 catch (...) {
227 // if an iteration fails during a time step, all previous iterations
228 // count as a failure as well
229 failureReport_ = report;
230 failureReport_ += model_->failureReport();
231 throw;
232 }
233 }
234 while ( (!converged && (iteration <= maxIter())) || (iteration <= minIter()));
235
236 if (!converged) {
237 failureReport_ = report;
238
239 std::string msg = "Solver convergence failure - Failed to complete a time step within " + std::to_string(maxIter()) + " iterations.";
240 OPM_THROW_NOLOG(TooManyIterations, msg);
241 }
242
243 // Do model-specific post-step actions.
244 report += model_->afterStep(timer);
245 report.converged = true;
246 return report;
247 }
248
251 { return failureReport_; }
252
254 int linearizations() const
255 { return linearizations_; }
256
259 { return nonlinearIterations_; }
260
263 { return linearIterations_; }
264
266 int wellIterations() const
267 { return wellIterations_; }
268
271 { return nonlinearIterationsLast_; }
272
275 { return linearIterationsLast_; }
276
279 { return wellIterationsLast_; }
280
281 std::vector<std::vector<double> >
282 computeFluidInPlace(const std::vector<int>& fipnum) const
283 { return model_->computeFluidInPlace(fipnum); }
284
286 const PhysicalModel& model() const
287 { return *model_; }
288
291 { return *model_; }
292
294 void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
295 const int it, bool& oscillate, bool& stagnate) const
296 {
297 detail::detectOscillations(residualHistory, it, model_->numPhases(),
298 this->relaxRelTol(), oscillate, stagnate);
299 }
300
301
304 template <class BVector>
305 void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
306 {
307 detail::stabilizeNonlinearUpdate(dx, dxOld, omega, this->relaxType());
308 }
309
311 double relaxMax() const
312 { return param_.relaxMax_; }
313
315 double relaxIncrement() const
316 { return param_.relaxIncrement_; }
317
319 NonlinearRelaxType relaxType() const
320 { return param_.relaxType_; }
321
323 double relaxRelTol() const
324 { return param_.relaxRelTol_; }
325
327 int maxIter() const
328 { return param_.maxIter_; }
329
331 int minIter() const
332 { return param_.minIter_; }
333
336 { param_ = param; }
337
338 private:
339 // --------- Data members ---------
340 SimulatorReportSingle failureReport_;
341 SolverParameters param_;
342 std::unique_ptr<PhysicalModel> model_;
343 int linearizations_;
344 int nonlinearIterations_;
345 int linearIterations_;
346 int wellIterations_;
347 int nonlinearIterationsLast_;
348 int linearIterationsLast_;
349 int wellIterationsLast_;
350 };
351
352} // namespace Opm
353
354#endif // OPM_NON_LINEAR_SOLVER_EBOS_HPP
Definition AquiferInterface.hpp:35
A nonlinear solver class suitable for general fully-implicit models, as well as pressure,...
Definition NonlinearSolverEbos.hpp:110
double relaxMax() const
The greatest relaxation factor (i.e. smallest factor) allowed.
Definition NonlinearSolverEbos.hpp:311
int maxIter() const
The maximum number of nonlinear iterations allowed.
Definition NonlinearSolverEbos.hpp:327
const SimulatorReportSingle & failureReport() const
return the statistics if the step() method failed
Definition NonlinearSolverEbos.hpp:250
PhysicalModel & model()
Mutable reference to physical model.
Definition NonlinearSolverEbos.hpp:290
const PhysicalModel & model() const
Reference to physical model.
Definition NonlinearSolverEbos.hpp:286
int wellIterationsLastStep() const
Number of well iterations used in all calls to step().
Definition NonlinearSolverEbos.hpp:278
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:294
double relaxIncrement() const
The step-change size for the relaxation factor.
Definition NonlinearSolverEbos.hpp:315
int linearIterationsLastStep() const
Number of linear solver iterations used in the last call to step().
Definition NonlinearSolverEbos.hpp:274
int linearIterations() const
Number of linear solver iterations used in all calls to step().
Definition NonlinearSolverEbos.hpp:262
int linearizations() const
Number of linearizations used in all calls to step().
Definition NonlinearSolverEbos.hpp:254
void stabilizeNonlinearUpdate(BVector &dx, BVector &dxOld, const double omega) const
Apply a stabilization to dx, depending on dxOld and relaxation parameters.
Definition NonlinearSolverEbos.hpp:305
int minIter() const
The minimum number of nonlinear iterations allowed.
Definition NonlinearSolverEbos.hpp:331
int nonlinearIterationsLastStep() const
Number of nonlinear solver iterations used in the last call to step().
Definition NonlinearSolverEbos.hpp:270
double relaxRelTol() const
The relaxation relative tolerance.
Definition NonlinearSolverEbos.hpp:323
int wellIterations() const
Number of well iterations used in all calls to step().
Definition NonlinearSolverEbos.hpp:266
NonlinearSolverEbos(const SolverParameters &param, std::unique_ptr< PhysicalModel > model)
Construct solver for a given model.
Definition NonlinearSolverEbos.hpp:178
void setParameters(const SolverParameters &param)
Set parameters to override those given at construction time.
Definition NonlinearSolverEbos.hpp:335
NonlinearRelaxType relaxType() const
The relaxation type (Dampen or SOR).
Definition NonlinearSolverEbos.hpp:319
int nonlinearIterations() const
Number of full nonlinear solver iterations used in all calls to step().
Definition NonlinearSolverEbos.hpp:258
Interface class for SimulatorTimer objects, to be improved.
Definition SimulatorTimerInterface.hpp:34
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:60
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
Definition NonlinearSolverEbos.hpp:116
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:34