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.

CMakeLists_files.cmake

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

opm/simulators/flow/FlowMain.hpp

@@ -51,7 +51,6 @@ 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; };
-struct Slave { static constexpr bool value = false; };
 
 } // namespace Opm::Parameters
 

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,
-    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<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([[maybe_unused]]int argc, char
         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<char *> slave_argv = getSlaveArgv(
@@ -134,6 +226,8 @@ void ReservoirCouplingMaster::spawnSlaveProcesses([[maybe_unused]]int argc, char
         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 @@ class ReservoirCouplingMaster {
 
     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