From a2bd8b58107f0635f4017f14c489ffbf155db29a Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 2 Jan 2025 14:18:56 +0100 Subject: [PATCH] split AdaptiveTimeStepping to use an impl file --- CMakeLists_files.cmake | 1 + .../timestepping/AdaptiveTimeStepping.hpp | 586 +---------------- .../AdaptiveTimeStepping_impl.hpp | 607 ++++++++++++++++++ 3 files changed, 643 insertions(+), 551 deletions(-) create mode 100644 opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 4588c2dd9..571550bca 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -928,6 +928,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/WriteSystemMatrixHelper.hpp opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp opm/simulators/timestepping/AdaptiveTimeStepping.hpp + opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp opm/simulators/timestepping/ConvergenceReport.hpp opm/simulators/timestepping/EclTimeSteppingParams.hpp opm/simulators/timestepping/TimeStepControl.hpp diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp index 8f75687cf..636f85b1f 100644 --- a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp @@ -6,19 +6,11 @@ #include #include -#include -#include #include -#include - -#include -#include - #include #include -#include #include #include @@ -28,18 +20,10 @@ #include #include -#include - -#include - -#include -#include #include #include #include #include -#include -#include #include #include @@ -67,6 +51,7 @@ struct MinTimeStepBasedOnNewtonIterations { static constexpr double value = 0.0; namespace Opm { +class UnitSystem; struct StepReport; namespace detail { @@ -114,30 +99,9 @@ void registerAdaptiveParameters(); AdaptiveTimeStepping() = default; //! \brief contructor taking parameter object - AdaptiveTimeStepping(const UnitSystem& unitSystem, - const double max_next_tstep = -1.0, - const bool terminalOutput = true) - : timeStepControl_() - , restartFactor_(Parameters::Get>()) // 0.33 - , growthFactor_(Parameters::Get>()) // 2.0 - , maxGrowth_(Parameters::Get>()) // 3.0 - , maxTimeStep_(Parameters::Get>() * 24 * 60 * 60) // 365.25 - , minTimeStep_(unitSystem.to_si(UnitSystem::measure::time, Parameters::Get>())) // 1e-12; - , ignoreConvergenceFailure_(Parameters::Get()) // false; - , solverRestartMax_(Parameters::Get()) // 10 - , solverVerbose_(Parameters::Get() > 0 && terminalOutput) // 2 - , timestepVerbose_(Parameters::Get() > 0 && terminalOutput) // 2 - , suggestedNextTimestep_((max_next_tstep <= 0 ? Parameters::Get() : max_next_tstep) * 24 * 60 * 60) // 1.0 - , fullTimestepInitially_(Parameters::Get()) // false - , timestepAfterEvent_(Parameters::Get>() * 24 * 60 * 60) // 1e30 - , useNewtonIteration_(false) - , minTimeStepBeforeShuttingProblematicWells_(Parameters::Get() * unit::day) - - { - init_(unitSystem); - } - - + explicit AdaptiveTimeStepping(const UnitSystem& unitSystem, + const double max_next_tstep = -1.0, + const bool terminalOutput = true); //! \brief contructor taking parameter object //! \param tuning Pointer to ecl TUNING keyword @@ -145,31 +109,9 @@ void registerAdaptiveParameters(); AdaptiveTimeStepping(double max_next_tstep, const Tuning& tuning, const UnitSystem& unitSystem, - 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()) // 10 - , solverVerbose_(Parameters::Get() > 0 && terminalOutput) // 2 - , timestepVerbose_(Parameters::Get() > 0 && terminalOutput) // 2 - , suggestedNextTimestep_(max_next_tstep <= 0 ? Parameters::Get() * 24 * 60 * 60 : max_next_tstep) // 1.0 - , fullTimestepInitially_(Parameters::Get()) // false - , timestepAfterEvent_(tuning.TMAXWC) // 1e30 - , useNewtonIteration_(false) - , minTimeStepBeforeShuttingProblematicWells_(Parameters::Get() * unit::day) - { - init_(unitSystem); - } + const bool terminalOutput = true); - static void registerParameters() - { - registerEclTimeSteppingParameters(); - detail::registerAdaptiveParameters(); - } + static void registerParameters(); /** \brief step method that acts like the solver::step method in a sub cycle of time steps @@ -180,294 +122,7 @@ void registerAdaptiveParameters(); SimulatorReport step(const SimulatorTimer& simulatorTimer, Solver& solver, const bool isEvent, - const std::function 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 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 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 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; - } + const std::function tuningUpdater); /** \brief Returns the simulator report for the failed substeps of the last * report step. @@ -478,142 +133,24 @@ void registerAdaptiveParameters(); void setSuggestedNextStep(const double x) { suggestedNextTimestep_ = x; } - 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 updateTUNING(double max_next_tstep, const Tuning& tuning); - 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; - } - } + void updateNEXTSTEP(double max_next_tstep); template - void serializeOp(Serializer& serializer) - { - serializer(timeStepControlType_); - switch (timeStepControlType_) { - case TimeStepControlType::HardCodedTimeStep: - allocAndSerialize(serializer); - break; - case TimeStepControlType::PIDAndIterationCount: - allocAndSerialize(serializer); - break; - case TimeStepControlType::SimpleIterationCount: - allocAndSerialize(serializer); - break; - case TimeStepControlType::PID: - allocAndSerialize(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_); - } + void serializeOp(Serializer& serializer); - static AdaptiveTimeStepping serializationTestObjectHardcoded() - { - return serializationTestObject_(); - } + static AdaptiveTimeStepping serializationTestObjectHardcoded(); + static AdaptiveTimeStepping serializationTestObjectPID(); + static AdaptiveTimeStepping serializationTestObjectPIDIt(); + static AdaptiveTimeStepping serializationTestObjectSimple(); - static AdaptiveTimeStepping serializationTestObjectPID() - { - return serializationTestObject_(); - } - - static AdaptiveTimeStepping serializationTestObjectPIDIt() - { - return serializationTestObject_(); - } - - static AdaptiveTimeStepping serializationTestObjectSimple() - { - return serializationTestObject_(); - } - - bool operator==(const AdaptiveTimeStepping& rhs) const - { - if (timeStepControlType_ != rhs.timeStepControlType_ || - (timeStepControl_ && !rhs.timeStepControl_) || - (!timeStepControl_ && rhs.timeStepControl_)) { - return false; - } - - bool result = false; - switch (timeStepControlType_) { - case TimeStepControlType::HardCodedTimeStep: - result = castAndComp(rhs); - break; - case TimeStepControlType::PIDAndIterationCount: - result = castAndComp(rhs); - break; - case TimeStepControlType::SimpleIterationCount: - result = castAndComp(rhs); - break; - case TimeStepControlType::PID: - result = castAndComp(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_; - } + bool operator==(const AdaptiveTimeStepping& rhs) const; private: template - static AdaptiveTimeStepping serializationTestObject_() - { - AdaptiveTimeStepping result; + static AdaptiveTimeStepping serializationTestObject_(); - 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::serializationTestObject()); - - return result; - } template void allocAndSerialize(Serializer& serializer) { @@ -632,82 +169,29 @@ void registerAdaptiveParameters(); } protected: - void init_(const UnitSystem& unitSystem) - { - // valid are "pid" and "pid+iteration" - std::string control = Parameters::Get(); // "pid" - - const double tol = Parameters::Get(); // 1e-1 - if (control == "pid") { - timeStepControl_ = std::make_unique(tol); - timeStepControlType_ = TimeStepControlType::PID; - } - else if (control == "pid+iteration") { - const int iterations = Parameters::Get(); // 30 - const double decayDampingFactor = Parameters::Get(); // 1.0 - const double growthDampingFactor = Parameters::Get(); // 3.2 - timeStepControl_ = std::make_unique(iterations, decayDampingFactor, growthDampingFactor, tol); - timeStepControlType_ = TimeStepControlType::PIDAndIterationCount; - } - else if (control == "pid+newtoniteration") { - const int iterations = Parameters::Get(); // 8 - const double decayDampingFactor = Parameters::Get(); // 1.0 - const double growthDampingFactor = Parameters::Get(); // 3.2 - const double nonDimensionalMinTimeStepIterations = Parameters::Get(); // 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(iterations, decayDampingFactor, - growthDampingFactor, tol, minTimeStepReducedByIterations); - timeStepControlType_ = TimeStepControlType::PIDAndIterationCount; - useNewtonIteration_ = true; - } - else if (control == "iterationcount") { - const int iterations = Parameters::Get(); // 30 - const double decayrate = Parameters::Get(); // 0.75 - const double growthrate = Parameters::Get(); // 1.25 - timeStepControl_ = std::make_unique(iterations, decayrate, growthrate); - timeStepControlType_ = TimeStepControlType::SimpleIterationCount; - } - else if (control == "newtoniterationcount") { - const int iterations = Parameters::Get(); // 8 - const double decayrate = Parameters::Get(); // 0.75 - const double growthrate = Parameters::Get(); // 1.25 - timeStepControl_ = std::make_unique(iterations, decayrate, growthrate); - useNewtonIteration_ = true; - timeStepControlType_ = TimeStepControlType::SimpleIterationCount; - } - else if (control == "hardcoded") { - const std::string filename = Parameters::Get(); // "timesteps" - timeStepControl_ = std::make_unique(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); - } + void init_(const UnitSystem& unitSystem); using TimeStepController = std::unique_ptr; - TimeStepControlType timeStepControlType_; //!< type of time step control object - TimeStepController timeStepControl_; //!< time step control object - 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 maxGrowth_; //!< factor that limits the maximum growth of a time step - double maxTimeStep_; //!< maximal allowed time step size in days - double minTimeStep_; //!< minimal allowed time step size before throwing - bool ignoreConvergenceFailure_; //!< continue instead of stop when minimum time step is reached - int solverRestartMax_; //!< how many restart of solver are allowed - bool solverVerbose_; //!< solver verbosity - bool timestepVerbose_; //!< timestep verbosity - double suggestedNextTimestep_; //!< suggested size of next timestep - bool fullTimestepInitially_; //!< beginning with the size of the time step from data file - double timestepAfterEvent_; //!< suggested size of timestep after an event - bool useNewtonIteration_; //!< use newton iteration count for adaptive time step control - double minTimeStepBeforeShuttingProblematicWells_; //! < shut problematic wells when time step size in days are less than this + TimeStepControlType timeStepControlType_{TimeStepControlType::PIDAndIterationCount}; //!< type of time step control object + TimeStepController timeStepControl_{}; //!< time step control object + 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 maxGrowth_{}; //!< factor that limits the maximum growth of a time step + double maxTimeStep_{}; //!< maximal allowed time step size in days + double minTimeStep_{}; //!< minimal allowed time step size before throwing + bool ignoreConvergenceFailure_{false}; //!< continue instead of stop when minimum time step is reached + int solverRestartMax_{}; //!< how many restart of solver are allowed + bool solverVerbose_{false}; //!< solver verbosity + bool timestepVerbose_{false}; //!< timestep verbosity + double suggestedNextTimestep_{}; //!< suggested size of next timestep + bool fullTimestepInitially_{false}; //!< beginning with the size of the time step from data file + double timestepAfterEvent_{}; //!< suggested size of timestep after an event + 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 }; } +#include + #endif // OPM_ADAPTIVE_TIME_STEPPING_HPP diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp new file mode 100644 index 000000000..2fbbc58ed --- /dev/null +++ b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp @@ -0,0 +1,607 @@ +/* +*/ +#ifndef OPM_ADAPTIVE_TIME_STEPPING_IMPL_HPP +#define OPM_ADAPTIVE_TIME_STEPPING_IMPL_HPP + +// Improve IDE experience +#ifndef OPM_ADAPTIVE_TIME_STEPPING_HPP +#include +#include +#endif + +#include +#include + +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include + +namespace Opm { + +template +AdaptiveTimeStepping:: +AdaptiveTimeStepping(const UnitSystem& unitSystem, + const double max_next_tstep, + const bool terminalOutput) + : timeStepControl_() + , restartFactor_(Parameters::Get>()) // 0.33 + , growthFactor_(Parameters::Get>()) // 2.0 + , maxGrowth_(Parameters::Get>()) // 3.0 + , maxTimeStep_(Parameters::Get>() * 24 * 60 * 60) // 365.25 + , minTimeStep_(unitSystem.to_si(UnitSystem::measure::time, Parameters::Get>())) // 1e-12; + , ignoreConvergenceFailure_(Parameters::Get()) // false; + , solverRestartMax_(Parameters::Get()) // 10 + , solverVerbose_(Parameters::Get() > 0 && terminalOutput) // 2 + , timestepVerbose_(Parameters::Get() > 0 && terminalOutput) // 2 + , suggestedNextTimestep_((max_next_tstep <= 0 ? Parameters::Get() : max_next_tstep) * 24 * 60 * 60) // 1.0 + , fullTimestepInitially_(Parameters::Get()) // false + , timestepAfterEvent_(Parameters::Get>() * 24 * 60 * 60) // 1e30 + , useNewtonIteration_(false) + , minTimeStepBeforeShuttingProblematicWells_(Parameters::Get() * unit::day) + +{ + init_(unitSystem); +} + +template +AdaptiveTimeStepping:: +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()) // 10 + , solverVerbose_(Parameters::Get() > 0 && terminalOutput) // 2 + , timestepVerbose_(Parameters::Get() > 0 && terminalOutput) // 2 + , suggestedNextTimestep_(max_next_tstep <= 0 ? Parameters::Get() * 24 * 60 * 60 : max_next_tstep) // 1.0 + , fullTimestepInitially_(Parameters::Get()) // false + , timestepAfterEvent_(tuning.TMAXWC) // 1e30 + , useNewtonIteration_(false) + , minTimeStepBeforeShuttingProblematicWells_(Parameters::Get() * unit::day) +{ + init_(unitSystem); +} + +template +void AdaptiveTimeStepping::registerParameters() +{ + registerEclTimeSteppingParameters(); + detail::registerAdaptiveParameters(); +} + +template +template +SimulatorReport +AdaptiveTimeStepping:: +step(const SimulatorTimer& simulatorTimer, + Solver& solver, + const bool isEvent, + const std::function 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 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 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 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; +} + +template +void AdaptiveTimeStepping:: +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 +void AdaptiveTimeStepping:: +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 +template +void AdaptiveTimeStepping:: +serializeOp(Serializer& serializer) +{ + serializer(timeStepControlType_); + switch (timeStepControlType_) { + case TimeStepControlType::HardCodedTimeStep: + allocAndSerialize(serializer); + break; + case TimeStepControlType::PIDAndIterationCount: + allocAndSerialize(serializer); + break; + case TimeStepControlType::SimpleIterationCount: + allocAndSerialize(serializer); + break; + case TimeStepControlType::PID: + allocAndSerialize(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 +AdaptiveTimeStepping +AdaptiveTimeStepping:: +serializationTestObjectHardcoded() +{ + return serializationTestObject_(); +} + +template +AdaptiveTimeStepping +AdaptiveTimeStepping:: +serializationTestObjectPID() +{ + return serializationTestObject_(); +} + +template +AdaptiveTimeStepping +AdaptiveTimeStepping:: +serializationTestObjectPIDIt() +{ + return serializationTestObject_(); +} + +template +AdaptiveTimeStepping +AdaptiveTimeStepping:: +serializationTestObjectSimple() +{ + return serializationTestObject_(); +} + +template +bool +AdaptiveTimeStepping:: +operator==(const AdaptiveTimeStepping& rhs) const +{ + if (timeStepControlType_ != rhs.timeStepControlType_ || + (timeStepControl_ && !rhs.timeStepControl_) || + (!timeStepControl_ && rhs.timeStepControl_)) { + return false; + } + + bool result = false; + switch (timeStepControlType_) { + case TimeStepControlType::HardCodedTimeStep: + result = castAndComp(rhs); + break; + case TimeStepControlType::PIDAndIterationCount: + result = castAndComp(rhs); + break; + case TimeStepControlType::SimpleIterationCount: + result = castAndComp(rhs); + break; + case TimeStepControlType::PID: + result = castAndComp(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 +template +AdaptiveTimeStepping +AdaptiveTimeStepping:: +serializationTestObject_() +{ + AdaptiveTimeStepping 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::serializationTestObject()); + + return result; +} + +template +void AdaptiveTimeStepping:: +init_(const UnitSystem& unitSystem) +{ + // valid are "pid" and "pid+iteration" + std::string control = Parameters::Get(); // "pid" + + const double tol = Parameters::Get(); // 1e-1 + if (control == "pid") { + timeStepControl_ = std::make_unique(tol); + timeStepControlType_ = TimeStepControlType::PID; + } + else if (control == "pid+iteration") { + const int iterations = Parameters::Get(); // 30 + const double decayDampingFactor = Parameters::Get(); // 1.0 + const double growthDampingFactor = Parameters::Get(); // 3.2 + timeStepControl_ = std::make_unique(iterations, decayDampingFactor, growthDampingFactor, tol); + timeStepControlType_ = TimeStepControlType::PIDAndIterationCount; + } + else if (control == "pid+newtoniteration") { + const int iterations = Parameters::Get(); // 8 + const double decayDampingFactor = Parameters::Get(); // 1.0 + const double growthDampingFactor = Parameters::Get(); // 3.2 + const double nonDimensionalMinTimeStepIterations = Parameters::Get(); // 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(iterations, decayDampingFactor, + growthDampingFactor, tol, minTimeStepReducedByIterations); + timeStepControlType_ = TimeStepControlType::PIDAndIterationCount; + useNewtonIteration_ = true; + } + else if (control == "iterationcount") { + const int iterations = Parameters::Get(); // 30 + const double decayrate = Parameters::Get(); // 0.75 + const double growthrate = Parameters::Get(); // 1.25 + timeStepControl_ = std::make_unique(iterations, decayrate, growthrate); + timeStepControlType_ = TimeStepControlType::SimpleIterationCount; + } + else if (control == "newtoniterationcount") { + const int iterations = Parameters::Get(); // 8 + const double decayrate = Parameters::Get(); // 0.75 + const double growthrate = Parameters::Get(); // 1.25 + timeStepControl_ = std::make_unique(iterations, decayrate, growthrate); + useNewtonIteration_ = true; + timeStepControlType_ = TimeStepControlType::SimpleIterationCount; + } + else if (control == "hardcoded") { + const std::string filename = Parameters::Get(); // "timesteps" + timeStepControl_ = std::make_unique(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); +} + +} // namespace Opm + +#endif