From d78cf2e32effe9906a7f9c6c2cfde7e807adf911 Mon Sep 17 00:00:00 2001 From: Erik Hide Saeternes Date: Mon, 10 Feb 2025 14:44:24 +0100 Subject: [PATCH 1/6] Implemented general third order controller for adaptive time stepping --- .../timestepping/AdaptiveTimeStepping.cpp | 8 ++ .../timestepping/AdaptiveTimeStepping.hpp | 3 +- .../AdaptiveTimeStepping_impl.hpp | 21 +++- .../timestepping/TimeStepControl.cpp | 102 ++++++++++++++++-- .../timestepping/TimeStepControl.hpp | 55 +++++++++- .../timestepping/TimeStepControlInterface.hpp | 4 +- 6 files changed, 176 insertions(+), 17 deletions(-) diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping.cpp b/opm/simulators/timestepping/AdaptiveTimeStepping.cpp index 330b0d6e2..df0f3f7f2 100644 --- a/opm/simulators/timestepping/AdaptiveTimeStepping.cpp +++ b/opm/simulators/timestepping/AdaptiveTimeStepping.cpp @@ -212,6 +212,14 @@ createController(const UnitSystem& unitSystem) true }; }}, + {"general3rdorder", + [tol]() { + return RetVal{ + TimeStepControlType::General3rdOrder, + std::make_unique(tol), + false + }; + }}, {"hardcoded", []() { const std::string filename = Parameters::Get(); // "timesteps" diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp index b11e38404..6bb51065a 100644 --- a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp @@ -173,7 +173,7 @@ private: const SimulatorTimer& simulatorTimer_() const; boost::posix_time::ptime startDateTime_() const; double timeStepControlComputeEstimate_( - const double dt, const int iterations, double elapsed) const; + const double dt, const int iterations, AdaptiveSimulatorTimer& substepTimer) const; bool timeStepVerbose_() const; void updateSuggestedNextStep_(); bool useNewtonIteration_() const; @@ -229,6 +229,7 @@ public: static AdaptiveTimeStepping serializationTestObjectPID(); static AdaptiveTimeStepping serializationTestObjectPIDIt(); static AdaptiveTimeStepping serializationTestObjectSimple(); + static AdaptiveTimeStepping serializationTestObject3rdOrder(); private: void maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step, diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp index 02b6294eb..d11d60bbb 100644 --- a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp @@ -24,6 +24,7 @@ #ifndef OPM_ADAPTIVE_TIME_STEPPING_HPP #include #include +#include #endif #include @@ -144,6 +145,9 @@ operator==(const AdaptiveTimeStepping& rhs) case TimeStepControlType::PID: result = castAndComp(rhs); break; + case TimeStepControlType::General3rdOrder: + result = castAndComp(rhs); + break; } return result && @@ -233,6 +237,9 @@ serializeOp(Serializer& serializer) case TimeStepControlType::PID: allocAndSerialize(serializer); break; + case TimeStepControlType::General3rdOrder: + allocAndSerialize(serializer); + break; } serializer(this->restart_factor_); serializer(this->growth_factor_); @@ -282,6 +289,14 @@ serializationTestObjectSimple() return serializationTestObject_(); } +template +AdaptiveTimeStepping +AdaptiveTimeStepping:: +serializationTestObject3rdOrder() +{ + return serializationTestObject_(); +} + template void @@ -771,7 +786,7 @@ run() const int iterations = getNumIterations_(substep_report); auto dt_estimate = timeStepControlComputeEstimate_( - dt, iterations, this->substep_timer_.simulationTimeElapsed()); + dt, iterations, this->substep_timer_); assert(dt_estimate > 0); dt_estimate = maybeRestrictTimeStepGrowth_(dt, dt_estimate, restarts); @@ -1223,12 +1238,12 @@ template template double AdaptiveTimeStepping::SubStepIteration:: -timeStepControlComputeEstimate_(const double dt, const int iterations, double elapsed) const +timeStepControlComputeEstimate_(const double dt, const int iterations, AdaptiveSimulatorTimer& substepTimer) const { // create object to compute the time error, simply forwards the call to the model SolutionTimeErrorSolverWrapper relative_change{solver_()}; return this->adaptive_time_stepping_.time_step_control_->computeTimeStepSize( - dt, iterations, relative_change, elapsed); + dt, iterations, relative_change, substepTimer); } template diff --git a/opm/simulators/timestepping/TimeStepControl.cpp b/opm/simulators/timestepping/TimeStepControl.cpp index 30aa2bc62..8d458e09a 100644 --- a/opm/simulators/timestepping/TimeStepControl.cpp +++ b/opm/simulators/timestepping/TimeStepControl.cpp @@ -33,6 +33,7 @@ #include #include #include +#include #include @@ -73,7 +74,7 @@ namespace Opm } double SimpleIterationCountTimeStepControl:: - computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& /* relativeChange */, const double /*simulationTimeElapsed */) const + computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& /* relativeChange */, const AdaptiveSimulatorTimer& substepTimer) const { double dtEstimate = dt ; @@ -134,10 +135,10 @@ namespace Opm } double HardcodedTimeStepControl:: - computeTimeStepSize( const double /*dt */, const int /*iterations */, const RelativeChangeInterface& /* relativeChange */ , const double simulationTimeElapsed) const + computeTimeStepSize( const double /*dt */, const int /*iterations */, const RelativeChangeInterface& /* relativeChange */ , const AdaptiveSimulatorTimer& substepTimer) const { - auto nextTime = std::upper_bound(subStepTime_.begin(), subStepTime_.end(), simulationTimeElapsed); - return (*nextTime - simulationTimeElapsed); + auto nextTime = std::upper_bound(subStepTime_.begin(), subStepTime_.end(), substepTimer.simulationTimeElapsed()); + return (*nextTime - substepTimer.simulationTimeElapsed()); } bool HardcodedTimeStepControl::operator==(const HardcodedTimeStepControl& ctrl) const @@ -170,7 +171,7 @@ namespace Opm } double PIDTimeStepControl:: - computeTimeStepSize( const double dt, const int /* iterations */, const RelativeChangeInterface& relChange, const double /*simulationTimeElapsed */) const + computeTimeStepSize( const double dt, const int /* iterations */, const RelativeChangeInterface& relChange, const AdaptiveSimulatorTimer& /* substepTimer */) const { // shift errors for( int i=0; i<2; ++i ) { @@ -249,9 +250,9 @@ namespace Opm } double PIDAndIterationCountTimeStepControl:: - computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relChange, const double simulationTimeElapsed ) const + computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relChange, const AdaptiveSimulatorTimer& substepTimer) const { - double dtEstimatePID = PIDTimeStepControl :: computeTimeStepSize( dt, iterations, relChange, simulationTimeElapsed); + double dtEstimatePID = PIDTimeStepControl :: computeTimeStepSize( dt, iterations, relChange, substepTimer); // adjust timesteps based on target iteration double dtEstimateIter; @@ -279,4 +280,91 @@ namespace Opm this->minTimeStepBasedOnIterations_ == ctrl.minTimeStepBasedOnIterations_; } + + + //////////////////////////////////////////////////////////// + // + // General3rdOrderController Implementation + // + //////////////////////////////////////////////////////////// + + General3rdOrderController::General3rdOrderController( const double tolerance, + const double safetyFactor, + const bool verbose) + : tolerance_( tolerance ) + , safetyFactor_( safetyFactor ) + , errors_( 3, tolerance_ ) + , timeSteps_ ( 3, 1.0 ) + , verbose_( verbose ) + {} + + General3rdOrderController + General3rdOrderController::serializationTestObject() + { + General3rdOrderController result(1.0, 0.8, true); + result.errors_ = {2.0, 3.0}; + + return result; + } + + double General3rdOrderController:: + computeTimeStepSize(const double dt, const int /*iterations */, const RelativeChangeInterface& relChange, const AdaptiveSimulatorTimer& substepTimer) const + { + // Shift errors and time steps + for( int i = 0; i < 2; ++i ) + { + errors_[i] = errors_[i+1]; + timeSteps_[i] = timeSteps_[i+1]; + } + + // Store new error and time step + const double error = relChange.relativeChange(); + errors_[2] = error; + timeSteps_[2] = dt; + for( int i = 0; i < 2; ++i ) + { + assert(std::isfinite(errors_[i])); + } + + if (errors_[0] == 0 || errors_[1] == 0 || errors_[2] == 0.) + { + if ( verbose_ ) + OpmLog::info("The solution between time steps does not change, there is no time step constraint from the controller."); + return std::numeric_limits::max(); + } + // Use an I controller after report time steps or chopped time steps + else if (substepTimer.currentStepNum() < 3 || substepTimer.lastStepFailed() || counterSinceFailure_ > 0) + { + if (substepTimer.lastStepFailed() || counterSinceFailure_ > 0) + counterSinceFailure_++; + if (counterSinceFailure_ > 1) + counterSinceFailure_ = 0; + const double newDt = dt * std::pow(safetyFactor_ * tolerance_ / errors_[2], 0.35); + return newDt; + } + // Use the general third order controller for all other time steps + else + { + const std::vector beta = { 0.125, 0.25, 0.125 }; + const std::vector alpha = { 0.375, 0.125 }; + const double newDt = dt * std::pow(safetyFactor_ * tolerance_ / errors_[2], beta[0]) * + std::pow(safetyFactor_ * tolerance_ / errors_[1], beta[1]) * + std::pow(safetyFactor_ * tolerance_ / errors_[0], beta[2]) * + std::pow(timeSteps_[2] / timeSteps_[1], -alpha[0]) * + std::pow(timeSteps_[1] / timeSteps_[0], -alpha[1]); + if( verbose_ ) + OpmLog::info(fmt::format("Computed step size (pow): {} days", unit::convert::to( newDt, unit::day ))); + return newDt; + } + } + + bool General3rdOrderController::operator==(const General3rdOrderController& ctrl) const + { + return this->tolerance_ == ctrl.tolerance_ && + this->safetyFactor_ == ctrl.safetyFactor_ && + this->errors_ == ctrl.errors_ && + this->timeSteps_ == ctrl.timeSteps_ && + this->verbose_ == ctrl.verbose_; + } + } // end namespace Opm diff --git a/opm/simulators/timestepping/TimeStepControl.hpp b/opm/simulators/timestepping/TimeStepControl.hpp index b6c465954..14109694d 100644 --- a/opm/simulators/timestepping/TimeStepControl.hpp +++ b/opm/simulators/timestepping/TimeStepControl.hpp @@ -22,6 +22,7 @@ #define OPM_TIMESTEPCONTROL_HEADER_INCLUDED #include +#include #include #include @@ -32,7 +33,8 @@ namespace Opm SimpleIterationCount, PID, PIDAndIterationCount, - HardCodedTimeStep + HardCodedTimeStep, + General3rdOrder }; /////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -62,7 +64,7 @@ namespace Opm double computeTimeStepSize(const double dt, const int iterations, const RelativeChangeInterface& /* relativeChange */, - const double /*simulationTimeElapsed */ ) const override; + const AdaptiveSimulatorTimer& /* substepTimer */ ) const override; template void serializeOp(Serializer& serializer) @@ -113,7 +115,7 @@ namespace Opm double computeTimeStepSize(const double dt, const int /* iterations */, const RelativeChangeInterface& relativeChange, - const double /*simulationTimeElapsed */ ) const override; + const AdaptiveSimulatorTimer& /* substepTimer */ ) const override; template void serializeOp(Serializer& serializer) @@ -162,7 +164,7 @@ namespace Opm double computeTimeStepSize(const double dt, const int iterations, const RelativeChangeInterface& relativeChange, - const double /*simulationTimeElapsed */ ) const override; + const AdaptiveSimulatorTimer& /* substepTimer */ ) const override; template void serializeOp(Serializer& serializer) @@ -183,6 +185,49 @@ namespace Opm const double minTimeStepBasedOnIterations_; }; + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// + /// General 3rd order controller + /// + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + class General3rdOrderController : public TimeStepControlInterface + { + public: + static constexpr TimeStepControlType Type = TimeStepControlType::General3rdOrder; + + General3rdOrderController( const double tolerance = 1e-3, + const double safetyFactor = 0.8, + const bool verbose = false ); + + static General3rdOrderController serializationTestObject(); + + double computeTimeStepSize(const double dt, + const int /*iterations */, + const RelativeChangeInterface& relativeChange, + const AdaptiveSimulatorTimer& substepTimer) const override; + + template + void serializeOp(Serializer& serializer) + { + serializer(tolerance_); + serializer(safetyFactor_); + serializer(errors_); + serializer(timeSteps_); + serializer(verbose_); + } + + bool operator==(const General3rdOrderController&) const; + + + protected: + const double tolerance_ = 1e-3; + const double safetyFactor_ = 0.8; + mutable std::vector errors_{}; + mutable std::vector timeSteps_{}; + mutable int counterSinceFailure_ = 0; + const bool verbose_ = false; + }; + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// /// HardcodedTimeStepControl @@ -208,7 +253,7 @@ namespace Opm double computeTimeStepSize(const double dt, const int /* iterations */, const RelativeChangeInterface& /*relativeChange */, - const double simulationTimeElapsed) const override; + const AdaptiveSimulatorTimer& substepTimer) const override; template void serializeOp(Serializer& serializer) diff --git a/opm/simulators/timestepping/TimeStepControlInterface.hpp b/opm/simulators/timestepping/TimeStepControlInterface.hpp index bdab5e3a4..3a04e22ac 100644 --- a/opm/simulators/timestepping/TimeStepControlInterface.hpp +++ b/opm/simulators/timestepping/TimeStepControlInterface.hpp @@ -19,6 +19,8 @@ #ifndef OPM_TIMESTEPCONTROLINTERFACE_HEADER_INCLUDED #define OPM_TIMESTEPCONTROLINTERFACE_HEADER_INCLUDED +#include + namespace Opm { @@ -56,7 +58,7 @@ namespace Opm /// \param timeError object to compute || u^n+1 - u^n || / || u^n+1 || /// /// \return suggested time step size for the next step - virtual double computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relativeChange , const double simulationTimeElapsed) const = 0; + virtual double computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relativeChange , const AdaptiveSimulatorTimer& substepTimer) const = 0; /// virtual destructor (empty) virtual ~TimeStepControlInterface () {} From 2be953f73d3e4b06bea6f58dfc81e0a8f410008c Mon Sep 17 00:00:00 2001 From: Erik Hide Saeternes Date: Tue, 11 Feb 2025 18:08:29 +0100 Subject: [PATCH 2/6] Added missing safety factor paramet for time step controller --- opm/simulators/timestepping/AdaptiveTimeStepping.cpp | 7 ++++++- opm/simulators/timestepping/AdaptiveTimeStepping.hpp | 1 + 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping.cpp b/opm/simulators/timestepping/AdaptiveTimeStepping.cpp index df0f3f7f2..12f9e98ab 100644 --- a/opm/simulators/timestepping/AdaptiveTimeStepping.cpp +++ b/opm/simulators/timestepping/AdaptiveTimeStepping.cpp @@ -136,6 +136,9 @@ void registerAdaptiveParameters() Parameters::Register ("The minimum time step size (in days for field and metric unit and hours for lab unit) " "can be reduced to based on newton iteration counts"); + Parameters::Register + ("Value to be multiplied with the time step control tolerance to ensure that the target " + "relative change is lower than the tolerance"); } std::tuple, bool> @@ -214,9 +217,11 @@ createController(const UnitSystem& unitSystem) }}, {"general3rdorder", [tol]() { + const double safetyFactor = Parameters::Get(); // 0.8 return RetVal{ TimeStepControlType::General3rdOrder, - std::make_unique(tol), + std::make_unique(tol, + safetyFactor), false }; }}, diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp index 6bb51065a..7eb6d8663 100644 --- a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp @@ -57,6 +57,7 @@ struct TimeStepControlGrowthDampingFactor { static constexpr double value = 3.2; struct TimeStepControlFileName { static constexpr auto value = "timesteps"; }; struct MinTimeStepBeforeShuttingProblematicWellsInDays { static constexpr double value = 0.01; }; struct MinTimeStepBasedOnNewtonIterations { static constexpr double value = 0.0; }; +struct TimeStepControlSafetyFactor { static constexpr double value = 0.8; }; } // namespace Opm::Parameters From 7ac57ae177096c144dada46716ff65e0fc86493d Mon Sep 17 00:00:00 2001 From: Erik Hide Saeternes Date: Thu, 13 Feb 2025 10:18:51 +0100 Subject: [PATCH 3/6] Added test --- tests/test_RestartSerialization.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_RestartSerialization.cpp b/tests/test_RestartSerialization.cpp index f1c48eead..e08671f6f 100644 --- a/tests/test_RestartSerialization.cpp +++ b/tests/test_RestartSerialization.cpp @@ -123,6 +123,7 @@ TEST_FOR_TYPE_NAMED_OBJ(ATS, AdaptiveTimeSteppingHardcoded, serializationTestObj TEST_FOR_TYPE_NAMED_OBJ(ATS, AdaptiveTimeSteppingPID, serializationTestObjectPID) TEST_FOR_TYPE_NAMED_OBJ(ATS, AdaptiveTimeSteppingPIDIt, serializationTestObjectPIDIt) TEST_FOR_TYPE_NAMED_OBJ(ATS, AdaptiveTimeSteppingSimple, serializationTestObjectSimple) +TEST_FOR_TYPE_NAMED_OBJ(ATS, AdaptiveTimeStepping3rdOrder, serializationTestObject3rdOrder) namespace Opm { using BPV = BlackOilPrimaryVariables; } TEST_FOR_TYPE_NAMED(BPV, BlackoilPrimaryVariables) From fb05bbf881a65e70eb21a907591c79a1850f8a88 Mon Sep 17 00:00:00 2001 From: Erik Hide Saeternes Date: Fri, 14 Feb 2025 13:56:59 +0100 Subject: [PATCH 4/6] Fixed a few details in the time step control --- opm/simulators/timestepping/TimeStepControl.cpp | 8 +++++--- opm/simulators/timestepping/TimeStepControl.hpp | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/opm/simulators/timestepping/TimeStepControl.cpp b/opm/simulators/timestepping/TimeStepControl.cpp index 8d458e09a..940d7ab7e 100644 --- a/opm/simulators/timestepping/TimeStepControl.cpp +++ b/opm/simulators/timestepping/TimeStepControl.cpp @@ -74,7 +74,7 @@ namespace Opm } double SimpleIterationCountTimeStepControl:: - computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& /* relativeChange */, const AdaptiveSimulatorTimer& substepTimer) const + computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& /* relativeChange */, const AdaptiveSimulatorTimer& /* substepTimer */) const { double dtEstimate = dt ; @@ -340,13 +340,15 @@ namespace Opm if (counterSinceFailure_ > 1) counterSinceFailure_ = 0; const double newDt = dt * std::pow(safetyFactor_ * tolerance_ / errors_[2], 0.35); + if( verbose_ ) + OpmLog::info(fmt::format("Computed step size (pow): {} days", unit::convert::to( newDt, unit::day ))); return newDt; } // Use the general third order controller for all other time steps else { - const std::vector beta = { 0.125, 0.25, 0.125 }; - const std::vector alpha = { 0.375, 0.125 }; + const std::array beta = { 0.125, 0.25, 0.125 }; + const std::array alpha = { 0.375, 0.125 }; const double newDt = dt * std::pow(safetyFactor_ * tolerance_ / errors_[2], beta[0]) * std::pow(safetyFactor_ * tolerance_ / errors_[1], beta[1]) * std::pow(safetyFactor_ * tolerance_ / errors_[0], beta[2]) * diff --git a/opm/simulators/timestepping/TimeStepControl.hpp b/opm/simulators/timestepping/TimeStepControl.hpp index 14109694d..2c59865c6 100644 --- a/opm/simulators/timestepping/TimeStepControl.hpp +++ b/opm/simulators/timestepping/TimeStepControl.hpp @@ -34,7 +34,7 @@ namespace Opm PID, PIDAndIterationCount, HardCodedTimeStep, - General3rdOrder + General3rdOrder, }; /////////////////////////////////////////////////////////////////////////////////////////////////////////////// From 229c8f77339b34d920ba378760a19add544dd1ae Mon Sep 17 00:00:00 2001 From: Erik Hide Saeternes Date: Fri, 14 Feb 2025 14:16:41 +0100 Subject: [PATCH 5/6] Fixed bug in std::array --- opm/simulators/timestepping/TimeStepControl.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/simulators/timestepping/TimeStepControl.cpp b/opm/simulators/timestepping/TimeStepControl.cpp index 940d7ab7e..3ce975262 100644 --- a/opm/simulators/timestepping/TimeStepControl.cpp +++ b/opm/simulators/timestepping/TimeStepControl.cpp @@ -347,8 +347,8 @@ namespace Opm // Use the general third order controller for all other time steps else { - const std::array beta = { 0.125, 0.25, 0.125 }; - const std::array alpha = { 0.375, 0.125 }; + const std::array beta = { 0.125, 0.25, 0.125 }; + const std::array alpha = { 0.375, 0.125 }; const double newDt = dt * std::pow(safetyFactor_ * tolerance_ / errors_[2], beta[0]) * std::pow(safetyFactor_ * tolerance_ / errors_[1], beta[1]) * std::pow(safetyFactor_ * tolerance_ / errors_[0], beta[2]) * From 565f5cffd8d171769f20ade78404f42c103b30ed Mon Sep 17 00:00:00 2001 From: Erik Hide Saeternes Date: Fri, 14 Feb 2025 16:29:19 +0100 Subject: [PATCH 6/6] Added #include --- opm/simulators/timestepping/TimeStepControl.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/simulators/timestepping/TimeStepControl.cpp b/opm/simulators/timestepping/TimeStepControl.cpp index 3ce975262..3d1513a3c 100644 --- a/opm/simulators/timestepping/TimeStepControl.cpp +++ b/opm/simulators/timestepping/TimeStepControl.cpp @@ -20,6 +20,7 @@ */ #include #include +#include #include #include #include