NonlinearSolver.hpp
Go to the documentation of this file.
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_HPP
22#define OPM_NON_LINEAR_SOLVER_HPP
23
24#include <dune/common/fmatrix.hh>
25#include <dune/istl/bcrsmatrix.hh>
26
27#include <opm/common/ErrorMacros.hpp>
28#include <opm/common/Exceptions.hpp>
29
32
35
39
40#include <memory>
41
42namespace Opm::Parameters {
43
44template<class Scalar>
45struct NewtonMaxRelax { static constexpr Scalar value = 0.5; };
46
47struct NewtonRelaxationType { static constexpr auto value = "dampen"; };
48
49} // namespace Opm::Parameters
50
51namespace Opm {
52
53// Available relaxation scheme types.
55 Dampen,
56 SOR
57};
58
59namespace detail {
60
62template<class Scalar>
63void detectOscillations(const std::vector<std::vector<Scalar>>& residualHistory,
64 const int it, const int numPhases, const Scalar relaxRelTol,
65 const int minimumOscillatingPhases,
66 bool& oscillate, bool& stagnate);
67
70template <class BVector, class Scalar>
71void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
72 const Scalar omega, NonlinearRelaxType relaxType);
73
74}
75
76// Solver parameters controlling nonlinear process.
77template<class Scalar>
79{
81 Scalar relaxMax_;
84
86
87 static void registerParameters();
88
89 void reset();
90};
91
94 template <class TypeTag, class PhysicalModel>
96 {
98
99 public:
101
102 // --------- Public methods ---------
103
112 std::unique_ptr<PhysicalModel> model)
113 : param_(param)
114 , model_(std::move(model))
115 , linearizations_(0)
116 , nonlinearIterations_(0)
117 , linearIterations_(0)
118 , wellIterations_(0)
119 , nonlinearIterationsLast_(0)
120 , linearIterationsLast_(0)
121 , wellIterationsLast_(0)
122 {
123 if (!model_) {
124 OPM_THROW(std::logic_error, "Must provide a non-null model argument for NonlinearSolver.");
125 }
126 }
127
128
130 {
132 report.global_time = timer.simulationTimeElapsed();
133 report.timestep_length = timer.currentStepLength();
134
135 // Do model-specific once-per-step calculations.
136 report += model_->prepareStep(timer);
137
138 int iteration = 0;
139
140 // Let the model do one nonlinear iteration.
141
142 // Set up for main solver loop.
143 bool converged = false;
144
145 // ---------- Main nonlinear solver loop ----------
146 do {
147 try {
148 // Do the nonlinear step. If we are in a converged state, the
149 // model will usually do an early return without an expensive
150 // solve, unless the newton_min_iter_ count has not been reached yet.
151 auto iterReport = model_->nonlinearIteration(iteration, timer, *this);
152 iterReport.global_time = timer.simulationTimeElapsed();
153 report += iterReport;
154 report.converged = iterReport.converged;
155
156 converged = report.converged;
157 iteration += 1;
158 }
159 catch (...) {
160 // if an iteration fails during a time step, all previous iterations
161 // count as a failure as well
162 failureReport_ = report;
163 failureReport_ += model_->failureReport();
164 throw;
165 }
166 }
167 while ( (!converged && (iteration <= this->model().param().newton_max_iter_)) ||
168 (iteration <= this->model().param().newton_min_iter_));
169
170 if (!converged) {
171 failureReport_ = report;
172
173 std::string msg = "Solver convergence failure - Failed to complete a time step within ";
174 msg += std::to_string(model().param().newton_max_iter_) + " iterations.";
175 OPM_THROW_NOLOG(TooManyIterations, msg);
176 }
177 auto relativeChange = model_->relativeChange();
178 if (timeStepControl && !timeStepControl->timeStepAccepted(relativeChange, timer.currentStepLength())) {
179 report.converged = false;
180 report.time_step_rejected = true;
181 failureReport_ = report;
182
183 std::string msg = "Relative change in solution for time step was " + std::to_string(relativeChange);
184 msg += ", which is larger than the tolerance accepted by the timestepping algorithm.";
185 OPM_THROW_NOLOG(TimeSteppingBreakdown, msg);
186 }
187
188 // Do model-specific post-step actions.
189 report += model_->afterStep(timer);
190 report.converged = true;
191 return report;
192 }
193
196 { return failureReport_; }
197
199 int linearizations() const
200 { return linearizations_; }
201
204 { return nonlinearIterations_; }
205
208 { return linearIterations_; }
209
211 int wellIterations() const
212 { return wellIterations_; }
213
216 { return nonlinearIterationsLast_; }
217
220 { return linearIterationsLast_; }
221
224 { return wellIterationsLast_; }
225
226 std::vector<std::vector<Scalar> >
227 computeFluidInPlace(const std::vector<int>& fipnum) const
228 { return model_->computeFluidInPlace(fipnum); }
229
231 const PhysicalModel& model() const
232 { return *model_; }
233
235 PhysicalModel& model()
236 { return *model_; }
237
239 void detectOscillations(const std::vector<std::vector<Scalar>>& residualHistory,
240 const int it, bool& oscillate, bool& stagnate) const
241 {
242 detail::detectOscillations(residualHistory, it, model_->numPhases(),
243 this->relaxRelTol(), 2, oscillate, stagnate);
244 }
245
246
249 template <class BVector>
250 void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const Scalar omega) const
251 {
252 detail::stabilizeNonlinearUpdate(dx, dxOld, omega, this->relaxType());
253 }
254
256 Scalar relaxMax() const
257 { return param_.relaxMax_; }
258
260 Scalar relaxIncrement() const
261 { return param_.relaxIncrement_; }
262
265 { return param_.relaxType_; }
266
268 Scalar relaxRelTol() const
269 { return param_.relaxRelTol_; }
270
273 { param_ = param; }
274
275 private:
276 // --------- Data members ---------
277 SimulatorReportSingle failureReport_;
278 SolverParameters param_;
279 std::unique_ptr<PhysicalModel> model_;
280 int linearizations_;
281 int nonlinearIterations_;
282 int linearIterations_;
283 int wellIterations_;
284 int nonlinearIterationsLast_;
285 int linearIterationsLast_;
286 int wellIterationsLast_;
287 };
288
289} // namespace Opm
290
291#endif // OPM_NON_LINEAR_SOLVER_HPP
Defines a type tags and some fundamental properties all models.
Definition: NonlinearSolver.hpp:96
int linearIterationsLastStep() const
Number of linear solver iterations used in the last call to step().
Definition: NonlinearSolver.hpp:219
void stabilizeNonlinearUpdate(BVector &dx, BVector &dxOld, const Scalar omega) const
Definition: NonlinearSolver.hpp:250
Scalar relaxRelTol() const
The relaxation relative tolerance.
Definition: NonlinearSolver.hpp:268
Scalar relaxMax() const
The greatest relaxation factor (i.e. smallest factor) allowed.
Definition: NonlinearSolver.hpp:256
int linearizations() const
Number of linearizations used in all calls to step().
Definition: NonlinearSolver.hpp:199
int wellIterationsLastStep() const
Number of well iterations used in all calls to step().
Definition: NonlinearSolver.hpp:223
int nonlinearIterations() const
Number of full nonlinear solver iterations used in all calls to step().
Definition: NonlinearSolver.hpp:203
void setParameters(const SolverParameters &param)
Set parameters to override those given at construction time.
Definition: NonlinearSolver.hpp:272
NonlinearRelaxType relaxType() const
The relaxation type (Dampen or SOR).
Definition: NonlinearSolver.hpp:264
int linearIterations() const
Number of linear solver iterations used in all calls to step().
Definition: NonlinearSolver.hpp:207
SimulatorReportSingle step(const SimulatorTimerInterface &timer, const TimeStepControlInterface *timeStepControl)
Definition: NonlinearSolver.hpp:129
int nonlinearIterationsLastStep() const
Number of nonlinear solver iterations used in the last call to step().
Definition: NonlinearSolver.hpp:215
const SimulatorReportSingle & failureReport() const
return the statistics if the step() method failed
Definition: NonlinearSolver.hpp:195
int wellIterations() const
Number of well iterations used in all calls to step().
Definition: NonlinearSolver.hpp:211
const PhysicalModel & model() const
Reference to physical model.
Definition: NonlinearSolver.hpp:231
void detectOscillations(const std::vector< std::vector< Scalar > > &residualHistory, const int it, bool &oscillate, bool &stagnate) const
Detect oscillation or stagnation in a given residual history.
Definition: NonlinearSolver.hpp:239
NonlinearSolverParameters< Scalar > SolverParameters
Definition: NonlinearSolver.hpp:100
NonlinearSolver(const SolverParameters &param, std::unique_ptr< PhysicalModel > model)
Definition: NonlinearSolver.hpp:111
Scalar relaxIncrement() const
The step-change size for the relaxation factor.
Definition: NonlinearSolver.hpp:260
std::vector< std::vector< Scalar > > computeFluidInPlace(const std::vector< int > &fipnum) const
Definition: NonlinearSolver.hpp:227
PhysicalModel & model()
Mutable reference to physical model.
Definition: NonlinearSolver.hpp:235
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:34
virtual double currentStepLength() const =0
virtual double simulationTimeElapsed() const =0
Definition: TimeStepControlInterface.hpp:51
virtual bool timeStepAccepted(const double error, const double timeStepJustCompleted) const =0
Definition: blackoilnewtonmethodparams.hpp:31
void stabilizeNonlinearUpdate(BVector &dx, BVector &dxOld, const Scalar omega, NonlinearRelaxType relaxType)
void detectOscillations(const std::vector< std::vector< Scalar > > &residualHistory, const int it, const int numPhases, const Scalar relaxRelTol, const int minimumOscillatingPhases, bool &oscillate, bool &stagnate)
Detect oscillation or stagnation in a given residual history.
Definition: blackoilboundaryratevector.hh:39
NonlinearRelaxType
Definition: NonlinearSolver.hpp:54
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:233
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
The Opm property system, traits with inheritance.
Definition: NonlinearSolver.hpp:79
Scalar relaxRelTol_
Definition: NonlinearSolver.hpp:83
Scalar relaxMax_
Definition: NonlinearSolver.hpp:81
NonlinearRelaxType relaxType_
Definition: NonlinearSolver.hpp:80
Scalar relaxIncrement_
Definition: NonlinearSolver.hpp:82
Definition: NonlinearSolver.hpp:45
static constexpr Scalar value
Definition: NonlinearSolver.hpp:45
Definition: NonlinearSolver.hpp:47
static constexpr auto value
Definition: NonlinearSolver.hpp:47
A struct for returning timing data from a simulator to its caller.
Definition: SimulatorReport.hpp:34
bool converged
Definition: SimulatorReport.hpp:55
double global_time
Definition: SimulatorReport.hpp:60
double timestep_length
Definition: SimulatorReport.hpp:61
bool time_step_rejected
Definition: SimulatorReport.hpp:56