23 #ifndef OPM_NEWTONSOLVER_IMPL_HEADER_INCLUDED
24 #define OPM_NEWTONSOLVER_IMPL_HEADER_INCLUDED
30 template <
class PhysicalModel>
32 std::unique_ptr<PhysicalModel> model)
34 model_(
std::move(model)),
40 template <
class PhysicalModel>
43 return newtonIterations_;
46 template <
class PhysicalModel>
49 return linearIterations_;
53 template <
class PhysicalModel>
61 model_->prepareStep(dt, reservoir_state, well_state);
65 std::vector<std::vector<double>> residual_norms_history;
68 model_->assemble(reservoir_state, well_state,
true);
69 residual_norms_history.push_back(model_->computeResidualNorms());
74 bool converged = model_->getConvergence(dt, iteration);
75 const int sizeNonLinear = model_->sizeNonLinear();
76 V dxOld = V::Zero(sizeNonLinear);
77 bool isOscillate =
false;
78 bool isStagnate =
false;
79 const enum RelaxType relaxtype = relaxType();
83 while ( (!converged && (iteration < maxIter())) || (minIter() > iteration)) {
85 V dx = model_->solveJacobianSystem();
88 linIters += model_->linearIterationsLastSolve();
91 detectNewtonOscillations(residual_norms_history, iteration, relaxRelTol(), isOscillate, isStagnate);
93 omega -= relaxIncrement();
94 omega = std::max(omega, relaxMax());
95 if (model_->terminalOutputEnabled()) {
96 std::cout <<
" Oscillating behavior detected: Relaxation set to " << omega << std::endl;
99 stabilizeNewton(dx, dxOld, omega, relaxtype);
103 model_->updateState(dx, reservoir_state, well_state);
106 model_->assemble(reservoir_state, well_state,
false);
107 residual_norms_history.push_back(model_->computeResidualNorms());
112 converged = model_->getConvergence(dt, iteration);
116 if (model_->terminalOutputEnabled()) {
117 std::cerr <<
"WARNING: Failed to compute converged solution in " << iteration <<
" iterations." << std::endl;
122 linearIterations_ += linIters;
123 newtonIterations_ += iteration;
124 linearIterationsLast_ = linIters;
125 newtonIterationsLast_ = iteration;
128 model_->afterStep(dt, reservoir_state, well_state);
135 template <
class PhysicalModel>
148 template <
class PhysicalModel>
156 template <
class PhysicalModel>
164 relax_max_ = param.getDefault(
"relax_max", relax_max_);
165 max_iter_ = param.getDefault(
"max_iter", max_iter_);
166 min_iter_ = param.getDefault(
"min_iter", min_iter_);
168 std::string relaxation_type = param.getDefault(
"relax_type", std::string(
"dampen"));
169 if (relaxation_type ==
"dampen") {
171 }
else if (relaxation_type ==
"sor") {
174 OPM_THROW(std::runtime_error,
"Unknown Relaxtion Type " << relaxation_type);
178 template <
class PhysicalModel>
181 const int it,
const double relaxRelTol_arg,
182 bool& oscillate,
bool& stagnate)
const
196 int oscillatePhase = 0;
197 const std::vector<double>& F0 = residual_history[it];
198 const std::vector<double>& F1 = residual_history[it - 1];
199 const std::vector<double>& F2 = residual_history[it - 2];
200 for (
int p= 0; p < model_->numPhases(); ++p){
201 const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
202 const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);
204 oscillatePhase += (d1 < relaxRelTol_arg) && (relaxRelTol_arg < d2);
208 stagnate = (stagnate && !(std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3));
211 oscillate = (oscillatePhase > 1);
215 template <
class PhysicalModel>
217 NewtonSolver<PhysicalModel>::stabilizeNewton(
V& dx,
V& dxOld,
const double omega,
223 const V tempDxOld = dxOld;
226 switch (relax_type) {
237 dx = dx*omega + (1.-omega)*tempDxOld;
240 OPM_THROW(std::runtime_error,
"Can only handle DAMPEN and SOR relaxation type.");
250 #endif // OPM_FULLYIMPLICITSOLVER_IMPL_HEADER_INCLUDED
int min_iter_
Definition: NewtonSolver.hpp:52
int max_iter_
Definition: NewtonSolver.hpp:51
PhysicalModel::ReservoirState ReservoirState
Definition: NewtonSolver.hpp:61
int step(const double dt, ReservoirState &reservoir_state, WellState &well_state)
Definition: NewtonSolver_impl.hpp:56
Definition: AdditionalObjectDeleter.hpp:22
void reset()
Definition: NewtonSolver_impl.hpp:137
PhysicalModel::WellState WellState
Definition: NewtonSolver.hpp:62
double relax_rel_tol_
Definition: NewtonSolver.hpp:50
unsigned int newtonIterations() const
Number of Newton iterations used in all calls to step().
Definition: NewtonSolver_impl.hpp:41
double relax_max_
Definition: NewtonSolver.hpp:48
NewtonSolver(const SolverParameters ¶m, std::unique_ptr< PhysicalModel > model)
Definition: NewtonSolver_impl.hpp:31
RelaxType
Definition: NewtonSolver.hpp:42
Definition: NewtonSolver.hpp:42
double relax_increment_
Definition: NewtonSolver.hpp:49
A Newton solver class suitable for general fully-implicit models.
Definition: NewtonSolver.hpp:33
SolverParameters()
Definition: NewtonSolver_impl.hpp:150
enum RelaxType relax_type_
Definition: NewtonSolver.hpp:47
unsigned int linearIterations() const
Number of linear solver iterations used in all calls to step().
Definition: NewtonSolver_impl.hpp:47
Definition: NewtonSolver.hpp:45
Definition: NewtonSolver.hpp:42
ADB::V V
Definition: NewtonSolver.hpp:38