diff --git a/opm/simulators/flow/ReservoirCoupling.cpp b/opm/simulators/flow/ReservoirCoupling.cpp index 69c4e3e62..1e25c7b59 100644 --- a/opm/simulators/flow/ReservoirCoupling.cpp +++ b/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 diff --git a/opm/simulators/flow/ReservoirCoupling.hpp b/opm/simulators/flow/ReservoirCoupling.hpp index 0eb312cbd..0d69c6a3a 100644 --- a/opm/simulators/flow/ReservoirCoupling.hpp +++ b/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 diff --git a/opm/simulators/flow/ReservoirCouplingMaster.cpp b/opm/simulators/flow/ReservoirCouplingMaster.cpp index 1bb618740..93e230312 100644 --- a/opm/simulators/flow/ReservoirCouplingMaster.cpp +++ b/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 diff --git a/opm/simulators/flow/ReservoirCouplingMaster.hpp b/opm/simulators/flow/ReservoirCouplingMaster.hpp index 9f11a1ed9..a50906863 100644 --- a/opm/simulators/flow/ReservoirCouplingMaster.hpp +++ b/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,