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

Merge pull request #5834 from akva2/split_adaptive_time_stepping

Browse Source 
split AdaptiveTimeStepping to use an impl file
This commit is contained in:
Bård Skaflestad 2025-01-03 10:32:21 +01:00 committed by GitHub
commit 1c6622e528
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
4 changed files with 698 additions and 551 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -928,6 +928,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/linalg/WriteSystemMatrixHelper.hpp opm/simulators/linalg/WriteSystemMatrixHelper.hpp
opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp
opm/simulators/timestepping/AdaptiveTimeStepping.hpp opm/simulators/timestepping/AdaptiveTimeStepping.hpp
opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp
opm/simulators/timestepping/ConvergenceReport.hpp opm/simulators/timestepping/ConvergenceReport.hpp
opm/simulators/timestepping/EclTimeSteppingParams.hpp opm/simulators/timestepping/EclTimeSteppingParams.hpp
opm/simulators/timestepping/TimeStepControl.hpp opm/simulators/timestepping/TimeStepControl.hpp

96
opm/simulators/timestepping/AdaptiveTimeStepping.cpp
View File

@ -138,5 +138,101 @@ void registerAdaptiveParameters()
"can be reduced to based on newton iteration counts"); "can be reduced to based on newton iteration counts");
} }
std::tuple<TimeStepControlType, std::unique_ptr<TimeStepControlInterface>, bool>
createController(const UnitSystem& unitSystem)
{
const double tol = Parameters::Get<Parameters::TimeStepControlTolerance>(); // 1e-1
using RetVal = std::tuple<TimeStepControlType, std::unique_ptr<TimeStepControlInterface>, bool>;
using Func = std::function<RetVal()>;
const auto creators = std::unordered_map<std::string, Func> {
{"pid",
[tol]() {
return RetVal{
TimeStepControlType::PID,
std::make_unique<PIDTimeStepControl>(tol),
false
};
}},
{"pid+iteration",
[tol]() {
const int iterations = Parameters::Get<Parameters::TimeStepControlTargetIterations>(); // 30
const double decayDampingFactor = Parameters::Get<Parameters::TimeStepControlDecayDampingFactor>(); // 1.0
const double growthDampingFactor = Parameters::Get<Parameters::TimeStepControlGrowthDampingFactor>(); // 3.2
return RetVal{
TimeStepControlType::PIDAndIterationCount,
std::make_unique<PIDAndIterationCountTimeStepControl>(iterations,
decayDampingFactor,
growthDampingFactor,
tol),
false
};
}},
{"pid+newtoniteration",
[tol, &unitSystem]() {
const int iterations = Parameters::Get<Parameters::TimeStepControlTargetNewtonIterations>(); // 8
const double decayDampingFactor = Parameters::Get<Parameters::TimeStepControlDecayDampingFactor>(); // 1.0
const double growthDampingFactor = Parameters::Get<Parameters::TimeStepControlGrowthDampingFactor>(); // 3.2
const double nonDimensionalMinTimeStepIterations = Parameters::Get<Parameters::MinTimeStepBasedOnNewtonIterations>(); // 0.0 by default
// the min time step can be reduced by the newton iteration numbers
double minTimeStepReducedByIterations = unitSystem.to_si(UnitSystem::measure::time,
nonDimensionalMinTimeStepIterations);
return RetVal{
TimeStepControlType::PIDAndIterationCount,
std::make_unique<PIDAndIterationCountTimeStepControl>(iterations,
decayDampingFactor,
growthDampingFactor,
tol,
minTimeStepReducedByIterations),
true
};
}},
{"iterationcount",
[]() {
const int iterations = Parameters::Get<Parameters::TimeStepControlTargetIterations>(); // 30
const double decayrate = Parameters::Get<Parameters::TimeStepControlDecayRate>(); // 0.75
const double growthrate = Parameters::Get<Parameters::TimeStepControlGrowthRate>(); // 1.25
return RetVal{
TimeStepControlType::SimpleIterationCount,
std::make_unique<SimpleIterationCountTimeStepControl>(iterations,
decayrate,
growthrate),
false
};
}},
{"newtoniterationcount",
[]() {
const int iterations = Parameters::Get<Parameters::TimeStepControlTargetNewtonIterations>(); // 8
const double decayrate = Parameters::Get<Parameters::TimeStepControlDecayRate>(); // 0.75
const double growthrate = Parameters::Get<Parameters::TimeStepControlGrowthRate>(); // 1.25
return RetVal{
TimeStepControlType::SimpleIterationCount,
std::make_unique<SimpleIterationCountTimeStepControl>(iterations,
decayrate,
growthrate),
true
};
}},
{"hardcoded",
[]() {
const std::string filename = Parameters::Get<Parameters::TimeStepControlFileName>(); // "timesteps"
return RetVal{
TimeStepControlType::HardCodedTimeStep,
std::make_unique<HardcodedTimeStepControl>(filename),
false
};
}},
};
const std::string control = Parameters::Get<Parameters::TimeStepControl>(); // "pid"
const auto it = creators.find(control);
if (it == creators.end()) {
OPM_THROW(std::runtime_error,
"Unsupported time step control selected " + control);
}
// invoke creator
return it->second();
}
} // namespace detail } // namespace detail
} // namespace Opm } // namespace Opm

590
opm/simulators/timestepping/AdaptiveTimeStepping.hpp
View File

@ -6,19 +6,11 @@
#include <dune/common/version.hh> #include <dune/common/version.hh>
#include <dune/istl/istlexception.hh> #include <dune/istl/istlexception.hh>
#include <opm/common/Exceptions.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/OpmLog/OpmLog.hpp> #include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/grid/utility/StopWatch.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/input/eclipse/Units/UnitSystem.hpp>
#include <opm/input/eclipse/Schedule/Tuning.hpp> #include <opm/input/eclipse/Schedule/Tuning.hpp>
#include <opm/models/utils/basicproperties.hh> #include <opm/models/utils/basicproperties.hh>
#include <opm/models/utils/parametersystem.hpp>
#include <opm/models/utils/propertysystem.hh> #include <opm/models/utils/propertysystem.hh>
#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp> #include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
@ -28,19 +20,12 @@
#include <opm/simulators/timestepping/TimeStepControl.hpp> #include <opm/simulators/timestepping/TimeStepControl.hpp>
#include <opm/simulators/timestepping/TimeStepControlInterface.hpp> #include <opm/simulators/timestepping/TimeStepControlInterface.hpp>
#include <opm/simulators/utils/phaseUsageFromDeck.hpp>
#include <fmt/format.h>
#include <algorithm>
#include <cassert>
#include <cmath> #include <cmath>
#include <functional> #include <functional>
#include <memory> #include <memory>
#include <set> #include <set>
#include <sstream>
#include <stdexcept>
#include <string> #include <string>
#include <tuple>
#include <vector> #include <vector>
namespace Opm::Parameters { namespace Opm::Parameters {
@ -67,6 +52,7 @@ struct MinTimeStepBasedOnNewtonIterations { static constexpr double value = 0.0;
namespace Opm { namespace Opm {
class UnitSystem;
struct StepReport; struct StepReport;
namespace detail { namespace detail {
@ -77,6 +63,11 @@ std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr
void registerAdaptiveParameters(); void registerAdaptiveParameters();
std::tuple<TimeStepControlType,
std::unique_ptr<TimeStepControlInterface>,
bool>
createController(const UnitSystem& unitSystem);
} }
// AdaptiveTimeStepping // AdaptiveTimeStepping
@ -114,30 +105,9 @@ void registerAdaptiveParameters();
AdaptiveTimeStepping() = default; AdaptiveTimeStepping() = default;
//! \brief contructor taking parameter object //! \brief contructor taking parameter object
AdaptiveTimeStepping(const UnitSystem& unitSystem, explicit AdaptiveTimeStepping(const UnitSystem& unitSystem,
const double max_next_tstep = -1.0, const double max_next_tstep = -1.0,
const bool terminalOutput = true) const bool terminalOutput = true);
: timeStepControl_()
, restartFactor_(Parameters::Get<Parameters::SolverRestartFactor<Scalar>>()) // 0.33
, growthFactor_(Parameters::Get<Parameters::SolverGrowthFactor<Scalar>>()) // 2.0
, maxGrowth_(Parameters::Get<Parameters::SolverMaxGrowth<Scalar>>()) // 3.0
, maxTimeStep_(Parameters::Get<Parameters::SolverMaxTimeStepInDays<Scalar>>() * 24 * 60 * 60) // 365.25
, minTimeStep_(unitSystem.to_si(UnitSystem::measure::time, Parameters::Get<Parameters::SolverMinTimeStep<Scalar>>())) // 1e-12;
, ignoreConvergenceFailure_(Parameters::Get<Parameters::SolverContinueOnConvergenceFailure>()) // false;
, solverRestartMax_(Parameters::Get<Parameters::SolverMaxRestarts>()) // 10
, solverVerbose_(Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminalOutput) // 2
, timestepVerbose_(Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminalOutput) // 2
, suggestedNextTimestep_((max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>() : max_next_tstep) * 24 * 60 * 60) // 1.0
, fullTimestepInitially_(Parameters::Get<Parameters::FullTimeStepInitially>()) // false
, timestepAfterEvent_(Parameters::Get<Parameters::TimeStepAfterEventInDays<Scalar>>() * 24 * 60 * 60) // 1e30
, useNewtonIteration_(false)
, minTimeStepBeforeShuttingProblematicWells_(Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
{
init_(unitSystem);
}
//! \brief contructor taking parameter object //! \brief contructor taking parameter object
//! \param tuning Pointer to ecl TUNING keyword //! \param tuning Pointer to ecl TUNING keyword
@ -145,31 +115,9 @@ void registerAdaptiveParameters();
AdaptiveTimeStepping(double max_next_tstep, AdaptiveTimeStepping(double max_next_tstep,
const Tuning& tuning, const Tuning& tuning,
const UnitSystem& unitSystem, const UnitSystem& unitSystem,
const bool terminalOutput = true) const bool terminalOutput = true);
: timeStepControl_()
, restartFactor_(tuning.TSFCNV)
, growthFactor_(tuning.TFDIFF)
, maxGrowth_(tuning.TSFMAX)
, maxTimeStep_(tuning.TSMAXZ) // 365.25
, minTimeStep_(tuning.TSFMIN) // 0.1;
, ignoreConvergenceFailure_(true)
, solverRestartMax_(Parameters::Get<Parameters::SolverMaxRestarts>()) // 10
, solverVerbose_(Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminalOutput) // 2
, timestepVerbose_(Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminalOutput) // 2
, suggestedNextTimestep_(max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>() * 24 * 60 * 60 : max_next_tstep) // 1.0
, fullTimestepInitially_(Parameters::Get<Parameters::FullTimeStepInitially>()) // false
, timestepAfterEvent_(tuning.TMAXWC) // 1e30
, useNewtonIteration_(false)
, minTimeStepBeforeShuttingProblematicWells_(Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
{
init_(unitSystem);
}
static void registerParameters() static void registerParameters();
{
registerEclTimeSteppingParameters<Scalar>();
detail::registerAdaptiveParameters();
}
/** \brief step method that acts like the solver::step method /** \brief step method that acts like the solver::step method
in a sub cycle of time steps in a sub cycle of time steps
@ -180,294 +128,7 @@ void registerAdaptiveParameters();
SimulatorReport step(const SimulatorTimer& simulatorTimer, SimulatorReport step(const SimulatorTimer& simulatorTimer,
Solver& solver, Solver& solver,
const bool isEvent, const bool isEvent,
const std::function<bool(const double, const double, const int)> tuningUpdater) const std::function<bool(const double, const double, const int)> tuningUpdater);
{
// Maybe update tuning
tuningUpdater(simulatorTimer.simulationTimeElapsed(), suggestedNextTimestep_, 0);
SimulatorReport report;
const double timestep = simulatorTimer.currentStepLength();
// init last time step as a fraction of the given time step
if (suggestedNextTimestep_ < 0) {
suggestedNextTimestep_ = restartFactor_ * timestep;
}
if (fullTimestepInitially_) {
suggestedNextTimestep_ = timestep;
}
// use seperate time step after event
if (isEvent && timestepAfterEvent_ > 0) {
suggestedNextTimestep_ = timestepAfterEvent_;
}
auto& simulator = solver.model().simulator();
auto& problem = simulator.problem();
// create adaptive step timer with previously used sub step size
AdaptiveSimulatorTimer substepTimer(simulatorTimer, suggestedNextTimestep_, maxTimeStep_);
// counter for solver restarts
int restarts = 0;
// sub step time loop
while (!substepTimer.done()) {
// Maybe update tuning
// get current delta t
auto oldValue = suggestedNextTimestep_;
if (tuningUpdater(substepTimer.simulationTimeElapsed(),
substepTimer.currentStepLength(),
substepTimer.currentStepNum())) {
// Use provideTimeStepEstimate to make we sure don't simulate longer than the report step is.
substepTimer.provideTimeStepEstimate(suggestedNextTimestep_);
suggestedNextTimestep_ = oldValue;
}
const double dt = substepTimer.currentStepLength();
if (timestepVerbose_) {
detail::logTimer(substepTimer);
}
SimulatorReportSingle substepReport;
std::string causeOfFailure;
try {
substepReport = solver.step(substepTimer);
if (solverVerbose_) {
// report number of linear iterations
OpmLog::debug("Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations));
}
}
catch (const TooManyIterations& e) {
substepReport = solver.failureReport();
causeOfFailure = "Solver convergence failure - Iteration limit reached";
logException_(e, solverVerbose_);
// since linearIterations is < 0 this will restart the solver
}
catch (const ConvergenceMonitorFailure& e) {
substepReport = solver.failureReport();
causeOfFailure = "Convergence monitor failure";
}
catch (const LinearSolverProblem& e) {
substepReport = solver.failureReport();
causeOfFailure = "Linear solver convergence failure";
logException_(e, solverVerbose_);
// since linearIterations is < 0 this will restart the solver
}
catch (const NumericalProblem& e) {
substepReport = solver.failureReport();
causeOfFailure = "Solver convergence failure - Numerical problem encountered";
logException_(e, solverVerbose_);
// since linearIterations is < 0 this will restart the solver
}
catch (const std::runtime_error& e) {
substepReport = solver.failureReport();
logException_(e, solverVerbose_);
// also catch linear solver not converged
}
catch (const Dune::ISTLError& e) {
substepReport = solver.failureReport();
logException_(e, solverVerbose_);
// also catch errors in ISTL AMG that occur when time step is too large
}
catch (const Dune::MatrixBlockError& e) {
substepReport = solver.failureReport();
logException_(e, solverVerbose_);
// this can be thrown by ISTL's ILU0 in block mode, yet is not an ISTLError
}
//Pass substep to eclwriter for summary output
simulator.problem().setSubStepReport(substepReport);
report += substepReport;
bool continue_on_uncoverged_solution = ignoreConvergenceFailure_ &&
!substepReport.converged &&
dt <= minTimeStep_;
if (continue_on_uncoverged_solution && solverVerbose_) {
const auto msg = fmt::format(
"Solver failed to converge but timestep "
"{} is smaller or equal to {}\n"
"which is the minimum threshold given "
"by option --solver-min-time-step\n",
dt, minTimeStep_
);
OpmLog::problem(msg);
}
if (substepReport.converged || continue_on_uncoverged_solution) {
// advance by current dt
++substepTimer;
// create object to compute the time error, simply forwards the call to the model
SolutionTimeErrorSolverWrapper<Solver> relativeChange(solver);
// compute new time step estimate
const int iterations = useNewtonIteration_ ? substepReport.total_newton_iterations
: substepReport.total_linear_iterations;
double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
substepTimer.simulationTimeElapsed());
assert(dtEstimate > 0);
// limit the growth of the timestep size by the growth factor
dtEstimate = std::min(dtEstimate, double(maxGrowth_ * dt));
assert(dtEstimate > 0);
// further restrict time step size growth after convergence problems
if (restarts > 0) {
dtEstimate = std::min(growthFactor_ * dt, dtEstimate);
// solver converged, reset restarts counter
restarts = 0;
}
if (timestepVerbose_) {
std::ostringstream ss;
substepReport.reportStep(ss);
OpmLog::info(ss.str());
}
// write data if outputWriter was provided
// if the time step is done we do not need
// to write it as this will be done by the simulator
// anyway.
if (!substepTimer.done()) {
time::StopWatch perfTimer;
perfTimer.start();
problem.writeOutput(true);
report.success.output_write_time += perfTimer.secsSinceStart();
}
// set new time step length
substepTimer.provideTimeStepEstimate(dtEstimate);
report.success.converged = substepTimer.done();
substepTimer.setLastStepFailed(false);
}
else { // in case of no convergence
substepTimer.setLastStepFailed(true);
// If we have restarted (i.e. cut the timestep) too
// many times, we have failed and throw an exception.
if (restarts >= solverRestartMax_) {
const auto msg = fmt::format(
"Solver failed to converge after cutting timestep {} times.",
restarts
);
if (solverVerbose_) {
OpmLog::error(msg);
}
// Use throw directly to prevent file and line
throw TimeSteppingBreakdown{msg};
}
// The new, chopped timestep.
const double newTimeStep = restartFactor_ * dt;
// If we have restarted (i.e. cut the timestep) too
// much, we have failed and throw an exception.
if (newTimeStep < minTimeStep_) {
const auto msg = fmt::format(
"Solver failed to converge after cutting timestep to {}\n"
"which is the minimum threshold given by option --solver-min-time-step\n",
minTimeStep_
);
if (solverVerbose_) {
OpmLog::error(msg);
}
// Use throw directly to prevent file and line
throw TimeSteppingBreakdown{msg};
}
// Define utility function for chopping timestep.
auto chopTimestep = [&]() {
substepTimer.provideTimeStepEstimate(newTimeStep);
if (solverVerbose_) {
const auto msg = fmt::format(
"{}\nTimestep chopped to {} days\n",
causeOfFailure,
std::to_string(unit::convert::to(substepTimer.currentStepLength(), unit::day))
);
OpmLog::problem(msg);
}
++restarts;
};
const double minimumChoppedTimestep = minTimeStepBeforeShuttingProblematicWells_;
if (newTimeStep > minimumChoppedTimestep) {
chopTimestep();
} else {
// We are below the threshold, and will check if there are any
// wells we should close rather than chopping again.
std::set<std::string> failing_wells = detail::consistentlyFailingWells(solver.model().stepReports());
if (failing_wells.empty()) {
// Found no wells to close, chop the timestep as above.
chopTimestep();
} else {
// Close all consistently failing wells that are not under group control
std::vector<std::string> shut_wells;
for (const auto& well : failing_wells) {
bool was_shut = solver.model().wellModel().forceShutWellByName(
well, substepTimer.simulationTimeElapsed(), /*dont_shut_grup_wells =*/ true);
if (was_shut) {
shut_wells.push_back(well);
}
}
// If no wells are closed we also try to shut wells under group control
if (shut_wells.empty()) {
for (const auto& well : failing_wells) {
bool was_shut = solver.model().wellModel().forceShutWellByName(
well, substepTimer.simulationTimeElapsed(), /*dont_shut_grup_wells =*/ false);
if (was_shut) {
shut_wells.push_back(well);
}
}
}
// If still no wells are closed we must fall back to chopping again
if (shut_wells.empty()) {
chopTimestep();
} else {
substepTimer.provideTimeStepEstimate(dt);
if (solverVerbose_) {
std::string msg;
msg = "\nProblematic well(s) were shut: ";
for (const auto& well : shut_wells) {
msg += well;
msg += " ";
}
msg += "(retrying timestep)\n";
OpmLog::problem(msg);
}
}
}
}
}
problem.setNextTimeStepSize(substepTimer.currentStepLength());
}
// store estimated time step for next reportStep
suggestedNextTimestep_ = substepTimer.currentStepLength();
if (timestepVerbose_) {
std::ostringstream ss;
substepTimer.report(ss);
ss << "Suggested next step size = " << unit::convert::to(suggestedNextTimestep_, unit::day) << " (days)" << std::endl;
OpmLog::debug(ss.str());
}
if (! std::isfinite(suggestedNextTimestep_)) { // check for NaN
suggestedNextTimestep_ = timestep;
}
return report;
}
/** \brief Returns the simulator report for the failed substeps of the last /** \brief Returns the simulator report for the failed substeps of the last
* report step. * report step.
@ -478,142 +139,24 @@ void registerAdaptiveParameters();
void setSuggestedNextStep(const double x) void setSuggestedNextStep(const double x)
{ suggestedNextTimestep_ = x; } { suggestedNextTimestep_ = x; }
void updateTUNING(double max_next_tstep, const Tuning& tuning) void updateTUNING(double max_next_tstep, const Tuning& tuning);
{
restartFactor_ = tuning.TSFCNV;
growthFactor_ = tuning.TFDIFF;
maxGrowth_ = tuning.TSFMAX;
maxTimeStep_ = tuning.TSMAXZ;
updateNEXTSTEP(max_next_tstep);
timestepAfterEvent_ = tuning.TMAXWC;
}
void updateNEXTSTEP(double max_next_tstep) void updateNEXTSTEP(double max_next_tstep);
{
// \Note Only update next suggested step if TSINIT was explicitly set in TUNING or NEXTSTEP is active.
if (max_next_tstep > 0) {
suggestedNextTimestep_ = max_next_tstep;
}
}
template<class Serializer> template<class Serializer>
void serializeOp(Serializer& serializer) void serializeOp(Serializer& serializer);
{
serializer(timeStepControlType_);
switch (timeStepControlType_) {
case TimeStepControlType::HardCodedTimeStep:
allocAndSerialize<HardcodedTimeStepControl>(serializer);
break;
case TimeStepControlType::PIDAndIterationCount:
allocAndSerialize<PIDAndIterationCountTimeStepControl>(serializer);
break;
case TimeStepControlType::SimpleIterationCount:
allocAndSerialize<SimpleIterationCountTimeStepControl>(serializer);
break;
case TimeStepControlType::PID:
allocAndSerialize<PIDTimeStepControl>(serializer);
break;
}
serializer(restartFactor_);
serializer(growthFactor_);
serializer(maxGrowth_);
serializer(maxTimeStep_);
serializer(minTimeStep_);
serializer(ignoreConvergenceFailure_);
serializer(solverRestartMax_);
serializer(solverVerbose_);
serializer(timestepVerbose_);
serializer(suggestedNextTimestep_);
serializer(fullTimestepInitially_);
serializer(timestepAfterEvent_);
serializer(useNewtonIteration_);
serializer(minTimeStepBeforeShuttingProblematicWells_);
}
static AdaptiveTimeStepping<TypeTag> serializationTestObjectHardcoded() static AdaptiveTimeStepping<TypeTag> serializationTestObjectHardcoded();
{ static AdaptiveTimeStepping<TypeTag> serializationTestObjectPID();
return serializationTestObject_<HardcodedTimeStepControl>(); static AdaptiveTimeStepping<TypeTag> serializationTestObjectPIDIt();
} static AdaptiveTimeStepping<TypeTag> serializationTestObjectSimple();
static AdaptiveTimeStepping<TypeTag> serializationTestObjectPID() bool operator==(const AdaptiveTimeStepping<TypeTag>& rhs) const;
{
return serializationTestObject_<PIDTimeStepControl>();
}
static AdaptiveTimeStepping<TypeTag> serializationTestObjectPIDIt()
{
return serializationTestObject_<PIDAndIterationCountTimeStepControl>();
}
static AdaptiveTimeStepping<TypeTag> serializationTestObjectSimple()
{
return serializationTestObject_<SimpleIterationCountTimeStepControl>();
}
bool operator==(const AdaptiveTimeStepping<TypeTag>& rhs) const
{
if (timeStepControlType_ != rhs.timeStepControlType_ ||
(timeStepControl_ && !rhs.timeStepControl_) ||
(!timeStepControl_ && rhs.timeStepControl_)) {
return false;
}
bool result = false;
switch (timeStepControlType_) {
case TimeStepControlType::HardCodedTimeStep:
result = castAndComp<HardcodedTimeStepControl>(rhs);
break;
case TimeStepControlType::PIDAndIterationCount:
result = castAndComp<PIDAndIterationCountTimeStepControl>(rhs);
break;
case TimeStepControlType::SimpleIterationCount:
result = castAndComp<SimpleIterationCountTimeStepControl>(rhs);
break;
case TimeStepControlType::PID:
result = castAndComp<PIDTimeStepControl>(rhs);
break;
}
return result &&
this->restartFactor_ == rhs.restartFactor_ &&
this->growthFactor_ == rhs.growthFactor_ &&
this->maxGrowth_ == rhs.maxGrowth_ &&
this->maxTimeStep_ == rhs.maxTimeStep_ &&
this->minTimeStep_ == rhs.minTimeStep_ &&
this->ignoreConvergenceFailure_ == rhs.ignoreConvergenceFailure_ &&
this->solverRestartMax_== rhs.solverRestartMax_ &&
this->solverVerbose_ == rhs.solverVerbose_ &&
this->fullTimestepInitially_ == rhs.fullTimestepInitially_ &&
this->timestepAfterEvent_ == rhs.timestepAfterEvent_ &&
this->useNewtonIteration_ == rhs.useNewtonIteration_ &&
this->minTimeStepBeforeShuttingProblematicWells_ ==
rhs.minTimeStepBeforeShuttingProblematicWells_;
}
private: private:
template<class Controller> template<class Controller>
static AdaptiveTimeStepping<TypeTag> serializationTestObject_() static AdaptiveTimeStepping<TypeTag> serializationTestObject_();
{
AdaptiveTimeStepping<TypeTag> result;
result.restartFactor_ = 1.0;
result.growthFactor_ = 2.0;
result.maxGrowth_ = 3.0;
result.maxTimeStep_ = 4.0;
result.minTimeStep_ = 5.0;
result.ignoreConvergenceFailure_ = true;
result.solverRestartMax_ = 6;
result.solverVerbose_ = true;
result.timestepVerbose_ = true;
result.suggestedNextTimestep_ = 7.0;
result.fullTimestepInitially_ = true;
result.useNewtonIteration_ = true;
result.minTimeStepBeforeShuttingProblematicWells_ = 9.0;
result.timeStepControlType_ = Controller::Type;
result.timeStepControl_ = std::make_unique<Controller>(Controller::serializationTestObject());
return result;
}
template<class T, class Serializer> template<class T, class Serializer>
void allocAndSerialize(Serializer& serializer) void allocAndSerialize(Serializer& serializer)
{ {
@ -632,82 +175,29 @@ void registerAdaptiveParameters();
} }
protected: protected:
void init_(const UnitSystem& unitSystem) void init_(const UnitSystem& unitSystem);
{
// valid are "pid" and "pid+iteration"
std::string control = Parameters::Get<Parameters::TimeStepControl>(); // "pid"
const double tol = Parameters::Get<Parameters::TimeStepControlTolerance>(); // 1e-1
if (control == "pid") {
timeStepControl_ = std::make_unique<PIDTimeStepControl>(tol);
timeStepControlType_ = TimeStepControlType::PID;
}
else if (control == "pid+iteration") {
const int iterations = Parameters::Get<Parameters::TimeStepControlTargetIterations>(); // 30
const double decayDampingFactor = Parameters::Get<Parameters::TimeStepControlDecayDampingFactor>(); // 1.0
const double growthDampingFactor = Parameters::Get<Parameters::TimeStepControlGrowthDampingFactor>(); // 3.2
timeStepControl_ = std::make_unique<PIDAndIterationCountTimeStepControl>(iterations, decayDampingFactor, growthDampingFactor, tol);
timeStepControlType_ = TimeStepControlType::PIDAndIterationCount;
}
else if (control == "pid+newtoniteration") {
const int iterations = Parameters::Get<Parameters::TimeStepControlTargetNewtonIterations>(); // 8
const double decayDampingFactor = Parameters::Get<Parameters::TimeStepControlDecayDampingFactor>(); // 1.0
const double growthDampingFactor = Parameters::Get<Parameters::TimeStepControlGrowthDampingFactor>(); // 3.2
const double nonDimensionalMinTimeStepIterations = Parameters::Get<Parameters::MinTimeStepBasedOnNewtonIterations>(); // 0.0 by default
// the min time step can be reduced by the newton iteration numbers
double minTimeStepReducedByIterations = unitSystem.to_si(UnitSystem::measure::time, nonDimensionalMinTimeStepIterations);
timeStepControl_ = std::make_unique<PIDAndIterationCountTimeStepControl>(iterations, decayDampingFactor,
growthDampingFactor, tol, minTimeStepReducedByIterations);
timeStepControlType_ = TimeStepControlType::PIDAndIterationCount;
useNewtonIteration_ = true;
}
else if (control == "iterationcount") {
const int iterations = Parameters::Get<Parameters::TimeStepControlTargetIterations>(); // 30
const double decayrate = Parameters::Get<Parameters::TimeStepControlDecayRate>(); // 0.75
const double growthrate = Parameters::Get<Parameters::TimeStepControlGrowthRate>(); // 1.25
timeStepControl_ = std::make_unique<SimpleIterationCountTimeStepControl>(iterations, decayrate, growthrate);
timeStepControlType_ = TimeStepControlType::SimpleIterationCount;
}
else if (control == "newtoniterationcount") {
const int iterations = Parameters::Get<Parameters::TimeStepControlTargetNewtonIterations>(); // 8
const double decayrate = Parameters::Get<Parameters::TimeStepControlDecayRate>(); // 0.75
const double growthrate = Parameters::Get<Parameters::TimeStepControlGrowthRate>(); // 1.25
timeStepControl_ = std::make_unique<SimpleIterationCountTimeStepControl>(iterations, decayrate, growthrate);
useNewtonIteration_ = true;
timeStepControlType_ = TimeStepControlType::SimpleIterationCount;
}
else if (control == "hardcoded") {
const std::string filename = Parameters::Get<Parameters::TimeStepControlFileName>(); // "timesteps"
timeStepControl_ = std::make_unique<HardcodedTimeStepControl>(filename);
timeStepControlType_ = TimeStepControlType::HardCodedTimeStep;
}
else
OPM_THROW(std::runtime_error,
"Unsupported time step control selected " + control);
// make sure growth factor is something reasonable
assert(growthFactor_ >= 1.0);
}
using TimeStepController = std::unique_ptr<TimeStepControlInterface>; using TimeStepController = std::unique_ptr<TimeStepControlInterface>;
TimeStepControlType timeStepControlType_; //!< type of time step control object TimeStepControlType timeStepControlType_{TimeStepControlType::PIDAndIterationCount}; //!< type of time step control object
TimeStepController timeStepControl_; //!< time step control object TimeStepController timeStepControl_{}; //!< time step control object
double restartFactor_; //!< factor to multiply time step with when solver fails to converge double restartFactor_{}; //!< factor to multiply time step with when solver fails to converge
double growthFactor_; //!< factor to multiply time step when solver recovered from failed convergence double growthFactor_{}; //!< factor to multiply time step when solver recovered from failed convergence
double maxGrowth_; //!< factor that limits the maximum growth of a time step double maxGrowth_{}; //!< factor that limits the maximum growth of a time step
double maxTimeStep_; //!< maximal allowed time step size in days double maxTimeStep_{}; //!< maximal allowed time step size in days
double minTimeStep_; //!< minimal allowed time step size before throwing double minTimeStep_{}; //!< minimal allowed time step size before throwing
bool ignoreConvergenceFailure_; //!< continue instead of stop when minimum time step is reached bool ignoreConvergenceFailure_{false}; //!< continue instead of stop when minimum time step is reached
int solverRestartMax_; //!< how many restart of solver are allowed int solverRestartMax_{}; //!< how many restart of solver are allowed
bool solverVerbose_; //!< solver verbosity bool solverVerbose_{false}; //!< solver verbosity
bool timestepVerbose_; //!< timestep verbosity bool timestepVerbose_{false}; //!< timestep verbosity
double suggestedNextTimestep_; //!< suggested size of next timestep double suggestedNextTimestep_{}; //!< suggested size of next timestep
bool fullTimestepInitially_; //!< beginning with the size of the time step from data file bool fullTimestepInitially_{false}; //!< beginning with the size of the time step from data file
double timestepAfterEvent_; //!< suggested size of timestep after an event double timestepAfterEvent_{}; //!< suggested size of timestep after an event
bool useNewtonIteration_; //!< use newton iteration count for adaptive time step control bool useNewtonIteration_{false}; //!< use newton iteration count for adaptive time step control
double minTimeStepBeforeShuttingProblematicWells_; //! < shut problematic wells when time step size in days are less than this double minTimeStepBeforeShuttingProblematicWells_{}; //! < shut problematic wells when time step size in days are less than this
}; };
} }
#include <opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp>
#endif // OPM_ADAPTIVE_TIME_STEPPING_HPP #endif // OPM_ADAPTIVE_TIME_STEPPING_HPP

560
opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp Normal file
View File

@ -0,0 +1,560 @@
/*
*/
#ifndef OPM_ADAPTIVE_TIME_STEPPING_IMPL_HPP
#define OPM_ADAPTIVE_TIME_STEPPING_IMPL_HPP
// Improve IDE experience
#ifndef OPM_ADAPTIVE_TIME_STEPPING_HPP
#include <config.h>
#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
#endif
#include <opm/common/Exceptions.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/grid/utility/StopWatch.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/input/eclipse/Units/UnitSystem.hpp>
#include <opm/models/utils/parametersystem.hpp>
#include <algorithm>
#include <cassert>
#include <cmath>
#include <sstream>
#include <stdexcept>
#include <fmt/format.h>
namespace Opm {
template<class TypeTag>
AdaptiveTimeStepping<TypeTag>::
AdaptiveTimeStepping(const UnitSystem& unitSystem,
const double max_next_tstep,
const bool terminalOutput)
: timeStepControl_()
, restartFactor_(Parameters::Get<Parameters::SolverRestartFactor<Scalar>>()) // 0.33
, growthFactor_(Parameters::Get<Parameters::SolverGrowthFactor<Scalar>>()) // 2.0
, maxGrowth_(Parameters::Get<Parameters::SolverMaxGrowth<Scalar>>()) // 3.0
, maxTimeStep_(Parameters::Get<Parameters::SolverMaxTimeStepInDays<Scalar>>() * 24 * 60 * 60) // 365.25
, minTimeStep_(unitSystem.to_si(UnitSystem::measure::time, Parameters::Get<Parameters::SolverMinTimeStep<Scalar>>())) // 1e-12;
, ignoreConvergenceFailure_(Parameters::Get<Parameters::SolverContinueOnConvergenceFailure>()) // false;
, solverRestartMax_(Parameters::Get<Parameters::SolverMaxRestarts>()) // 10
, solverVerbose_(Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminalOutput) // 2
, timestepVerbose_(Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminalOutput) // 2
, suggestedNextTimestep_((max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>() : max_next_tstep) * 24 * 60 * 60) // 1.0
, fullTimestepInitially_(Parameters::Get<Parameters::FullTimeStepInitially>()) // false
, timestepAfterEvent_(Parameters::Get<Parameters::TimeStepAfterEventInDays<Scalar>>() * 24 * 60 * 60) // 1e30
, useNewtonIteration_(false)
, minTimeStepBeforeShuttingProblematicWells_(Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
{
init_(unitSystem);
}
template<class TypeTag>
AdaptiveTimeStepping<TypeTag>::
AdaptiveTimeStepping(double max_next_tstep,
const Tuning& tuning,
const UnitSystem& unitSystem,
const bool terminalOutput)
: timeStepControl_()
, restartFactor_(tuning.TSFCNV)
, growthFactor_(tuning.TFDIFF)
, maxGrowth_(tuning.TSFMAX)
, maxTimeStep_(tuning.TSMAXZ) // 365.25
, minTimeStep_(tuning.TSFMIN) // 0.1;
, ignoreConvergenceFailure_(true)
, solverRestartMax_(Parameters::Get<Parameters::SolverMaxRestarts>()) // 10
, solverVerbose_(Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminalOutput) // 2
, timestepVerbose_(Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminalOutput) // 2
, suggestedNextTimestep_(max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>() * 24 * 60 * 60 : max_next_tstep) // 1.0
, fullTimestepInitially_(Parameters::Get<Parameters::FullTimeStepInitially>()) // false
, timestepAfterEvent_(tuning.TMAXWC) // 1e30
, useNewtonIteration_(false)
, minTimeStepBeforeShuttingProblematicWells_(Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
{
init_(unitSystem);
}
template<class TypeTag>
void AdaptiveTimeStepping<TypeTag>::registerParameters()
{
registerEclTimeSteppingParameters<Scalar>();
detail::registerAdaptiveParameters();
}
template<class TypeTag>
template<class Solver>
SimulatorReport
AdaptiveTimeStepping<TypeTag>::
step(const SimulatorTimer& simulatorTimer,
Solver& solver,
const bool isEvent,
const std::function<bool(const double, const double, const int)> tuningUpdater)
{
// Maybe update tuning
tuningUpdater(simulatorTimer.simulationTimeElapsed(), suggestedNextTimestep_, 0);
SimulatorReport report;
const double timestep = simulatorTimer.currentStepLength();
// init last time step as a fraction of the given time step
if (suggestedNextTimestep_ < 0) {
suggestedNextTimestep_ = restartFactor_ * timestep;
}
if (fullTimestepInitially_) {
suggestedNextTimestep_ = timestep;
}
// use seperate time step after event
if (isEvent && timestepAfterEvent_ > 0) {
suggestedNextTimestep_ = timestepAfterEvent_;
}
auto& simulator = solver.model().simulator();
auto& problem = simulator.problem();
// create adaptive step timer with previously used sub step size
AdaptiveSimulatorTimer substepTimer(simulatorTimer, suggestedNextTimestep_, maxTimeStep_);
// counter for solver restarts
int restarts = 0;
// sub step time loop
while (!substepTimer.done()) {
// Maybe update tuning
// get current delta t
auto oldValue = suggestedNextTimestep_;
if (tuningUpdater(substepTimer.simulationTimeElapsed(),
substepTimer.currentStepLength(),
substepTimer.currentStepNum())) {
// Use provideTimeStepEstimate to make we sure don't simulate longer than the report step is.
substepTimer.provideTimeStepEstimate(suggestedNextTimestep_);
suggestedNextTimestep_ = oldValue;
}
const double dt = substepTimer.currentStepLength();
if (timestepVerbose_) {
detail::logTimer(substepTimer);
}
SimulatorReportSingle substepReport;
std::string causeOfFailure;
try {
substepReport = solver.step(substepTimer);
if (solverVerbose_) {
// report number of linear iterations
OpmLog::debug("Overall linear iterations used: " +
std::to_string(substepReport.total_linear_iterations));
}
}
catch (const TooManyIterations& e) {
substepReport = solver.failureReport();
causeOfFailure = "Solver convergence failure - Iteration limit reached";
logException_(e, solverVerbose_);
// since linearIterations is < 0 this will restart the solver
}
catch (const ConvergenceMonitorFailure& e) {
substepReport = solver.failureReport();
causeOfFailure = "Convergence monitor failure";
}
catch (const LinearSolverProblem& e) {
substepReport = solver.failureReport();
causeOfFailure = "Linear solver convergence failure";
logException_(e, solverVerbose_);
// since linearIterations is < 0 this will restart the solver
}
catch (const NumericalProblem& e) {
substepReport = solver.failureReport();
causeOfFailure = "Solver convergence failure - Numerical problem encountered";
logException_(e, solverVerbose_);
// since linearIterations is < 0 this will restart the solver
}
catch (const std::runtime_error& e) {
substepReport = solver.failureReport();
logException_(e, solverVerbose_);
// also catch linear solver not converged
}
catch (const Dune::ISTLError& e) {
substepReport = solver.failureReport();
logException_(e, solverVerbose_);
// also catch errors in ISTL AMG that occur when time step is too large
}
catch (const Dune::MatrixBlockError& e) {
substepReport = solver.failureReport();
logException_(e, solverVerbose_);
// this can be thrown by ISTL's ILU0 in block mode, yet is not an ISTLError
}
//Pass substep to eclwriter for summary output
simulator.problem().setSubStepReport(substepReport);
report += substepReport;
bool continue_on_uncoverged_solution = ignoreConvergenceFailure_ &&
!substepReport.converged &&
dt <= minTimeStep_;
if (continue_on_uncoverged_solution && solverVerbose_) {
const auto msg = fmt::format(
"Solver failed to converge but timestep "
"{} is smaller or equal to {}\n"
"which is the minimum threshold given "
"by option --solver-min-time-step\n",
dt, minTimeStep_
);
OpmLog::problem(msg);
}
if (substepReport.converged || continue_on_uncoverged_solution) {
// advance by current dt
++substepTimer;
// create object to compute the time error, simply forwards the call to the model
SolutionTimeErrorSolverWrapper<Solver> relativeChange(solver);
// compute new time step estimate
const int iterations = useNewtonIteration_ ? substepReport.total_newton_iterations
: substepReport.total_linear_iterations;
double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
substepTimer.simulationTimeElapsed());
assert(dtEstimate > 0);
// limit the growth of the timestep size by the growth factor
dtEstimate = std::min(dtEstimate, double(maxGrowth_ * dt));
assert(dtEstimate > 0);
// further restrict time step size growth after convergence problems
if (restarts > 0) {
dtEstimate = std::min(growthFactor_ * dt, dtEstimate);
// solver converged, reset restarts counter
restarts = 0;
}
if (timestepVerbose_) {
std::ostringstream ss;
substepReport.reportStep(ss);
OpmLog::info(ss.str());
}
// write data if outputWriter was provided
// if the time step is done we do not need
// to write it as this will be done by the simulator
// anyway.
if (!substepTimer.done()) {
time::StopWatch perfTimer;
perfTimer.start();
problem.writeOutput(true);
report.success.output_write_time += perfTimer.secsSinceStart();
}
// set new time step length
substepTimer.provideTimeStepEstimate(dtEstimate);
report.success.converged = substepTimer.done();
substepTimer.setLastStepFailed(false);
}
else { // in case of no convergence
substepTimer.setLastStepFailed(true);
// If we have restarted (i.e. cut the timestep) too
// many times, we have failed and throw an exception.
if (restarts >= solverRestartMax_) {
const auto msg = fmt::format(
"Solver failed to converge after cutting timestep {} times.",
restarts
);
if (solverVerbose_) {
OpmLog::error(msg);
}
// Use throw directly to prevent file and line
throw TimeSteppingBreakdown{msg};
}
// The new, chopped timestep.
const double newTimeStep = restartFactor_ * dt;
// If we have restarted (i.e. cut the timestep) too
// much, we have failed and throw an exception.
if (newTimeStep < minTimeStep_) {
const auto msg = fmt::format(
"Solver failed to converge after cutting timestep to {}\n"
"which is the minimum threshold given by option --solver-min-time-step\n",
minTimeStep_
);
if (solverVerbose_) {
OpmLog::error(msg);
}
// Use throw directly to prevent file and line
throw TimeSteppingBreakdown{msg};
}
// Define utility function for chopping timestep.
auto chopTimestep = [&]() {
substepTimer.provideTimeStepEstimate(newTimeStep);
if (solverVerbose_) {
const auto msg = fmt::format(
"{}\nTimestep chopped to {} days\n",
causeOfFailure,
unit::convert::to(substepTimer.currentStepLength(), unit::day)
);
OpmLog::problem(msg);
}
++restarts;
};
const double minimumChoppedTimestep = minTimeStepBeforeShuttingProblematicWells_;
if (newTimeStep > minimumChoppedTimestep) {
chopTimestep();
} else {
// We are below the threshold, and will check if there are any
// wells we should close rather than chopping again.
std::set<std::string> failing_wells = detail::consistentlyFailingWells(solver.model().stepReports());
if (failing_wells.empty()) {
// Found no wells to close, chop the timestep as above.
chopTimestep();
} else {
// Close all consistently failing wells that are not under group control
std::vector<std::string> shut_wells;
for (const auto& well : failing_wells) {
bool was_shut = solver.model().wellModel().forceShutWellByName(
well, substepTimer.simulationTimeElapsed(), /*dont_shut_grup_wells =*/ true);
if (was_shut) {
shut_wells.push_back(well);
}
}
// If no wells are closed we also try to shut wells under group control
if (shut_wells.empty()) {
for (const auto& well : failing_wells) {
bool was_shut = solver.model().wellModel().forceShutWellByName(
well, substepTimer.simulationTimeElapsed(), /*dont_shut_grup_wells =*/ false);
if (was_shut) {
shut_wells.push_back(well);
}
}
}
// If still no wells are closed we must fall back to chopping again
if (shut_wells.empty()) {
chopTimestep();
} else {
substepTimer.provideTimeStepEstimate(dt);
if (solverVerbose_) {
const std::string msg =
fmt::format("\nProblematic well(s) were shut: {}"
"(retrying timestep)\n",
fmt::join(shut_wells, " "));
OpmLog::problem(msg);
}
}
}
}
}
problem.setNextTimeStepSize(substepTimer.currentStepLength());
}
// store estimated time step for next reportStep
suggestedNextTimestep_ = substepTimer.currentStepLength();
if (timestepVerbose_) {
std::ostringstream ss;
substepTimer.report(ss);
ss << "Suggested next step size = " << unit::convert::to(suggestedNextTimestep_, unit::day) << " (days)" << std::endl;
OpmLog::debug(ss.str());
}
if (! std::isfinite(suggestedNextTimestep_)) { // check for NaN
suggestedNextTimestep_ = timestep;
}
return report;
}
template<class TypeTag>
void AdaptiveTimeStepping<TypeTag>::
updateTUNING(double max_next_tstep, const Tuning& tuning)
{
restartFactor_ = tuning.TSFCNV;
growthFactor_ = tuning.TFDIFF;
maxGrowth_ = tuning.TSFMAX;
maxTimeStep_ = tuning.TSMAXZ;
updateNEXTSTEP(max_next_tstep);
timestepAfterEvent_ = tuning.TMAXWC;
}
template<class TypeTag>
void AdaptiveTimeStepping<TypeTag>::
updateNEXTSTEP(double max_next_tstep)
{
// \Note Only update next suggested step if TSINIT was explicitly set in TUNING or NEXTSTEP is active.
if (max_next_tstep > 0) {
suggestedNextTimestep_ = max_next_tstep;
}
}
template<class TypeTag>
template<class Serializer>
void AdaptiveTimeStepping<TypeTag>::
serializeOp(Serializer& serializer)
{
serializer(timeStepControlType_);
switch (timeStepControlType_) {
case TimeStepControlType::HardCodedTimeStep:
allocAndSerialize<HardcodedTimeStepControl>(serializer);
break;
case TimeStepControlType::PIDAndIterationCount:
allocAndSerialize<PIDAndIterationCountTimeStepControl>(serializer);
break;
case TimeStepControlType::SimpleIterationCount:
allocAndSerialize<SimpleIterationCountTimeStepControl>(serializer);
break;
case TimeStepControlType::PID:
allocAndSerialize<PIDTimeStepControl>(serializer);
break;
}
serializer(restartFactor_);
serializer(growthFactor_);
serializer(maxGrowth_);
serializer(maxTimeStep_);
serializer(minTimeStep_);
serializer(ignoreConvergenceFailure_);
serializer(solverRestartMax_);
serializer(solverVerbose_);
serializer(timestepVerbose_);
serializer(suggestedNextTimestep_);
serializer(fullTimestepInitially_);
serializer(timestepAfterEvent_);
serializer(useNewtonIteration_);
serializer(minTimeStepBeforeShuttingProblematicWells_);
}
template<class TypeTag>
AdaptiveTimeStepping<TypeTag>
AdaptiveTimeStepping<TypeTag>::
serializationTestObjectHardcoded()
{
return serializationTestObject_<HardcodedTimeStepControl>();
}
template<class TypeTag>
AdaptiveTimeStepping<TypeTag>
AdaptiveTimeStepping<TypeTag>::
serializationTestObjectPID()
{
return serializationTestObject_<PIDTimeStepControl>();
}
template<class TypeTag>
AdaptiveTimeStepping<TypeTag>
AdaptiveTimeStepping<TypeTag>::
serializationTestObjectPIDIt()
{
return serializationTestObject_<PIDAndIterationCountTimeStepControl>();
}
template<class TypeTag>
AdaptiveTimeStepping<TypeTag>
AdaptiveTimeStepping<TypeTag>::
serializationTestObjectSimple()
{
return serializationTestObject_<SimpleIterationCountTimeStepControl>();
}
template<class TypeTag>
bool
AdaptiveTimeStepping<TypeTag>::
operator==(const AdaptiveTimeStepping<TypeTag>& rhs) const
{
if (timeStepControlType_ != rhs.timeStepControlType_ ||
(timeStepControl_ && !rhs.timeStepControl_) ||
(!timeStepControl_ && rhs.timeStepControl_)) {
return false;
}
bool result = false;
switch (timeStepControlType_) {
case TimeStepControlType::HardCodedTimeStep:
result = castAndComp<HardcodedTimeStepControl>(rhs);
break;
case TimeStepControlType::PIDAndIterationCount:
result = castAndComp<PIDAndIterationCountTimeStepControl>(rhs);
break;
case TimeStepControlType::SimpleIterationCount:
result = castAndComp<SimpleIterationCountTimeStepControl>(rhs);
break;
case TimeStepControlType::PID:
result = castAndComp<PIDTimeStepControl>(rhs);
break;
}
return result
&& this->restartFactor_ == rhs.restartFactor_
&& this->growthFactor_ == rhs.growthFactor_
&& this->maxGrowth_ == rhs.maxGrowth_
&& this->maxTimeStep_ == rhs.maxTimeStep_
&& this->minTimeStep_ == rhs.minTimeStep_
&& this->ignoreConvergenceFailure_ == rhs.ignoreConvergenceFailure_
&& this->solverRestartMax_== rhs.solverRestartMax_
&& this->solverVerbose_ == rhs.solverVerbose_
&& this->fullTimestepInitially_ == rhs.fullTimestepInitially_
&& this->timestepAfterEvent_ == rhs.timestepAfterEvent_
&& this->useNewtonIteration_ == rhs.useNewtonIteration_
&& this->minTimeStepBeforeShuttingProblematicWells_ ==
rhs.minTimeStepBeforeShuttingProblematicWells_;
}
template<class TypeTag>
template<class Controller>
AdaptiveTimeStepping<TypeTag>
AdaptiveTimeStepping<TypeTag>::
serializationTestObject_()
{
AdaptiveTimeStepping<TypeTag> result;
result.restartFactor_ = 1.0;
result.growthFactor_ = 2.0;
result.maxGrowth_ = 3.0;
result.maxTimeStep_ = 4.0;
result.minTimeStep_ = 5.0;
result.ignoreConvergenceFailure_ = true;
result.solverRestartMax_ = 6;
result.solverVerbose_ = true;
result.timestepVerbose_ = true;
result.suggestedNextTimestep_ = 7.0;
result.fullTimestepInitially_ = true;
result.useNewtonIteration_ = true;
result.minTimeStepBeforeShuttingProblematicWells_ = 9.0;
result.timeStepControlType_ = Controller::Type;
result.timeStepControl_ = std::make_unique<Controller>(Controller::serializationTestObject());
return result;
}
template<class TypeTag>
void AdaptiveTimeStepping<TypeTag>::
init_(const UnitSystem& unitSystem)
{
std::tie(timeStepControlType_,
timeStepControl_,
useNewtonIteration_) = detail::createController(unitSystem);
// make sure growth factor is something reasonable
if (growthFactor_ < 1.0) {
OPM_THROW(std::runtime_error,
"Growth factor cannot be less than 1.");
}
}
} // namespace Opm
#endif