Eliminate TimePoint class

Make fuzzy comparison between two dates explicit.

opm/simulators/flow/ReservoirCoupling.cpp

@@ -74,5 +74,35 @@ void setErrhandler(MPI_Comm comm, bool is_master)
     MPI_Comm_set_errhandler(comm, errhandler);
 }
 
+bool Seconds::compare_eq(double a, double b)
+{
+    // Are a and b equal?
+    return std::abs(a - b) < std::max(abstol, reltol * std::max(std::abs(a), std::abs(b)));
+}
+
+bool Seconds::compare_gt_or_eq(double a, double b)
+{
+    // Is a greater than or equal to b?
+    if (compare_eq(a, b)) {
+        return true;
+    }
+    return a > b;
+}
+
+bool Seconds::compare_gt(double a, double b)
+{
+    // Is a greater than b?
+    return !compare_eq(a, b) && a > b;
+}
+
+bool Seconds::compare_lt_or_eq(double a, double b)
+{
+    // Is a less than or equal to b?
+    if (compare_eq(a, b)) {
+        return true;
+    }
+    return a < b;
+}
+
 } // namespace ReservoirCoupling
 } // namespace Opm

opm/simulators/flow/ReservoirCoupling.hpp

@@ -37,158 +37,39 @@ enum class MessageTag : int {
     MasterGroupNamesSize,
 };
 
-// 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;
-    }
-};
-
 // Helper functions
 void custom_error_handler_(MPI_Comm* comm, int* err, const std::string &msg);
 void setErrhandler(MPI_Comm comm, bool is_master);
 
+// Utility class for comparing double values representing epoch dates (seconds since
+// unix epoch) or elapsed time (seconds since the start of the simulation).
+// NOTE: It is important that when comparing against start of a report step or similar, that
+//   that we do not miss these due to numerical issues. This is because communication between
+//   master and slave processes are based on these points in time.
+// NOTE: Epoch values in this century (2000-2100) lies in the range of [1e9,4e9], and a double variable cannot
+//  represent such large values with high precision. For example, the date 01-01-2020 is equal
+//  to 1.5778368e9 seconds and adding 1e-7 seconds to this value will not change the value.
+//  So microseconds (1e-6) is approximately the smallest time unit we can represent for such a number.
+// NOTE: Report steps seems to have a maximum resolution of whole seconds, see stepLength() in
+//   Schedule.cpp in opm-common, which returns the step length in seconds.
+struct Seconds {
+    static constexpr double abstol = 1e-15;
+    static constexpr double reltol = 1e-15;
+    // We will will use the following expression to determine if two values a and b are equal:
+    //   |a - b| <= tol = abstol + reltol * max(|a|, |b|)
+    // For example, assume abstol = reltol = 1e-15, then the following holds:
+    // - If |a| and |b| are below 1, then the absolute tolerance applies.
+    // - If a and b are above 1, then the relative tolerance applies.
+    // For example, for dates in the range 01-01-2000 to 01-01-2100, epoch values will be in the range
+    //  [1e9, 4e9]. And we have 1e-15 * 1e9 = 1e-6, so numbers differing below one microsecond will
+    // be considered equal.
+    // NOTE: The above is not true for numbers close to zero, but we do not expect to compare such numbers.
+    static bool compare_eq(double a, double b);
+    static bool compare_gt(double a, double b);
+    static bool compare_gt_or_eq(double a, double b);
+    static bool compare_lt_or_eq(double a, double b);
+};
+
 } // namespace ReservoirCoupling
 } // namespace Opm
 

opm/simulators/flow/ReservoirCouplingMaster.cpp

@@ -81,21 +81,21 @@ 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 = this->schedule_.getStartTime();
-    TimePoint step_start_date{start_date + elapsed_time};
-    TimePoint step_end_date{step_start_date + suggested_timestep_original};
-    TimePoint suggested_timestep{suggested_timestep_original};
+    double step_start_date{start_date + elapsed_time};
+    double step_end_date{step_start_date + suggested_timestep_original};
+    double suggested_timestep{suggested_timestep_original};
     auto num_slaves = this->numSlavesStarted();
     for (std::size_t i = 0; i < num_slaves; i++) {
         double slave_start_date = 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) {
+        double slave_next_report_date{this->slave_next_report_time_offsets_[i] + slave_start_date};
+        if (Seconds::compare_gt_or_eq(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) {
+        double slave_elapsed_time;
+        if (Seconds::compare_lt_or_eq(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) {
+            if (Seconds::compare_gt(slave_next_report_date, step_end_date)) {
                 // The slave process will not report during this timestep
                 continue;
             }
@@ -109,7 +109,7 @@ maybeChopSubStep(double suggested_timestep_original, double elapsed_time) const
         suggested_timestep = slave_elapsed_time;
         step_end_date = step_start_date + suggested_timestep;
     }
-    return suggested_timestep.getTime();
+    return suggested_timestep;
 }
 
 void

opm/simulators/flow/ReservoirCouplingMaster.hpp

@@ -35,7 +35,7 @@ namespace Opm {
 class ReservoirCouplingMaster {
 public:
     using MessageTag = ReservoirCoupling::MessageTag;
-    using TimePoint = ReservoirCoupling::TimePoint;
+    using Seconds = ReservoirCoupling::Seconds;
     ReservoirCouplingMaster(
         const Parallel::Communication &comm,
         const Schedule &schedule,