From 85c69f7d4410de784c73ec99a2ac43dfcd7ceede Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 14 Aug 2025 11:02:24 +0200 Subject: [PATCH] Add End-of-Episode Predicate to Flow Problem Classes This is arguably a hack. The default implementation in terms of Simulator<>::episodeWillBeOver() is known to have non-trivial false negatives in certain regression tests. These are in turn possibly related to using the simulation start time point (often a large number of seconds) and subtracting values of similar size. The custom implementation in FlowProblemBlackoil<> uses a floating-point comparison involving Schedule::seconds() instead, thereby reducing the size of the numbers involved. We're still vulnerable to "catastrophic cancellation" errors, but potentially less so. A complete fix for the problem may involve switching to representing the simulated time variable as an integer instead of a floating-point type. --- opm/simulators/flow/FlowProblem.hpp | 14 ++++-- opm/simulators/flow/FlowProblemBlackoil.hpp | 47 ++++++++++++++++----- opm/simulators/flow/FlowProblemComp.hpp | 20 ++++----- 3 files changed, 58 insertions(+), 23 deletions(-) diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index d1a8759647..3ac570d239 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/opm/simulators/flow/FlowProblem.hpp @@ -473,10 +473,11 @@ public: virtual void writeOutput(bool verbose) { OPM_TIMEBLOCK(problemWriteOutput); - // use the generic code to prepare the output fields and to - // write the desired VTK files. + if (Parameters::Get() || - this->simulator().episodeWillBeOver()) { + this->episodeWillBeOver()) + { + // Create VTK output as needed. ParentType::writeOutput(verbose); } } @@ -1719,6 +1720,13 @@ protected: BCData bcindex_; bool nonTrivialBoundaryConditions_ = false; bool first_step_ = true; + + /// Whether or not the current episode will end at the end of the + /// current time step. + virtual bool episodeWillBeOver() const + { + return this->simulator().episodeWillBeOver(); + } }; } // namespace Opm diff --git a/opm/simulators/flow/FlowProblemBlackoil.hpp b/opm/simulators/flow/FlowProblemBlackoil.hpp index f37aa13db1..e82e701f09 100644 --- a/opm/simulators/flow/FlowProblemBlackoil.hpp +++ b/opm/simulators/flow/FlowProblemBlackoil.hpp @@ -62,6 +62,7 @@ #include #include #include +#include #include #include #include @@ -457,9 +458,7 @@ public: using GridType = std::remove_cv_t>; constexpr bool isCpGrid = std::is_same_v; if (!isCpGrid || (grid.maxLevel() == 0)) { - const bool isSubStep = !this->simulator().episodeWillBeOver(); - - this->eclWriter_->evalSummaryState(isSubStep); + this->eclWriter_->evalSummaryState(!this->episodeWillBeOver()); } { @@ -526,22 +525,24 @@ public: * \brief Write the requested quantities of the current solution into the output * files. */ - void writeOutput(bool verbose) override + void writeOutput(const bool verbose) override { FlowProblemType::writeOutput(verbose); - bool isSubStep = !this->simulator().episodeWillBeOver(); + const auto isSubStep = !this->episodeWillBeOver(); + + auto localCellData = data::Solution {}; - data::Solution localCellData = {}; #if HAVE_DAMARIS // N.B. the Damaris output has to be done before the ECL output as the ECL one // does all kinds of std::move() relocation of data - if (enableDamarisOutput_) { - damarisWriter_->writeOutput(localCellData, isSubStep) ; + if (this->enableDamarisOutput_ && (this->damarisWriter_ != nullptr)) { + this->damarisWriter_->writeOutput(localCellData, isSubStep); } #endif - if (enableEclOutput_){ - eclWriter_->writeOutput(std::move(localCellData), isSubStep); + + if (this->enableEclOutput_ && (this->eclWriter_ != nullptr)) { + this->eclWriter_->writeOutput(std::move(localCellData), isSubStep); } } @@ -1679,6 +1680,32 @@ protected: ActionHandler actionHandler_; ModuleParams moduleParams_; + +private: + /// Whether or not the current epsiode will end at the end of the + /// current time step. + /// + /// Custom implementation for black-oil cases. + /// + /// [2025-08-14] This is arguably a hack. We intentionally do not use + /// the 'simulator()'s episodeWillBeOver() function, as that has been + /// shown to have non-trivial false negatives in the context of certain + /// regression tests. It is likely that we will need to revisit this + /// predicate in the future. + bool episodeWillBeOver() const override + { + const auto currTime = this->simulator().time() + + this->simulator().timeStepSize(); + + const auto nextReportStep = + this->simulator().vanguard().schedule() + .seconds(this->simulator().episodeIndex() + 1); + + const auto isSubStep = (nextReportStep - currTime) + > (2 * std::numeric_limits::epsilon()) * nextReportStep; + + return !isSubStep; + } }; } // namespace Opm diff --git a/opm/simulators/flow/FlowProblemComp.hpp b/opm/simulators/flow/FlowProblemComp.hpp index 796d53b042..bf5d054555 100644 --- a/opm/simulators/flow/FlowProblemComp.hpp +++ b/opm/simulators/flow/FlowProblemComp.hpp @@ -256,8 +256,6 @@ public: { FlowProblemType::endTimeStep(); - const bool isSubStep = !this->simulator().episodeWillBeOver(); - // after the solution is updated, the values in output module also needs to be updated this->eclWriter_->mutableOutputModule().invalidateLocalData(); @@ -267,9 +265,8 @@ public: using GridType = std::remove_cv_t>; constexpr bool isCpGrid = std::is_same_v; if (!isCpGrid || (grid.maxLevel() == 0)) { - this->eclWriter_->evalSummaryState(isSubStep); + this->eclWriter_->evalSummaryState(! this->episodeWillBeOver()); } - } void writeReports(const SimulatorTimer& timer) { @@ -286,13 +283,16 @@ public: { FlowProblemType::writeOutput(verbose); - const bool isSubStep = !this->simulator().episodeWillBeOver(); + if (! this->enableEclOutput_) { + return; + } - data::Solution localCellData = {}; - if (enableEclOutput_) { - if (Parameters::Get() || !isSubStep) { - eclWriter_->writeOutput(std::move(localCellData), isSubStep); - } + const auto isSubStep = !this->episodeWillBeOver(); + + if (!isSubStep || Parameters::Get()) { + auto localCellData = data::Solution {}; + + this->eclWriter_->writeOutput(std::move(localCellData), isSubStep); } }