Loading...
Searching...
No Matches
NonlinearSolver.hpp
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
Copyright 2015 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_NON_LINEAR_SOLVER_HPP
#define OPM_NON_LINEAR_SOLVER_HPP
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/models/nonlinear/newtonmethodparams.hpp>
#include <opm/models/nonlinear/newtonmethodproperties.hh>
#include <opm/models/utils/propertysystem.hh>
#include <opm/models/utils/basicproperties.hh>
#include <opm/simulators/timestepping/SimulatorReport.hpp>
#include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
#include <opm/simulators/timestepping/TimeStepControl.hpp>
#include <memory>
namespace Opm::Parameters {
template<class Scalar>
struct NewtonMaxRelax { static constexpr Scalar value = 0.5; };
struct NewtonRelaxationType { static constexpr auto value = "dampen"; };
} // namespace Opm::Parameters
namespace Opm {
// Available relaxation scheme types.
enum class NonlinearRelaxType {
Dampen,
SOR
};
namespace detail {
template<class Scalar>
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);
template <class BVector, class Scalar>
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
const Scalar omega, NonlinearRelaxType relaxType);
}
// Solver parameters controlling nonlinear process.
template<class Scalar>
struct NonlinearSolverParameters
{
NonlinearRelaxType relaxType_;
Scalar relaxMax_;
Scalar relaxIncrement_;
Scalar relaxRelTol_;
NonlinearSolverParameters();
static void registerParameters();
void reset();
};
template <class TypeTag, class PhysicalModel>
class NonlinearSolver
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using SolverParameters = NonlinearSolverParameters<Scalar>;
// --------- Public methods ---------
NonlinearSolver(const SolverParameters& param,
std::unique_ptr<PhysicalModel> model)
: param_(param)
, model_(std::move(model))
, linearizations_(0)
, nonlinearIterations_(0)
, linearIterations_(0)
, wellIterations_(0)
, nonlinearIterationsLast_(0)
, linearIterationsLast_(0)
, wellIterationsLast_(0)
{
if (!model_) {
OPM_THROW(std::logic_error, "Must provide a non-null model argument for NonlinearSolver.");
}
}
SimulatorReportSingle step(const SimulatorTimerInterface& timer, const TimeStepControlInterface *timeStepControl)
{
SimulatorReportSingle report;
report.global_time = timer.simulationTimeElapsed();
report.timestep_length = timer.currentStepLength();
// Do model-specific once-per-step calculations.
report += model_->prepareStep(timer);
int iteration = 0;
// Let the model do one nonlinear iteration.
// Set up for main solver loop.
bool converged = false;
// ---------- Main nonlinear solver loop ----------
do {
try {
// Do the nonlinear step. If we are in a converged state, the
// model will usually do an early return without an expensive
// solve, unless the newton_min_iter_ count has not been reached yet.
auto iterReport = model_->nonlinearIteration(iteration, timer, *this);
iterReport.global_time = timer.simulationTimeElapsed();
report += iterReport;
report.converged = iterReport.converged;
converged = report.converged;
iteration += 1;
}
catch (...) {
// if an iteration fails during a time step, all previous iterations
// count as a failure as well
failureReport_ = report;
failureReport_ += model_->failureReport();
throw;
}
}
while ( (!converged && (iteration <= this->model().param().newton_max_iter_)) ||
(iteration <= this->model().param().newton_min_iter_));
if (!converged) {
failureReport_ = report;
std::string msg = "Solver convergence failure - Failed to complete a time step within ";
msg += std::to_string(model().param().newton_max_iter_) + " iterations.";
OPM_THROW_NOLOG(TooManyIterations, msg);
}
auto relativeChange = model_->relativeChange();
if (timeStepControl && !timeStepControl->timeStepAccepted(relativeChange, timer.currentStepLength())) {
report.converged = false;
report.time_step_rejected = true;
failureReport_ = report;
std::string msg = "Relative change in solution for time step was " + std::to_string(relativeChange);
msg += ", which is larger than the tolerance accepted by the timestepping algorithm.";
OPM_THROW_NOLOG(TimeSteppingBreakdown, msg);
}
// Do model-specific post-step actions.
report += model_->afterStep(timer);
report.converged = true;
return report;
}
const SimulatorReportSingle& failureReport() const
{ return failureReport_; }
int linearizations() const
{ return linearizations_; }
int nonlinearIterations() const
{ return nonlinearIterations_; }
int linearIterations() const
{ return linearIterations_; }
int wellIterations() const
{ return wellIterations_; }
int nonlinearIterationsLastStep() const
{ return nonlinearIterationsLast_; }
int linearIterationsLastStep() const
{ return linearIterationsLast_; }
int wellIterationsLastStep() const
{ return wellIterationsLast_; }
std::vector<std::vector<Scalar> >
computeFluidInPlace(const std::vector<int>& fipnum) const
{ return model_->computeFluidInPlace(fipnum); }
const PhysicalModel& model() const
{ return *model_; }
PhysicalModel& model()
{ return *model_; }
void detectOscillations(const std::vector<std::vector<Scalar>>& residualHistory,
const int it, bool& oscillate, bool& stagnate) const
{
detail::detectOscillations(residualHistory, it, model_->numPhases(),
this->relaxRelTol(), 2, oscillate, stagnate);
}
template <class BVector>
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const Scalar omega) const
{
detail::stabilizeNonlinearUpdate(dx, dxOld, omega, this->relaxType());
}
Scalar relaxMax() const
{ return param_.relaxMax_; }
Scalar relaxIncrement() const
{ return param_.relaxIncrement_; }
NonlinearRelaxType relaxType() const
{ return param_.relaxType_; }
Scalar relaxRelTol() const
{ return param_.relaxRelTol_; }
void setParameters(const SolverParameters& param)
{ param_ = param; }
private:
// --------- Data members ---------
SimulatorReportSingle failureReport_;
SolverParameters param_;
std::unique_ptr<PhysicalModel> model_;
int linearizations_;
int nonlinearIterations_;
int linearIterations_;
int wellIterations_;
int nonlinearIterationsLast_;
int linearIterationsLast_;
int wellIterationsLast_;
};
} // namespace Opm
#endif // OPM_NON_LINEAR_SOLVER_HPP
basicproperties.hh
Defines a type tags and some fundamental properties all models.
Opm::NonlinearSolver::linearIterationsLastStep
int linearIterationsLastStep() const
Number of linear solver iterations used in the last call to step().
Definition NonlinearSolver.hpp:219
Opm::NonlinearSolver::stabilizeNonlinearUpdate
void stabilizeNonlinearUpdate(BVector &dx, BVector &dxOld, const Scalar omega) const
Apply a stabilization to dx, depending on dxOld and relaxation parameters.
Definition NonlinearSolver.hpp:250
Opm::NonlinearSolver::relaxRelTol
Scalar relaxRelTol() const
The relaxation relative tolerance.
Definition NonlinearSolver.hpp:268
Opm::NonlinearSolver::relaxMax
Scalar relaxMax() const
The greatest relaxation factor (i.e. smallest factor) allowed.
Definition NonlinearSolver.hpp:256
Opm::NonlinearSolver::linearizations
int linearizations() const
Number of linearizations used in all calls to step().
Definition NonlinearSolver.hpp:199
Opm::NonlinearSolver::wellIterationsLastStep
int wellIterationsLastStep() const
Number of well iterations used in all calls to step().
Definition NonlinearSolver.hpp:223
Opm::NonlinearSolver::nonlinearIterations
int nonlinearIterations() const
Number of full nonlinear solver iterations used in all calls to step().
Definition NonlinearSolver.hpp:203
Opm::NonlinearSolver::setParameters
void setParameters(const SolverParameters &param)
Set parameters to override those given at construction time.
Definition NonlinearSolver.hpp:272
Opm::NonlinearSolver::relaxType
NonlinearRelaxType relaxType() const
The relaxation type (Dampen or SOR).
Definition NonlinearSolver.hpp:264
Opm::NonlinearSolver::linearIterations
int linearIterations() const
Number of linear solver iterations used in all calls to step().
Definition NonlinearSolver.hpp:207
Opm::NonlinearSolver::nonlinearIterationsLastStep
int nonlinearIterationsLastStep() const
Number of nonlinear solver iterations used in the last call to step().
Definition NonlinearSolver.hpp:215
Opm::NonlinearSolver::failureReport
const SimulatorReportSingle & failureReport() const
return the statistics if the step() method failed
Definition NonlinearSolver.hpp:195
Opm::NonlinearSolver::wellIterations
int wellIterations() const
Number of well iterations used in all calls to step().
Definition NonlinearSolver.hpp:211
Opm::NonlinearSolver< TypeTag, Model >::model
const Model & model() const
Definition NonlinearSolver.hpp:231
Opm::NonlinearSolver::detectOscillations
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
Opm::NonlinearSolver::NonlinearSolver
NonlinearSolver(const SolverParameters &param, std::unique_ptr< PhysicalModel > model)
Construct solver for a given model.
Definition NonlinearSolver.hpp:111
Opm::NonlinearSolver::relaxIncrement
Scalar relaxIncrement() const
The step-change size for the relaxation factor.
Definition NonlinearSolver.hpp:260
Opm::NonlinearSolver::model
PhysicalModel & model()
Mutable reference to physical model.
Definition NonlinearSolver.hpp:235
Opm::SimulatorTimerInterface
Interface class for SimulatorTimer objects, to be improved.
Definition SimulatorTimerInterface.hpp:34
Opm::SimulatorTimerInterface::currentStepLength
virtual double currentStepLength() const =0
Current step length.
Opm::SimulatorTimerInterface::simulationTimeElapsed
virtual double simulationTimeElapsed() const =0
Time elapsed since the start of the simulation until the beginning of the current time step [s].
Opm::TimeStepControlInterface
TimeStepControlInterface.
Definition TimeStepControlInterface.hpp:51
Opm::TimeStepControlInterface::timeStepAccepted
virtual bool timeStepAccepted(const double error, const double timeStepJustCompleted) const =0
For the general 3rd order controller, the internal shifting of errors and time steps happens here,...
Opm
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilbioeffectsmodules.hh:43
Opm::GetPropType
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
propertysystem.hh
The Opm property system, traits with inheritance.
Opm::NonlinearSolverParameters
Definition NonlinearSolver.hpp:79
Opm::Parameters::NewtonMaxRelax
Definition NonlinearSolver.hpp:45
Opm::Parameters::NewtonRelaxationType
Definition NonlinearSolver.hpp:47
Opm::SimulatorReportSingle
A struct for returning timing data from a simulator to its caller.
Definition SimulatorReport.hpp:34