NonlinearSolver_impl.hpp
1 /*
2  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
3  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
4  Copyright 2015 NTNU
5  Copyright 2015 IRIS AS
6 
7  This file is part of the Open Porous Media project (OPM).
8 
9  OPM is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation, either version 3 of the License, or
12  (at your option) any later version.
13 
14  OPM is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  along with OPM. If not, see <http://www.gnu.org/licenses/>.
21 */
22 
23 #ifndef OPM_NONLINEARSOLVER_IMPL_HEADER_INCLUDED
24 #define OPM_NONLINEARSOLVER_IMPL_HEADER_INCLUDED
25 
26 #include <opm/autodiff/NonlinearSolver.hpp>
27 #include <opm/common/Exceptions.hpp>
28 #include <opm/common/ErrorMacros.hpp>
29 
30 namespace Opm
31 {
32  template <class PhysicalModel>
33  NonlinearSolver<PhysicalModel>::NonlinearSolver(const SolverParameters& param,
34  std::unique_ptr<PhysicalModel> model_arg)
35  : param_(param),
36  model_(std::move(model_arg)),
37  linearizations_(0),
38  nonlinearIterations_(0),
39  linearIterations_(0),
40  wellIterations_(0),
41  nonlinearIterationsLast_(0),
42  linearIterationsLast_(0),
43  wellIterationsLast_(0)
44  {
45  if (!model_) {
46  OPM_THROW(std::logic_error, "Must provide a non-null model argument for NonlinearSolver.");
47  }
48  }
49 
50  template <class PhysicalModel>
52  {
53  return linearizations_;
54  }
55 
56  template <class PhysicalModel>
58  {
59  return nonlinearIterations_;
60  }
61 
62  template <class PhysicalModel>
64  {
65  return linearIterations_;
66  }
67 
68  template <class PhysicalModel>
70  {
71  return wellIterations_;
72  }
73 
74  template <class PhysicalModel>
75  const PhysicalModel& NonlinearSolver<PhysicalModel>::model() const
76  {
77  return *model_;
78  }
79 
80  template <class PhysicalModel>
82  {
83  return *model_;
84  }
85 
86  template <class PhysicalModel>
88  {
89  return nonlinearIterationsLast_;
90  }
91 
92  template <class PhysicalModel>
94  {
95  return linearIterationsLast_;
96  }
97 
98  template <class PhysicalModel>
100  {
101  return wellIterationsLast_;
102  }
103 
104  template <class PhysicalModel>
105  SimulatorReport
108  ReservoirState& reservoir_state,
109  WellState& well_state)
110  {
111  return step(timer, reservoir_state, well_state, reservoir_state, well_state);
112  }
113 
114 
115 
116  template <class PhysicalModel>
117  SimulatorReport
120  const ReservoirState& initial_reservoir_state,
121  const WellState& initial_well_state,
122  ReservoirState& reservoir_state,
123  WellState& well_state)
124  {
125  SimulatorReport iterReport;
126  SimulatorReport report;
127  failureReport_ = SimulatorReport();
128 
129  // Do model-specific once-per-step calculations.
130  model_->prepareStep(timer, initial_reservoir_state, initial_well_state);
131 
132  int iteration = 0;
133 
134  // Let the model do one nonlinear iteration.
135 
136  // Set up for main solver loop.
137  bool converged = false;
138 
139  // ---------- Main nonlinear solver loop ----------
140  do {
141  try {
142  // Do the nonlinear step. If we are in a converged state, the
143  // model will usually do an early return without an expensive
144  // solve, unless the minIter() count has not been reached yet.
145  iterReport = model_->nonlinearIteration(iteration, timer, *this, reservoir_state, well_state);
146 
147  report += iterReport;
148  report.converged = iterReport.converged;
149 
150  converged = report.converged;
151  iteration += 1;
152  }
153  catch (...) {
154  // if an iteration fails during a time step, all previous iterations
155  // count as a failure as well
156  failureReport_ += report;
157  failureReport_ += model_->failureReport();
158  throw;
159  }
160  } while ( (!converged && (iteration <= maxIter())) || (iteration <= minIter()));
161 
162  if (!converged) {
163  failureReport_ += report;
164 
165  std::string msg = "Solver convergence failure - Failed to complete a time step within " + std::to_string(maxIter()) + " iterations.";
166  OPM_THROW_NOLOG(Opm::TooManyIterations, msg);
167  }
168 
169  // Do model-specific post-step actions.
170  model_->afterStep(timer, reservoir_state, well_state);
171  report.converged = true;
172 
173  return report;
174  }
175 
176 
177 
178  template <class PhysicalModel>
180  reset()
181  {
182  // default values for the solver parameters
183  relax_type_ = DAMPEN;
184  relax_max_ = 0.5;
185  relax_increment_ = 0.1;
186  relax_rel_tol_ = 0.2;
187  max_iter_ = 10;
188  min_iter_ = 1;
189  }
190 
191  template <class PhysicalModel>
192  NonlinearSolver<PhysicalModel>::SolverParameters::
193  SolverParameters()
194  {
195  // set default values
196  reset();
197  }
198 
199  template <class PhysicalModel>
200  NonlinearSolver<PhysicalModel>::SolverParameters::
201  SolverParameters( const ParameterGroup& param )
202  {
203  // set default values
204  reset();
205 
206  // overload with given parameters
207  relax_max_ = param.getDefault("relax_max", relax_max_);
208  max_iter_ = param.getDefault("max_iter", max_iter_);
209  min_iter_ = param.getDefault("min_iter", min_iter_);
210 
211  std::string relaxation_type = param.getDefault("relax_type", std::string("dampen"));
212  if (relaxation_type == "dampen") {
213  relax_type_ = DAMPEN;
214  } else if (relaxation_type == "sor") {
215  relax_type_ = SOR;
216  } else {
217  OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxation_type);
218  }
219  }
220 
221  template <class PhysicalModel>
222  void
223  NonlinearSolver<PhysicalModel>::detectOscillations(const std::vector<std::vector<double>>& residual_history,
224  const int it,
225  bool& oscillate, bool& stagnate) const
226  {
227  // The detection of oscillation in two primary variable results in the report of the detection
228  // of oscillation for the solver.
229  // Only the saturations are used for oscillation detection for the black oil model.
230  // Stagnate is not used for any treatment here.
231 
232  if ( it < 2 ) {
233  oscillate = false;
234  stagnate = false;
235  return;
236  }
237 
238  stagnate = true;
239  int oscillatePhase = 0;
240  const std::vector<double>& F0 = residual_history[it];
241  const std::vector<double>& F1 = residual_history[it - 1];
242  const std::vector<double>& F2 = residual_history[it - 2];
243  for (int p= 0; p < model_->numPhases(); ++p){
244  const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
245  const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);
246 
247  oscillatePhase += (d1 < relaxRelTol()) && (relaxRelTol() < d2);
248 
249  // Process is 'stagnate' unless at least one phase
250  // exhibits significant residual change.
251  stagnate = (stagnate && !(std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3));
252  }
253 
254  oscillate = (oscillatePhase > 1);
255  }
256 
257 
258  template <class PhysicalModel>
259  template <class BVector>
260  void
261  NonlinearSolver<PhysicalModel>::stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
262  {
263  // The dxOld is updated with dx.
264  // If omega is equal to 1., no relaxtion will be appiled.
265 
266  BVector tempDxOld = dxOld;
267  dxOld = dx;
268 
269  switch (relaxType()) {
270  case DAMPEN: {
271  if (omega == 1.) {
272  return;
273  }
274  auto i = dx.size();
275  for (i = 0; i < dx.size(); ++i) {
276  dx[i] *= omega;
277  }
278  return;
279  }
280  case SOR: {
281  if (omega == 1.) {
282  return;
283  }
284  auto i = dx.size();
285  for (i = 0; i < dx.size(); ++i) {
286  dx[i] *= omega;
287  tempDxOld[i] *= (1.-omega);
288  dx[i] += tempDxOld[i];
289  }
290  return;
291  }
292  default:
293  OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type.");
294  }
295 
296  return;
297  }
298 } // namespace Opm
299 
300 
301 #endif // OPM_FULLYIMPLICITSOLVER_IMPL_HEADER_INCLUDED
A nonlinear solver class suitable for general fully-implicit models, as well as pressure, transport and sequential models.
Definition: NonlinearSolver.hpp:37
NonlinearSolver(const SolverParameters &param, std::unique_ptr< PhysicalModel > model)
Construct solver for a given model.
Definition: NonlinearSolver_impl.hpp:33
Definition: AutoDiff.hpp:297
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:35