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.
This commit is contained in:
Bård Skaflestad committed 2025-08-15 17:47:45 +02:00
1 parent ab75376e01
commit 85c69f7d44
3 files changed
+58 -23

No files matched your search

+11 -3
View File
@@ -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<Parameters::EnableWriteAllSolutions>() ||
this->simulator().episodeWillBeOver()) {
this->episodeWillBeOver())
{
// Create VTK output as needed.
ParentType::writeOutput(verbose);
}
}
@@ -1719,6 +1720,13 @@ protected:
BCData<int> 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
+37 -10
View File
@@ -62,6 +62,7 @@
#include <algorithm>
#include <cstddef>
#include <functional>
#include <limits>
#include <memory>
#include <stdexcept>
#include <string>
@@ -457,9 +458,7 @@ public:
using GridType = std::remove_cv_t<std::remove_reference_t<decltype(grid)>>;
constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
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<Scalar> 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<float>::epsilon()) * nextReportStep;
return !isSubStep;
}
};
} // namespace Opm
+10 -10
View File
@@ -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<std::remove_reference_t<decltype(grid)>>;
constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
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<Parameters::EnableWriteAllSolutions>() || !isSubStep) {
eclWriter_->writeOutput(std::move(localCellData), isSubStep);
}
const auto isSubStep = !this->episodeWillBeOver();
if (!isSubStep || Parameters::Get<Parameters::EnableWriteAllSolutions>()) {
auto localCellData = data::Solution {};
this->eclWriter_->writeOutput(std::move(localCellData), isSubStep);
}
}