Timestepping for reservoir coupling

Implement adaptive time stepping for master and slave procesess
when using reservoir coupling. The original adaptive time stepping method
is refactored at the same time.

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 {

opm/simulators/flow/ReservoirCoupling.hpp

@@ -20,15 +20,18 @@
 #ifndef OPM_RESERVOIR_COUPLING_HPP
 #define OPM_RESERVOIR_COUPLING_HPP
 #include <mpi.h>
+#include <cmath>
+#include <iostream>
 #include <memory>
 
 namespace Opm {
 namespace ReservoirCoupling {
 
 enum class MessageTag : int {
-    SimulationStartDate,
+    SlaveSimulationStartDate,
-    SimulationEndDate,
+    SlaveProcessTermination,
-    SlaveProcessTermination
+    SlaveNextReportDate,
+    SlaveNextTimeStep,
 };
 
 // Custom deleter for MPI_Comm
@@ -43,6 +46,154 @@ struct MPI_Comm_Deleter {
 
 using MPI_Comm_Ptr = std::unique_ptr<MPI_Comm, MPI_Comm_Deleter>;
 
+// 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
 

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<double>(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<double>(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(
+                &timestep,
+                /*count=*/1,
+                /*datatype=*/MPI_DOUBLE,
+                /*dest_rank=*/0,
+                /*tag=*/static_cast<int>(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<int>(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<std::time_t, long>::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<int>(MessageTag::SimulationStartDate),
+                /*tag=*/static_cast<int>(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 dir_path{directory_path};
-        std::filesystem::path data_file(data_file_name);
+        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<char *> 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_);
 }
 
 // ------------------

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<char *> 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<MPI_Comm_Ptr> master_slave_comm_; // MPI communicators for the slave processes
     std::vector<std::string> slave_names_;
     std::vector<std::time_t> slave_start_dates_;
+    std::vector<double> slave_next_report_time_offsets_;  // Elapsed time from the beginning of the simulation
 };
 
 } // namespace Opm

opm/simulators/flow/ReservoirCouplingSlave.cpp

@@ -29,15 +29,18 @@
 #include <opm/simulators/utils/ParallelCommunication.hpp>
 
 #include <vector>
+#include <fmt/format.h>
 
 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(
+            &timestep,
+            /*count=*/1,
+            /*datatype=*/MPI_DOUBLE,
+            /*source_rank=*/0,
+            /*tag=*/static_cast<int>(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(&timestep, 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<int>(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<std::time_t, long>::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<int>(MessageTag::SimulationStartDate),
+            /*tag=*/static_cast<int>(MessageTag::SlaveSimulationStartDate),
             *this->slave_master_comm_
         );
         OpmLog::info("Sent simulation start date to master process from rank 0");

opm/simulators/flow/ReservoirCouplingSlave.hpp

@@ -23,6 +23,7 @@
 #include <opm/simulators/flow/ReservoirCoupling.hpp>
 #include <opm/input/eclipse/Schedule/Schedule.hpp>
 #include <opm/simulators/utils/ParallelCommunication.hpp>
+#include <opm/simulators/timestepping/SimulatorTimer.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
 
 #include <mpi.h>
@@ -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

opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp

@@ -206,11 +206,12 @@ public:
     {
 #if HAVE_MPI
         auto slave_mode = Parameters::Get<Parameters::Slave>();
+        FlowGenericVanguard::comm().barrier();
         if (slave_mode) {
             this->reservoirCouplingSlave_ =
                 std::make_unique<ReservoirCouplingSlave>(
                     FlowGenericVanguard::comm(),
-                    this->schedule()
+                    this->schedule(), timer
                 );
             this->reservoirCouplingSlave_->sendSimulationStartDateToMasterProcess();
         }
@@ -252,7 +253,14 @@ public:
             else {
                 adaptiveTimeStepping_ = std::make_unique<TimeStepper>(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

opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp

@@ -36,25 +36,28 @@
 namespace Opm
 {
     AdaptiveSimulatorTimer::
-    AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer,
+    AdaptiveSimulatorTimer( const boost::posix_time::ptime simulation_start_time,
-                            const double lastStepTaken,
+                            const double step_length,
-                            const double maxTimeStep )
+                            const double elapsed_time,
-        : start_date_time_( std::make_shared<boost::posix_time::ptime>(timer.startDateTime()) )
+                            const double last_step_taken,
-        , start_time_( timer.simulationTimeElapsed() )
+                            const int report_step,
-        , total_time_( start_time_ + timer.currentStepLength() )
+                            const double max_time_step )
-        , report_step_( timer.reportStepNum() )
+        : start_date_time_{ std::make_shared<boost::posix_time::ptime>(simulation_start_time) }
-        , max_time_step_( maxTimeStep )
+        , start_time_{ elapsed_time }
-        , current_time_( start_time_ )
+        , total_time_{ start_time_ + step_length }
-        , dt_( 0.0 )
+        , report_step_{ report_step }
-        , current_step_( 0 )
+        , max_time_step_{ max_time_step }
-        , steps_()
+        , current_time_{ start_time_ }
-        , lastStepFailed_( false )
+        , 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

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,
+        AdaptiveSimulatorTimer( const boost::posix_time::ptime simulation_start_time,
-                                const double lastStepTaken,
+                                const double step_length,
-                                const double maxTimeStep = std::numeric_limits<double>::max() );
+                                const double elapsed_time,
+                                const double last_step_taken,
+                                const int report_step,
+                                const double max_time_step = std::numeric_limits<double>::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_;
 
     };
 

opm/simulators/timestepping/AdaptiveTimeStepping.hpp

@@ -19,12 +19,23 @@
 #include <opm/simulators/timestepping/TimeStepControl.hpp>
 #include <opm/simulators/timestepping/TimeStepControlInterface.hpp>
 
-#include <cmath>
+#if HAVE_MPI
+#define RESERVOIR_COUPLING_ENABLED
+#endif
+#ifdef RESERVOIR_COUPLING_ENABLED
+#include <opm/simulators/flow/ReservoirCoupling.hpp>
+#include <opm/simulators/flow/ReservoirCouplingMaster.hpp>
+#include <opm/simulators/flow/ReservoirCouplingSlave.hpp>
+#endif
+
 #include <functional>
 #include <memory>
 #include <set>
+#include <sstream>
+#include <stdexcept>
 #include <string>
 #include <tuple>
+#include <utility>
 #include <vector>
 
 namespace Opm::Parameters {
@@ -55,148 +66,208 @@ class UnitSystem;
 struct StepReport;
 
 namespace detail {
-
+    void logTimer(const AdaptiveSimulatorTimer& substep_timer);
-void logTimer(const AdaptiveSimulatorTimer& substepTimer);
+    std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr);
-
+    void registerAdaptiveParameters();
-std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr);
-
-void registerAdaptiveParameters();
-
-std::tuple<TimeStepControlType,
-           std::unique_ptr<TimeStepControlInterface>,
-           bool>
-createController(const UnitSystem& unitSystem);
-
 }
 
-    // AdaptiveTimeStepping
+template<class TypeTag>
-    //---------------------
+class AdaptiveTimeStepping
-    template<class TypeTag>
+{
-    class AdaptiveTimeStepping
+private:
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    template <class Solver>
+    class SolutionTimeErrorSolverWrapper : public RelativeChangeInterface
     {
-        using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-        template <class Solver>
-        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<class E>
-        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
+    // Forward declaration of SubStepIteration
-        explicit AdaptiveTimeStepping(const UnitSystem& unitSystem,
+    // TODO: This forward declaration is not enough for GCC version 9.4.0 (released June 2021),
-                                      const double max_next_tstep = -1.0,
+    //  but it works fine with GCC version 13.2.0 (released July 2023).
-                                      const bool terminalOutput = true);
+    template <class Solver> class SubStepIteration;
 
-        //! \brief contructor taking parameter object
+    template <class Solver>
-        //! \param tuning Pointer to ecl TUNING keyword
+    class SubStepper {
-        //! \param timeStep current report step
+    public:
-        AdaptiveTimeStepping(double max_next_tstep,
+        SubStepper(
-                             const Tuning& tuning,
+            AdaptiveTimeStepping<TypeTag>& adaptive_time_stepping,
-                             const UnitSystem& unitSystem,
+            const SimulatorTimer& simulator_timer,
-                             const bool terminalOutput = true);
+            Solver& solver,
+            const bool is_event,
+            const std::function<bool(const double, const double, const int)>& tuning_updater
+        );
 
-        static void registerParameters();
+        AdaptiveTimeStepping<TypeTag>& getAdaptiveTimerStepper();
-
+        SimulatorReport run();
-        /** \brief  step method that acts like the solver::step method
+        friend class AdaptiveTimeStepping<TypeTag>::template SubStepIteration<Solver>;
-                    in a sub cycle of time steps
-            \param tuningUpdater Function used to update TUNING parameters before each
-                                 time step. ACTIONX might change tuning.
-        */
-        template <class Solver>
-        SimulatorReport step(const SimulatorTimer& simulatorTimer,
-                             Solver& solver,
-                             const bool isEvent,
-                             const std::function<bool(const double, const double, const int)> 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<class Serializer>
-        void serializeOp(Serializer& serializer);
-
-        static AdaptiveTimeStepping<TypeTag> serializationTestObjectHardcoded();
-        static AdaptiveTimeStepping<TypeTag> serializationTestObjectPID();
-        static AdaptiveTimeStepping<TypeTag> serializationTestObjectPIDIt();
-        static AdaptiveTimeStepping<TypeTag> serializationTestObjectSimple();
-
-        bool operator==(const AdaptiveTimeStepping<TypeTag>& rhs) const;
 
     private:
-        template<class Controller>
+        bool isReservoirCouplingMaster_() const;
-        static AdaptiveTimeStepping<TypeTag> serializationTestObject_();
+        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<class T, class Serializer>
+        AdaptiveTimeStepping<TypeTag>& adaptive_time_stepping_;
-        void allocAndSerialize(Serializer& serializer)
+        const SimulatorTimer& simulator_timer_;
-        {
+        Solver& solver_;
-            if (!serializer.isSerializing()) {
+        const bool is_event_;
-                timeStepControl_ = std::make_unique<T>();
+        const std::function<bool(double elapsed, double dt, int sub_step_number)>& tuning_updater_;
-            }
-            serializer(*static_cast<T*>(timeStepControl_.get()));
-        }
-
-        template<class T>
-        bool castAndComp(const AdaptiveTimeStepping<TypeTag>& Rhs) const
-        {
-            const T* lhs = static_cast<const T*>(timeStepControl_.get());
-            const T* rhs = static_cast<const T*>(Rhs.timeStepControl_.get());
-            return *lhs == *rhs;
-        }
-
-    protected:
-        void init_(const UnitSystem& unitSystem);
-
-        using TimeStepController = std::unique_ptr<TimeStepControlInterface>;
-
-        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
     };
-}
+
+    template <class Solver>
+    class SubStepIteration {
+    public:
+        SubStepIteration(
+            SubStepper<Solver>& 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<Solver>& substepper_;
+        AdaptiveSimulatorTimer& substep_timer_;
+        const double original_time_step_;
+        const bool final_step_;
+        std::string cause_of_failure_;
+        AdaptiveTimeStepping<TypeTag>& 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<TypeTag>& 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 <class Solver>
+    SimulatorReport step(const SimulatorTimer& simulator_timer,
+                         Solver& solver,
+                         const bool is_event,
+                         const std::function<bool(const double, const double, const int)>
+                            tuning_updater);
+
+    void updateTUNING(double max_next_tstep, const Tuning& tuning);
+    void updateNEXTSTEP(double max_next_tstep);
+
+    template<class Serializer>
+    void serializeOp(Serializer& serializer);
+
+    static AdaptiveTimeStepping<TypeTag> serializationTestObjectHardcoded();
+    static AdaptiveTimeStepping<TypeTag> serializationTestObjectPID();
+    static AdaptiveTimeStepping<TypeTag> serializationTestObjectPIDIt();
+    static AdaptiveTimeStepping<TypeTag> serializationTestObjectSimple();
+
+private:
+    void maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step,
+                                                              const bool is_event);
+
+    template<class Controller>
+    static AdaptiveTimeStepping<TypeTag> serializationTestObject_();
+
+    template<class T, class Serializer>
+    void allocAndSerialize(Serializer& serializer);
+
+    template<class T>
+    bool castAndComp(const AdaptiveTimeStepping<TypeTag>& Rhs) const;
+
+protected:
+    void init_(const UnitSystem& unitSystem);
+
+    using TimeStepController = std::unique_ptr<TimeStepControlInterface>;
+
+    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 <opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp>
-
 #endif // OPM_ADAPTIVE_TIME_STEPPING_HPP

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
+        // Pure virtual functions to be implemented by derived classes
-        /// after initialization but before the simulation has started
+        // -----------------------------------------------------------
-        /// is timestep number zero.
-        virtual int currentStepNum() const = 0;
 
-        /// Current report step number. This might differ from currentStepNum in case of sub stepping
+        /// advance time by currentStepLength
-        virtual int reportStepNum() const { return currentStepNum(); }
+        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
+        /// Previous report step length. This is the length of the step that
-        virtual bool lastStepFailed() const = 0;
+        /// 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;
     };
 
 

opm/simulators/utils/readDeck.cpp

@@ -504,7 +504,6 @@ Opm::setupLogging(Parallel::Communication& comm,
     }
     logFileStream << ".PRT";
     debugFileStream << ".DBG";
-
     FileOutputMode output;
     {
         static std::map<std::string, FileOutputMode> stringToOutputMode =
@@ -567,7 +566,6 @@ void Opm::readDeck(Opm::Parallel::Communication    comm,
                    const bool                      slaveMode)
 {
     auto errorGuard = std::make_unique<ErrorGuard>();
-
     int parseSuccess = 1; // > 0 is success
     std::string failureMessage;