diff --git a/opm/simulators/flow/FlowMain.hpp b/opm/simulators/flow/FlowMain.hpp index db48ee6b3..77a1f140e 100644 --- a/opm/simulators/flow/FlowMain.hpp +++ b/opm/simulators/flow/FlowMain.hpp @@ -51,6 +51,7 @@ namespace Opm::Parameters { // Do not merge parallel output files or warn about them struct EnableLoggingFalloutWarning { static constexpr bool value = false; }; struct OutputInterval { static constexpr int value = 1; }; + } // namespace Opm::Parameters namespace Opm { diff --git a/opm/simulators/flow/ReservoirCoupling.hpp b/opm/simulators/flow/ReservoirCoupling.hpp index 9b3f777fa..75cffd73e 100644 --- a/opm/simulators/flow/ReservoirCoupling.hpp +++ b/opm/simulators/flow/ReservoirCoupling.hpp @@ -20,15 +20,18 @@ #ifndef OPM_RESERVOIR_COUPLING_HPP #define OPM_RESERVOIR_COUPLING_HPP #include +#include +#include #include namespace Opm { namespace ReservoirCoupling { enum class MessageTag : int { - SimulationStartDate, - SimulationEndDate, - SlaveProcessTermination + SlaveSimulationStartDate, + SlaveProcessTermination, + SlaveNextReportDate, + SlaveNextTimeStep, }; // Custom deleter for MPI_Comm @@ -43,6 +46,154 @@ struct MPI_Comm_Deleter { using MPI_Comm_Ptr = std::unique_ptr; +// This class represents a time point. +// It is currently used to represent an epoch time (a double value in seconds since the epoch), +// or an elapsed time (a double value in seconds since the start of the simulation). +// To avoid numerical issues when adding or subtracting time points and then later comparing +// for equality with for example a given report date, we use a tolerance value. +class TimePoint { +private: + double time; + // TODO: Epoch values often lies in the range of [1e9,1e11], so a tolerance value of 1e-10 + // might be a little too small. However, for elapsed time values, the range is often + // in the range of [0, 1e8], so a tolerance value of 1e-10 should be sufficient. + // NOTE: 1 nano-second = 1e-9 seconds + static constexpr double tol = 1e-10; // Tolerance value + +public: + TimePoint() : time(0.0) {} + explicit TimePoint(double t) : time(t) {} + TimePoint(const TimePoint& other) : time(other.time) {} + + // Assignment operator for double + TimePoint& operator=(double t) { + time = t; + return *this; + } + + // Copy assignment operator + TimePoint& operator=(const TimePoint& other) { + if (this != &other) { + time = other.time; + } + return *this; + } + + double getTime() const { return time; } + + // Equality operator + bool operator==(const TimePoint& other) const { + return std::abs(time - other.time) < tol; + } + + // Inequality operator + bool operator!=(const TimePoint& other) const { + return !(*this == other); + } + + // Less than operator + bool operator<(const TimePoint& other) const { + return (time < other.time) && !(*this == other); + } + + // Comparison operator: double < TimePoint + friend bool operator<(double lhs, const TimePoint& rhs) { + return lhs < rhs.time; + } + + // Comparison operator: TimePoint < double + bool operator<(double rhs) const { + return time < rhs; + } + + // Less than or equal to operator + bool operator<=(const TimePoint& other) const { + return (time < other.time) || (*this == other); + } + + // Comparison operator: double <= TimePoint + friend bool operator<=(double lhs, const TimePoint& rhs) { + return lhs <= rhs.time; + } + + // Comparison operator: TimePoint <= double + bool operator<=(double rhs) const { + return time <= rhs; + } + + // Greater than operator + bool operator>(const TimePoint& other) const { + return (time > other.time) && !(*this == other); + } + + // Comparison operator: double > TimePoint + friend bool operator>(double lhs, const TimePoint& rhs) { + return lhs > rhs.time; + } + + // Comparison operator: TimePoint > double + bool operator>(double rhs) const { + return time > rhs; + } + + // Greater than or equal to operator + bool operator>=(const TimePoint& other) const { + return (time > other.time) || (*this == other); + } + + // Comparison operator: TimePoint >= double + bool operator>=(double rhs) const { + return time >= rhs; + } + + // Comparison operator: double >= TimePoint + friend bool operator>=(double lhs, const TimePoint& rhs) { + return lhs >= rhs.time; + } + + // Addition operator: TimePoint + TimePoint (summing their times) + TimePoint operator+(const TimePoint& other) const { + return TimePoint(time + other.time); + } + + // Addition operator: TimePoint + double + TimePoint operator+(double delta) const { + return TimePoint(time + delta); + } + + // Friend addition operator: double + TimePoint + friend TimePoint operator+(double lhs, const TimePoint& rhs) { + return TimePoint(lhs + rhs.time); + } + + // Overload += operator for adding a double + TimePoint& operator+=(double delta) { + time += delta; + return *this; + } + + // Subtraction operator: TimePoint - TimePoint (resulting in a new TimePoint) + TimePoint operator-(const TimePoint& other) const { + return TimePoint(time - other.time); + } + + // Subtraction operator: TimePoint - double + TimePoint operator-(double delta) const { + return TimePoint(time - delta); + } + + // Friend subtraction operator: double - TimePoint + friend TimePoint operator-(double lhs, const TimePoint& rhs) { + return TimePoint(lhs - rhs.time); + } + + // Stream insertion operator for easy printing + friend std::ostream& operator<<(std::ostream& os, const TimePoint& tp) { + os << tp.time; + return os; + } +}; + } // namespace ReservoirCoupling } // namespace Opm diff --git a/opm/simulators/flow/ReservoirCouplingMaster.cpp b/opm/simulators/flow/ReservoirCouplingMaster.cpp index 3113dd039..57bf674fc 100644 --- a/opm/simulators/flow/ReservoirCouplingMaster.cpp +++ b/opm/simulators/flow/ReservoirCouplingMaster.cpp @@ -44,41 +44,133 @@ ReservoirCouplingMaster::ReservoirCouplingMaster( ) : comm_{comm}, schedule_{schedule} -{ } +{ + this->start_date_ = this->schedule_.getStartTime(); +} // ------------------ // Public methods // ------------------ +double ReservoirCouplingMaster::maybeChopSubStep( + double suggested_timestep_original, double elapsed_time +) const +{ + // Check if the suggested timestep needs to be adjusted based on the slave processes' + // next report step, or if the slave process has not started yet: the start of a slave process. + double start_date = static_cast(this->start_date_); + TimePoint step_start_date{start_date + elapsed_time}; + TimePoint step_end_date{step_start_date + suggested_timestep_original}; + TimePoint suggested_timestep{suggested_timestep_original}; + for (unsigned int i = 0; i < this->num_slaves_; i++) { + double slave_start_date = static_cast(this->slave_start_dates_[i]); + TimePoint slave_next_report_date{this->slave_next_report_time_offsets_[i] + slave_start_date}; + if (slave_start_date > step_end_date) { + // The slave process has not started yet, and will not start during this timestep + continue; + } + TimePoint slave_elapsed_time; + if (slave_start_date <= step_start_date) { + // The slave process has already started, and will continue during this timestep + if (slave_next_report_date > step_end_date) { + // The slave process will not report during this timestep + continue; + } + // The slave process will report during this timestep + slave_elapsed_time = slave_next_report_date - step_start_date; + } + else { + // The slave process will start during the timestep, but not at the beginning + slave_elapsed_time = slave_start_date - step_start_date; + } + suggested_timestep = slave_elapsed_time; + step_end_date = step_start_date + suggested_timestep; + } + return suggested_timestep.getTime(); +} + +void ReservoirCouplingMaster::sendNextTimeStepToSlaves(double timestep) { + if (this->comm_.rank() == 0) { + for (unsigned int i = 0; i < this->master_slave_comm_.size(); i++) { + MPI_Send( + ×tep, + /*count=*/1, + /*datatype=*/MPI_DOUBLE, + /*dest_rank=*/0, + /*tag=*/static_cast(MessageTag::SlaveNextTimeStep), + *this->master_slave_comm_[i].get() + ); + OpmLog::info(fmt::format( + "Sent next time step {} from master process rank 0 to slave process " + "rank 0 with name: {}", timestep, this->slave_names_[i]) + ); + } + } +} + +void ReservoirCouplingMaster::receiveNextReportDateFromSlaves() { + if (this->comm_.rank() == 0) { + for (unsigned int i = 0; i < this->master_slave_comm_.size(); i++) { + double slave_next_report_time_offset; // Elapsed time from the beginning of the simulation + int result = MPI_Recv( + &slave_next_report_time_offset, + /*count=*/1, + /*datatype=*/MPI_DOUBLE, + /*source_rank=*/0, + /*tag=*/static_cast(MessageTag::SlaveNextReportDate), + *this->master_slave_comm_[i].get(), + MPI_STATUS_IGNORE + ); + if (result != MPI_SUCCESS) { + OPM_THROW(std::runtime_error, "Failed to receive next report date from slave process"); + } + this->slave_next_report_time_offsets_[i] = slave_next_report_time_offset; + OpmLog::info( + fmt::format( + "Received simulation slave next report date from slave process with name: {}. " + "Next report date: {}", this->slave_names_[i], slave_next_report_time_offset + ) + ); + } + } + this->comm_.broadcast( + this->slave_next_report_time_offsets_.data(), this->num_slaves_, /*emitter_rank=*/0 + ); + OpmLog::info("Broadcasted slave next report dates to all ranks"); +} + void ReservoirCouplingMaster::receiveSimulationStartDateFromSlaves() { - this->slave_start_dates_.resize(this->num_slaves_); if (this->comm_.rank() == 0) { // Ensure that std::time_t is of type long since we are sending it over MPI with MPI_LONG static_assert(std::is_same::value, "std::time_t is not of type long"); for (unsigned int i = 0; i < this->master_slave_comm_.size(); i++) { - std::time_t start_date; + std::time_t slave_start_date; int result = MPI_Recv( - &start_date, + &slave_start_date, /*count=*/1, /*datatype=*/MPI_LONG, /*source_rank=*/0, - /*tag=*/static_cast(MessageTag::SimulationStartDate), + /*tag=*/static_cast(MessageTag::SlaveSimulationStartDate), *this->master_slave_comm_[i].get(), MPI_STATUS_IGNORE ); if (result != MPI_SUCCESS) { OPM_THROW(std::runtime_error, "Failed to receive simulation start date from slave process"); } - this->slave_start_dates_[i] = start_date; + if (slave_start_date < this->start_date_) { + OPM_THROW(std::runtime_error, "Slave process start date is before master process start date"); + } + this->slave_start_dates_[i] = slave_start_date; OpmLog::info( fmt::format( "Received simulation start date from slave process with name: {}. " - "Start date: {}", this->slave_names_[i], start_date + "Start date: {}", this->slave_names_[i], slave_start_date ) ); } } this->comm_.broadcast(this->slave_start_dates_.data(), this->num_slaves_, /*emitter_rank=*/0); + OpmLog::info("Broadcasted slave start dates to all ranks"); } // NOTE: This functions is executed for all ranks, but only rank 0 will spawn @@ -91,8 +183,8 @@ void ReservoirCouplingMaster::spawnSlaveProcesses(int argc, char **argv) { const auto& data_file_name = slave.dataFilename(); const auto& directory_path = slave.directoryPath(); // Concatenate the directory path and the data file name to get the full path - std::filesystem::path dir_path(directory_path); - std::filesystem::path data_file(data_file_name); + std::filesystem::path dir_path{directory_path}; + std::filesystem::path data_file{data_file_name}; std::filesystem::path full_path = dir_path / data_file; std::string log_filename; // the getSlaveArgv() function will set this std::vector slave_argv = getSlaveArgv( @@ -134,6 +226,8 @@ void ReservoirCouplingMaster::spawnSlaveProcesses(int argc, char **argv) { this->slave_names_.push_back(slave_name); this->num_slaves_++; } + this->slave_start_dates_.resize(this->num_slaves_); + this->slave_next_report_time_offsets_.resize(this->num_slaves_); } // ------------------ diff --git a/opm/simulators/flow/ReservoirCouplingMaster.hpp b/opm/simulators/flow/ReservoirCouplingMaster.hpp index 6d1ed57b1..182ec0f01 100644 --- a/opm/simulators/flow/ReservoirCouplingMaster.hpp +++ b/opm/simulators/flow/ReservoirCouplingMaster.hpp @@ -36,11 +36,14 @@ class ReservoirCouplingMaster { public: using MPI_Comm_Ptr = ReservoirCoupling::MPI_Comm_Ptr; using MessageTag = ReservoirCoupling::MessageTag; - + using TimePoint = ReservoirCoupling::TimePoint; ReservoirCouplingMaster(const Parallel::Communication &comm, const Schedule &schedule); + double maybeChopSubStep(double suggested_timestep, double current_time) const; void spawnSlaveProcesses(int argc, char **argv); void receiveSimulationStartDateFromSlaves(); + void receiveNextReportDateFromSlaves(); + void sendNextTimeStepToSlaves(double next_time_step); private: std::vector getSlaveArgv( @@ -53,10 +56,12 @@ private: const Parallel::Communication &comm_; const Schedule& schedule_; + std::time_t start_date_; // Master process' simulation start date std::size_t num_slaves_ = 0; // Initially zero, will be updated in spawnSlaveProcesses() std::vector master_slave_comm_; // MPI communicators for the slave processes std::vector slave_names_; std::vector slave_start_dates_; + std::vector slave_next_report_time_offsets_; // Elapsed time from the beginning of the simulation }; } // namespace Opm diff --git a/opm/simulators/flow/ReservoirCouplingSlave.cpp b/opm/simulators/flow/ReservoirCouplingSlave.cpp index f369b0bf1..5f567feb6 100644 --- a/opm/simulators/flow/ReservoirCouplingSlave.cpp +++ b/opm/simulators/flow/ReservoirCouplingSlave.cpp @@ -29,15 +29,18 @@ #include #include +#include namespace Opm { ReservoirCouplingSlave::ReservoirCouplingSlave( const Parallel::Communication &comm, - const Schedule &schedule + const Schedule &schedule, + const SimulatorTimer &timer ) : comm_{comm}, - schedule_{schedule} + schedule_{schedule}, + timer_{timer} { this->slave_master_comm_ = MPI_Comm_Ptr(new MPI_Comm(MPI_COMM_NULL)); MPI_Comm_get_parent(this->slave_master_comm_.get()); @@ -46,8 +49,51 @@ ReservoirCouplingSlave::ReservoirCouplingSlave( } } -void ReservoirCouplingSlave::sendSimulationStartDateToMasterProcess() { +double ReservoirCouplingSlave::receiveNextTimeStepFromMaster() { + double timestep; + if (this->comm_.rank() == 0) { + int result = MPI_Recv( + ×tep, + /*count=*/1, + /*datatype=*/MPI_DOUBLE, + /*source_rank=*/0, + /*tag=*/static_cast(MessageTag::SlaveNextTimeStep), + *this->slave_master_comm_, + MPI_STATUS_IGNORE + ); + if (result != MPI_SUCCESS) { + OPM_THROW(std::runtime_error, "Failed to receive next time step from master"); + } + OpmLog::info( + fmt::format("Slave rank 0 received next timestep {} from master.", timestep) + ); + } + this->comm_.broadcast(×tep, 1, /*emitter_rank=*/0); + OpmLog::info("Broadcasted slave next time step to all ranks"); + return timestep; +} + +void ReservoirCouplingSlave::sendNextReportDateToMasterProcess() { + if (this->comm_.rank() == 0) { + double elapsed_time = this->timer_.simulationTimeElapsed(); + double current_step_length = this->timer_.currentStepLength(); + // NOTE: This is an offset in seconds from the start date, so it will be 0 if the next report + // would be the start date. In general, it should be a positive number. + double next_report_time_offset = elapsed_time + current_step_length; + MPI_Send( + &next_report_time_offset, + /*count=*/1, + /*datatype=*/MPI_DOUBLE, + /*dest_rank=*/0, + /*tag=*/static_cast(MessageTag::SlaveNextReportDate), + *this->slave_master_comm_ + ); + OpmLog::info("Sent next report date to master process from rank 0"); + } +} + +void ReservoirCouplingSlave::sendSimulationStartDateToMasterProcess() { if (this->comm_.rank() == 0) { // Ensure that std::time_t is of type long since we are sending it over MPI with MPI_LONG static_assert(std::is_same::value, "std::time_t is not of type long"); @@ -57,7 +103,7 @@ void ReservoirCouplingSlave::sendSimulationStartDateToMasterProcess() { /*count=*/1, /*datatype=*/MPI_LONG, /*dest_rank=*/0, - /*tag=*/static_cast(MessageTag::SimulationStartDate), + /*tag=*/static_cast(MessageTag::SlaveSimulationStartDate), *this->slave_master_comm_ ); OpmLog::info("Sent simulation start date to master process from rank 0"); diff --git a/opm/simulators/flow/ReservoirCouplingSlave.hpp b/opm/simulators/flow/ReservoirCouplingSlave.hpp index e34b8cfeb..b7579a890 100644 --- a/opm/simulators/flow/ReservoirCouplingSlave.hpp +++ b/opm/simulators/flow/ReservoirCouplingSlave.hpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include @@ -36,15 +37,19 @@ public: using MPI_Comm_Ptr = ReservoirCoupling::MPI_Comm_Ptr; using MessageTag = ReservoirCoupling::MessageTag; - ReservoirCouplingSlave(const Parallel::Communication &comm, const Schedule &schedule); + ReservoirCouplingSlave( + const Parallel::Communication &comm, const Schedule &schedule, const SimulatorTimer &timer + ); void sendSimulationStartDateToMasterProcess(); + void sendNextReportDateToMasterProcess(); + double receiveNextTimeStepFromMaster(); private: const Parallel::Communication &comm_; const Schedule& schedule_; + const SimulatorTimer &timer_; // MPI parent communicator for a slave process MPI_Comm_Ptr slave_master_comm_{nullptr}; - }; } // namespace Opm diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp index 1a8dbfb40..07f7bcb19 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp @@ -206,11 +206,12 @@ public: { #if HAVE_MPI auto slave_mode = Parameters::Get(); + FlowGenericVanguard::comm().barrier(); if (slave_mode) { this->reservoirCouplingSlave_ = std::make_unique( FlowGenericVanguard::comm(), - this->schedule() + this->schedule(), timer ); this->reservoirCouplingSlave_->sendSimulationStartDateToMasterProcess(); } @@ -252,7 +253,14 @@ public: else { adaptiveTimeStepping_ = std::make_unique(unitSystem, max_next_tstep, terminalOutput_); } - +#if HAVE_MPI + if (slave_mode) { + adaptiveTimeStepping_->setReservoirCouplingSlave(this->reservoirCouplingSlave_.get()); + } + else if (this->schedule()[0].rescoup().masterMode()) { + adaptiveTimeStepping_->setReservoirCouplingMaster(this->reservoirCouplingMaster_.get()); + } +#endif if (isRestart()) { // For restarts the simulator may have gotten some information // about the next timestep size from the OPMEXTRA field diff --git a/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp b/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp index 2c36efcf5..4fe45e5ab 100644 --- a/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp +++ b/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp @@ -36,25 +36,28 @@ namespace Opm { AdaptiveSimulatorTimer:: - AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer, - const double lastStepTaken, - const double maxTimeStep ) - : start_date_time_( std::make_shared(timer.startDateTime()) ) - , start_time_( timer.simulationTimeElapsed() ) - , total_time_( start_time_ + timer.currentStepLength() ) - , report_step_( timer.reportStepNum() ) - , max_time_step_( maxTimeStep ) - , current_time_( start_time_ ) - , dt_( 0.0 ) - , current_step_( 0 ) - , steps_() - , lastStepFailed_( false ) + AdaptiveSimulatorTimer( const boost::posix_time::ptime simulation_start_time, + const double step_length, + const double elapsed_time, + const double last_step_taken, + const int report_step, + const double max_time_step ) + : start_date_time_{ std::make_shared(simulation_start_time) } + , start_time_{ elapsed_time } + , total_time_{ start_time_ + step_length } + , report_step_{ report_step } + , max_time_step_{ max_time_step } + , current_time_{ start_time_ } + , dt_{ 0.0 } + , current_step_{ 0 } + , steps_{} + , last_step_failed_{ false } { // reserve memory for sub steps steps_.reserve( 10 ); // set appropriate value for dt_ - provideTimeStepEstimate( lastStepTaken ); + provideTimeStepEstimate( last_step_taken ); } bool AdaptiveSimulatorTimer::initialStep () const diff --git a/opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp b/opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp index 3dffd4908..31ae49f4c 100644 --- a/opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp +++ b/opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp @@ -44,9 +44,12 @@ namespace Opm /// \param timer in case of sub stepping this is the outer timer /// \param lastStepTaken last suggested time step /// \param maxTimeStep maximum time step allowed - AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer, - const double lastStepTaken, - const double maxTimeStep = std::numeric_limits::max() ); + AdaptiveSimulatorTimer( const boost::posix_time::ptime simulation_start_time, + const double step_length, + const double elapsed_time, + const double last_step_taken, + const int report_step, + const double max_time_step = std::numeric_limits::max() ); /// \brief advance time by currentStepLength AdaptiveSimulatorTimer& operator++ (); @@ -101,10 +104,10 @@ namespace Opm boost::posix_time::ptime startDateTime() const; /// \brief Return true if last time step failed - bool lastStepFailed() const {return lastStepFailed_;} + bool lastStepFailed() const {return last_step_failed_;} /// \brief tell the timestepper whether timestep failed or not - void setLastStepFailed(bool lastStepFailed) {lastStepFailed_ = lastStepFailed;} + void setLastStepFailed(bool last_step_failed) {last_step_failed_ = last_step_failed;} /// return copy of object virtual std::unique_ptr< SimulatorTimerInterface > clone() const; @@ -121,7 +124,7 @@ namespace Opm int current_step_; std::vector< double > steps_; - bool lastStepFailed_; + bool last_step_failed_; }; diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp index 822d5ee8a..daf0f4553 100644 --- a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp @@ -19,12 +19,23 @@ #include #include -#include +#if HAVE_MPI +#define RESERVOIR_COUPLING_ENABLED +#endif +#ifdef RESERVOIR_COUPLING_ENABLED +#include +#include +#include +#endif + #include #include #include +#include +#include #include #include +#include #include namespace Opm::Parameters { @@ -55,148 +66,208 @@ class UnitSystem; struct StepReport; namespace detail { - -void logTimer(const AdaptiveSimulatorTimer& substepTimer); - -std::set consistentlyFailingWells(const std::vector& sr); - -void registerAdaptiveParameters(); - -std::tuple, - bool> -createController(const UnitSystem& unitSystem); - + void logTimer(const AdaptiveSimulatorTimer& substep_timer); + std::set consistentlyFailingWells(const std::vector& sr); + void registerAdaptiveParameters(); } - // AdaptiveTimeStepping - //--------------------- - template - class AdaptiveTimeStepping +template +class AdaptiveTimeStepping +{ +private: + using Scalar = GetPropType; + template + class SolutionTimeErrorSolverWrapper : public RelativeChangeInterface { - using Scalar = GetPropType; - template - class SolutionTimeErrorSolverWrapper : public RelativeChangeInterface - { - const Solver& solver_; - public: - SolutionTimeErrorSolverWrapper(const Solver& solver) - : solver_(solver) - {} - - /// return || u^n+1 - u^n || / || u^n+1 || - double relativeChange() const - { return solver_.model().relativeChange(); } - }; - - template - void logException_(const E& exception, bool verbose) - { - if (verbose) { - std::string message; - message = "Caught Exception: "; - message += exception.what(); - OpmLog::debug(message); - } - } - public: - AdaptiveTimeStepping() = default; + SolutionTimeErrorSolverWrapper(const Solver& solver); + double relativeChange() const; + private: + const Solver& solver_; + }; - //! \brief contructor taking parameter object - explicit AdaptiveTimeStepping(const UnitSystem& unitSystem, - const double max_next_tstep = -1.0, - const bool terminalOutput = true); + // Forward declaration of SubStepIteration + // TODO: This forward declaration is not enough for GCC version 9.4.0 (released June 2021), + // but it works fine with GCC version 13.2.0 (released July 2023). + template class SubStepIteration; - //! \brief contructor taking parameter object - //! \param tuning Pointer to ecl TUNING keyword - //! \param timeStep current report step - AdaptiveTimeStepping(double max_next_tstep, - const Tuning& tuning, - const UnitSystem& unitSystem, - const bool terminalOutput = true); + template + class SubStepper { + public: + SubStepper( + AdaptiveTimeStepping& adaptive_time_stepping, + const SimulatorTimer& simulator_timer, + Solver& solver, + const bool is_event, + const std::function& tuning_updater + ); - static void registerParameters(); - - /** \brief step method that acts like the solver::step method - in a sub cycle of time steps - \param tuningUpdater Function used to update TUNING parameters before each - time step. ACTIONX might change tuning. - */ - template - SimulatorReport step(const SimulatorTimer& simulatorTimer, - Solver& solver, - const bool isEvent, - const std::function tuningUpdater); - - /** \brief Returns the simulator report for the failed substeps of the last - * report step. - */ - double suggestedNextStep() const - { return suggestedNextTimestep_; } - - void setSuggestedNextStep(const double x) - { suggestedNextTimestep_ = x; } - - void updateTUNING(double max_next_tstep, const Tuning& tuning); - - void updateNEXTSTEP(double max_next_tstep); - - template - void serializeOp(Serializer& serializer); - - static AdaptiveTimeStepping serializationTestObjectHardcoded(); - static AdaptiveTimeStepping serializationTestObjectPID(); - static AdaptiveTimeStepping serializationTestObjectPIDIt(); - static AdaptiveTimeStepping serializationTestObjectSimple(); - - bool operator==(const AdaptiveTimeStepping& rhs) const; + AdaptiveTimeStepping& getAdaptiveTimerStepper(); + SimulatorReport run(); + friend class AdaptiveTimeStepping::template SubStepIteration; private: - template - static AdaptiveTimeStepping serializationTestObject_(); + bool isReservoirCouplingMaster_() const; + bool isReservoirCouplingSlave_() const; + void maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double originalTimeStep); + bool maybeUpdateTuning_(double elapsed, double dt, int sub_step_number) const; + double maxTimeStep_() const; + SimulatorReport runStepOriginal_(); +#ifdef RESERVOIR_COUPLING_ENABLED + ReservoirCouplingMaster& reservoirCouplingMaster_(); + ReservoirCouplingSlave& reservoirCouplingSlave_(); + SimulatorReport runStepReservoirCouplingMaster_(); + SimulatorReport runStepReservoirCouplingSlave_(); +#endif + double suggestedNextTimestep_() const; - template - void allocAndSerialize(Serializer& serializer) - { - if (!serializer.isSerializing()) { - timeStepControl_ = std::make_unique(); - } - serializer(*static_cast(timeStepControl_.get())); - } - - template - bool castAndComp(const AdaptiveTimeStepping& Rhs) const - { - const T* lhs = static_cast(timeStepControl_.get()); - const T* rhs = static_cast(Rhs.timeStepControl_.get()); - return *lhs == *rhs; - } - - protected: - void init_(const UnitSystem& unitSystem); - - using TimeStepController = std::unique_ptr; - - 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 + AdaptiveTimeStepping& adaptive_time_stepping_; + const SimulatorTimer& simulator_timer_; + Solver& solver_; + const bool is_event_; + const std::function& tuning_updater_; }; -} + + template + class SubStepIteration { + public: + SubStepIteration( + SubStepper& substepper, + AdaptiveSimulatorTimer& substep_timer, + const double original_time_step, + const bool final_step + ); + + SimulatorReport run(); + + private: + bool checkContinueOnUnconvergedSolution_(double dt) const; + void checkTimeStepMaxRestartLimit_(const int restarts) const; + void checkTimeStepMinLimit_(const int new_time_step) const; + void chopTimeStep_(const double new_time_step); + bool chopTimeStepOrCloseFailingWells_(const int new_time_step); + boost::posix_time::ptime currentDateTime_() const; + int getNumIterations_(const SimulatorReportSingle &substep_report) const; + double growthFactor_() const; + bool ignoreConvergenceFailure_() const; + void maybeReportSubStep_(SimulatorReportSingle substep_report) const; + double maybeRestrictTimeStepGrowth_( + const double dt, double dt_estimate, const int restarts) const; + void maybeUpdateTuningAndTimeStep_(); + double maxGrowth_() const; + double minTimeStepBeforeClosingWells_() const; + double minTimeStep_() const; + double restartFactor_() const; + SimulatorReportSingle runSubStep_(); + int solverRestartMax_() const; + double suggestedNextTimestep_() const; + void setSuggestedNextStep_(double step); + void setTimeStep_(double dt_estimate); + Solver& solver_() const; + bool solverVerbose_() const; + const SimulatorTimer& simulatorTimer_() const; + boost::posix_time::ptime startDateTime_() const; + double timeStepControlComputeEstimate_( + const double dt, const int iterations, double elapsed) const; + bool timeStepVerbose_() const; + void updateSuggestedNextStep_(); + bool useNewtonIteration_() const; + double writeOutput_() const; + + SubStepper& substepper_; + AdaptiveSimulatorTimer& substep_timer_; + const double original_time_step_; + const bool final_step_; + std::string cause_of_failure_; + AdaptiveTimeStepping& adaptive_time_stepping_; + }; + +public: + AdaptiveTimeStepping() = default; + + AdaptiveTimeStepping( + const UnitSystem& unitSystem, + const double max_next_tstep = -1.0, + const bool terminalOutput = true + ); + + AdaptiveTimeStepping( + double max_next_tstep, + const Tuning& tuning, + const UnitSystem& unitSystem, + const bool terminalOutput = true + ); + bool operator==(const AdaptiveTimeStepping& rhs); + + static void registerParameters(); +#ifdef RESERVOIR_COUPLING_ENABLED + void setReservoirCouplingMaster(ReservoirCouplingMaster *reservoir_coupling_master); + void setReservoirCouplingSlave(ReservoirCouplingSlave *reservoir_coupling_slave); +#endif + void setSuggestedNextStep(const double x); + double suggestedNextStep() const; + + template + SimulatorReport step(const SimulatorTimer& simulator_timer, + Solver& solver, + const bool is_event, + const std::function + tuning_updater); + + void updateTUNING(double max_next_tstep, const Tuning& tuning); + void updateNEXTSTEP(double max_next_tstep); + + template + void serializeOp(Serializer& serializer); + + static AdaptiveTimeStepping serializationTestObjectHardcoded(); + static AdaptiveTimeStepping serializationTestObjectPID(); + static AdaptiveTimeStepping serializationTestObjectPIDIt(); + static AdaptiveTimeStepping serializationTestObjectSimple(); + +private: + void maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step, + const bool is_event); + + template + static AdaptiveTimeStepping serializationTestObject_(); + + template + void allocAndSerialize(Serializer& serializer); + + template + bool castAndComp(const AdaptiveTimeStepping& Rhs) const; + +protected: + void init_(const UnitSystem& unitSystem); + + using TimeStepController = std::unique_ptr; + + TimeStepControlType time_step_control_type_; //!< type of time step control object + TimeStepController time_step_control_; //!< time step control object + double restart_factor_; //!< factor to multiply time step with when solver fails to converge + double growth_factor_; //!< factor to multiply time step when solver recovered from failed convergence + double max_growth_; //!< factor that limits the maximum growth of a time step + double max_time_step_; //!< maximal allowed time step size in days + double min_time_step_; //!< minimal allowed time step size before throwing + bool ignore_convergence_failure_; //!< continue instead of stop when minimum time step is reached + int solver_restart_max_; //!< how many restart of solver are allowed + bool solver_verbose_; //!< solver verbosity + bool timestep_verbose_; //!< timestep verbosity + double suggested_next_timestep_; //!< suggested size of next timestep + bool full_timestep_initially_; //!< beginning with the size of the time step from data file + double timestep_after_event_; //!< suggested size of timestep after an event + bool use_newton_iteration_; //!< use newton iteration count for adaptive time step control + + //! < shut problematic wells when time step size in days are less than this + double min_time_step_before_shutting_problematic_wells_; +#ifdef RESERVOIR_COUPLING_ENABLED + ReservoirCouplingMaster *reservoir_coupling_master_ = nullptr; + ReservoirCouplingSlave *reservoir_coupling_slave_ = nullptr; +#endif +}; + +} // namespace Opm #include - #endif // OPM_ADAPTIVE_TIME_STEPPING_HPP diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp index 26343d47f..74c6c2c23 100644 --- a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp @@ -1,11 +1,24 @@ /* + Copyright 2024 Equinor ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . */ + #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 @@ -31,387 +44,178 @@ #include namespace Opm { +/********************************************* + * Public methods of AdaptiveTimeStepping + * ******************************************/ + +//! \brief contructor taking parameter object 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) +AdaptiveTimeStepping::AdaptiveTimeStepping( + const UnitSystem& unit_system, + const double max_next_tstep, + const bool terminal_output +) + : time_step_control_{} + , restart_factor_{Parameters::Get>()} // 0.33 + , growth_factor_{Parameters::Get>()} // 2.0 + , max_growth_{Parameters::Get>()} // 3.0 + , max_time_step_{ + Parameters::Get>() * 24 * 60 * 60} // 365.25 + , min_time_step_{ + unit_system.to_si(UnitSystem::measure::time, + Parameters::Get>())} // 1e-12; + , ignore_convergence_failure_{ + Parameters::Get()} // false; + , solver_restart_max_{Parameters::Get()} // 10 + , solver_verbose_{Parameters::Get() > 0 && terminal_output} // 2 + , timestep_verbose_{Parameters::Get() > 0 && terminal_output} // 2 + , suggested_next_timestep_{ + (max_next_tstep <= 0 ? Parameters::Get() + : max_next_tstep) * 24 * 60 * 60} // 1.0 + , full_timestep_initially_{Parameters::Get()} // false + , timestep_after_event_{ + Parameters::Get>() * 24 * 60 * 60} // 1e30 + , use_newton_iteration_{false} + , min_time_step_before_shutting_problematic_wells_{ + Parameters::Get() * unit::day} { - init_(unitSystem); + init_(unit_system); +} + +//! \brief contructor +//! \param tuning Pointer to ecl TUNING keyword +template +AdaptiveTimeStepping::AdaptiveTimeStepping(double max_next_tstep, + const Tuning& tuning, + const UnitSystem& unit_system, + const bool terminal_output +) + : time_step_control_{} + , restart_factor_{tuning.TSFCNV} + , growth_factor_{tuning.TFDIFF} + , max_growth_{tuning.TSFMAX} + , max_time_step_{tuning.TSMAXZ} // 365.25 + , min_time_step_{tuning.TSFMIN} // 0.1; + , ignore_convergence_failure_{true} + , solver_restart_max_{Parameters::Get()} // 10 + , solver_verbose_{Parameters::Get() > 0 && terminal_output} // 2 + , timestep_verbose_{Parameters::Get() > 0 && terminal_output} // 2 + , suggested_next_timestep_{ + max_next_tstep <= 0 ? Parameters::Get() * 24 * 60 * 60 + : max_next_tstep} // 1.0 + , full_timestep_initially_{Parameters::Get()} // false + , timestep_after_event_{tuning.TMAXWC} // 1e30 + , use_newton_iteration_{false} + , min_time_step_before_shutting_problematic_wells_{ + Parameters::Get() * unit::day} +{ + init_(unit_system); } template +bool 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) +operator==(const AdaptiveTimeStepping& rhs) { - init_(unitSystem); + if (this->time_step_control_type_ != rhs.time_step_control_type_ || + (this->time_step_control_ && !rhs.time_step_control_) || + (!this->time_step_control_ && rhs.time_step_control_)) { + return false; + } + + bool result = false; + switch (this->time_step_control_type_) { + 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->restart_factor_ == rhs.restart_factor_ && + this->growth_factor_ == rhs.growth_factor_ && + this->max_growth_ == rhs.max_growth_ && + this->max_time_step_ == rhs.max_time_step_ && + this->min_time_step_ == rhs.min_time_step_ && + this->ignore_convergence_failure_ == rhs.ignore_convergence_failure_ && + this->solver_restart_max_== rhs.solver_restart_max_ && + this->solver_verbose_ == rhs.solver_verbose_ && + this->full_timestep_initially_ == rhs.full_timestep_initially_ && + this->timestep_after_event_ == rhs.timestep_after_event_ && + this->use_newton_iteration_ == rhs.use_newton_iteration_ && + this->min_time_step_before_shutting_problematic_wells_ == + rhs.min_time_step_before_shutting_problematic_wells_; } template -void AdaptiveTimeStepping::registerParameters() +void +AdaptiveTimeStepping:: +registerParameters() { registerEclTimeSteppingParameters(); detail::registerAdaptiveParameters(); } +#ifdef RESERVOIR_COUPLING_ENABLED template -template +void +AdaptiveTimeStepping:: +setReservoirCouplingMaster(ReservoirCouplingMaster *reservoir_coupling_master) +{ + this->reservoir_coupling_master_ = reservoir_coupling_master; +} + +template +void +AdaptiveTimeStepping:: +setReservoirCouplingSlave(ReservoirCouplingSlave *reservoir_coupling_slave) +{ + this->reservoir_coupling_slave_ = reservoir_coupling_slave; +} +#endif + +/** \brief step method that acts like the solver::step method + in a sub cycle of time steps + \param tuningUpdater Function used to update TUNING parameters before each + time step. ACTIONX might change tuning. +*/ +template +template SimulatorReport AdaptiveTimeStepping:: -step(const SimulatorTimer& simulatorTimer, +step(const SimulatorTimer& simulator_timer, Solver& solver, - const bool isEvent, - const std::function tuningUpdater) + const bool is_event, + const std::function tuning_updater +) { - // 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, - 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_) { - 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 -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; - } + SubStepper sub_stepper{ + *this, simulator_timer, solver, is_event, tuning_updater, + }; + return sub_stepper.run(); } template template -void AdaptiveTimeStepping:: +void +AdaptiveTimeStepping:: serializeOp(Serializer& serializer) { - serializer(timeStepControlType_); - switch (timeStepControlType_) { + serializer(this->time_step_control_type_); + switch (this->time_step_control_type_) { case TimeStepControlType::HardCodedTimeStep: allocAndSerialize(serializer); break; @@ -425,20 +229,20 @@ serializeOp(Serializer& serializer) 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_); + serializer(this->restart_factor_); + serializer(this->growth_factor_); + serializer(this->max_growth_); + serializer(this->max_time_step_); + serializer(this->min_time_step_); + serializer(this->ignore_convergence_failure_); + serializer(this->solver_restart_max_); + serializer(this->solver_verbose_); + serializer(this->timestep_verbose_); + serializer(this->suggested_next_timestep_); + serializer(this->full_timestep_initially_); + serializer(this->timestep_after_event_); + serializer(this->use_newton_iteration_); + serializer(this->min_time_step_before_shutting_problematic_wells_); } template @@ -473,47 +277,95 @@ serializationTestObjectSimple() return serializationTestObject_(); } + template +void +AdaptiveTimeStepping:: +setSuggestedNextStep(const double x) +{ + this->suggested_next_timestep_ = x; +} + +template +double +AdaptiveTimeStepping:: +suggestedNextStep() const +{ + return this->suggested_next_timestep_; +} + + +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) { + this->suggested_next_timestep_ = max_next_tstep; + } +} + +template +void +AdaptiveTimeStepping:: +updateTUNING(double max_next_tstep, const Tuning& tuning) +{ + this->restart_factor_ = tuning.TSFCNV; + this->growth_factor_ = tuning.TFDIFF; + this->max_growth_ = tuning.TSFMAX; + this->max_time_step_ = tuning.TSMAXZ; + updateNEXTSTEP(max_next_tstep); + this->timestep_after_event_ = tuning.TMAXWC; +} + +/********************************************* + * Private methods of AdaptiveTimeStepping + * ******************************************/ + +template +template +void +AdaptiveTimeStepping:: +allocAndSerialize(Serializer& serializer) +{ + if (!serializer.isSerializing()) { + this->time_step_control_ = std::make_unique(); + } + serializer(*static_cast(this->time_step_control_.get())); +} + +template +template bool AdaptiveTimeStepping:: -operator==(const AdaptiveTimeStepping& rhs) const +castAndComp(const AdaptiveTimeStepping& Rhs) const { - if (timeStepControlType_ != rhs.timeStepControlType_ || - (timeStepControl_ && !rhs.timeStepControl_) || - (!timeStepControl_ && rhs.timeStepControl_)) { - return false; + const T* lhs = static_cast(this->time_step_control_.get()); + const T* rhs = static_cast(Rhs.time_step_control_.get()); + return *lhs == *rhs; +} + +template +void +AdaptiveTimeStepping:: +maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step, + bool is_event) +{ + // init last time step as a fraction of the given time step + if (this->suggested_next_timestep_ < 0) { + this->suggested_next_timestep_ = this->restart_factor_ * original_time_step; } - 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; + if (this->full_timestep_initially_) { + this->suggested_next_timestep_ = original_time_step; } - 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_; + // use seperate time step after event + if (is_event && this->timestep_after_event_ > 0) { + this->suggested_next_timestep_ = this->timestep_after_event_; + } } template @@ -524,33 +376,37 @@ 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()); + result.restart_factor_ = 1.0; + result.growth_factor_ = 2.0; + result.max_growth_ = 3.0; + result.max_time_step_ = 4.0; + result.min_time_step_ = 5.0; + result.ignore_convergence_failure_ = true; + result.solver_restart_max_ = 6; + result.solver_verbose_ = true; + result.timestep_verbose_ = true; + result.suggested_next_timestep_ = 7.0; + result.full_timestep_initially_ = true; + result.use_newton_iteration_ = true; + result.min_time_step_before_shutting_problematic_wells_ = 9.0; + result.time_step_control_type_ = Controller::Type; + result.time_step_control_ = + std::make_unique(Controller::serializationTestObject()); return result; } +/********************************************* + * Protected methods of AdaptiveTimeStepping + * ******************************************/ + template void AdaptiveTimeStepping:: init_(const UnitSystem& unitSystem) { - std::tie(timeStepControlType_, - timeStepControl_, - useNewtonIteration_) = detail::createController(unitSystem); - + std::tie(time_step_control_type_, + time_step_control_, + use_newton_iteration_) = detail::createController(unitSystem); // make sure growth factor is something reasonable if (growthFactor_ < 1.0) { OPM_THROW(std::runtime_error, @@ -558,6 +414,855 @@ init_(const UnitSystem& unitSystem) } } -} // namespace Opm + + +/************************************************ + * Private class SubStepper public methods + ************************************************/ + +template +template +AdaptiveTimeStepping::SubStepper:: +SubStepper( + AdaptiveTimeStepping& adaptive_time_stepping, + const SimulatorTimer& simulator_timer, + Solver& solver, + const bool is_event, + const std::function& tuning_updater + +) + : adaptive_time_stepping_{adaptive_time_stepping} + , simulator_timer_{simulator_timer} + , solver_{solver} + , is_event_{is_event} + , tuning_updater_{tuning_updater} +{ +} + +template +template +AdaptiveTimeStepping& +AdaptiveTimeStepping::SubStepper:: +getAdaptiveTimerStepper() +{ + return adaptive_time_stepping_; +} + +template +template +SimulatorReport +AdaptiveTimeStepping::SubStepper:: +run() +{ +#ifdef RESERVOIR_COUPLING_ENABLED + if (isReservoirCouplingSlave_() && reservoirCouplingSlave_().activated()) { + return runStepReservoirCouplingSlave_(); + } + else if (isReservoirCouplingMaster_() && reservoirCouplingMaster_().activated()) { + return runStepReservoirCouplingMaster_(); + } + else { + return runStepOriginal_(); + } +#else + return runStepOriginal_(); +#endif +} + +/************************************************ + * Private class SubStepper private methods + ************************************************/ + + +template +template +bool +AdaptiveTimeStepping::SubStepper:: +isReservoirCouplingMaster_() const +{ + return this->adaptive_time_stepping_.reservoir_coupling_master_ != nullptr; +} + +template +template +bool +AdaptiveTimeStepping::SubStepper:: +isReservoirCouplingSlave_() const +{ + return this->adaptive_time_stepping_.reservoir_coupling_slave_ != nullptr; +} + +template +template +void +AdaptiveTimeStepping::SubStepper:: +maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step) +{ + this->adaptive_time_stepping_.maybeModifySuggestedTimeStepAtBeginningOfReportStep_( + original_time_step, this->is_event_ + ); +} + +// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicitBlackoil::runStep() +// It has to be called for each substep since TUNING might have been changed for next sub step due +// to ACTIONX (via NEXTSTEP) or WCYCLE keywords. +template +template +bool +AdaptiveTimeStepping::SubStepper:: +maybeUpdateTuning_(double elapsed, double dt, int sub_step_number) const +{ + return this->tuning_updater_(elapsed, dt, sub_step_number); +} + +template +template +double +AdaptiveTimeStepping::SubStepper:: +maxTimeStep_() const +{ + return this->adaptive_time_stepping_.max_time_step_; +} + +template +template +SimulatorReport +AdaptiveTimeStepping::SubStepper:: +runStepOriginal_() +{ + auto elapsed = this->simulator_timer_.simulationTimeElapsed(); + auto original_time_step = this->simulator_timer_.currentStepLength(); + auto report_step = this->simulator_timer_.reportStepNum(); + maybeUpdateTuning_(elapsed, original_time_step, report_step); + maybeModifySuggestedTimeStepAtBeginningOfReportStep_(original_time_step); + + AdaptiveSimulatorTimer substep_timer{ + this->simulator_timer_.startDateTime(), + original_time_step, + elapsed, + suggestedNextTimestep_(), + report_step, + maxTimeStep_() + }; + SubStepIteration substepIteration{*this, substep_timer, original_time_step, /*final_step=*/true}; + return substepIteration.run(); +} + +#ifdef RESERVOIR_COUPLING_ENABLED +template +template +ReservoirCouplingMaster& +AdaptiveTimeStepping::SubStepper:: +reservoirCouplingMaster_() +{ + return *adaptive_time_stepping_.reservoir_coupling_master_; +} +#endif + +#ifdef RESERVOIR_COUPLING_ENABLED +template +template +ReservoirCouplingSlave& +AdaptiveTimeStepping::SubStepper:: +reservoirCouplingSlave_() +{ + return *this->adaptive_time_stepping_.reservoir_coupling_slave_; +} +#endif + +#ifdef RESERVOIR_COUPLING_ENABLED +template +template +SimulatorReport +AdaptiveTimeStepping::SubStepper:: +runStepReservoirCouplingMaster_() +{ + bool substep_done = false; + int iteration = 0; + const double original_time_step = this->simulator_timer_.currentStepLength(); + double current_time{this->simulator_timer_.simulationTimeElapsed()}; + double step_end_time = current_time + original_time_step; + auto current_step_length = original_time_step; + SimulatorReport report; + while(!substep_done) { + reservoirCouplingMaster_().receiveNextReportDateFromSlaves(); + if (iteration == 0) { + maybeUpdateTuning_(current_time, current_step_length, /*substep=*/0); + } + current_step_length = reservoirCouplingMaster_().maybeChopSubStep( + current_step_length, current_time); + reservoirCouplingMaster_().sendNextTimeStepToSlaves(current_step_length); + if (iteration == 0) { + maybeModifySuggestedTimeStepAtBeginningOfReportStep_(current_step_length); + } + AdaptiveSimulatorTimer substep_timer{ + this->simulator_timer_.startDateTime(), + /*stepLength=*/current_step_length, + /*elapsedTime=*/current_time, + /*timeStepEstimate=*/suggestedNextTimestep_(), + /*reportStep=*/this->simulator_timer_.reportStepNum(), + maxTimeStep_() + }; + bool final_step = ReservoirCoupling::Seconds::compare_gt_or_eq( + current_time + current_step_length, step_end_time + ); + SubStepIteration substepIteration{*this, substep_timer, current_step_length, final_step}; + auto sub_steps_report = substepIteration.run(); + report += sub_steps_report; + current_time += current_step_length; + if (final_step) { + break; + } + iteration++; + } + return report; +} +#endif + +#ifdef RESERVOIR_COUPLING_ENABLED +template +template +SimulatorReport +AdaptiveTimeStepping::SubStepper:: +runStepReservoirCouplingSlave_() +{ + bool substep_done = false; + int iteration = 0; + const double original_time_step = this->simulator_timer_.currentStepLength(); + double current_time{this->simulator_timer_.simulationTimeElapsed()}; + double step_end_time = current_time + original_time_step; + SimulatorReport report; + while(!substep_done) { + reservoirCouplingSlave_().sendNextReportDateToMasterProcess(); + auto timestep = reservoirCouplingSlave_().receiveNextTimeStepFromMaster(); + if (iteration == 0) { + maybeUpdateTuning_(current_time, original_time_step, /*substep=*/0); + maybeModifySuggestedTimeStepAtBeginningOfReportStep_(timestep); + } + AdaptiveSimulatorTimer substep_timer{ + this->simulator_timer_.startDateTime(), + /*step_length=*/timestep, + /*elapsed_time=*/current_time, + /*time_step_estimate=*/suggestedNextTimestep_(), + this->simulator_timer_.reportStepNum(), + maxTimeStep_() + }; + bool final_step = ReservoirCoupling::Seconds::compare_gt_or_eq( + current_time + timestep, step_end_time + ); + SubStepIteration substepIteration{*this, substep_timer, timestep, final_step}; + auto sub_steps_report = substepIteration.run(); + report += sub_steps_report; + current_time += timestep; + if (final_step) { + substep_done = true; + break; + } + iteration++; + } + return report; +} #endif + +template +template +double +AdaptiveTimeStepping::SubStepper:: +suggestedNextTimestep_() const +{ + return this->adaptive_time_stepping_.suggestedNextStep(); +} + + + +/************************************************ + * Private class SubStepIteration public methods + ************************************************/ + +template +template +AdaptiveTimeStepping::SubStepIteration:: +SubStepIteration( + SubStepper& substepper, + AdaptiveSimulatorTimer& substep_timer, + const double original_time_step, + bool final_step +) + : substepper_{substepper} + , substep_timer_{substep_timer} + , original_time_step_{original_time_step} + , final_step_{final_step} + , adaptive_time_stepping_{substepper.getAdaptiveTimerStepper()} +{ +} + +template +template +SimulatorReport +AdaptiveTimeStepping::SubStepIteration:: +run() +{ + auto& simulator = solver_().model().simulator(); + auto& problem = simulator.problem(); + // counter for solver restarts + int restarts = 0; + SimulatorReport report; + + // sub step time loop + while (!this->substep_timer_.done()) { + maybeUpdateTuningAndTimeStep_(); + const double dt = this->substep_timer_.currentStepLength(); + if (timeStepVerbose_()) { + detail::logTimer(this->substep_timer_); + } + + auto substep_report = runSubStep_(); + + //Pass substep to eclwriter for summary output + problem.setSubStepReport(substep_report); + + report += substep_report; + + if (substep_report.converged || checkContinueOnUnconvergedSolution_(dt)) { + ++this->substep_timer_; // advance by current dt + + const int iterations = getNumIterations_(substep_report); + auto dt_estimate = timeStepControlComputeEstimate_( + dt, iterations, this->substep_timer_.simulationTimeElapsed()); + + assert(dt_estimate > 0); + dt_estimate = maybeRestrictTimeStepGrowth_(dt, dt_estimate, restarts); + restarts = 0; // solver converged, reset restarts counter + + maybeReportSubStep_(substep_report); + if (this->final_step_ && this->substep_timer_.done()) { + // if the time step is done we do not need to write it as this will be done + // by the simulator anyway. + } + else { + report.success.output_write_time += writeOutput_(); + } + + // set new time step length + setTimeStep_(dt_estimate); + + report.success.converged = this->substep_timer_.done(); + this->substep_timer_.setLastStepFailed(false); + } + else { // in case of no convergence + this->substep_timer_.setLastStepFailed(true); + checkTimeStepMaxRestartLimit_(restarts); + + const double new_time_step = restartFactor_() * dt; + checkTimeStepMinLimit_(new_time_step); + bool wells_shut = false; + if (new_time_step > minTimeStepBeforeClosingWells_()) { + chopTimeStep_(new_time_step); + } else { + wells_shut = chopTimeStepOrCloseFailingWells_(new_time_step); + } + if (wells_shut) { + setTimeStep_(dt); // retry the old timestep + } + else { + restarts++; // only increase if no wells were shut + } + } + problem.setNextTimeStepSize(this->substep_timer_.currentStepLength()); + } + updateSuggestedNextStep_(); + return report; +} + + +/************************************************ + * Private class SubStepIteration private methods + ************************************************/ + + +template +template +bool +AdaptiveTimeStepping::SubStepIteration:: +checkContinueOnUnconvergedSolution_(double dt) const +{ + bool continue_on_uncoverged_solution = ignoreConvergenceFailure_() && dt <= minTimeStep_(); + if (continue_on_uncoverged_solution && solverVerbose_()) { + // NOTE: This method is only called if the solver failed to converge. + 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); + } + return continue_on_uncoverged_solution; +} + +template +template +void +AdaptiveTimeStepping::SubStepIteration:: +checkTimeStepMaxRestartLimit_(const int restarts) const +{ + // 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}; + } +} + +template +template +void +AdaptiveTimeStepping::SubStepIteration:: +checkTimeStepMinLimit_(const int new_time_step) const +{ + // If we have restarted (i.e. cut the timestep) too + // much, we have failed and throw an exception. + if (new_time_step < 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}; + } +} + +template +template +void +AdaptiveTimeStepping::SubStepIteration:: +chopTimeStep_(const double new_time_step) +{ + setTimeStep_(new_time_step); + if (solverVerbose_()) { + const auto msg = fmt::format("{}\nTimestep chopped to {} days\n", + this->cause_of_failure_, + unit::convert::to(this->substep_timer_.currentStepLength(), unit::day)); + OpmLog::problem(msg); + } +} + +template +template +bool +AdaptiveTimeStepping::SubStepIteration:: +chopTimeStepOrCloseFailingWells_(const int new_time_step) +{ + bool wells_shut = false; + // 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 + chopTimeStep_(new_time_step); + } 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, this->substep_timer_.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, this->substep_timer_.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_(new_time_step); + } else { + wells_shut = true; + if (solverVerbose_()) { + const std::string msg = + fmt::format("\nProblematic well(s) were shut: {}" + "(retrying timestep)\n", + fmt::join(shut_wells, " ")); + OpmLog::problem(msg); + } + } + } + return wells_shut; +} + +template +template +boost::posix_time::ptime +AdaptiveTimeStepping::SubStepIteration:: +currentDateTime_() const +{ + return simulatorTimer_().currentDateTime(); +} + +template +template +int +AdaptiveTimeStepping::SubStepIteration:: +getNumIterations_(const SimulatorReportSingle &substep_report) const +{ + if (useNewtonIteration_()) { + return substep_report.total_newton_iterations; + } + else { + return substep_report.total_linear_iterations; + } +} + +template +template +double +AdaptiveTimeStepping::SubStepIteration:: +growthFactor_() const +{ + return this->adaptive_time_stepping_.growth_factor_; +} + +template +template +bool +AdaptiveTimeStepping::SubStepIteration:: +ignoreConvergenceFailure_() const +{ + return adaptive_time_stepping_.ignore_convergence_failure_; +} + +template +template +double +AdaptiveTimeStepping::SubStepIteration:: +maxGrowth_() const +{ + return this->adaptive_time_stepping_.max_growth_; +} + +template +template +void +AdaptiveTimeStepping::SubStepIteration:: +maybeReportSubStep_(SimulatorReportSingle substep_report) const +{ + if (timeStepVerbose_()) { + std::ostringstream ss; + substep_report.reportStep(ss); + OpmLog::info(ss.str()); + } +} + +template +template +double +AdaptiveTimeStepping::SubStepIteration:: +maybeRestrictTimeStepGrowth_(const double dt, double dt_estimate, const int restarts) const +{ + // limit the growth of the timestep size by the growth factor + dt_estimate = std::min(dt_estimate, double(maxGrowth_() * dt)); + assert(dt_estimate > 0); + // further restrict time step size growth after convergence problems + if (restarts > 0) { + dt_estimate = std::min(growthFactor_() * dt, dt_estimate); + } + + return dt_estimate; +} + +// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicitBlackoil::runStep() +// It has to be called for each substep since TUNING might have been changed for next sub step due +// to ACTIONX (via NEXTSTEP) or WCYCLE keywords. +template +template +void +AdaptiveTimeStepping::SubStepIteration:: +maybeUpdateTuningAndTimeStep_() +{ + // TODO: This function is currently only called if NEXTSTEP is activated from ACTIONX or + // if the WCYCLE keyword needs to modify the current timestep. So this method should rather + // be named maybeUpdateTimeStep_() or similar, since it should not update the tuning. However, + // the current definition of the maybeUpdateTuning_() callback is actually calling + // adaptiveTimeStepping_->updateTUNING(max_next_tstep, tuning) which is updating the tuning + // see SimulatorFullyImplicitBlackoil::runStep() for more details. + auto old_value = suggestedNextTimestep_(); + if (this->substepper_.maybeUpdateTuning_(this->substep_timer_.simulationTimeElapsed(), + this->substep_timer_.currentStepLength(), + this->substep_timer_.currentStepNum())) + { + // Either NEXTSTEP and WCYCLE wants to change the current time step, but they cannot + // change the current time step directly. Instead, they change the suggested next time step + // by calling updateNEXTSTEP() via the maybeUpdateTuning() callback. We now need to update + // the current time step to the new suggested time step and reset the suggested time step + // to the old value. + setTimeStep_(suggestedNextTimestep_()); + setSuggestedNextStep_(old_value); + } +} + +template +template +double +AdaptiveTimeStepping::SubStepIteration:: +minTimeStepBeforeClosingWells_() const +{ + return this->adaptive_time_stepping_.min_time_step_before_shutting_problematic_wells_; +} + +template +template +double +AdaptiveTimeStepping::SubStepIteration:: +minTimeStep_() const +{ + return this->adaptive_time_stepping_.min_time_step_; +} + +template +template +double +AdaptiveTimeStepping::SubStepIteration:: +restartFactor_() const +{ + return this->adaptive_time_stepping_.restart_factor_; +} + +template +template +SimulatorReportSingle +AdaptiveTimeStepping::SubStepIteration:: +runSubStep_() +{ + SimulatorReportSingle substep_report; + + auto handleFailure = [this, &substep_report] + (const std::string& failure_reason, const std::exception& e, bool log_exception = true) + { + substep_report = solver_().failureReport(); + this->cause_of_failure_ = failure_reason; + if (log_exception && solverVerbose_()) { + std::string message; + message = "Caught Exception: "; + message += e.what(); + OpmLog::debug(message); + } + }; + + try { + substep_report = solver_().step(this->substep_timer_); + if (solverVerbose_()) { + // report number of linear iterations + OpmLog::debug("Overall linear iterations used: " + + std::to_string(substep_report.total_linear_iterations)); + } + } + catch (const TooManyIterations& e) { + handleFailure("Solver convergence failure - Iteration limit reached", e); + } + catch (const ConvergenceMonitorFailure& e) { + handleFailure("Convergence monitor failure", e, /*log_exception=*/false); + } + catch (const LinearSolverProblem& e) { + handleFailure("Linear solver convergence failure", e); + } + catch (const NumericalProblem& e) { + handleFailure("Solver convergence failure - Numerical problem encountered", e); + } + catch (const std::runtime_error& e) { + handleFailure("Runtime error encountered", e); + } + catch (const Dune::ISTLError& e) { + handleFailure("ISTL error - Time step too large", e); + } + catch (const Dune::MatrixBlockError& e) { + handleFailure("Matrix block error", e); + } + + return substep_report; +} + +template +template +void +AdaptiveTimeStepping::SubStepIteration:: +setTimeStep_(double dt_estimate) +{ + this->substep_timer_.provideTimeStepEstimate(dt_estimate); +} + +template +template +Solver& +AdaptiveTimeStepping::SubStepIteration:: +solver_() const +{ + return this->substepper_.solver_; +} + + +template +template +int +AdaptiveTimeStepping::SubStepIteration:: +solverRestartMax_() const +{ + return this->adaptive_time_stepping_.solver_restart_max_; +} + +template +template +void +AdaptiveTimeStepping::SubStepIteration:: +setSuggestedNextStep_(double step) +{ + this->adaptive_time_stepping_.setSuggestedNextStep(step); +} + +template +template +const SimulatorTimer& +AdaptiveTimeStepping::SubStepIteration:: +simulatorTimer_() const +{ + return this->substepper_.simulator_timer_; +} + +template +template +bool +AdaptiveTimeStepping::SubStepIteration:: +solverVerbose_() const +{ + return this->adaptive_time_stepping_.solver_verbose_; +} + +template +template +boost::posix_time::ptime +AdaptiveTimeStepping::SubStepIteration:: +startDateTime_() const +{ + return simulatorTimer_().startDateTime(); +} + +template +template +double +AdaptiveTimeStepping::SubStepIteration:: +suggestedNextTimestep_() const +{ + return this->adaptive_time_stepping_.suggestedNextStep(); +} + +template +template +double +AdaptiveTimeStepping::SubStepIteration:: +timeStepControlComputeEstimate_(const double dt, const int iterations, double elapsed) 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); +} + +template +template +bool +AdaptiveTimeStepping::SubStepIteration:: +timeStepVerbose_() const +{ + return this->adaptive_time_stepping_.timestep_verbose_; +} + +// The suggested time step is the stepsize that will be used as a first try for +// the next sub step. It is updated at the end of each substep. It can also be +// updated by the TUNING or NEXTSTEP keywords at the beginning of each report step or +// at the beginning of each substep by the ACTIONX keyword (via NEXTSTEP), this is +// done by the maybeUpdateTuning_() method which is called at the beginning of each substep +// (and the begginning of each report step). Note that the WCYCLE keyword can also update the +// suggested time step via the maybeUpdateTuning_() method. +template +template +void +AdaptiveTimeStepping::SubStepIteration:: +updateSuggestedNextStep_() +{ + auto suggested_next_step = this->substep_timer_.currentStepLength(); + if (! std::isfinite(suggested_next_step)) { // check for NaN + suggested_next_step = this->original_time_step_; + } + if (timeStepVerbose_()) { + std::ostringstream ss; + this->substep_timer_.report(ss); + ss << "Suggested next step size = " + << unit::convert::to(suggested_next_step, unit::day) << " (days)" << std::endl; + OpmLog::debug(ss.str()); + } + setSuggestedNextStep_(suggested_next_step); +} + +template +template +bool +AdaptiveTimeStepping::SubStepIteration:: +useNewtonIteration_() const +{ + return this->adaptive_time_stepping_.use_newton_iteration_; +} + +template +template +double +AdaptiveTimeStepping::SubStepIteration:: +writeOutput_() const +{ + time::StopWatch perf_timer; + perf_timer.start(); + auto& problem = solver_().model().simulator().problem(); + problem.writeOutput(true); + return perf_timer.secsSinceStart(); +} + +/************************************************ + * Private class SolutionTimeErrorSolverWrapper + * **********************************************/ + +template +template +AdaptiveTimeStepping::SolutionTimeErrorSolverWrapper::SolutionTimeErrorSolverWrapper( + const Solver& solver +) + : solver_{solver} +{} + +template +template +double AdaptiveTimeStepping::SolutionTimeErrorSolverWrapper::relativeChange() const +{ + // returns: || u^n+1 - u^n || / || u^n+1 || + return solver_.model().relativeChange(); +} + +} // namespace Opm diff --git a/opm/simulators/timestepping/SimulatorTimerInterface.hpp b/opm/simulators/timestepping/SimulatorTimerInterface.hpp index 6875db871..3948609bf 100644 --- a/opm/simulators/timestepping/SimulatorTimerInterface.hpp +++ b/opm/simulators/timestepping/SimulatorTimerInterface.hpp @@ -40,14 +40,15 @@ namespace Opm /// destructor virtual ~SimulatorTimerInterface() {} - /// Current step number. This is the number of timesteps that - /// has been completed from the start of the run. The time - /// after initialization but before the simulation has started - /// is timestep number zero. - virtual int currentStepNum() const = 0; + // ----------------------------------------------------------- + // Pure virtual functions to be implemented by derived classes + // ----------------------------------------------------------- - /// Current report step number. This might differ from currentStepNum in case of sub stepping - virtual int reportStepNum() const { return currentStepNum(); } + /// advance time by currentStepLength + virtual void advance() = 0 ; + + /// return copy of current timer instance + virtual std::unique_ptr< SimulatorTimerInterface > clone () const = 0; /// Current step length. This is the length of the step /// the simulator will take in the next iteration. @@ -55,6 +56,25 @@ namespace Opm /// @note if done(), it is an error to call currentStepLength(). virtual double currentStepLength() const = 0; + /// Current step number. This is the number of timesteps that + /// has been completed from the start of the run. The time + /// after initialization but before the simulation has started + /// is timestep number zero. + virtual int currentStepNum() const = 0; + + /// Return true if timer indicates that simulation of timer interval is finished + virtual bool done() const = 0; + + /// Whether the current step is the first step. + virtual bool initialStep() const = 0; + + /// Return true if last time step failed + virtual bool lastStepFailed() const = 0; + + /// Time elapsed since the start of the simulation until the + /// beginning of the current time step [s]. + virtual double simulationTimeElapsed() const = 0; + /// Previous step length. This is the length of the step that /// was taken to arrive at this time. /// @@ -63,30 +83,13 @@ namespace Opm /// it is an error to call stepLengthTaken(). virtual double stepLengthTaken () const = 0; - /// Previous report step length. This is the length of the step that - /// was taken to arrive at this report step time. - /// - /// @note if no increments have been done (i.e. the timer is - /// still in its constructed state and reportStepNum() == 0), - /// it is an error to call stepLengthTaken(). - virtual double reportStepLengthTaken () const { return stepLengthTaken(); } - - /// Time elapsed since the start of the simulation until the - /// beginning of the current time step [s]. - virtual double simulationTimeElapsed() const = 0; - - /// advance time by currentStepLength - virtual void advance() = 0 ; - - /// Return true if timer indicates that simulation of timer interval is finished - virtual bool done() const = 0; - - /// Whether the current step is the first step. - virtual bool initialStep() const = 0; - /// Return start date of simulation virtual boost::posix_time::ptime startDateTime() const = 0; + // ----------------------------------------------------------------- + // Virtual functions (not pure) to allow for default implementations + // ----------------------------------------------------------------- + /// Return the current time as a posix time object. virtual boost::posix_time::ptime currentDateTime() const; @@ -94,11 +97,17 @@ namespace Opm /// 1970) until the current time step begins [s]. virtual time_t currentPosixTime() const; - /// Return true if last time step failed - virtual bool lastStepFailed() const = 0; + /// Previous report step length. This is the length of the step that + /// was taken to arrive at this report step time. + /// + /// @note if no increments have been done (i.e. the timer is + /// still in its constructed state and reportStepNum() == 0), + /// it is an error to call stepLengthTaken(). + virtual double reportStepLengthTaken () const { return stepLengthTaken(); } + + /// Current report step number. This might differ from currentStepNum in case of sub stepping + virtual int reportStepNum() const { return currentStepNum(); } - /// return copy of current timer instance - virtual std::unique_ptr< SimulatorTimerInterface > clone () const = 0; }; diff --git a/opm/simulators/utils/readDeck.cpp b/opm/simulators/utils/readDeck.cpp index 102a9f89e..aae46fcd0 100644 --- a/opm/simulators/utils/readDeck.cpp +++ b/opm/simulators/utils/readDeck.cpp @@ -504,7 +504,6 @@ Opm::setupLogging(Parallel::Communication& comm, } logFileStream << ".PRT"; debugFileStream << ".DBG"; - FileOutputMode output; { static std::map stringToOutputMode = @@ -567,7 +566,6 @@ void Opm::readDeck(Opm::Parallel::Communication comm, const bool slaveMode) { auto errorGuard = std::make_unique(); - int parseSuccess = 1; // > 0 is success std::string failureMessage;