OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Make nonlinearIteration() the only interface to the model.

Browse Source 
This means that details such as calling assemble(), solveJacobianSystem(),
updateState() etc. are now left to the model class. This will make it easier
to create new model classes with different behaviour (such as sequential models).
This commit is contained in:
Atgeirr Flø Rasmussen
2015-10-02 13:51:40 +02:00
parent 7c21a630e5
commit 29a1a891d2
4 changed files with 113 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
opm/autodiff/BlackoilModelBase.hpp
View File

@@ -76,6 +76,17 @@ namespace Opm {
/// Class used for reporting the outcome of a nonlinearIteration() call.
struct IterationReport
{
bool failed;
bool converged;
int linear_iterations;
};
/// Traits to encapsulate the types used by classes using or
/// extending this model. Forward declared here, must be
/// specialised for each concrete model class.
@@ -155,6 +166,22 @@ namespace Opm {
ReservoirState& reservoir_state,
WellState& well_state);
/// Called once per nonlinear iteration.
/// This model will perform a Newton-Raphson update, changing reservoir_state
/// and well_state. It will also use the nonlinear_solver to do relaxation of
/// updates if necessary.
/// \param[in] iteration should be 0 for the first call of a new timestep
/// \param[in] dt time step size
/// \param[in] nonlinear_solver nonlinear solver used (for oscillation/relaxation control)
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
template <class NonlinearSolverType>
IterationReport nonlinearIteration(const int iteration,
const double dt,
NonlinearSolverType& nonlinear_solver,
ReservoirState& reservoir_state,
WellState& well_state);
/// Called once after each time step.
/// In this class, this function does nothing.
/// \param[in] dt time step size
@@ -290,6 +317,9 @@ namespace Opm {
std::vector<int> primalVariable_;
V pvdt_;
std::vector<std::string> material_name_;
std::vector<std::vector<double>> residual_norms_history_;
double current_relaxation_;
V dx_old_;
// --------- Protected methods ---------

51
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@@ -180,6 +180,7 @@ namespace detail {
{ 1.1169, 1.0031, 0.0031 }} ) // the default magic numbers
, terminal_output_ (terminal_output)
, material_name_{ "Water", "Oil", "Gas" }
, current_relaxation_(1.0)
{
assert(numMaterials() == 3); // Due to the material_name_ init above.
#if HAVE_MPI
@@ -227,6 +228,56 @@ namespace detail {
template <class Grid, class Implementation>
template <class NonlinearSolverType>
IterationReport
BlackoilModelBase<Grid, Implementation>::
nonlinearIteration(const int iteration,
const double dt,
NonlinearSolverType& nonlinear_solver,
ReservoirState& reservoir_state,
WellState& well_state)
{
if (iteration == 0) {
// For each iteration we store in a vector the norms of the residual of
// the mass balance for each active phase, the well flux and the well equations.
residual_norms_history_.clear();
current_relaxation_ = 1.0;
dx_old_ = V::Zero(sizeNonLinear());
}
assemble(reservoir_state, well_state, iteration == 0);
residual_norms_history_.push_back(computeResidualNorms());
const bool converged = getConvergence(dt, iteration);
if (!converged) {
V dx = solveJacobianSystem();
// Stabilize the nonlinear update.
bool isOscillate = false;
bool isStagnate = false;
nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate);
if (isOscillate) {
current_relaxation_ -= nonlinear_solver.relaxIncrement();
current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
if (terminalOutputEnabled()) {
std::cout << " Oscillating behavior detected: Relaxation set to "
<< current_relaxation_ << std::endl;
}
}
nonlinear_solver.stabilizeNonlinearUpdate(dx, dx_old_, current_relaxation_);
// Apply the update, applying model-dependent
// limitations and chopping of the update.
updateState(dx, reservoir_state, well_state);
}
const bool failed = false; // Not needed in this model.
const int linear_iters = converged ? 0 : linearIterationsLastSolve();
return IterationReport{ failed, converged, linear_iters };
}
template <class Grid, class Implementation>
void
BlackoilModelBase<Grid, Implementation>::

21
opm/autodiff/NonlinearSolver.hpp
View File

@@ -97,9 +97,22 @@ namespace Opm {
/// Number of linear solver iterations used in the last call to step().
unsigned int linearIterationsLastStep() const;
/// return reference to physical model
/// Reference to physical model.
const PhysicalModel& model() const;
/// Detect oscillation or stagnation in a given residual history.
void detectOscillations(const std::vector<std::vector<double>>& residual_history,
const int it, bool& oscillate, bool& stagnate) const;
/// Apply a stabilization to dx, depending on dxOld and relaxation parameters.
void stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega) const;
/// The greatest relaxation factor (i.e. smallest factor) allowed.
double relaxMax() const { return param_.relax_max_; }
/// The step-change size for the relaxation factor.
double relaxIncrement() const { return param_.relax_increment_; }
private:
// --------- Data members ---------
SolverParameters param_;
@@ -111,15 +124,9 @@ namespace Opm {
// --------- Private methods ---------
enum RelaxType relaxType() const { return param_.relax_type_; }
double relaxMax() const { return param_.relax_max_; }
double relaxIncrement() const { return param_.relax_increment_; }
double relaxRelTol() const { return param_.relax_rel_tol_; }
double maxIter() const { return param_.max_iter_; }
double minIter() const { return param_.min_iter_; }
void detectOscillations(const std::vector<std::vector<double>>& residual_history,
const int it, const double relaxRelTol,
bool& oscillate, bool& stagnate) const;
void stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const;
};
} // namespace Opm

67
opm/autodiff/NonlinearSolver_impl.hpp
View File

@@ -81,57 +81,27 @@ namespace Opm
// Do model-specific once-per-step calculations.
model_->prepareStep(dt, reservoir_state, well_state);
// For each iteration we store in a vector the norms of the residual of
// the mass balance for each active phase, the well flux and the well equations.
std::vector<std::vector<double>> residual_norms_history;
int iteration = 0;
// Assemble residual and Jacobian, store residual norms.
model_->assemble(reservoir_state, well_state, true);
residual_norms_history.push_back(model_->computeResidualNorms());
// Let the model do one nonlinear iteration.
// Set up for main solver loop.
double omega = 1.0;
int iteration = 0;
bool converged = model_->getConvergence(dt, iteration);
const int sizeNonLinear = model_->sizeNonLinear();
V dxOld = V::Zero(sizeNonLinear);
bool isOscillate = false;
bool isStagnate = false;
const enum RelaxType relaxtype = relaxType();
int linIters = 0;
bool converged = false;
// ---------- Main nonlinear solver loop ----------
while ( (!converged && (iteration < maxIter())) || (minIter() > iteration)) {
// Compute the update to the primary variables.
V dx = model_->solveJacobianSystem();
// Store number of linear iterations used.
linIters += model_->linearIterationsLastSolve();
// Stabilize the nonlinear update.
detectOscillations(residual_norms_history, iteration, relaxRelTol(), isOscillate, isStagnate);
if (isOscillate) {
omega -= relaxIncrement();
omega = std::max(omega, relaxMax());
if (model_->terminalOutputEnabled()) {
std::cout << " Oscillating behavior detected: Relaxation set to " << omega << std::endl;
}
do {
IterationReport report = model_->nonlinearIteration(iteration, dt, *this, reservoir_state, well_state);
if (report.failed) {
OPM_THROW(Opm::NumericalProblem, "Failed to complete a nonlinear iteration.");
}
stabilizeNonlinearUpdate(dx, dxOld, omega, relaxtype);
// Apply the update, the model may apply model-dependent
// limitations and chopping of the update.
model_->updateState(dx, reservoir_state, well_state);
// Assemble residual and Jacobian, store residual norms.
model_->assemble(reservoir_state, well_state, false);
residual_norms_history.push_back(model_->computeResidualNorms());
// increase iteration counter
if (report.converged) {
assert(report.linear_iterations == 0);
converged = true;
}
linIters += report.linear_iterations;
++iteration;
converged = model_->getConvergence(dt, iteration);
}
} while ( (!converged && (iteration <= maxIter())) || (minIter() > iteration));
if (!converged) {
if (model_->terminalOutputEnabled()) {
@@ -141,7 +111,7 @@ namespace Opm
}
linearIterations_ += linIters;
nonlinearIterations_ += iteration;
nonlinearIterations_ += iteration - 1; // Since the last one will always be trivial.
linearIterationsLast_ = linIters;
nonlinearIterationsLast_ = iteration;
@@ -199,7 +169,7 @@ namespace Opm
template <class PhysicalModel>
void
NonlinearSolver<PhysicalModel>::detectOscillations(const std::vector<std::vector<double>>& residual_history,
const int it, const double relaxRelTol_arg,
const int it,
bool& oscillate, bool& stagnate) const
{
// The detection of oscillation in two primary variable results in the report of the detection
@@ -222,7 +192,7 @@ namespace Opm
const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);
oscillatePhase += (d1 < relaxRelTol_arg) && (relaxRelTol_arg < d2);
oscillatePhase += (d1 < relaxRelTol()) && (relaxRelTol() < d2);
// Process is 'stagnate' unless at least one phase
// exhibits significant residual change.
@@ -235,8 +205,7 @@ namespace Opm
template <class PhysicalModel>
void
NonlinearSolver<PhysicalModel>::stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega,
const RelaxType relax_type) const
NonlinearSolver<PhysicalModel>::stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega) const
{
// The dxOld is updated with dx.
// If omega is equal to 1., no relaxtion will be appiled.
@@ -244,7 +213,7 @@ namespace Opm
const V tempDxOld = dxOld;
dxOld = dx;
switch (relax_type) {
switch (relaxType()) {
case DAMPEN:
if (omega == 1.) {
return;