mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Merge pull request #536 from atgeirr/refactor-nonlinear
Refactor nonlinear solver
This commit is contained in:
commit
022bd6b8fe
@ -126,8 +126,8 @@ list (APPEND PUBLIC_HEADER_FILES
|
|||||||
opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
|
opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
|
||||||
opm/autodiff/NewtonIterationBlackoilSimple.hpp
|
opm/autodiff/NewtonIterationBlackoilSimple.hpp
|
||||||
opm/autodiff/NewtonIterationUtilities.hpp
|
opm/autodiff/NewtonIterationUtilities.hpp
|
||||||
opm/autodiff/NewtonSolver.hpp
|
opm/autodiff/NonlinearSolver.hpp
|
||||||
opm/autodiff/NewtonSolver_impl.hpp
|
opm/autodiff/NonlinearSolver_impl.hpp
|
||||||
opm/autodiff/LinearisedBlackoilResidual.hpp
|
opm/autodiff/LinearisedBlackoilResidual.hpp
|
||||||
opm/autodiff/ParallelDebugOutput.hpp
|
opm/autodiff/ParallelDebugOutput.hpp
|
||||||
opm/autodiff/RateConverter.hpp
|
opm/autodiff/RateConverter.hpp
|
||||||
|
@ -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
|
/// Traits to encapsulate the types used by classes using or
|
||||||
/// extending this model. Forward declared here, must be
|
/// extending this model. Forward declared here, must be
|
||||||
/// specialised for each concrete model class.
|
/// specialised for each concrete model class.
|
||||||
@ -155,6 +166,22 @@ namespace Opm {
|
|||||||
ReservoirState& reservoir_state,
|
ReservoirState& reservoir_state,
|
||||||
WellState& well_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.
|
/// Called once after each time step.
|
||||||
/// In this class, this function does nothing.
|
/// In this class, this function does nothing.
|
||||||
/// \param[in] dt time step size
|
/// \param[in] dt time step size
|
||||||
@ -290,6 +317,9 @@ namespace Opm {
|
|||||||
std::vector<int> primalVariable_;
|
std::vector<int> primalVariable_;
|
||||||
V pvdt_;
|
V pvdt_;
|
||||||
std::vector<std::string> material_name_;
|
std::vector<std::string> material_name_;
|
||||||
|
std::vector<std::vector<double>> residual_norms_history_;
|
||||||
|
double current_relaxation_;
|
||||||
|
V dx_old_;
|
||||||
|
|
||||||
// --------- Protected methods ---------
|
// --------- Protected methods ---------
|
||||||
|
|
||||||
|
@ -180,6 +180,7 @@ namespace detail {
|
|||||||
{ 1.1169, 1.0031, 0.0031 }} ) // the default magic numbers
|
{ 1.1169, 1.0031, 0.0031 }} ) // the default magic numbers
|
||||||
, terminal_output_ (terminal_output)
|
, terminal_output_ (terminal_output)
|
||||||
, material_name_{ "Water", "Oil", "Gas" }
|
, material_name_{ "Water", "Oil", "Gas" }
|
||||||
|
, current_relaxation_(1.0)
|
||||||
{
|
{
|
||||||
assert(numMaterials() == 3); // Due to the material_name_ init above.
|
assert(numMaterials() == 3); // Due to the material_name_ init above.
|
||||||
#if HAVE_MPI
|
#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>
|
template <class Grid, class Implementation>
|
||||||
void
|
void
|
||||||
BlackoilModelBase<Grid, Implementation>::
|
BlackoilModelBase<Grid, Implementation>::
|
||||||
|
@ -18,8 +18,8 @@
|
|||||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
*/
|
*/
|
||||||
|
|
||||||
#ifndef OPM_NEWTONSOLVER_HEADER_INCLUDED
|
#ifndef OPM_NONLINEARSOLVER_HEADER_INCLUDED
|
||||||
#define OPM_NEWTONSOLVER_HEADER_INCLUDED
|
#define OPM_NONLINEARSOLVER_HEADER_INCLUDED
|
||||||
|
|
||||||
#include <opm/autodiff/AutoDiffBlock.hpp>
|
#include <opm/autodiff/AutoDiffBlock.hpp>
|
||||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||||
@ -28,9 +28,10 @@
|
|||||||
namespace Opm {
|
namespace Opm {
|
||||||
|
|
||||||
|
|
||||||
/// A Newton solver class suitable for general fully-implicit models.
|
/// A nonlinear solver class suitable for general fully-implicit models,
|
||||||
|
/// as well as pressure, transport and sequential models.
|
||||||
template <class PhysicalModel>
|
template <class PhysicalModel>
|
||||||
class NewtonSolver
|
class NonlinearSolver
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
// --------- Types and enums ---------
|
// --------- Types and enums ---------
|
||||||
@ -38,18 +39,18 @@ namespace Opm {
|
|||||||
typedef ADB::V V;
|
typedef ADB::V V;
|
||||||
typedef ADB::M M;
|
typedef ADB::M M;
|
||||||
|
|
||||||
// The Newton relaxation scheme type
|
// Available relaxation scheme types.
|
||||||
enum RelaxType { DAMPEN, SOR };
|
enum RelaxType { DAMPEN, SOR };
|
||||||
|
|
||||||
// Solver parameters controlling nonlinear Newton process.
|
// Solver parameters controlling nonlinear process.
|
||||||
struct SolverParameters
|
struct SolverParameters
|
||||||
{
|
{
|
||||||
enum RelaxType relax_type_;
|
enum RelaxType relax_type_;
|
||||||
double relax_max_;
|
double relax_max_;
|
||||||
double relax_increment_;
|
double relax_increment_;
|
||||||
double relax_rel_tol_;
|
double relax_rel_tol_;
|
||||||
int max_iter_; // max newton iterations
|
int max_iter_; // max nonlinear iterations
|
||||||
int min_iter_; // min newton iterations
|
int min_iter_; // min nonlinear iterations
|
||||||
|
|
||||||
explicit SolverParameters( const parameter::ParameterGroup& param );
|
explicit SolverParameters( const parameter::ParameterGroup& param );
|
||||||
SolverParameters();
|
SolverParameters();
|
||||||
@ -66,11 +67,11 @@ namespace Opm {
|
|||||||
/// Construct solver for a given model.
|
/// Construct solver for a given model.
|
||||||
///
|
///
|
||||||
/// The model is a std::unique_ptr because the object to which model points to is
|
/// The model is a std::unique_ptr because the object to which model points to is
|
||||||
/// not allowed to be deleted as long as the NewtonSolver object exists.
|
/// not allowed to be deleted as long as the NonlinearSolver object exists.
|
||||||
///
|
///
|
||||||
/// \param[in] param parameters controlling nonlinear Newton process
|
/// \param[in] param parameters controlling nonlinear process
|
||||||
/// \param[in, out] model physical simulation model.
|
/// \param[in, out] model physical simulation model.
|
||||||
explicit NewtonSolver(const SolverParameters& param,
|
explicit NonlinearSolver(const SolverParameters& param,
|
||||||
std::unique_ptr<PhysicalModel> model);
|
std::unique_ptr<PhysicalModel> model);
|
||||||
|
|
||||||
/// Take a single forward step, after which the states will be modified
|
/// Take a single forward step, after which the states will be modified
|
||||||
@ -84,44 +85,51 @@ namespace Opm {
|
|||||||
ReservoirState& reservoir_state,
|
ReservoirState& reservoir_state,
|
||||||
WellState& well_state);
|
WellState& well_state);
|
||||||
|
|
||||||
/// Number of Newton iterations used in all calls to step().
|
/// Number of nonlinear solver iterations used in all calls to step().
|
||||||
unsigned int newtonIterations() const;
|
unsigned int nonlinearIterations() const;
|
||||||
|
|
||||||
/// Number of linear solver iterations used in all calls to step().
|
/// Number of linear solver iterations used in all calls to step().
|
||||||
unsigned int linearIterations() const;
|
unsigned int linearIterations() const;
|
||||||
|
|
||||||
/// Number of linear solver iterations used in the last call to step().
|
/// Number of nonlinear solver iterations used in the last call to step().
|
||||||
unsigned int newtonIterationsLastStep() const;
|
unsigned int nonlinearIterationsLastStep() const;
|
||||||
|
|
||||||
/// Number of linear solver iterations used in the last call to step().
|
/// Number of linear solver iterations used in the last call to step().
|
||||||
unsigned int linearIterationsLastStep() const;
|
unsigned int linearIterationsLastStep() const;
|
||||||
|
|
||||||
/// return reference to physical model
|
/// Reference to physical model.
|
||||||
const PhysicalModel& model() const;
|
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:
|
private:
|
||||||
// --------- Data members ---------
|
// --------- Data members ---------
|
||||||
SolverParameters param_;
|
SolverParameters param_;
|
||||||
std::unique_ptr<PhysicalModel> model_;
|
std::unique_ptr<PhysicalModel> model_;
|
||||||
unsigned int newtonIterations_;
|
unsigned int nonlinearIterations_;
|
||||||
unsigned int linearIterations_;
|
unsigned int linearIterations_;
|
||||||
unsigned int newtonIterationsLast_;
|
unsigned int nonlinearIterationsLast_;
|
||||||
unsigned int linearIterationsLast_;
|
unsigned int linearIterationsLast_;
|
||||||
|
|
||||||
// --------- Private methods ---------
|
// --------- Private methods ---------
|
||||||
enum RelaxType relaxType() const { return param_.relax_type_; }
|
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 relaxRelTol() const { return param_.relax_rel_tol_; }
|
||||||
double maxIter() const { return param_.max_iter_; }
|
double maxIter() const { return param_.max_iter_; }
|
||||||
double minIter() const { return param_.min_iter_; }
|
double minIter() const { return param_.min_iter_; }
|
||||||
void detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
|
|
||||||
const int it, const double relaxRelTol,
|
|
||||||
bool& oscillate, bool& stagnate) const;
|
|
||||||
void stabilizeNewton(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const;
|
|
||||||
};
|
};
|
||||||
} // namespace Opm
|
} // namespace Opm
|
||||||
|
|
||||||
#include "NewtonSolver_impl.hpp"
|
#include "NonlinearSolver_impl.hpp"
|
||||||
|
|
||||||
#endif // OPM_NEWTONSOLVER_HEADER_INCLUDED
|
#endif // OPM_NONLINEARSOLVER_HEADER_INCLUDED
|
@ -20,45 +20,62 @@
|
|||||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
*/
|
*/
|
||||||
|
|
||||||
#ifndef OPM_NEWTONSOLVER_IMPL_HEADER_INCLUDED
|
#ifndef OPM_NONLINEARSOLVER_IMPL_HEADER_INCLUDED
|
||||||
#define OPM_NEWTONSOLVER_IMPL_HEADER_INCLUDED
|
#define OPM_NONLINEARSOLVER_IMPL_HEADER_INCLUDED
|
||||||
|
|
||||||
#include <opm/autodiff/NewtonSolver.hpp>
|
#include <opm/autodiff/NonlinearSolver.hpp>
|
||||||
|
|
||||||
namespace Opm
|
namespace Opm
|
||||||
{
|
{
|
||||||
template <class PhysicalModel>
|
template <class PhysicalModel>
|
||||||
NewtonSolver<PhysicalModel>::NewtonSolver(const SolverParameters& param,
|
NonlinearSolver<PhysicalModel>::NonlinearSolver(const SolverParameters& param,
|
||||||
std::unique_ptr<PhysicalModel> model)
|
std::unique_ptr<PhysicalModel> model)
|
||||||
: param_(param),
|
: param_(param),
|
||||||
model_(std::move(model)),
|
model_(std::move(model)),
|
||||||
newtonIterations_(0),
|
nonlinearIterations_(0),
|
||||||
linearIterations_(0)
|
linearIterations_(0),
|
||||||
|
nonlinearIterationsLast_(0),
|
||||||
|
linearIterationsLast_(0)
|
||||||
{
|
{
|
||||||
|
if (!model_) {
|
||||||
|
OPM_THROW(std::logic_error, "Must provide a non-null model argument for NonlinearSolver.");
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
template <class PhysicalModel>
|
template <class PhysicalModel>
|
||||||
unsigned int NewtonSolver<PhysicalModel>::newtonIterations () const
|
unsigned int NonlinearSolver<PhysicalModel>::nonlinearIterations() const
|
||||||
{
|
{
|
||||||
return newtonIterations_;
|
return nonlinearIterations_;
|
||||||
}
|
}
|
||||||
|
|
||||||
template <class PhysicalModel>
|
template <class PhysicalModel>
|
||||||
unsigned int NewtonSolver<PhysicalModel>::linearIterations () const
|
unsigned int NonlinearSolver<PhysicalModel>::linearIterations() const
|
||||||
{
|
{
|
||||||
return linearIterations_;
|
return linearIterations_;
|
||||||
}
|
}
|
||||||
|
|
||||||
template <class PhysicalModel>
|
template <class PhysicalModel>
|
||||||
const PhysicalModel& NewtonSolver<PhysicalModel>::model() const
|
const PhysicalModel& NonlinearSolver<PhysicalModel>::model() const
|
||||||
{
|
{
|
||||||
assert( model_ );
|
|
||||||
return *model_;
|
return *model_;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
template <class PhysicalModel>
|
||||||
|
unsigned int NonlinearSolver<PhysicalModel>::nonlinearIterationsLastStep() const
|
||||||
|
{
|
||||||
|
return nonlinearIterationsLast_;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <class PhysicalModel>
|
||||||
|
unsigned int NonlinearSolver<PhysicalModel>::linearIterationsLastStep() const
|
||||||
|
{
|
||||||
|
return linearIterationsLast_;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
template <class PhysicalModel>
|
template <class PhysicalModel>
|
||||||
int
|
int
|
||||||
NewtonSolver<PhysicalModel>::
|
NonlinearSolver<PhysicalModel>::
|
||||||
step(const double dt,
|
step(const double dt,
|
||||||
ReservoirState& reservoir_state,
|
ReservoirState& reservoir_state,
|
||||||
WellState& well_state)
|
WellState& well_state)
|
||||||
@ -66,57 +83,27 @@ namespace Opm
|
|||||||
// Do model-specific once-per-step calculations.
|
// Do model-specific once-per-step calculations.
|
||||||
model_->prepareStep(dt, reservoir_state, well_state);
|
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;
|
|
||||||
|
|
||||||
// Assemble residual and Jacobian, store residual norms.
|
|
||||||
model_->assemble(reservoir_state, well_state, true);
|
|
||||||
residual_norms_history.push_back(model_->computeResidualNorms());
|
|
||||||
|
|
||||||
// Set up for main Newton loop.
|
|
||||||
double omega = 1.0;
|
|
||||||
int iteration = 0;
|
int iteration = 0;
|
||||||
bool converged = model_->getConvergence(dt, iteration);
|
|
||||||
const int sizeNonLinear = model_->sizeNonLinear();
|
// Let the model do one nonlinear iteration.
|
||||||
V dxOld = V::Zero(sizeNonLinear);
|
|
||||||
bool isOscillate = false;
|
// Set up for main solver loop.
|
||||||
bool isStagnate = false;
|
|
||||||
const enum RelaxType relaxtype = relaxType();
|
|
||||||
int linIters = 0;
|
int linIters = 0;
|
||||||
|
bool converged = false;
|
||||||
|
|
||||||
// ---------- Main Newton loop ----------
|
// ---------- Main nonlinear solver loop ----------
|
||||||
while ( (!converged && (iteration < maxIter())) || (minIter() > iteration)) {
|
do {
|
||||||
// Compute the Newton update to the primary variables.
|
IterationReport report = model_->nonlinearIteration(iteration, dt, *this, reservoir_state, well_state);
|
||||||
V dx = model_->solveJacobianSystem();
|
if (report.failed) {
|
||||||
|
OPM_THROW(Opm::NumericalProblem, "Failed to complete a nonlinear iteration.");
|
||||||
// Store number of linear iterations used.
|
|
||||||
linIters += model_->linearIterationsLastSolve();
|
|
||||||
|
|
||||||
// Stabilize the Newton update.
|
|
||||||
detectNewtonOscillations(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;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
stabilizeNewton(dx, dxOld, omega, relaxtype);
|
if (report.converged) {
|
||||||
|
assert(report.linear_iterations == 0);
|
||||||
// Apply the update, the model may apply model-dependent
|
converged = true;
|
||||||
// limitations and chopping of the update.
|
}
|
||||||
model_->updateState(dx, reservoir_state, well_state);
|
linIters += report.linear_iterations;
|
||||||
|
|
||||||
// Assemble residual and Jacobian, store residual norms.
|
|
||||||
model_->assemble(reservoir_state, well_state, false);
|
|
||||||
residual_norms_history.push_back(model_->computeResidualNorms());
|
|
||||||
|
|
||||||
// increase iteration counter
|
|
||||||
++iteration;
|
++iteration;
|
||||||
|
} while ( (!converged && (iteration <= maxIter())) || (minIter() > iteration));
|
||||||
converged = model_->getConvergence(dt, iteration);
|
|
||||||
}
|
|
||||||
|
|
||||||
if (!converged) {
|
if (!converged) {
|
||||||
if (model_->terminalOutputEnabled()) {
|
if (model_->terminalOutputEnabled()) {
|
||||||
@ -126,9 +113,9 @@ namespace Opm
|
|||||||
}
|
}
|
||||||
|
|
||||||
linearIterations_ += linIters;
|
linearIterations_ += linIters;
|
||||||
newtonIterations_ += iteration;
|
nonlinearIterations_ += iteration - 1; // Since the last one will always be trivial.
|
||||||
linearIterationsLast_ = linIters;
|
linearIterationsLast_ = linIters;
|
||||||
newtonIterationsLast_ = iteration;
|
nonlinearIterationsLast_ = iteration;
|
||||||
|
|
||||||
// Do model-specific post-step actions.
|
// Do model-specific post-step actions.
|
||||||
model_->afterStep(dt, reservoir_state, well_state);
|
model_->afterStep(dt, reservoir_state, well_state);
|
||||||
@ -139,7 +126,7 @@ namespace Opm
|
|||||||
|
|
||||||
|
|
||||||
template <class PhysicalModel>
|
template <class PhysicalModel>
|
||||||
void NewtonSolver<PhysicalModel>::SolverParameters::
|
void NonlinearSolver<PhysicalModel>::SolverParameters::
|
||||||
reset()
|
reset()
|
||||||
{
|
{
|
||||||
// default values for the solver parameters
|
// default values for the solver parameters
|
||||||
@ -152,7 +139,7 @@ namespace Opm
|
|||||||
}
|
}
|
||||||
|
|
||||||
template <class PhysicalModel>
|
template <class PhysicalModel>
|
||||||
NewtonSolver<PhysicalModel>::SolverParameters::
|
NonlinearSolver<PhysicalModel>::SolverParameters::
|
||||||
SolverParameters()
|
SolverParameters()
|
||||||
{
|
{
|
||||||
// set default values
|
// set default values
|
||||||
@ -160,7 +147,7 @@ namespace Opm
|
|||||||
}
|
}
|
||||||
|
|
||||||
template <class PhysicalModel>
|
template <class PhysicalModel>
|
||||||
NewtonSolver<PhysicalModel>::SolverParameters::
|
NonlinearSolver<PhysicalModel>::SolverParameters::
|
||||||
SolverParameters( const parameter::ParameterGroup& param )
|
SolverParameters( const parameter::ParameterGroup& param )
|
||||||
{
|
{
|
||||||
// set default values
|
// set default values
|
||||||
@ -183,9 +170,9 @@ namespace Opm
|
|||||||
|
|
||||||
template <class PhysicalModel>
|
template <class PhysicalModel>
|
||||||
void
|
void
|
||||||
NewtonSolver<PhysicalModel>::detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
|
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
|
bool& oscillate, bool& stagnate) const
|
||||||
{
|
{
|
||||||
// The detection of oscillation in two primary variable results in the report of the detection
|
// The detection of oscillation in two primary variable results in the report of the detection
|
||||||
// of oscillation for the solver.
|
// of oscillation for the solver.
|
||||||
@ -207,7 +194,7 @@ namespace Opm
|
|||||||
const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
|
const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
|
||||||
const double d2 = std::abs((F0[p] - F1[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
|
// Process is 'stagnate' unless at least one phase
|
||||||
// exhibits significant residual change.
|
// exhibits significant residual change.
|
||||||
@ -220,8 +207,7 @@ namespace Opm
|
|||||||
|
|
||||||
template <class PhysicalModel>
|
template <class PhysicalModel>
|
||||||
void
|
void
|
||||||
NewtonSolver<PhysicalModel>::stabilizeNewton(V& dx, V& dxOld, const double omega,
|
NonlinearSolver<PhysicalModel>::stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega) const
|
||||||
const RelaxType relax_type) const
|
|
||||||
{
|
{
|
||||||
// The dxOld is updated with dx.
|
// The dxOld is updated with dx.
|
||||||
// If omega is equal to 1., no relaxtion will be appiled.
|
// If omega is equal to 1., no relaxtion will be appiled.
|
||||||
@ -229,7 +215,7 @@ namespace Opm
|
|||||||
const V tempDxOld = dxOld;
|
const V tempDxOld = dxOld;
|
||||||
dxOld = dx;
|
dxOld = dx;
|
||||||
|
|
||||||
switch (relax_type) {
|
switch (relaxType()) {
|
||||||
case DAMPEN:
|
case DAMPEN:
|
||||||
if (omega == 1.) {
|
if (omega == 1.) {
|
||||||
return;
|
return;
|
@ -26,7 +26,6 @@
|
|||||||
#include <opm/common/ErrorMacros.hpp>
|
#include <opm/common/ErrorMacros.hpp>
|
||||||
|
|
||||||
#include <opm/autodiff/GeoProps.hpp>
|
#include <opm/autodiff/GeoProps.hpp>
|
||||||
#include <opm/autodiff/NewtonSolver.hpp>
|
|
||||||
#include <opm/autodiff/BlackoilModel.hpp>
|
#include <opm/autodiff/BlackoilModel.hpp>
|
||||||
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
|
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
|
||||||
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
|
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
|
||||||
|
@ -106,7 +106,7 @@ namespace Opm
|
|||||||
output_writer_.restore( timer, state, prev_well_state, restorefilename, desiredRestoreStep );
|
output_writer_.restore( timer, state, prev_well_state, restorefilename, desiredRestoreStep );
|
||||||
}
|
}
|
||||||
|
|
||||||
unsigned int totalNewtonIterations = 0;
|
unsigned int totalNonlinearIterations = 0;
|
||||||
unsigned int totalLinearIterations = 0;
|
unsigned int totalLinearIterations = 0;
|
||||||
|
|
||||||
// Main simulation loop.
|
// Main simulation loop.
|
||||||
@ -167,8 +167,8 @@ namespace Opm
|
|||||||
// take time that was used to solve system for this reportStep
|
// take time that was used to solve system for this reportStep
|
||||||
solver_timer.stop();
|
solver_timer.stop();
|
||||||
|
|
||||||
// accumulate the number of Newton and Linear Iterations
|
// accumulate the number of nonlinear and linear Iterations
|
||||||
totalNewtonIterations += solver->newtonIterations();
|
totalNonlinearIterations += solver->nonlinearIterations();
|
||||||
totalLinearIterations += solver->linearIterations();
|
totalLinearIterations += solver->linearIterations();
|
||||||
|
|
||||||
// Report timing.
|
// Report timing.
|
||||||
@ -201,7 +201,7 @@ namespace Opm
|
|||||||
report.pressure_time = stime;
|
report.pressure_time = stime;
|
||||||
report.transport_time = 0.0;
|
report.transport_time = 0.0;
|
||||||
report.total_time = total_timer.secsSinceStart();
|
report.total_time = total_timer.secsSinceStart();
|
||||||
report.total_newton_iterations = totalNewtonIterations;
|
report.total_newton_iterations = totalNonlinearIterations;
|
||||||
report.total_linear_iterations = totalLinearIterations;
|
report.total_linear_iterations = totalLinearIterations;
|
||||||
return report;
|
return report;
|
||||||
}
|
}
|
||||||
|
@ -1,5 +1,5 @@
|
|||||||
/*
|
/*
|
||||||
Copyright 2013 SINTEF ICT, Applied Mathematics.
|
Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
|
||||||
Copyright 2015 Andreas Lauser
|
Copyright 2015 Andreas Lauser
|
||||||
|
|
||||||
This file is part of the Open Porous Media project (OPM).
|
This file is part of the Open Porous Media project (OPM).
|
||||||
@ -21,9 +21,8 @@
|
|||||||
#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
|
#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
|
||||||
#define OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
|
#define OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
|
||||||
|
|
||||||
#include "SimulatorBase.hpp"
|
#include <opm/autodiff/SimulatorBase.hpp>
|
||||||
|
#include <opm/autodiff/NonlinearSolver.hpp>
|
||||||
#include "NewtonSolver.hpp"
|
|
||||||
|
|
||||||
namespace Opm {
|
namespace Opm {
|
||||||
|
|
||||||
@ -38,7 +37,7 @@ struct SimulatorTraits<SimulatorFullyImplicitBlackoil<GridT> >
|
|||||||
typedef BlackoilOutputWriter OutputWriter;
|
typedef BlackoilOutputWriter OutputWriter;
|
||||||
typedef GridT Grid;
|
typedef GridT Grid;
|
||||||
typedef BlackoilModel<Grid> Model;
|
typedef BlackoilModel<Grid> Model;
|
||||||
typedef NewtonSolver<Model> Solver;
|
typedef NonlinearSolver<Model> Solver;
|
||||||
};
|
};
|
||||||
|
|
||||||
/// a simulator for the blackoil model
|
/// a simulator for the blackoil model
|
||||||
|
@ -32,7 +32,7 @@
|
|||||||
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
|
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
|
||||||
#include <opm/autodiff/SolventPropsAdFromDeck.hpp>
|
#include <opm/autodiff/SolventPropsAdFromDeck.hpp>
|
||||||
#include <opm/autodiff/RateConverter.hpp>
|
#include <opm/autodiff/RateConverter.hpp>
|
||||||
#include <opm/autodiff/NewtonSolver.hpp>
|
#include <opm/autodiff/NonlinearSolver.hpp>
|
||||||
#include <opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp>
|
#include <opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp>
|
||||||
|
|
||||||
#include <opm/core/grid.h>
|
#include <opm/core/grid.h>
|
||||||
@ -88,7 +88,7 @@ namespace Opm
|
|||||||
typedef BlackoilOutputWriter OutputWriter;
|
typedef BlackoilOutputWriter OutputWriter;
|
||||||
typedef GridT Grid;
|
typedef GridT Grid;
|
||||||
typedef BlackoilSolventModel<Grid> Model;
|
typedef BlackoilSolventModel<Grid> Model;
|
||||||
typedef NewtonSolver<Model> Solver;
|
typedef NonlinearSolver<Model> Solver;
|
||||||
};
|
};
|
||||||
|
|
||||||
/// Class collecting all necessary components for a blackoil simulation with polymer
|
/// Class collecting all necessary components for a blackoil simulation with polymer
|
||||||
|
Loading…
Reference in New Issue
Block a user