diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index cda462e0fa5..4d8cd0c39da 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -1172,7 +1172,20 @@ if(dune-alugrid_FOUND)
examples/fracture_discretefracture.cpp
)
endif()
-
+if(USE_MPI)
+ list (APPEND MAIN_SOURCE_FILES
+ opm/simulators/flow/ReservoirCoupling.cpp
+ opm/simulators/flow/ReservoirCouplingMaster.cpp
+ opm/simulators/flow/ReservoirCouplingSlave.cpp
+ opm/simulators/flow/ReservoirCouplingSpawnSlaves.cpp
+ )
+ list (APPEND PUBLIC_HEADER_FILES
+ opm/simulators/flow/ReservoirCoupling.hpp
+ opm/simulators/flow/ReservoirCouplingMaster.hpp
+ opm/simulators/flow/ReservoirCouplingSlave.hpp
+ opm/simulators/flow/ReservoirCouplingSpawnSlaves.hpp
+ )
+endif()
if(HYPRE_FOUND)
list(APPEND PUBLIC_HEADER_FILES
opm/simulators/linalg/HyprePreconditioner.hpp
diff --git a/opm/simulators/flow/FlowGenericVanguard.cpp b/opm/simulators/flow/FlowGenericVanguard.cpp
index 1622411b98f..1f6e0af49e5 100644
--- a/opm/simulators/flow/FlowGenericVanguard.cpp
+++ b/opm/simulators/flow/FlowGenericVanguard.cpp
@@ -177,7 +177,7 @@ void FlowGenericVanguard::readDeck(const std::string& filename)
modelParams_.actionState_,
modelParams_.wtestState_,
modelParams_.eclSummaryConfig_,
- nullptr, "normal", "normal", "100", false, false, false, {});
+ nullptr, "normal", "normal", "100", false, false, false, {}, /*slaveMode=*/false);
modelParams_.setupTime_ = setupTimer.stop();
}
diff --git a/opm/simulators/flow/FlowMain.hpp b/opm/simulators/flow/FlowMain.hpp
index 24a888c637a..610b735e5fe 100644
--- a/opm/simulators/flow/FlowMain.hpp
+++ b/opm/simulators/flow/FlowMain.hpp
@@ -32,6 +32,9 @@
#include <opm/simulators/flow/FlowUtils.hpp>
#include <opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp>
+#if HAVE_MPI
+#define RESERVOIR_COUPLING_ENABLED
+#endif
#if HAVE_DUNE_FEM
#include <dune/fem/misc/mpimanager.hh>
#else
@@ -50,7 +53,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; };
} // namespace Opm::Parameters
@@ -366,7 +368,11 @@ namespace Opm {
// Callback that will be called from runSimulatorInitOrRun_().
int runSimulatorRunCallback_()
{
+#ifdef RESERVOIR_COUPLING_ENABLED
+ SimulatorReport report = simulator_->run(*simtimer_, this->argc_, this->argv_);
+#else
SimulatorReport report = simulator_->run(*simtimer_);
+#endif
runSimulatorAfterSim_(report);
return report.success.exit_status;
}
@@ -374,7 +380,11 @@ namespace Opm {
// Callback that will be called from runSimulatorInitOrRun_().
int runSimulatorInitCallback_()
{
+#ifdef RESERVOIR_COUPLING_ENABLED
+ simulator_->init(*simtimer_, this->argc_, this->argv_);
+#else
simulator_->init(*simtimer_);
+#endif
return EXIT_SUCCESS;
}
diff --git a/opm/simulators/flow/Main.cpp b/opm/simulators/flow/Main.cpp
index 3abc4c07da7..9ce461cce0e 100644
--- a/opm/simulators/flow/Main.cpp
+++ b/opm/simulators/flow/Main.cpp
@@ -44,11 +44,21 @@
#include <amgx_c.h>
#endif
+#include <iostream>
+// NOTE: There is no C++ header replacement for these C posix headers (as of C++17)
+#include <fcntl.h> // for open()
+#include <unistd.h> // for dup2(), close()
+
+#include <iostream>
+
namespace Opm {
Main::Main(int argc, char** argv, bool ownMPI)
: argc_(argc), argv_(argv), ownMPI_(ownMPI)
{
+#if HAVE_MPI
+ maybeSaveReservoirCouplingSlaveLogFilename_();
+#endif
if (ownMPI_) {
initMPI();
}
@@ -132,6 +142,53 @@ Main::~Main()
#endif
}
+#if HAVE_MPI
+void Main::maybeSaveReservoirCouplingSlaveLogFilename_()
+{
+ // If first command line argument is "--slave-log-file=<filename>",
+ // then redirect stdout and stderr to the specified file.
+ if (this->argc_ >= 2) {
+ std::string_view arg = this->argv_[1];
+ if (arg.substr(0, 17) == "--slave-log-file=") {
+ // For now, just save the basename of the filename and we will append the rank
+ // to the filename after having called MPI_Init() to avoid all ranks outputting
+ // to the same file.
+ this->reservoirCouplingSlaveOutputFilename_ = arg.substr(17);
+ this->argc_ -= 1;
+ char *program_name = this->argv_[0];
+ this->argv_ += 1;
+ // We assume the "argv" array pointers remain valid (not freed) for the lifetime
+ // of this program, so the following assignment is safe.
+ // Overwrite the first argument with the program name, i.e. pretend the program
+ // was called with the same arguments, but without the "--slave-log-file" argument.
+ this->argv_[0] = program_name;
+ }
+ }
+}
+#endif
+#if HAVE_MPI
+void Main::maybeRedirectReservoirCouplingSlaveOutput_() {
+ if (!this->reservoirCouplingSlaveOutputFilename_.empty()) {
+ std::string filename = this->reservoirCouplingSlaveOutputFilename_
+ + "." + std::to_string(FlowGenericVanguard::comm().rank()) + ".log";
+ int fd = open(filename.c_str(), O_WRONLY | O_CREAT | O_TRUNC, 0644);
+ if (fd == -1) {
+ std::string error_msg = "Slave: Failed to open stdout+stderr file" + filename;
+ perror(error_msg.c_str());
+ MPI_Abort(MPI_COMM_WORLD, /*status=*/1);
+ }
+ // Redirect stdout and stderr to the file.
+ if (dup2(fd, fileno(stdout)) == -1 || dup2(fileno(stdout), fileno(stderr)) == -1) {
+ std::string error_msg = "Slave: Failed to redirect stdout+stderr to " + filename;
+ perror(error_msg.c_str());
+ close(fd);
+ MPI_Abort(MPI_COMM_WORLD, /*status=*/1);
+ }
+ close(fd);
+ }
+}
+#endif
+
void Main::setArgvArgc_(const std::string& filename)
{
this->deckFilename_ = filename;
@@ -166,6 +223,7 @@ void Main::initMPI()
handleTestSplitCommunicatorCmdLine_();
#if HAVE_MPI
+ maybeRedirectReservoirCouplingSlaveOutput_();
if (test_split_comm_ && FlowGenericVanguard::comm().size() > 1) {
int world_rank = FlowGenericVanguard::comm().rank();
int color = (world_rank == 0);
@@ -230,6 +288,7 @@ void Main::readDeck(const std::string& deckFilename,
const bool keepKeywords,
const std::size_t numThreads,
const int output_param,
+ const bool slaveMode,
const std::string& parameters,
std::string_view moduleVersion,
std::string_view compileTimestamp)
@@ -265,7 +324,8 @@ void Main::readDeck(const std::string& deckFilename,
init_from_restart_file,
outputCout_,
keepKeywords,
- outputInterval);
+ outputInterval,
+ slaveMode);
verifyValidCellGeometry(FlowGenericVanguard::comm(), *this->eclipseState_);
@@ -294,7 +354,7 @@ void Main::setupDamaris(const std::string& outputDir )
//const auto find_replace_map = Opm::DamarisOutput::DamarisKeywords<PreTypeTag>(EclGenericVanguard::comm(), outputDir);
std::map<std::string, std::string> find_replace_map;
find_replace_map = DamarisOutput::getDamarisKeywords(FlowGenericVanguard::comm(), outputDir);
-
+
// By default EnableDamarisOutputCollective is true so all simulation results will
// be written into one single file for each iteration using Parallel HDF5.
// If set to false, FilePerCore mode is used in Damaris, then simulation results in each
diff --git a/opm/simulators/flow/Main.hpp b/opm/simulators/flow/Main.hpp
index 9a4e0f47c60..6a9e1a7b8e8 100644
--- a/opm/simulators/flow/Main.hpp
+++ b/opm/simulators/flow/Main.hpp
@@ -150,7 +150,8 @@ class Main
~Main();
void setArgvArgc_(const std::string& filename);
-
+ void maybeSaveReservoirCouplingSlaveLogFilename_();
+ void maybeRedirectReservoirCouplingSlaveOutput_();
void initMPI();
int runDynamic()
@@ -413,6 +414,7 @@ class Main
keepKeywords,
getNumThreads(),
Parameters::Get<Parameters::EclOutputInterval>(),
+ Parameters::Get<Parameters::Slave>(),
cmdline_params,
Opm::moduleVersion(),
Opm::compileTimestamp());
@@ -697,6 +699,7 @@ class Main
const bool keepKeywords,
const std::size_t numThreads,
const int output_param,
+ const bool slaveMode,
const std::string& parameters,
std::string_view moduleVersion,
std::string_view compileTimestamp);
@@ -740,6 +743,9 @@ class Main
// To demonstrate run with non_world_comm
bool test_split_comm_ = false;
bool isSimulationRank_ = true;
+#if HAVE_MPI
+ std::string reservoirCouplingSlaveOutputFilename_{};
+#endif
#if HAVE_DAMARIS
bool enableDamarisOutput_ = false;
#endif
diff --git a/opm/simulators/flow/ReservoirCoupling.cpp b/opm/simulators/flow/ReservoirCoupling.cpp
new file mode 100644
index 00000000000..8c4499a87ba
--- /dev/null
+++ b/opm/simulators/flow/ReservoirCoupling.cpp
@@ -0,0 +1,116 @@
+/*
+ Copyright 2024 Equinor ASA
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#include <opm/simulators/flow/ReservoirCoupling.hpp>
+#include <opm/simulators/utils/ParallelCommunication.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+
+#include <fmt/format.h>
+
+namespace Opm {
+namespace ReservoirCoupling {
+
+void custom_error_handler_(MPI_Comm* comm, int* err, const std::string &msg)
+{
+ // It can be useful to have a custom error handler for debugging purposes.
+ // If not, MPI will use the default error handler which aborts the program, and
+ // it can be difficult to determine exactly where the error occurred. With a custom
+ // error handler you can at least set a breakpoint here to get a backtrace.
+ int rank;
+ MPI_Comm_rank(*comm, &rank);
+ char err_string[MPI_MAX_ERROR_STRING];
+ int len;
+ MPI_Error_string(*err, err_string, &len);
+ const std::string error_msg = fmt::format(
+ "Reservoir coupling MPI error handler {} rank {}: {}", msg, rank, err_string);
+ // NOTE: The output to terminal vie stderr or stdout has been redirected to files for
+ // the slaves, see Main.cpp. So the following output will not be visible in the terminal
+ // if we are called from a slave process.
+ // std::cerr << error_msg << std::endl;
+ OpmLog::error(error_msg); // Output to log file, also shows the message on stdout for master.
+ MPI_Abort(*comm, *err);
+}
+
+void custom_error_handler_slave_(MPI_Comm* comm, int* err, ...)
+{
+ custom_error_handler_(comm, err, "slave");
+}
+
+void custom_error_handler_master_(MPI_Comm* comm, int* err, ...)
+{
+ custom_error_handler_(comm, err, "master");
+}
+
+void setErrhandler(MPI_Comm comm, bool is_master)
+{
+ MPI_Errhandler errhandler;
+ // NOTE: Lambdas with captures cannot be used as C function pointers, also
+ // converting lambdas that use ellipsis "..." as last argument to a C function pointer
+ // is not currently possible, so we need to use static functions instead.
+ if (is_master) {
+ MPI_Comm_create_errhandler(custom_error_handler_master_, &errhandler);
+ }
+ else {
+ MPI_Comm_create_errhandler(custom_error_handler_slave_, &errhandler);
+ }
+ // NOTE: The errhandler is a handle (an integer) that is associated with the communicator
+ // that is why we pass this by value below. And it is safe to free the errhandler after it has
+ // been associated with the communicator.
+ MPI_Comm_set_errhandler(comm, errhandler);
+ // Mark the error handler for deletion. According to the documentation: "The error handler will
+ // be deallocated after all the objects associated with it (communicator, window, or file) have
+ // been deallocated." So the error handler will still be in effect until the communicator is
+ // deallocated.
+ MPI_Errhandler_free(&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
new file mode 100644
index 00000000000..2308d4ce454
--- /dev/null
+++ b/opm/simulators/flow/ReservoirCoupling.hpp
@@ -0,0 +1,113 @@
+/*
+ Copyright 2024 Equinor ASA
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#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 {
+ SlaveSimulationStartDate,
+ SlaveActivationDate,
+ SlaveProcessTermination,
+ SlaveNextReportDate,
+ SlaveNextTimeStep,
+ MasterGroupNames,
+ MasterGroupNamesSize,
+};
+
+// Helper functions
+void custom_error_handler_(MPI_Comm* comm, int* err, const std::string &msg);
+void setErrhandler(MPI_Comm comm, bool is_master);
+
+/// \brief Utility class for comparing double values representing epoch dates or elapsed time.
+///
+/// This class is used to compare double values that represent:
+/// - Epoch dates (seconds since the Unix epoch)
+/// - Elapsed time (seconds since the start of the simulation)
+///
+/// \note When comparing against the start of a report step or similar, it is important not to miss
+/// these points due to numerical issues. Communication between master and slave processes is based
+/// on these specific points in time.
+///
+/// \note Epoch values in this century (2000-2100) fall within the range [1e9, 4e9]. A double variable
+/// cannot represent such large values with high precision. For example, the date 01-01-2020 corresponds
+/// to 1.5778368e9 seconds, and adding 1e-7 seconds to this value does not change it. Microseconds (1e-6)
+/// are approximately the smallest time units that can be represented accurately for such numbers.
+///
+/// \note Report steps appear to have a maximum resolution of whole seconds. See the `stepLength()`
+/// function in `Schedule.cpp` in the `opm-common` module, which returns the step length in seconds.
+struct Seconds {
+ /// \brief Absolute tolerance used for comparisons.
+ static constexpr double abstol = 1e-15;
+
+ /// \brief Relative tolerance used for comparisons.
+ static constexpr double reltol = 1e-15;
+
+ /// \brief Determines if two double values are equal within a specified tolerance.
+ ///
+ /// Two values \a a and \a b are considered equal if:
+ /// \f[ |a - b| \leq \text{tol} = \text{abstol} + \text{reltol} \times \max(|a|, |b|) \f]
+ ///
+ /// For example, if \a abstol = \a reltol = 1e-15:
+ /// - If \a |a| and \a |b| are below 1, the absolute tolerance applies.
+ /// - If \a a and \a b are above 1, the relative tolerance applies.
+ ///
+ /// \note For epoch dates between 01-01-2000 and 01-01-2100, epoch values are in the range [1e9, 4e9].
+ /// Therefore, \f$ 10^{-15} \times 10^9 = 10^{-6} \f$, meaning differences smaller than one microsecond
+ /// will be considered equal.
+ ///
+ /// \note This approach is not accurate for numbers close to zero, but such comparisons are not expected.
+ ///
+ /// \param a First double value.
+ /// \param b Second double value.
+ /// \return True if \a a and \a b are considered equal within the tolerance.
+ static bool compare_eq(double a, double b);
+
+ /// \brief Determines if \a a is greater than \a b within the specified tolerance.
+ ///
+ /// \param a First double value.
+ /// \param b Second double value.
+ /// \return True if \a a is greater than \a b.
+ static bool compare_gt(double a, double b);
+
+ /// \brief Determines if \a a is greater than \a b within the specified tolerance.
+ ///
+ /// \param a First double value.
+ /// \param b Second double value.
+ /// \return True if \a a is greater than \a b.
+ static bool compare_gt_or_eq(double a, double b);
+
+ /// \brief Determines if \a a is less than or equal to \a b within the specified tolerance.
+ ///
+ /// \param a First double value.
+ /// \param b Second double value.
+ /// \return True if \a a is less than or equal to \a b.
+ static bool compare_lt_or_eq(double a, double b);
+};
+
+} // namespace ReservoirCoupling
+} // namespace Opm
+
+#endif // OPM_RESERVOIR_COUPLING_HPP
diff --git a/opm/simulators/flow/ReservoirCouplingMaster.cpp b/opm/simulators/flow/ReservoirCouplingMaster.cpp
new file mode 100644
index 00000000000..c0e82fe3279
--- /dev/null
+++ b/opm/simulators/flow/ReservoirCouplingMaster.cpp
@@ -0,0 +1,217 @@
+/*
+ Copyright 2024 Equinor ASA
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#include <opm/simulators/flow/ReservoirCouplingMaster.hpp>
+#include <opm/simulators/flow/ReservoirCouplingSpawnSlaves.hpp>
+
+#include <opm/input/eclipse/Schedule/ResCoup/ReservoirCouplingInfo.hpp>
+#include <opm/input/eclipse/Schedule/ResCoup/MasterGroup.hpp>
+#include <opm/input/eclipse/Schedule/ResCoup/Slaves.hpp>
+#include <opm/common/ErrorMacros.hpp>
+
+#include <opm/simulators/utils/ParallelCommunication.hpp>
+
+#include <opm/input/eclipse/Schedule/Schedule.hpp>
+
+
+#include <filesystem>
+#include <vector>
+
+#include <fmt/format.h>
+
+namespace Opm {
+
+ReservoirCouplingMaster::
+ReservoirCouplingMaster(
+ const Parallel::Communication &comm,
+ const Schedule &schedule,
+ int argc, char **argv
+) :
+ comm_{comm},
+ schedule_{schedule},
+ argc_{argc},
+ argv_{argv}
+{
+ this->activation_date_ = this->getMasterActivationDate_();
+}
+
+// ------------------
+// Public methods
+// ------------------
+
+void
+ReservoirCouplingMaster::
+maybeSpawnSlaveProcesses(int report_step)
+{
+ if (this->numSlavesStarted() > 0) { // We have already spawned the slave processes
+ return;
+ }
+ const auto& rescoup = this->schedule_[report_step].rescoup();
+ auto slave_count = rescoup.slaveCount();
+ auto master_group_count = rescoup.masterGroupCount();
+ if (slave_count > 0 && master_group_count > 0) {
+ ReservoirCouplingSpawnSlaves spawn_slaves{*this, rescoup, report_step};
+ spawn_slaves.spawn();
+ }
+}
+
+
+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.
+ // NOTE: getStartTime() returns a std::time_t value, which is typically a long integer. It should
+ // be possible to represent reasonable epoch values within a double. See comment for
+ // getMasterActivationDate_() for more information.
+ double start_date = this->schedule_.getStartTime();
+ 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();
+ // Determine the minimum step_end_date and the corresponding suggested_timestep such that no
+ // slave process will report or start during the timestep [step_start_date, step_end_date]
+ // where suggested_timestep = step_end_date - step_start_date
+ for (std::size_t i = 0; i < num_slaves; i++) {
+ double slave_start_date = this->slave_start_dates_[i];
+ 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;
+ }
+ 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 (Seconds::compare_gt(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;
+}
+
+void
+ReservoirCouplingMaster::
+sendNextTimeStepToSlaves(double timestep)
+{
+ if (this->comm_.rank() == 0) {
+ for (unsigned int i = 0; i < this->master_slave_comm_.size(); i++) {
+ MPI_Send(
+ ×tep,
+ /*count=*/1,
+ /*datatype=*/MPI_DOUBLE,
+ /*dest_rank=*/0,
+ /*tag=*/static_cast<int>(MessageTag::SlaveNextTimeStep),
+ this->getSlaveComm(i)
+ );
+ 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()
+{
+ auto num_slaves = this->numSlavesStarted();
+ OpmLog::info("Receiving next report dates from slave processes");
+ if (this->comm_.rank() == 0) {
+ for (unsigned int i = 0; i < num_slaves; i++) {
+ double slave_next_report_time_offset; // Elapsed time from the beginning of the simulation
+ // NOTE: All slave-master communicators have set a custom error handler, which eventually
+ // will call MPI_Abort() so there is no need to check the return value of any MPI_Recv()
+ // or MPI_Send() calls.
+ MPI_Recv(
+ &slave_next_report_time_offset,
+ /*count=*/1,
+ /*datatype=*/MPI_DOUBLE,
+ /*source_rank=*/0,
+ /*tag=*/static_cast<int>(MessageTag::SlaveNextReportDate),
+ this->getSlaveComm(i),
+ MPI_STATUS_IGNORE
+ );
+ 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(), /*count=*/num_slaves, /*emitter_rank=*/0
+ );
+ OpmLog::info("Broadcasted slave next report dates to all ranks");
+}
+
+
+std::size_t
+ReservoirCouplingMaster::
+numSlavesStarted() const
+{
+ return this->slave_names_.size();
+}
+
+// ------------------
+// Private methods
+// ------------------
+
+double
+ReservoirCouplingMaster::
+getMasterActivationDate_() const
+{
+ // Assume master mode is activated when the first SLAVES keyword is encountered in the schedule
+ // NOTE: getStartTime() returns a std::time_t value, which is typically a long integer representing
+ // the number of seconds since the epoch (1970-01-01 00:00:00 UTC)
+ // The maximum integer that can be represented by a double is 2^53 - 1, which is approximately
+ // 9e15. This corresponds to a date in the year 2.85e8 or 285 million years into the future.
+ // So we should be able to represent reasonable epoch values within a double.
+ double start_date = this->schedule_.getStartTime();
+ for (std::size_t report_step = 0; report_step < this->schedule_.size(); ++report_step) {
+ auto rescoup = this->schedule_[report_step].rescoup();
+ if (rescoup.slaveCount() > 0) {
+ return start_date + this->schedule_.seconds(report_step);
+ }
+ }
+ // NOTE: Consistency between SLAVES and GRUPMAST keywords has already been checked in
+ // init() in SimulatorFullyImplicitBlackoil.hpp
+ OPM_THROW(std::runtime_error, "Reservoir coupling: Failed to find master activation time: "
+ "No SLAVES keyword found in schedule");
+}
+
+
+} // namespace Opm
+
diff --git a/opm/simulators/flow/ReservoirCouplingMaster.hpp b/opm/simulators/flow/ReservoirCouplingMaster.hpp
new file mode 100644
index 00000000000..c8de9316ffc
--- /dev/null
+++ b/opm/simulators/flow/ReservoirCouplingMaster.hpp
@@ -0,0 +1,91 @@
+/*
+ Copyright 2024 Equinor ASA
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_RESERVOIR_COUPLING_MASTER_HPP
+#define OPM_RESERVOIR_COUPLING_MASTER_HPP
+
+#include <opm/simulators/utils/ParallelCommunication.hpp>
+#include <opm/simulators/flow/ReservoirCoupling.hpp>
+#include <opm/input/eclipse/Schedule/Schedule.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+
+#include <mpi.h>
+
+#include <filesystem>
+#include <vector>
+
+namespace Opm {
+
+class ReservoirCouplingMaster {
+public:
+ using MessageTag = ReservoirCoupling::MessageTag;
+ using Seconds = ReservoirCoupling::Seconds;
+ ReservoirCouplingMaster(
+ const Parallel::Communication &comm,
+ const Schedule &schedule,
+ int argc, char **argv
+ );
+
+ bool activated() { return this->numSlavesStarted() > 0; }
+ void addSlaveCommunicator(MPI_Comm comm) {
+ this->master_slave_comm_.push_back(comm);
+ }
+ void addSlaveName(const std::string &name) { this->slave_names_.push_back(name); }
+ void addSlaveStartDate(std::time_t date) { this->slave_start_dates_.push_back(date); }
+ double getActivationDate() const { return this->activation_date_; }
+ int getArgc() const { return this->argc_; }
+ char *getArgv(int index) const { return this->argv_[index]; }
+ char **getArgv() const { return this->argv_; }
+ const Parallel::Communication &getComm() const { return this->comm_; }
+ double getSimulationStartDate() const { return this->schedule_.getStartTime(); }
+ MPI_Comm getSlaveComm(int index) const { return this->master_slave_comm_[index]; }
+ const std::string &getSlaveName(int index) const { return this->slave_names_[index]; }
+ const double *getSlaveStartDates() { return this->slave_start_dates_.data(); }
+ double maybeChopSubStep(double suggested_timestep, double current_time) const;
+ void maybeSpawnSlaveProcesses(int report_step);
+ std::size_t numSlavesStarted() const;
+ void receiveNextReportDateFromSlaves();
+ void resizeSlaveStartDates(int size) { this->slave_start_dates_.resize(size); }
+ void resizeNextReportDates(int size) { this->slave_next_report_time_offsets_.resize(size); }
+ void sendNextTimeStepToSlaves(double next_time_step);
+
+private:
+ double getMasterActivationDate_() const;
+
+ const Parallel::Communication &comm_;
+ const Schedule& schedule_;
+ int argc_;
+ char **argv_;
+ // NOTE: MPI_Comm is just an integer handle, so we can just copy it into the vector
+ std::vector<MPI_Comm> master_slave_comm_; // MPI communicators for the slave processes
+ std::vector<std::string> slave_names_;
+ // The start dates are in whole seconds since the epoch. We use a double to store the value
+ // since both schedule_.getStartTime() and schedule_.stepLength(report_step) returns
+ // a double value representing whole seconds.
+ // However, note that schedule_[report_step].start_time() returns a time_point
+ // which can include milliseconds. The double values are also convenient when we need to
+ // to add fractions of seconds for sub steps to the start date.
+ std::vector<double> slave_start_dates_;
+ // Elapsed time from the beginning of the simulation
+ std::vector<double> slave_next_report_time_offsets_;
+ double activation_date_{0.0}; // The date when SLAVES is encountered in the schedule
+};
+
+} // namespace Opm
+#endif // OPM_RESERVOIR_COUPLING_MASTER_HPP
diff --git a/opm/simulators/flow/ReservoirCouplingSlave.cpp b/opm/simulators/flow/ReservoirCouplingSlave.cpp
new file mode 100644
index 00000000000..eb16256af37
--- /dev/null
+++ b/opm/simulators/flow/ReservoirCouplingSlave.cpp
@@ -0,0 +1,276 @@
+/*
+ Copyright 2024 Equinor ASA
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+#include <opm/simulators/flow/ReservoirCoupling.hpp>
+#include <opm/simulators/flow/ReservoirCouplingSlave.hpp>
+
+#include <opm/input/eclipse/Schedule/ResCoup/ReservoirCouplingInfo.hpp>
+#include <opm/input/eclipse/Schedule/ResCoup/MasterGroup.hpp>
+#include <opm/input/eclipse/Schedule/ResCoup/Slaves.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/simulators/utils/ParallelCommunication.hpp>
+
+#include <dune/common/parallel/mpitraits.hh>
+
+#include <vector>
+#include <fmt/format.h>
+
+namespace Opm {
+
+ReservoirCouplingSlave::
+ReservoirCouplingSlave(
+ const Parallel::Communication &comm,
+ const Schedule &schedule,
+ const SimulatorTimer &timer
+) :
+ comm_{comm},
+ schedule_{schedule},
+ timer_{timer}
+{
+ this->slave_master_comm_ = MPI_COMM_NULL;
+ MPI_Comm_get_parent(&this->slave_master_comm_);
+ if (this->slave_master_comm_ == MPI_COMM_NULL) {
+ OPM_THROW(std::runtime_error, "Slave process is not spawned by a master process");
+ }
+ // NOTE: By installing a custom error handler for all slave-master communicators, which
+ // eventually will call MPI_Abort(), there is no need to check the return value of any
+ // MPI_Recv() or MPI_Send() calls as errors will be caught by the error handler.
+ ReservoirCoupling::setErrhandler(this->slave_master_comm_, /*is_master=*/false);
+}
+
+double
+ReservoirCouplingSlave::
+receiveNextTimeStepFromMaster() {
+ double timestep;
+ if (this->comm_.rank() == 0) {
+ // NOTE: All slave-master communicators have set a custom error handler, which eventually
+ // will call MPI_Abort() so there is no need to check the return value of any MPI_Recv()
+ // or MPI_Send() calls.
+ MPI_Recv(
+ ×tep,
+ /*count=*/1,
+ /*datatype=*/MPI_DOUBLE,
+ /*source_rank=*/0,
+ /*tag=*/static_cast<int>(MessageTag::SlaveNextTimeStep),
+ this->slave_master_comm_,
+ MPI_STATUS_IGNORE
+ );
+ OpmLog::info(
+ fmt::format("Slave rank 0 received next timestep {} from master.", timestep)
+ );
+ }
+ this->comm_.broadcast(×tep, /*count=*/1, /*emitter_rank=*/0);
+ OpmLog::info("Broadcasted slave next time step to all ranks");
+ return timestep;
+}
+
+void
+ReservoirCouplingSlave::
+receiveMasterGroupNamesFromMasterProcess() {
+ std::size_t size;
+ std::vector<char> group_names;
+ if (this->comm_.rank() == 0) {
+ auto MPI_SIZE_T_TYPE = Dune::MPITraits<std::size_t>::getType();
+ // NOTE: All slave-master communicators have set a custom error handler, which eventually
+ // will call MPI_Abort() so there is no need to check the return value of any MPI_Recv()
+ // or MPI_Send() calls.
+ MPI_Recv(
+ &size,
+ /*count=*/1,
+ /*datatype=*/MPI_SIZE_T_TYPE,
+ /*source_rank=*/0,
+ /*tag=*/static_cast<int>(MessageTag::MasterGroupNamesSize),
+ this->slave_master_comm_,
+ MPI_STATUS_IGNORE
+ );
+ OpmLog::info("Received master group names size from master process rank 0");
+ group_names.resize(size);
+ MPI_Recv(
+ group_names.data(),
+ /*count=*/size,
+ /*datatype=*/MPI_CHAR,
+ /*source_rank=*/0,
+ /*tag=*/static_cast<int>(MessageTag::MasterGroupNames),
+ this->slave_master_comm_,
+ MPI_STATUS_IGNORE
+ );
+ OpmLog::info("Received master group names from master process rank 0");
+ }
+ this->comm_.broadcast(&size, /*count=*/1, /*emitter_rank=*/0);
+ if (this->comm_.rank() != 0) {
+ group_names.resize(size);
+ }
+ this->comm_.broadcast(group_names.data(), /*count=*/size, /*emitter_rank=*/0);
+ this->saveMasterGroupNamesAsMap_(group_names);
+ this->checkGrupSlavGroupNames_();
+}
+
+void
+ReservoirCouplingSlave::
+sendNextReportDateToMasterProcess() const
+{
+ if (this->comm_.rank() == 0) {
+ double elapsed_time = this->timer_.simulationTimeElapsed();
+ double current_step_length = this->timer_.currentStepLength();
+ // NOTE: This is an offset in seconds from the start date, so it will be 0 if the next report
+ // would be the start date. In general, it should be a positive number.
+ double next_report_time_offset = elapsed_time + current_step_length;
+ // NOTE: All slave-master communicators have set a custom error handler, which eventually
+ // will call MPI_Abort() so there is no need to check the return value of any MPI_Recv()
+ // or MPI_Send() calls.
+ MPI_Send(
+ &next_report_time_offset,
+ /*count=*/1,
+ /*datatype=*/MPI_DOUBLE,
+ /*dest_rank=*/0,
+ /*tag=*/static_cast<int>(MessageTag::SlaveNextReportDate),
+ this->slave_master_comm_
+ );
+ OpmLog::info("Sent next report date to master process from rank 0");
+ }
+}
+
+void
+ReservoirCouplingSlave::
+sendActivationDateToMasterProcess() const
+{
+ if (this->comm_.rank() == 0) {
+ // NOTE: The master process needs the s
+ double activation_date = this->getGrupSlavActivationDate_();
+ // NOTE: All slave-master communicators have set a custom error handler, which eventually
+ // will call MPI_Abort() so there is no need to check the return value of any MPI_Recv()
+ // or MPI_Send() calls.
+ MPI_Send(
+ &activation_date,
+ /*count=*/1,
+ /*datatype=*/MPI_DOUBLE,
+ /*dest_rank=*/0,
+ /*tag=*/static_cast<int>(MessageTag::SlaveActivationDate),
+ this->slave_master_comm_
+ );
+ OpmLog::info("Sent simulation activation date to master process from rank 0");
+ }
+}
+
+void
+ReservoirCouplingSlave::
+sendSimulationStartDateToMasterProcess() const
+{
+ if (this->comm_.rank() == 0) {
+ // NOTE: The master process needs the s
+ double start_date = this->schedule_.getStartTime();
+ // NOTE: All slave-master communicators have set a custom error handler, which eventually
+ // will call MPI_Abort() so there is no need to check the return value of any MPI_Recv()
+ // or MPI_Send() calls.
+ MPI_Send(
+ &start_date,
+ /*count=*/1,
+ /*datatype=*/MPI_DOUBLE,
+ /*dest_rank=*/0,
+ /*tag=*/static_cast<int>(MessageTag::SlaveSimulationStartDate),
+ this->slave_master_comm_
+ );
+ OpmLog::info("Sent simulation start date to master process from rank 0");
+ }
+}
+
+// ------------------
+// Private methods
+// ------------------
+
+void
+ReservoirCouplingSlave::
+checkGrupSlavGroupNames_()
+{
+ // Validate that each slave group name has a corresponding master group name
+ bool grup_slav_found = false;
+ for (std::size_t report_step = 0; report_step < this->schedule_.size(); ++report_step) {
+ auto rescoup = this->schedule_[report_step].rescoup();
+ if (rescoup.grupSlavCount() > 0) {
+ grup_slav_found = true;
+ auto grup_slavs = rescoup.grupSlavs();
+ for (const auto& [slave_group_name, grup_slav] : grup_slavs) {
+ auto master_group_name_it = this->master_group_names_.find(slave_group_name);
+ if (master_group_name_it == this->master_group_names_.end()) {
+ OPM_THROW(std::runtime_error,
+ "Reservoir coupling: Failed to find master group name for slave group: "
+ + slave_group_name);
+ }
+ else {
+ auto master_group_name = master_group_name_it->second;
+ if (grup_slav.masterGroupName() != master_group_name) {
+ OPM_THROW(std::runtime_error,
+ "Reservoir coupling: Inconsistent master group name for slave group: "
+ + slave_group_name);
+ }
+ }
+ }
+ }
+ }
+ if (!grup_slav_found) {
+ OPM_THROW(std::runtime_error, "Reservoir coupling: Failed to find slave group names: "
+ "No GRUPSLAV keyword found in schedule");
+ }
+}
+
+double
+ReservoirCouplingSlave::
+getGrupSlavActivationDate_() const
+{
+ double start_date = this->schedule_.getStartTime();
+ for (std::size_t report_step = 0; report_step < this->schedule_.size(); ++report_step) {
+ auto rescoup = this->schedule_[report_step].rescoup();
+ if (rescoup.grupSlavCount() > 0) {
+ return start_date + this->schedule_.seconds(report_step);
+ }
+ }
+ OPM_THROW(std::runtime_error, "Reservoir coupling: Failed to find slave activation time: "
+ "No GRUPSLAV keyword found in schedule");
+}
+
+void
+ReservoirCouplingSlave::
+maybeActivate(int report_step) {
+ if (!this->activated()) {
+ auto rescoup = this->schedule_[report_step].rescoup();
+ if (rescoup.grupSlavCount() > 0) {
+ this->activated_ = true;
+ }
+ }
+}
+
+void
+ReservoirCouplingSlave::
+saveMasterGroupNamesAsMap_(const std::vector<char>& group_names) {
+ // Deserialize the group names vector into a map of slavegroup names -> mastergroup names
+ auto total_size = group_names.size();
+ std::size_t offset = 0;
+ while (offset < total_size) {
+ std::string master_group{group_names.data() + offset};
+ offset += master_group.size() + 1;
+ assert(offset < total_size);
+ std::string slave_group{group_names.data() + offset};
+ offset += slave_group.size() + 1;
+ this->master_group_names_[slave_group] = master_group;
+ }
+}
+
+
+} // namespace Opm
diff --git a/opm/simulators/flow/ReservoirCouplingSlave.hpp b/opm/simulators/flow/ReservoirCouplingSlave.hpp
new file mode 100644
index 00000000000..61a8d3c0cd8
--- /dev/null
+++ b/opm/simulators/flow/ReservoirCouplingSlave.hpp
@@ -0,0 +1,65 @@
+/*
+ Copyright 2024 Equinor ASA
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_RESERVOIR_COUPLING_SLAVE_HPP
+#define OPM_RESERVOIR_COUPLING_SLAVE_HPP
+
+#include <opm/simulators/flow/ReservoirCoupling.hpp>
+#include <opm/input/eclipse/Schedule/Schedule.hpp>
+#include <opm/simulators/utils/ParallelCommunication.hpp>
+#include <opm/simulators/timestepping/SimulatorTimer.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+
+#include <mpi.h>
+
+#include <vector>
+
+namespace Opm {
+
+class ReservoirCouplingSlave {
+public:
+ using MessageTag = ReservoirCoupling::MessageTag;
+
+ ReservoirCouplingSlave(
+ const Parallel::Communication &comm, const Schedule &schedule, const SimulatorTimer &timer
+ );
+ bool activated() const { return activated_; }
+ void maybeActivate(int report_step);
+ void sendActivationDateToMasterProcess() const;
+ void sendNextReportDateToMasterProcess() const;
+ void sendSimulationStartDateToMasterProcess() const;
+ void receiveMasterGroupNamesFromMasterProcess();
+ double receiveNextTimeStepFromMaster();
+
+private:
+ void checkGrupSlavGroupNames_();
+ double getGrupSlavActivationDate_() const;
+ void saveMasterGroupNamesAsMap_(const std::vector<char>& group_names);
+
+ const Parallel::Communication &comm_;
+ const Schedule& schedule_;
+ const SimulatorTimer &timer_;
+ // MPI parent communicator for a slave process
+ MPI_Comm slave_master_comm_{MPI_COMM_NULL};
+ std::map<std::string, std::string> master_group_names_;
+ bool activated_{false};
+};
+
+} // namespace Opm
+#endif // OPM_RESERVOIR_COUPLING_SLAVE_HPP
diff --git a/opm/simulators/flow/ReservoirCouplingSpawnSlaves.cpp b/opm/simulators/flow/ReservoirCouplingSpawnSlaves.cpp
new file mode 100644
index 00000000000..741bbb224aa
--- /dev/null
+++ b/opm/simulators/flow/ReservoirCouplingSpawnSlaves.cpp
@@ -0,0 +1,335 @@
+/*
+ Copyright 2024 Equinor ASA
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#include <opm/simulators/flow/ReservoirCouplingSpawnSlaves.hpp>
+
+#include <opm/input/eclipse/Schedule/ResCoup/ReservoirCouplingInfo.hpp>
+#include <opm/input/eclipse/Schedule/ResCoup/MasterGroup.hpp>
+#include <opm/input/eclipse/Schedule/ResCoup/Slaves.hpp>
+#include <opm/common/ErrorMacros.hpp>
+
+#include <opm/simulators/utils/ParallelCommunication.hpp>
+
+#include <opm/input/eclipse/Schedule/Schedule.hpp>
+
+#include <dune/common/parallel/mpitraits.hh>
+
+#include <filesystem>
+#include <vector>
+
+#include <fmt/format.h>
+
+namespace Opm {
+
+// ------------------
+// Public methods
+// ------------------
+
+ReservoirCouplingSpawnSlaves::
+ReservoirCouplingSpawnSlaves(
+ ReservoirCouplingMaster &master,
+ const ReservoirCoupling::CouplingInfo &rescoup,
+ const int report_step
+) :
+ master_{master},
+ rescoup_{rescoup},
+ report_step_{report_step},
+ comm_{master.getComm()}
+{
+}
+
+void
+ReservoirCouplingSpawnSlaves::
+spawn()
+{
+ // spawn the slave processes and get the simulation start date from the slaves,
+ // and finally send the master group names to the slaves
+ this->spawnSlaveProcesses_();
+ this->receiveActivationDateFromSlaves_();
+ this->receiveSimulationStartDateFromSlaves_();
+ this->sendMasterGroupNamesToSlaves_();
+ this->prepareTimeStepping_();
+ OpmLog::info("Reservoir coupling slave processes was spawned successfully");
+}
+
+
+// ------------------
+// Private methods
+// ------------------
+
+std::vector<char *> ReservoirCouplingSpawnSlaves::
+getSlaveArgv_(
+ const std::filesystem::path &data_file,
+ const std::string &slave_name,
+ std::string &log_filename) const
+{
+ // Calculate the size of the slave_argv vector like this:
+ // - We will not use the first argument in argv, as it is the program name
+ // - Replace the data file name in argv with the data_file path
+ // - Insert as first argument --slave-log-file=<slave_name>.log
+ // - Also add the argument "--slave=true" to the argv
+ // - Add a nullptr at the end of the argv
+ // So the size of the slave_argv vector will be argc + 2
+ //
+ // Assume: All parameters will be on the form --parameter=value (as a string without spaces)
+ //
+ // Important: The returned vector will have pointers to argv pointers,
+ // data_file string buffer, and slave_name string buffer. So the caller
+ // must ensure that these buffers are not deallocated before the slave_argv has
+ // been used.
+ auto argc = this->master_.getArgc();
+ auto argv = this->master_.getArgv();
+ std::vector<char *> slave_argv(argc + 2);
+ log_filename = "--slave-log-file=" + slave_name; // .log extension will be added by the slave process
+ slave_argv[0] = const_cast<char*>(log_filename.c_str());
+ for (int i = 1; i < argc; i++) {
+ // Check if the argument starts with "--", if not, we will assume it is a positional argument
+ // and we will replace it with the data file path
+ if (std::string(argv[i]).substr(0, 2) == "--") {
+ slave_argv[i] = argv[i];
+ } else {
+ slave_argv[i] = const_cast<char*>(data_file.c_str());
+ }
+ }
+ slave_argv[argc] = const_cast<char *>("--slave=true");
+ slave_argv[argc+1] = nullptr;
+ return slave_argv;
+}
+
+std::pair<std::vector<char>, std::size_t>
+ReservoirCouplingSpawnSlaves::
+getMasterGroupNamesForSlave_(const std::string &slave_name) const
+{
+ // For the given slave name, get all pairs of master group names and slave group names
+ // Serialize the data such that it can be sent over MPI in one chunk
+ auto master_groups = this->rescoup_.masterGroups();
+ std::vector<std::string> data;
+ std::vector<std::string> master_group_names;
+ for (const auto& [group_name, master_group] : master_groups) {
+ if (master_group.slaveName() == slave_name) {
+ data.push_back(group_name);
+ data.push_back(master_group.slaveGroupName());
+ }
+ }
+ assert(data.size() % 2 == 0);
+ assert(data.size() > 0);
+ return this->serializeStrings_(std::move(data));
+}
+
+void
+ReservoirCouplingSpawnSlaves::
+prepareTimeStepping_()
+{
+ // Prepare the time stepping for the master process
+ // This is done after the slave processes have been spawned
+ // and the master group names have been sent to the slaves
+ auto num_slaves = this->master_.numSlavesStarted();
+ this->master_.resizeNextReportDates(num_slaves);
+}
+
+void
+ReservoirCouplingSpawnSlaves::
+receiveActivationDateFromSlaves_()
+{
+ // Currently, we only use the activation date to check that no slave process
+ // starts before the master process.
+ auto num_slaves = this->master_.numSlavesStarted();
+ if (this->comm_.rank() == 0) {
+ for (unsigned int i = 0; i < num_slaves; i++) {
+ double slave_activation_date;
+ // NOTE: All slave-master communicators have set a custom error handler, which eventually
+ // will call MPI_Abort() so there is no need to check the return value of any MPI_Recv()
+ // or MPI_Send() calls.
+ MPI_Recv(
+ &slave_activation_date,
+ /*count=*/1,
+ /*datatype=*/MPI_DOUBLE,
+ /*source_rank=*/0,
+ /*tag=*/static_cast<int>(MessageTag::SlaveActivationDate),
+ this->master_.getSlaveComm(i),
+ MPI_STATUS_IGNORE
+ );
+ if (slave_activation_date < this->master_.getActivationDate()) {
+ OPM_THROW(std::runtime_error, "Slave process start date is earlier than "
+ "the master process' activation date");
+ }
+ OpmLog::info(
+ fmt::format(
+ "Received activation date from slave process with name: {}. "
+ "Activation date: {}", this->master_.getSlaveName(i), slave_activation_date
+ )
+ );
+ }
+ }
+}
+
+void
+ReservoirCouplingSpawnSlaves::
+receiveSimulationStartDateFromSlaves_()
+{
+ auto num_slaves = this->master_.numSlavesStarted();
+ if (this->comm_.rank() == 0) {
+ for (unsigned int i = 0; i < num_slaves; i++) {
+ double slave_start_date;
+ // NOTE: All slave-master communicators have set a custom error handler, which eventually
+ // will call MPI_Abort() so there is no need to check the return value of any MPI_Recv()
+ // or MPI_Send() calls.
+ MPI_Recv(
+ &slave_start_date,
+ /*count=*/1,
+ /*datatype=*/MPI_DOUBLE,
+ /*source_rank=*/0,
+ /*tag=*/static_cast<int>(MessageTag::SlaveSimulationStartDate),
+ this->master_.getSlaveComm(i),
+ MPI_STATUS_IGNORE
+ );
+ this->master_.addSlaveStartDate(slave_start_date);
+ OpmLog::info(
+ fmt::format(
+ "Received start date from slave process with name: {}. "
+ "Start date: {}", this->master_.getSlaveName(i), slave_start_date
+ )
+ );
+ }
+ }
+ if (this->comm_.rank() != 0) {
+ // Ensure that all ranks have the same number of slave activation dates
+ this->master_.resizeSlaveStartDates(num_slaves);
+ }
+ const double* data = this->master_.getSlaveStartDates();
+ this->comm_.broadcast(const_cast<double *>(data), /*count=*/num_slaves, /*emitter_rank=*/0);
+ OpmLog::info("Broadcasted slave start dates to all ranks");
+}
+
+void
+ReservoirCouplingSpawnSlaves::
+sendMasterGroupNamesToSlaves_()
+{
+ if (this->comm_.rank() == 0) {
+ auto num_slaves = this->master_.numSlavesStarted();
+ for (unsigned int i = 0; i < num_slaves; i++) {
+ auto slave_name = this->master_.getSlaveName(i);
+ auto [group_names, size] = this->getMasterGroupNamesForSlave_(slave_name);
+ // NOTE: All slave-master communicators have set a custom error handler, which eventually
+ // will call MPI_Abort() so there is no need to check the return value of any MPI_Recv()
+ // or MPI_Send() calls.
+ static_assert(std::is_same_v<decltype(size), std::size_t>, "size must be of type std::size_t");
+ auto MPI_SIZE_T_TYPE = Dune::MPITraits<std::size_t>::getType();
+ MPI_Send(
+ &size,
+ /*count=*/1,
+ /*datatype=*/MPI_SIZE_T_TYPE,
+ /*dest_rank=*/0,
+ /*tag=*/static_cast<int>(MessageTag::MasterGroupNamesSize),
+ this->master_.getSlaveComm(i)
+ );
+ MPI_Send(
+ group_names.data(),
+ /*count=*/size,
+ /*datatype=*/MPI_CHAR,
+ /*dest_rank=*/0,
+ /*tag=*/static_cast<int>(MessageTag::MasterGroupNames),
+ this->master_.getSlaveComm(i)
+ );
+ OpmLog::info(fmt::format(
+ "Sent master group names to slave process rank 0 with name: {}", slave_name)
+ );
+ }
+ }
+}
+
+std::pair<std::vector<char>, std::size_t>
+ReservoirCouplingSpawnSlaves::
+serializeStrings_(std::vector<std::string> data) const
+{
+ std::size_t total_size = 0;
+ std::vector<char> serialized_data;
+ for (const auto& str: data) {
+ for (const auto& c: str) {
+ serialized_data.push_back(c);
+ }
+ serialized_data.push_back('\0');
+ total_size += str.size() + 1;
+ }
+ return {serialized_data, total_size};
+}
+
+// NOTE: This functions is executed for all ranks, but only rank 0 will spawn
+// the slave processes
+void
+ReservoirCouplingSpawnSlaves::
+spawnSlaveProcesses_()
+{
+ char *flow_program_name = this->master_.getArgv(0);
+ for (const auto& [slave_name, slave] : this->rescoup_.slaves()) {
+ MPI_Comm master_slave_comm = MPI_COMM_NULL;
+ 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 full_path = dir_path / data_file;
+ std::string log_filename; // the getSlaveArgv() function will set this
+ std::vector<char *> slave_argv = this->getSlaveArgv_(
+ full_path, slave_name, log_filename
+ );
+ auto num_procs = slave.numprocs();
+ std::vector<int> errcodes(num_procs);
+ // TODO: We need to decide how to handle the output from the slave processes..
+ // As far as I can tell, open MPI does not support redirecting the output
+ // to a file, so we might need to implement a custom solution for this
+ int spawn_result = MPI_Comm_spawn(
+ flow_program_name,
+ slave_argv.data(),
+ /*maxprocs=*/num_procs,
+ /*info=*/MPI_INFO_NULL,
+ /*root=*/0, // Rank 0 spawns the slave processes
+ /*comm=*/this->comm_,
+ /*intercomm=*/&master_slave_comm,
+ /*array_of_errcodes=*/errcodes.data()
+ );
+ if (spawn_result != MPI_SUCCESS || (master_slave_comm == MPI_COMM_NULL)) {
+ for (unsigned int i = 0; i < num_procs; i++) {
+ if (errcodes[i] != MPI_SUCCESS) {
+ char error_string[MPI_MAX_ERROR_STRING];
+ int length_of_error_string;
+ MPI_Error_string(errcodes[i], error_string, &length_of_error_string);
+ OpmLog::info(fmt::format("Error spawning process {}: {}", i, error_string));
+ }
+ }
+ OPM_THROW(std::runtime_error, "Failed to spawn slave process");
+ }
+ // NOTE: By installing a custom error handler for all slave-master communicators, which
+ // eventually will call MPI_Abort(), there is no need to check the return value of any
+ // MPI_Recv() or MPI_Send() calls as errors will be caught by the error handler.
+ ReservoirCoupling::setErrhandler(master_slave_comm, /*is_master=*/true);
+ OpmLog::info(
+ fmt::format(
+ "Spawned reservoir coupling slave simulation for slave with name: "
+ "{}. Standard output logfile name: {}.log", slave_name, slave_name
+ )
+ );
+ this->master_.addSlaveCommunicator(master_slave_comm);
+ this->master_.addSlaveName(slave_name);
+ }
+}
+
+} // namespace Opm
diff --git a/opm/simulators/flow/ReservoirCouplingSpawnSlaves.hpp b/opm/simulators/flow/ReservoirCouplingSpawnSlaves.hpp
new file mode 100644
index 00000000000..054f9939b23
--- /dev/null
+++ b/opm/simulators/flow/ReservoirCouplingSpawnSlaves.hpp
@@ -0,0 +1,72 @@
+/*
+ Copyright 2024 Equinor ASA
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_RESERVOIR_COUPLING_SPAWN_SLAVES_HPP
+#define OPM_RESERVOIR_COUPLING_SPAWN_SLAVES_HPP
+
+#include <opm/simulators/utils/ParallelCommunication.hpp>
+#include <opm/simulators/flow/ReservoirCoupling.hpp>
+#include <opm/simulators/flow/ReservoirCouplingMaster.hpp>
+
+#include <opm/input/eclipse/Schedule/ResCoup/ReservoirCouplingInfo.hpp>
+#include <opm/input/eclipse/Schedule/Schedule.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+
+#include <mpi.h>
+
+#include <filesystem>
+#include <vector>
+
+namespace Opm {
+
+class ReservoirCouplingSpawnSlaves {
+public:
+ using MessageTag = ReservoirCoupling::MessageTag;
+
+ ReservoirCouplingSpawnSlaves(
+ ReservoirCouplingMaster &master,
+ const ReservoirCoupling::CouplingInfo &rescoup,
+ const int report_step
+ );
+
+ void spawn();
+
+private:
+ std::pair<std::vector<char>, std::size_t>
+ getMasterGroupNamesForSlave_(const std::string &slave_name) const;
+ std::vector<char *> getSlaveArgv_(
+ const std::filesystem::path &data_file,
+ const std::string &slave_name,
+ std::string &log_filename) const;
+ void prepareTimeStepping_();
+ void receiveActivationDateFromSlaves_();
+ void receiveSimulationStartDateFromSlaves_();
+ void sendMasterGroupNamesToSlaves_();
+ std::pair<std::vector<char>, std::size_t>
+ serializeStrings_(std::vector<std::string> data) const;
+ void spawnSlaveProcesses_();
+
+ ReservoirCouplingMaster &master_;
+ const ReservoirCoupling::CouplingInfo &rescoup_;
+ const int report_step_;
+ const Parallel::Communication &comm_;
+};
+
+} // namespace Opm
+#endif // OPM_RESERVOIR_COUPLING_SPAWN_SLAVES_HPP
diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.cpp b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.cpp
index aac95efe23c..e4f2a8a796c 100644
--- a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.cpp
+++ b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.cpp
@@ -58,6 +58,10 @@ void registerSimulatorParameters()
("FileName for .OPMRST file used to load serialized state. "
"If empty, CASENAME.OPMRST is used.");
Parameters::Hide<Parameters::LoadFile>();
+ Parameters::Register<Parameters::Slave>
+ ("Specify if the simulation is a slave simulation in a master-slave simulation");
+ Parameters::Hide<Parameters::Slave>();
+
}
} // namespace Opm::detail
diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp
index 42e18dd98c3..4198cc3777c 100644
--- a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp
+++ b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp
@@ -24,6 +24,18 @@
#include <opm/common/ErrorMacros.hpp>
+#if HAVE_MPI
+#define RESERVOIR_COUPLING_ENABLED
+#endif
+#ifdef RESERVOIR_COUPLING_ENABLED
+#include <opm/input/eclipse/Schedule/ResCoup/ReservoirCouplingInfo.hpp>
+#include <opm/input/eclipse/Schedule/ResCoup/MasterGroup.hpp>
+#include <opm/input/eclipse/Schedule/ResCoup/Slaves.hpp>
+#include <opm/simulators/flow/ReservoirCouplingMaster.hpp>
+#include <opm/simulators/flow/ReservoirCouplingSlave.hpp>
+#include <opm/common/Exceptions.hpp>
+#endif
+
#include <opm/input/eclipse/Units/UnitSystem.hpp>
#include <opm/grid/utility/StopWatch.hpp>
@@ -63,6 +75,7 @@ struct SaveStep { static constexpr auto* value = ""; };
struct SaveFile { static constexpr auto* value = ""; };
struct LoadFile { static constexpr auto* value = ""; };
struct LoadStep { static constexpr int value = -1; };
+struct Slave { static constexpr bool value = false; };
} // namespace Opm::Parameters
@@ -78,6 +91,8 @@ namespace Opm {
template<class TypeTag>
class SimulatorFullyImplicitBlackoil : private SerializableSim
{
+protected:
+ struct MPI_Comm_Deleter;
public:
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Grid = GetPropType<TypeTag, Properties::Grid>;
@@ -171,9 +186,15 @@ class SimulatorFullyImplicitBlackoil : private SerializableSim
/// \param[in,out] timer governs the requested reporting timesteps
/// \param[in,out] state state of reservoir: pressure, fluxes
/// \return simulation report, with timing data
+#ifdef RESERVOIR_COUPLING_ENABLED
+ SimulatorReport run(SimulatorTimer& timer, int argc, char** argv)
+ {
+ init(timer, argc, argv);
+#else
SimulatorReport run(SimulatorTimer& timer)
{
init(timer);
+#endif
// Make cache up to date. No need for updating it in elementCtx.
// NB! Need to be at the correct step in case of restart
simulator_.setEpisodeIndex(timer.currentStepNum());
@@ -188,8 +209,69 @@ class SimulatorFullyImplicitBlackoil : private SerializableSim
return finalize();
}
+#ifdef RESERVOIR_COUPLING_ENABLED
+ // This method should only be called if slave mode (i.e. Parameters::Get<Parameters::Slave>())
+ // is false. We try to determine if this is a normal flow simulation or a reservoir
+ // coupling master. It is a normal flow simulation if the schedule does not contain
+ // any SLAVES and GRUPMAST keywords.
+ bool checkRunningAsReservoirCouplingMaster()
+ {
+ for (std::size_t report_step = 0; report_step < this->schedule().size(); ++report_step) {
+ auto rescoup = this->schedule()[report_step].rescoup();
+ auto slave_count = rescoup.slaveCount();
+ auto master_group_count = rescoup.masterGroupCount();
+ // - GRUPMAST and SLAVES keywords need to be specified at the same report step
+ // - They can only occur once in the schedule
+ if (slave_count > 0 && master_group_count > 0) {
+ return true;
+ }
+ else if (slave_count > 0 && master_group_count == 0) {
+ throw ReservoirCouplingError(
+ "Inconsistent reservoir coupling master schedule: "
+ "Slave count is greater than 0 but master group count is 0"
+ );
+ }
+ else if (slave_count == 0 && master_group_count > 0) {
+ throw ReservoirCouplingError(
+ "Inconsistent reservoir coupling master schedule: "
+ "Master group count is greater than 0 but slave count is 0"
+ );
+ }
+ }
+ return false;
+ }
+#endif
+
+#ifdef RESERVOIR_COUPLING_ENABLED
+ // NOTE: The argc and argv will be used when launching a slave process
+ void init(SimulatorTimer &timer, int argc, char** argv)
+ {
+ auto slave_mode = Parameters::Get<Parameters::Slave>();
+ if (slave_mode) {
+ this->reservoirCouplingSlave_ =
+ std::make_unique<ReservoirCouplingSlave>(
+ FlowGenericVanguard::comm(),
+ this->schedule(), timer
+ );
+ this->reservoirCouplingSlave_->sendActivationDateToMasterProcess();
+ this->reservoirCouplingSlave_->sendSimulationStartDateToMasterProcess();
+ this->reservoirCouplingSlave_->receiveMasterGroupNamesFromMasterProcess();
+ }
+ else {
+ auto master_mode = checkRunningAsReservoirCouplingMaster();
+ if (master_mode) {
+ this->reservoirCouplingMaster_ =
+ std::make_unique<ReservoirCouplingMaster>(
+ FlowGenericVanguard::comm(),
+ this->schedule(),
+ argc, argv
+ );
+ }
+ }
+#else
void init(SimulatorTimer &timer)
{
+#endif
simulator_.setEpisodeIndex(-1);
// Create timers and file for writing timing info.
@@ -212,7 +294,14 @@ class SimulatorFullyImplicitBlackoil : private SerializableSim
else {
adaptiveTimeStepping_ = std::make_unique<TimeStepper>(unitSystem, max_next_tstep, terminalOutput_);
}
-
+#ifdef RESERVOIR_COUPLING_ENABLED
+ if (this->reservoirCouplingSlave_) {
+ adaptiveTimeStepping_->setReservoirCouplingSlave(this->reservoirCouplingSlave_.get());
+ }
+ else if (this->reservoirCouplingMaster_) {
+ adaptiveTimeStepping_->setReservoirCouplingMaster(this->reservoirCouplingMaster_.get());
+ }
+#endif
if (isRestart()) {
// For restarts the simulator may have gotten some information
// about the next timestep size from the OPMEXTRA field
@@ -358,6 +447,14 @@ class SimulatorFullyImplicitBlackoil : private SerializableSim
tuningUpdater(timer.simulationTimeElapsed(),
this->adaptiveTimeStepping_->suggestedNextStep(), 0);
+#ifdef RESERVOIR_COUPLING_ENABLED
+ if (this->reservoirCouplingMaster_) {
+ this->reservoirCouplingMaster_->maybeSpawnSlaveProcesses(timer.currentStepNum());
+ }
+ else if (this->reservoirCouplingSlave_) {
+ this->reservoirCouplingSlave_->maybeActivate(timer.currentStepNum());
+ }
+#endif
const auto& events = schedule()[timer.currentStepNum()].events();
bool event = events.hasEvent(ScheduleEvents::NEW_WELL) ||
events.hasEvent(ScheduleEvents::INJECTION_TYPE_CHANGED) ||
@@ -555,6 +652,12 @@ class SimulatorFullyImplicitBlackoil : private SerializableSim
SimulatorConvergenceOutput convergence_output_;
+#ifdef RESERVOIR_COUPLING_ENABLED
+ bool slaveMode_{false};
+ std::unique_ptr<ReservoirCouplingMaster> reservoirCouplingMaster_{nullptr};
+ std::unique_ptr<ReservoirCouplingSlave> reservoirCouplingSlave_{nullptr};
+#endif
+
SimulatorSerializer serializer_;
};
diff --git a/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp b/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp
index 2c36efcf56d..4fe45e5ab35 100644
--- a/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp
+++ b/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp
@@ -36,25 +36,28 @@
namespace Opm
{
AdaptiveSimulatorTimer::
- AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer,
- const double lastStepTaken,
- const double maxTimeStep )
- : start_date_time_( std::make_shared<boost::posix_time::ptime>(timer.startDateTime()) )
- , start_time_( timer.simulationTimeElapsed() )
- , total_time_( start_time_ + timer.currentStepLength() )
- , report_step_( timer.reportStepNum() )
- , max_time_step_( maxTimeStep )
- , current_time_( start_time_ )
- , dt_( 0.0 )
- , current_step_( 0 )
- , steps_()
- , lastStepFailed_( false )
+ AdaptiveSimulatorTimer( const boost::posix_time::ptime simulation_start_time,
+ const double step_length,
+ const double elapsed_time,
+ const double last_step_taken,
+ const int report_step,
+ const double max_time_step )
+ : start_date_time_{ std::make_shared<boost::posix_time::ptime>(simulation_start_time) }
+ , start_time_{ elapsed_time }
+ , total_time_{ start_time_ + step_length }
+ , report_step_{ report_step }
+ , max_time_step_{ max_time_step }
+ , current_time_{ start_time_ }
+ , dt_{ 0.0 }
+ , current_step_{ 0 }
+ , steps_{}
+ , last_step_failed_{ false }
{
// reserve memory for sub steps
steps_.reserve( 10 );
// set appropriate value for dt_
- provideTimeStepEstimate( lastStepTaken );
+ provideTimeStepEstimate( last_step_taken );
}
bool AdaptiveSimulatorTimer::initialStep () const
diff --git a/opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp b/opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp
index 3dffd49084d..31ae49f4ce8 100644
--- a/opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp
+++ b/opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp
@@ -44,9 +44,12 @@ namespace Opm
/// \param timer in case of sub stepping this is the outer timer
/// \param lastStepTaken last suggested time step
/// \param maxTimeStep maximum time step allowed
- AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer,
- const double lastStepTaken,
- const double maxTimeStep = std::numeric_limits<double>::max() );
+ AdaptiveSimulatorTimer( const boost::posix_time::ptime simulation_start_time,
+ const double step_length,
+ const double elapsed_time,
+ const double last_step_taken,
+ const int report_step,
+ const double max_time_step = std::numeric_limits<double>::max() );
/// \brief advance time by currentStepLength
AdaptiveSimulatorTimer& operator++ ();
@@ -101,10 +104,10 @@ namespace Opm
boost::posix_time::ptime startDateTime() const;
/// \brief Return true if last time step failed
- bool lastStepFailed() const {return lastStepFailed_;}
+ bool lastStepFailed() const {return last_step_failed_;}
/// \brief tell the timestepper whether timestep failed or not
- void setLastStepFailed(bool lastStepFailed) {lastStepFailed_ = lastStepFailed;}
+ void setLastStepFailed(bool last_step_failed) {last_step_failed_ = last_step_failed;}
/// return copy of object
virtual std::unique_ptr< SimulatorTimerInterface > clone() const;
@@ -121,7 +124,7 @@ namespace Opm
int current_step_;
std::vector< double > steps_;
- bool lastStepFailed_;
+ bool last_step_failed_;
};
diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp
index 822d5ee8a57..a90aefd8e70 100644
--- a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp
+++ b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp
@@ -19,12 +19,23 @@
#include <opm/simulators/timestepping/TimeStepControl.hpp>
#include <opm/simulators/timestepping/TimeStepControlInterface.hpp>
-#include <cmath>
+#if HAVE_MPI
+#define RESERVOIR_COUPLING_ENABLED
+#endif
+#ifdef RESERVOIR_COUPLING_ENABLED
+#include <opm/simulators/flow/ReservoirCoupling.hpp>
+#include <opm/simulators/flow/ReservoirCouplingMaster.hpp>
+#include <opm/simulators/flow/ReservoirCouplingSlave.hpp>
+#endif
+
#include <functional>
#include <memory>
#include <set>
+#include <sstream>
+#include <stdexcept>
#include <string>
#include <tuple>
+#include <utility>
#include <vector>
namespace Opm::Parameters {
@@ -55,148 +66,210 @@ class UnitSystem;
struct StepReport;
namespace detail {
-
-void logTimer(const AdaptiveSimulatorTimer& substepTimer);
-
-std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr);
-
-void registerAdaptiveParameters();
-
-std::tuple<TimeStepControlType,
- std::unique_ptr<TimeStepControlInterface>,
- bool>
-createController(const UnitSystem& unitSystem);
-
+ void logTimer(const AdaptiveSimulatorTimer& substep_timer);
+ std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr);
+ void registerAdaptiveParameters();
+ std::tuple<TimeStepControlType, std::unique_ptr<TimeStepControlInterface>, bool>
+ createController(const UnitSystem& unitSystem);
}
- // AdaptiveTimeStepping
- //---------------------
- template<class TypeTag>
- class AdaptiveTimeStepping
+template<class TypeTag>
+class AdaptiveTimeStepping
+{
+private:
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ template <class Solver>
+ class SolutionTimeErrorSolverWrapper : public RelativeChangeInterface
{
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- template <class Solver>
- class SolutionTimeErrorSolverWrapper : public RelativeChangeInterface
- {
- const Solver& solver_;
- public:
- SolutionTimeErrorSolverWrapper(const Solver& solver)
- : solver_(solver)
- {}
-
- /// return || u^n+1 - u^n || / || u^n+1 ||
- double relativeChange() const
- { return solver_.model().relativeChange(); }
- };
-
- template<class E>
- void logException_(const E& exception, bool verbose)
- {
- if (verbose) {
- std::string message;
- message = "Caught Exception: ";
- message += exception.what();
- OpmLog::debug(message);
- }
- }
-
public:
- AdaptiveTimeStepping() = default;
-
- //! \brief contructor taking parameter object
- explicit AdaptiveTimeStepping(const UnitSystem& unitSystem,
- const double max_next_tstep = -1.0,
- const bool terminalOutput = true);
-
- //! \brief contructor taking parameter object
- //! \param tuning Pointer to ecl TUNING keyword
- //! \param timeStep current report step
- AdaptiveTimeStepping(double max_next_tstep,
- const Tuning& tuning,
- const UnitSystem& unitSystem,
- const bool terminalOutput = true);
-
- static void registerParameters();
-
- /** \brief step method that acts like the solver::step method
- in a sub cycle of time steps
- \param tuningUpdater Function used to update TUNING parameters before each
- time step. ACTIONX might change tuning.
- */
- template <class Solver>
- SimulatorReport step(const SimulatorTimer& simulatorTimer,
- Solver& solver,
- const bool isEvent,
- const std::function<bool(const double, const double, const int)> tuningUpdater);
-
- /** \brief Returns the simulator report for the failed substeps of the last
- * report step.
- */
- double suggestedNextStep() const
- { return suggestedNextTimestep_; }
-
- void setSuggestedNextStep(const double x)
- { suggestedNextTimestep_ = x; }
+ SolutionTimeErrorSolverWrapper(const Solver& solver);
+ double relativeChange() const;
+ private:
+ const Solver& solver_;
+ };
- void updateTUNING(double max_next_tstep, const Tuning& tuning);
+ // Forward declaration of SubStepIteration
+ // TODO: This forward declaration is not enough for GCC version 9.4.0 (released June 2021),
+ // but it works fine with GCC version 13.2.0 (released July 2023).
+ template <class Solver> class SubStepIteration;
- void updateNEXTSTEP(double max_next_tstep);
+ template <class Solver>
+ class SubStepper {
+ public:
+ SubStepper(
+ AdaptiveTimeStepping<TypeTag>& adaptive_time_stepping,
+ const SimulatorTimer& simulator_timer,
+ Solver& solver,
+ const bool is_event,
+ const std::function<bool(const double, const double, const int)>& tuning_updater
+ );
+
+ AdaptiveTimeStepping<TypeTag>& getAdaptiveTimerStepper();
+ SimulatorReport run();
+ friend class AdaptiveTimeStepping<TypeTag>::template SubStepIteration<Solver>;
- template<class Serializer>
- void serializeOp(Serializer& serializer);
+ private:
+ bool isReservoirCouplingMaster_() const;
+ bool isReservoirCouplingSlave_() const;
+ void maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double originalTimeStep);
+ bool maybeUpdateTuning_(double elapsed, double dt, int sub_step_number) const;
+ double maxTimeStep_() const;
+ SimulatorReport runStepOriginal_();
+#ifdef RESERVOIR_COUPLING_ENABLED
+ ReservoirCouplingMaster& reservoirCouplingMaster_();
+ ReservoirCouplingSlave& reservoirCouplingSlave_();
+ SimulatorReport runStepReservoirCouplingMaster_();
+ SimulatorReport runStepReservoirCouplingSlave_();
+#endif
+ double suggestedNextTimestep_() const;
+
+ AdaptiveTimeStepping<TypeTag>& adaptive_time_stepping_;
+ const SimulatorTimer& simulator_timer_;
+ Solver& solver_;
+ const bool is_event_;
+ const std::function<bool(double elapsed, double dt, int sub_step_number)>& tuning_updater_;
+ };
- static AdaptiveTimeStepping<TypeTag> serializationTestObjectHardcoded();
- static AdaptiveTimeStepping<TypeTag> serializationTestObjectPID();
- static AdaptiveTimeStepping<TypeTag> serializationTestObjectPIDIt();
- static AdaptiveTimeStepping<TypeTag> serializationTestObjectSimple();
+ template <class Solver>
+ class SubStepIteration {
+ public:
+ SubStepIteration(
+ SubStepper<Solver>& substepper,
+ AdaptiveSimulatorTimer& substep_timer,
+ const double original_time_step,
+ const bool final_step
+ );
- bool operator==(const AdaptiveTimeStepping<TypeTag>& rhs) const;
+ SimulatorReport run();
private:
- template<class Controller>
- static AdaptiveTimeStepping<TypeTag> serializationTestObject_();
-
- template<class T, class Serializer>
- void allocAndSerialize(Serializer& serializer)
- {
- if (!serializer.isSerializing()) {
- timeStepControl_ = std::make_unique<T>();
- }
- serializer(*static_cast<T*>(timeStepControl_.get()));
- }
-
- template<class T>
- bool castAndComp(const AdaptiveTimeStepping<TypeTag>& Rhs) const
- {
- const T* lhs = static_cast<const T*>(timeStepControl_.get());
- const T* rhs = static_cast<const T*>(Rhs.timeStepControl_.get());
- return *lhs == *rhs;
- }
-
- protected:
- void init_(const UnitSystem& unitSystem);
-
- using TimeStepController = std::unique_ptr<TimeStepControlInterface>;
-
- TimeStepControlType timeStepControlType_{TimeStepControlType::PIDAndIterationCount}; //!< type of time step control object
- TimeStepController timeStepControl_{}; //!< time step control object
- double restartFactor_{}; //!< factor to multiply time step with when solver fails to converge
- double growthFactor_{}; //!< factor to multiply time step when solver recovered from failed convergence
- double maxGrowth_{}; //!< factor that limits the maximum growth of a time step
- double maxTimeStep_{}; //!< maximal allowed time step size in days
- double minTimeStep_{}; //!< minimal allowed time step size before throwing
- bool ignoreConvergenceFailure_{false}; //!< continue instead of stop when minimum time step is reached
- int solverRestartMax_{}; //!< how many restart of solver are allowed
- bool solverVerbose_{false}; //!< solver verbosity
- bool timestepVerbose_{false}; //!< timestep verbosity
- double suggestedNextTimestep_{}; //!< suggested size of next timestep
- bool fullTimestepInitially_{false}; //!< beginning with the size of the time step from data file
- double timestepAfterEvent_{}; //!< suggested size of timestep after an event
- bool useNewtonIteration_{false}; //!< use newton iteration count for adaptive time step control
- double minTimeStepBeforeShuttingProblematicWells_{}; //! < shut problematic wells when time step size in days are less than this
+ bool checkContinueOnUnconvergedSolution_(double dt) const;
+ void checkTimeStepMaxRestartLimit_(const int restarts) const;
+ void checkTimeStepMinLimit_(const int new_time_step) const;
+ void chopTimeStep_(const double new_time_step);
+ bool chopTimeStepOrCloseFailingWells_(const int new_time_step);
+ boost::posix_time::ptime currentDateTime_() const;
+ int getNumIterations_(const SimulatorReportSingle &substep_report) const;
+ double growthFactor_() const;
+ bool ignoreConvergenceFailure_() const;
+ void maybeReportSubStep_(SimulatorReportSingle substep_report) const;
+ double maybeRestrictTimeStepGrowth_(
+ const double dt, double dt_estimate, const int restarts) const;
+ void maybeUpdateTuningAndTimeStep_();
+ double maxGrowth_() const;
+ double minTimeStepBeforeClosingWells_() const;
+ double minTimeStep_() const;
+ double restartFactor_() const;
+ SimulatorReportSingle runSubStep_();
+ int solverRestartMax_() const;
+ double suggestedNextTimestep_() const;
+ void setSuggestedNextStep_(double step);
+ void setTimeStep_(double dt_estimate);
+ Solver& solver_() const;
+ bool solverVerbose_() const;
+ const SimulatorTimer& simulatorTimer_() const;
+ boost::posix_time::ptime startDateTime_() const;
+ double timeStepControlComputeEstimate_(
+ const double dt, const int iterations, double elapsed) const;
+ bool timeStepVerbose_() const;
+ void updateSuggestedNextStep_();
+ bool useNewtonIteration_() const;
+ double writeOutput_() const;
+
+ SubStepper<Solver>& substepper_;
+ AdaptiveSimulatorTimer& substep_timer_;
+ const double original_time_step_;
+ const bool final_step_;
+ std::string cause_of_failure_;
+ AdaptiveTimeStepping<TypeTag>& adaptive_time_stepping_;
};
-}
-#include <opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp>
+public:
+ AdaptiveTimeStepping() = default;
+
+ AdaptiveTimeStepping(
+ const UnitSystem& unitSystem,
+ const double max_next_tstep = -1.0,
+ const bool terminalOutput = true
+ );
+
+ AdaptiveTimeStepping(
+ double max_next_tstep,
+ const Tuning& tuning,
+ const UnitSystem& unitSystem,
+ const bool terminalOutput = true
+ );
+ bool operator==(const AdaptiveTimeStepping<TypeTag>& rhs);
+
+ static void registerParameters();
+#ifdef RESERVOIR_COUPLING_ENABLED
+ void setReservoirCouplingMaster(ReservoirCouplingMaster *reservoir_coupling_master);
+ void setReservoirCouplingSlave(ReservoirCouplingSlave *reservoir_coupling_slave);
+#endif
+ void setSuggestedNextStep(const double x);
+ double suggestedNextStep() const;
+
+ template <class Solver>
+ SimulatorReport step(const SimulatorTimer& simulator_timer,
+ Solver& solver,
+ const bool is_event,
+ const std::function<bool(const double, const double, const int)>
+ tuning_updater);
+
+ void updateTUNING(double max_next_tstep, const Tuning& tuning);
+ void updateNEXTSTEP(double max_next_tstep);
+
+ template<class Serializer>
+ void serializeOp(Serializer& serializer);
+
+ static AdaptiveTimeStepping<TypeTag> serializationTestObjectHardcoded();
+ static AdaptiveTimeStepping<TypeTag> serializationTestObjectPID();
+ static AdaptiveTimeStepping<TypeTag> serializationTestObjectPIDIt();
+ static AdaptiveTimeStepping<TypeTag> serializationTestObjectSimple();
+
+private:
+ void maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step,
+ const bool is_event);
+
+ template<class Controller>
+ static AdaptiveTimeStepping<TypeTag> serializationTestObject_();
+
+ template<class T, class Serializer>
+ void allocAndSerialize(Serializer& serializer);
+
+ template<class T>
+ bool castAndComp(const AdaptiveTimeStepping<TypeTag>& Rhs) const;
+
+protected:
+ void init_(const UnitSystem& unitSystem);
+
+ using TimeStepController = std::unique_ptr<TimeStepControlInterface>;
+
+ TimeStepControlType time_step_control_type_; //!< type of time step control object
+ TimeStepController time_step_control_; //!< time step control object
+ double restart_factor_; //!< factor to multiply time step with when solver fails to converge
+ double growth_factor_; //!< factor to multiply time step when solver recovered from failed convergence
+ double max_growth_; //!< factor that limits the maximum growth of a time step
+ double max_time_step_; //!< maximal allowed time step size in days
+ double min_time_step_; //!< minimal allowed time step size before throwing
+ bool ignore_convergence_failure_; //!< continue instead of stop when minimum time step is reached
+ int solver_restart_max_; //!< how many restart of solver are allowed
+ bool solver_verbose_; //!< solver verbosity
+ bool timestep_verbose_; //!< timestep verbosity
+ double suggested_next_timestep_; //!< suggested size of next timestep
+ bool full_timestep_initially_; //!< beginning with the size of the time step from data file
+ double timestep_after_event_; //!< suggested size of timestep after an event
+ bool use_newton_iteration_; //!< use newton iteration count for adaptive time step control
+
+ //! < shut problematic wells when time step size in days are less than this
+ double min_time_step_before_shutting_problematic_wells_;
+#ifdef RESERVOIR_COUPLING_ENABLED
+ ReservoirCouplingMaster *reservoir_coupling_master_ = nullptr;
+ ReservoirCouplingSlave *reservoir_coupling_slave_ = nullptr;
+#endif
+};
+
+} // namespace Opm
+#include <opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp>
#endif // OPM_ADAPTIVE_TIME_STEPPING_HPP
diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp
index 26343d47ff0..4e113fe5e93 100644
--- a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp
+++ b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp
@@ -1,5 +1,22 @@
/*
+ Copyright 2024 Equinor ASA.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
+
#ifndef OPM_ADAPTIVE_TIME_STEPPING_IMPL_HPP
#define OPM_ADAPTIVE_TIME_STEPPING_IMPL_HPP
@@ -27,537 +44,1262 @@
#include <sstream>
#include <stdexcept>
+#include <boost/date_time/posix_time/posix_time.hpp>
#include <fmt/format.h>
#include <fmt/ranges.h>
namespace Opm {
+/*********************************************
+ * Public methods of AdaptiveTimeStepping
+ * ******************************************/
+
+
+//! \brief contructor taking parameter object
+template<class TypeTag>
+AdaptiveTimeStepping<TypeTag>::AdaptiveTimeStepping(
+ const UnitSystem& unit_system,
+ const double max_next_tstep,
+ const bool terminal_output
+)
+ : time_step_control_{}
+ , restart_factor_{Parameters::Get<Parameters::SolverRestartFactor<Scalar>>()} // 0.33
+ , growth_factor_{Parameters::Get<Parameters::SolverGrowthFactor<Scalar>>()} // 2.0
+ , max_growth_{Parameters::Get<Parameters::SolverMaxGrowth<Scalar>>()} // 3.0
+ , max_time_step_{
+ Parameters::Get<Parameters::SolverMaxTimeStepInDays<Scalar>>() * 24 * 60 * 60} // 365.25
+ , min_time_step_{
+ unit_system.to_si(UnitSystem::measure::time,
+ Parameters::Get<Parameters::SolverMinTimeStep<Scalar>>())} // 1e-12;
+ , ignore_convergence_failure_{
+ Parameters::Get<Parameters::SolverContinueOnConvergenceFailure>()} // false;
+ , solver_restart_max_{Parameters::Get<Parameters::SolverMaxRestarts>()} // 10
+ , solver_verbose_{Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminal_output} // 2
+ , timestep_verbose_{Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminal_output} // 2
+ , suggested_next_timestep_{
+ (max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>()
+ : max_next_tstep) * 24 * 60 * 60} // 1.0
+ , full_timestep_initially_{Parameters::Get<Parameters::FullTimeStepInitially>()} // false
+ , timestep_after_event_{
+ Parameters::Get<Parameters::TimeStepAfterEventInDays<Scalar>>() * 24 * 60 * 60} // 1e30
+ , use_newton_iteration_{false}
+ , min_time_step_before_shutting_problematic_wells_{
+ Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day}
+
+{
+ init_(unit_system);
+}
+
+//! \brief contructor
+//! \param tuning Pointer to ecl TUNING keyword
+template<class TypeTag>
+AdaptiveTimeStepping<TypeTag>::AdaptiveTimeStepping(double max_next_tstep,
+ const Tuning& tuning,
+ const UnitSystem& unit_system,
+ const bool terminal_output
+)
+ : time_step_control_{}
+ , restart_factor_{tuning.TSFCNV}
+ , growth_factor_{tuning.TFDIFF}
+ , max_growth_{tuning.TSFMAX}
+ , max_time_step_{tuning.TSMAXZ} // 365.25
+ , min_time_step_{tuning.TSFMIN} // 0.1;
+ , ignore_convergence_failure_{true}
+ , solver_restart_max_{Parameters::Get<Parameters::SolverMaxRestarts>()} // 10
+ , solver_verbose_{Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminal_output} // 2
+ , timestep_verbose_{Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminal_output} // 2
+ , suggested_next_timestep_{
+ max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>() * 24 * 60 * 60
+ : max_next_tstep} // 1.0
+ , full_timestep_initially_{Parameters::Get<Parameters::FullTimeStepInitially>()} // false
+ , timestep_after_event_{tuning.TMAXWC} // 1e30
+ , use_newton_iteration_{false}
+ , min_time_step_before_shutting_problematic_wells_{
+ Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day}
+{
+ init_(unit_system);
+}
template<class TypeTag>
+bool
AdaptiveTimeStepping<TypeTag>::
-AdaptiveTimeStepping(const UnitSystem& unitSystem,
- const double max_next_tstep,
- const bool terminalOutput)
- : timeStepControl_()
- , restartFactor_(Parameters::Get<Parameters::SolverRestartFactor<Scalar>>()) // 0.33
- , growthFactor_(Parameters::Get<Parameters::SolverGrowthFactor<Scalar>>()) // 2.0
- , maxGrowth_(Parameters::Get<Parameters::SolverMaxGrowth<Scalar>>()) // 3.0
- , maxTimeStep_(Parameters::Get<Parameters::SolverMaxTimeStepInDays<Scalar>>() * 24 * 60 * 60) // 365.25
- , minTimeStep_(unitSystem.to_si(UnitSystem::measure::time, Parameters::Get<Parameters::SolverMinTimeStep<Scalar>>())) // 1e-12;
- , ignoreConvergenceFailure_(Parameters::Get<Parameters::SolverContinueOnConvergenceFailure>()) // false;
- , solverRestartMax_(Parameters::Get<Parameters::SolverMaxRestarts>()) // 10
- , solverVerbose_(Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminalOutput) // 2
- , timestepVerbose_(Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminalOutput) // 2
- , suggestedNextTimestep_((max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>() : max_next_tstep) * 24 * 60 * 60) // 1.0
- , fullTimestepInitially_(Parameters::Get<Parameters::FullTimeStepInitially>()) // false
- , timestepAfterEvent_(Parameters::Get<Parameters::TimeStepAfterEventInDays<Scalar>>() * 24 * 60 * 60) // 1e30
- , useNewtonIteration_(false)
- , minTimeStepBeforeShuttingProblematicWells_(Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
-
-{
- init_(unitSystem);
+operator==(const AdaptiveTimeStepping<TypeTag>& rhs)
+{
+ if (this->time_step_control_type_ != rhs.time_step_control_type_ ||
+ (this->time_step_control_ && !rhs.time_step_control_) ||
+ (!this->time_step_control_ && rhs.time_step_control_)) {
+ return false;
+ }
+
+ bool result = false;
+ switch (this->time_step_control_type_) {
+ case TimeStepControlType::HardCodedTimeStep:
+ result = castAndComp<HardcodedTimeStepControl>(rhs);
+ break;
+ case TimeStepControlType::PIDAndIterationCount:
+ result = castAndComp<PIDAndIterationCountTimeStepControl>(rhs);
+ break;
+ case TimeStepControlType::SimpleIterationCount:
+ result = castAndComp<SimpleIterationCountTimeStepControl>(rhs);
+ break;
+ case TimeStepControlType::PID:
+ result = castAndComp<PIDTimeStepControl>(rhs);
+ break;
+ }
+
+ return result &&
+ this->restart_factor_ == rhs.restart_factor_ &&
+ this->growth_factor_ == rhs.growth_factor_ &&
+ this->max_growth_ == rhs.max_growth_ &&
+ this->max_time_step_ == rhs.max_time_step_ &&
+ this->min_time_step_ == rhs.min_time_step_ &&
+ this->ignore_convergence_failure_ == rhs.ignore_convergence_failure_ &&
+ this->solver_restart_max_== rhs.solver_restart_max_ &&
+ this->solver_verbose_ == rhs.solver_verbose_ &&
+ this->full_timestep_initially_ == rhs.full_timestep_initially_ &&
+ this->timestep_after_event_ == rhs.timestep_after_event_ &&
+ this->use_newton_iteration_ == rhs.use_newton_iteration_ &&
+ this->min_time_step_before_shutting_problematic_wells_ ==
+ rhs.min_time_step_before_shutting_problematic_wells_;
}
template<class TypeTag>
+void
AdaptiveTimeStepping<TypeTag>::
-AdaptiveTimeStepping(double max_next_tstep,
- const Tuning& tuning,
- const UnitSystem& unitSystem,
- const bool terminalOutput)
- : timeStepControl_()
- , restartFactor_(tuning.TSFCNV)
- , growthFactor_(tuning.TFDIFF)
- , maxGrowth_(tuning.TSFMAX)
- , maxTimeStep_(tuning.TSMAXZ) // 365.25
- , minTimeStep_(tuning.TSFMIN) // 0.1;
- , ignoreConvergenceFailure_(true)
- , solverRestartMax_(Parameters::Get<Parameters::SolverMaxRestarts>()) // 10
- , solverVerbose_(Parameters::Get<Parameters::SolverVerbosity>() > 0 && terminalOutput) // 2
- , timestepVerbose_(Parameters::Get<Parameters::TimeStepVerbosity>() > 0 && terminalOutput) // 2
- , suggestedNextTimestep_(max_next_tstep <= 0 ? Parameters::Get<Parameters::InitialTimeStepInDays>() * 24 * 60 * 60 : max_next_tstep) // 1.0
- , fullTimestepInitially_(Parameters::Get<Parameters::FullTimeStepInitially>()) // false
- , timestepAfterEvent_(tuning.TMAXWC) // 1e30
- , useNewtonIteration_(false)
- , minTimeStepBeforeShuttingProblematicWells_(Parameters::Get<Parameters::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
-{
- init_(unitSystem);
-}
-
-template<class TypeTag>
-void AdaptiveTimeStepping<TypeTag>::registerParameters()
+registerParameters()
{
registerEclTimeSteppingParameters<Scalar>();
detail::registerAdaptiveParameters();
}
+#ifdef RESERVOIR_COUPLING_ENABLED
template<class TypeTag>
-template<class Solver>
+void
+AdaptiveTimeStepping<TypeTag>::
+setReservoirCouplingMaster(ReservoirCouplingMaster *reservoir_coupling_master)
+{
+ this->reservoir_coupling_master_ = reservoir_coupling_master;
+}
+
+template<class TypeTag>
+void
+AdaptiveTimeStepping<TypeTag>::
+setReservoirCouplingSlave(ReservoirCouplingSlave *reservoir_coupling_slave)
+{
+ this->reservoir_coupling_slave_ = reservoir_coupling_slave;
+}
+#endif
+
+/** \brief step method that acts like the solver::step method
+ in a sub cycle of time steps
+ \param tuningUpdater Function used to update TUNING parameters before each
+ time step. ACTIONX might change tuning.
+*/
+template<class TypeTag>
+template <class Solver>
SimulatorReport
AdaptiveTimeStepping<TypeTag>::
-step(const SimulatorTimer& simulatorTimer,
+step(const SimulatorTimer& simulator_timer,
Solver& solver,
- const bool isEvent,
- const std::function<bool(const double, const double, const int)> tuningUpdater)
+ const bool is_event,
+ const std::function<bool(const double /*current_time*/,
+ const double /*dt*/,
+ const int /*substep_number*/
+ )> tuning_updater
+)
{
- // Maybe update tuning
- tuningUpdater(simulatorTimer.simulationTimeElapsed(), suggestedNextTimestep_, 0);
- SimulatorReport report;
- const double timestep = simulatorTimer.currentStepLength();
+ SubStepper sub_stepper{
+ *this, simulator_timer, solver, is_event, tuning_updater,
+ };
+ return sub_stepper.run();
+}
+
+template<class TypeTag>
+template<class Serializer>
+void
+AdaptiveTimeStepping<TypeTag>::
+serializeOp(Serializer& serializer)
+{
+ serializer(this->time_step_control_type_);
+ switch (this->time_step_control_type_) {
+ case TimeStepControlType::HardCodedTimeStep:
+ allocAndSerialize<HardcodedTimeStepControl>(serializer);
+ break;
+ case TimeStepControlType::PIDAndIterationCount:
+ allocAndSerialize<PIDAndIterationCountTimeStepControl>(serializer);
+ break;
+ case TimeStepControlType::SimpleIterationCount:
+ allocAndSerialize<SimpleIterationCountTimeStepControl>(serializer);
+ break;
+ case TimeStepControlType::PID:
+ allocAndSerialize<PIDTimeStepControl>(serializer);
+ break;
+ }
+ serializer(this->restart_factor_);
+ serializer(this->growth_factor_);
+ serializer(this->max_growth_);
+ serializer(this->max_time_step_);
+ serializer(this->min_time_step_);
+ serializer(this->ignore_convergence_failure_);
+ serializer(this->solver_restart_max_);
+ serializer(this->solver_verbose_);
+ serializer(this->timestep_verbose_);
+ serializer(this->suggested_next_timestep_);
+ serializer(this->full_timestep_initially_);
+ serializer(this->timestep_after_event_);
+ serializer(this->use_newton_iteration_);
+ serializer(this->min_time_step_before_shutting_problematic_wells_);
+}
+
+template<class TypeTag>
+AdaptiveTimeStepping<TypeTag>
+AdaptiveTimeStepping<TypeTag>::
+serializationTestObjectHardcoded()
+{
+ return serializationTestObject_<HardcodedTimeStepControl>();
+}
+
+template<class TypeTag>
+AdaptiveTimeStepping<TypeTag>
+AdaptiveTimeStepping<TypeTag>::
+serializationTestObjectPID()
+{
+ return serializationTestObject_<PIDTimeStepControl>();
+}
+
+template<class TypeTag>
+AdaptiveTimeStepping<TypeTag>
+AdaptiveTimeStepping<TypeTag>::
+serializationTestObjectPIDIt()
+{
+ return serializationTestObject_<PIDAndIterationCountTimeStepControl>();
+}
+
+template<class TypeTag>
+AdaptiveTimeStepping<TypeTag>
+AdaptiveTimeStepping<TypeTag>::
+serializationTestObjectSimple()
+{
+ return serializationTestObject_<SimpleIterationCountTimeStepControl>();
+}
+
+
+template<class TypeTag>
+void
+AdaptiveTimeStepping<TypeTag>::
+setSuggestedNextStep(const double x)
+{
+ this->suggested_next_timestep_ = x;
+}
+
+template<class TypeTag>
+double
+AdaptiveTimeStepping<TypeTag>::
+suggestedNextStep() const
+{
+ return this->suggested_next_timestep_;
+}
+
+
+template<class TypeTag>
+void
+AdaptiveTimeStepping<TypeTag>::
+updateNEXTSTEP(double max_next_tstep)
+{
+ // \Note Only update next suggested step if TSINIT was explicitly
+ // set in TUNING or NEXTSTEP is active.
+ if (max_next_tstep > 0) {
+ this->suggested_next_timestep_ = max_next_tstep;
+ }
+}
+template<class TypeTag>
+void
+AdaptiveTimeStepping<TypeTag>::
+updateTUNING(double max_next_tstep, const Tuning& tuning)
+{
+ this->restart_factor_ = tuning.TSFCNV;
+ this->growth_factor_ = tuning.TFDIFF;
+ this->max_growth_ = tuning.TSFMAX;
+ this->max_time_step_ = tuning.TSMAXZ;
+ updateNEXTSTEP(max_next_tstep);
+ this->timestep_after_event_ = tuning.TMAXWC;
+}
+
+/*********************************************
+ * Private methods of AdaptiveTimeStepping
+ * ******************************************/
+
+template<class TypeTag>
+template<class T, class Serializer>
+void
+AdaptiveTimeStepping<TypeTag>::
+allocAndSerialize(Serializer& serializer)
+{
+ if (!serializer.isSerializing()) {
+ this->time_step_control_ = std::make_unique<T>();
+ }
+ serializer(*static_cast<T*>(this->time_step_control_.get()));
+}
+
+template<class TypeTag>
+template<class T>
+bool
+AdaptiveTimeStepping<TypeTag>::
+castAndComp(const AdaptiveTimeStepping<TypeTag>& Rhs) const
+{
+ const T* lhs = static_cast<const T*>(this->time_step_control_.get());
+ const T* rhs = static_cast<const T*>(Rhs.time_step_control_.get());
+ return *lhs == *rhs;
+}
+
+template<class TypeTag>
+void
+AdaptiveTimeStepping<TypeTag>::
+maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step,
+ bool is_event)
+{
// init last time step as a fraction of the given time step
- if (suggestedNextTimestep_ < 0) {
- suggestedNextTimestep_ = restartFactor_ * timestep;
+ if (this->suggested_next_timestep_ < 0) {
+ this->suggested_next_timestep_ = this->restart_factor_ * original_time_step;
}
- if (fullTimestepInitially_) {
- suggestedNextTimestep_ = timestep;
+ if (this->full_timestep_initially_) {
+ this->suggested_next_timestep_ = original_time_step;
}
// use seperate time step after event
- if (isEvent && timestepAfterEvent_ > 0) {
- suggestedNextTimestep_ = timestepAfterEvent_;
+ if (is_event && this->timestep_after_event_ > 0) {
+ this->suggested_next_timestep_ = this->timestep_after_event_;
}
+}
- auto& simulator = solver.model().simulator();
- auto& problem = simulator.problem();
+template<class TypeTag>
+template<class Controller>
+AdaptiveTimeStepping<TypeTag>
+AdaptiveTimeStepping<TypeTag>::
+serializationTestObject_()
+{
+ AdaptiveTimeStepping<TypeTag> result;
- // create adaptive step timer with previously used sub step size
- AdaptiveSimulatorTimer substepTimer(simulatorTimer, suggestedNextTimestep_, maxTimeStep_);
+ result.restart_factor_ = 1.0;
+ result.growth_factor_ = 2.0;
+ result.max_growth_ = 3.0;
+ result.max_time_step_ = 4.0;
+ result.min_time_step_ = 5.0;
+ result.ignore_convergence_failure_ = true;
+ result.solver_restart_max_ = 6;
+ result.solver_verbose_ = true;
+ result.timestep_verbose_ = true;
+ result.suggested_next_timestep_ = 7.0;
+ result.full_timestep_initially_ = true;
+ result.use_newton_iteration_ = true;
+ result.min_time_step_before_shutting_problematic_wells_ = 9.0;
+ result.time_step_control_type_ = Controller::Type;
+ result.time_step_control_ =
+ std::make_unique<Controller>(Controller::serializationTestObject());
- // counter for solver restarts
- int restarts = 0;
+ return result;
+}
- // sub step time loop
- while (!substepTimer.done()) {
- // Maybe update tuning
- // get current delta t
- auto oldValue = suggestedNextTimestep_;
- if (tuningUpdater(substepTimer.simulationTimeElapsed(),
- substepTimer.currentStepLength(),
- substepTimer.currentStepNum())) {
- // Use provideTimeStepEstimate to make we sure don't simulate longer than the report step is.
- substepTimer.provideTimeStepEstimate(suggestedNextTimestep_);
- suggestedNextTimestep_ = oldValue;
- }
- const double dt = substepTimer.currentStepLength();
- if (timestepVerbose_) {
- detail::logTimer(substepTimer);
- }
+/*********************************************
+ * Protected methods of AdaptiveTimeStepping
+ * ******************************************/
+
+template<class TypeTag>
+void AdaptiveTimeStepping<TypeTag>::
+init_(const UnitSystem& unitSystem)
+{
+ std::tie(time_step_control_type_,
+ time_step_control_,
+ use_newton_iteration_) = detail::createController(unitSystem);
+ // make sure growth factor is something reasonable
+ if (this->growth_factor_ < 1.0) {
+ OPM_THROW(std::runtime_error,
+ "Growth factor cannot be less than 1.");
+ }
+}
- SimulatorReportSingle substepReport;
- std::string causeOfFailure;
- try {
- substepReport = solver.step(substepTimer);
- if (solverVerbose_) {
- // report number of linear iterations
- OpmLog::debug("Overall linear iterations used: " +
- std::to_string(substepReport.total_linear_iterations));
- }
- }
- catch (const TooManyIterations& e) {
- substepReport = solver.failureReport();
- causeOfFailure = "Solver convergence failure - Iteration limit reached";
- logException_(e, solverVerbose_);
- // since linearIterations is < 0 this will restart the solver
+/************************************************
+ * Private class SubStepper public methods
+ ************************************************/
+
+template<class TypeTag>
+template<class Solver>
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+SubStepper(
+ AdaptiveTimeStepping<TypeTag>& adaptive_time_stepping,
+ const SimulatorTimer& simulator_timer,
+ Solver& solver,
+ const bool is_event,
+ const std::function<bool(const double /*current_time*/,
+ const double /*dt*/,
+ const int /*substep_number*/
+ )>& tuning_updater
+
+)
+ : adaptive_time_stepping_{adaptive_time_stepping}
+ , simulator_timer_{simulator_timer}
+ , solver_{solver}
+ , is_event_{is_event}
+ , tuning_updater_{tuning_updater}
+{
+}
+
+template<class TypeTag>
+template<class Solver>
+AdaptiveTimeStepping<TypeTag>&
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+getAdaptiveTimerStepper()
+{
+ return adaptive_time_stepping_;
+}
+
+template<class TypeTag>
+template<class Solver>
+SimulatorReport
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+run()
+{
+#ifdef RESERVOIR_COUPLING_ENABLED
+ if (isReservoirCouplingSlave_() && reservoirCouplingSlave_().activated()) {
+ return runStepReservoirCouplingSlave_();
+ }
+ else if (isReservoirCouplingMaster_() && reservoirCouplingMaster_().activated()) {
+ return runStepReservoirCouplingMaster_();
+ }
+ else {
+ return runStepOriginal_();
+ }
+#else
+ return runStepOriginal_();
+#endif
+}
+
+/************************************************
+ * Private class SubStepper private methods
+ ************************************************/
+
+
+template<class TypeTag>
+template<class Solver>
+bool
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+isReservoirCouplingMaster_() const
+{
+ return this->adaptive_time_stepping_.reservoir_coupling_master_ != nullptr;
+}
+
+template<class TypeTag>
+template<class Solver>
+bool
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+isReservoirCouplingSlave_() const
+{
+ return this->adaptive_time_stepping_.reservoir_coupling_slave_ != nullptr;
+}
+
+template<class TypeTag>
+template<class Solver>
+void
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step)
+{
+ this->adaptive_time_stepping_.maybeModifySuggestedTimeStepAtBeginningOfReportStep_(
+ original_time_step, this->is_event_
+ );
+}
+
+// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicitBlackoil::runStep()
+// It has to be called for each substep since TUNING might have been changed for next sub step due
+// to ACTIONX (via NEXTSTEP) or WCYCLE keywords.
+template<class TypeTag>
+template<class Solver>
+bool
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+maybeUpdateTuning_(double elapsed, double dt, int sub_step_number) const
+{
+ return this->tuning_updater_(elapsed, dt, sub_step_number);
+}
+
+template<class TypeTag>
+template<class Solver>
+double
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+maxTimeStep_() const
+{
+ return this->adaptive_time_stepping_.max_time_step_;
+}
+
+template <class TypeTag>
+template <class Solver>
+SimulatorReport
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+runStepOriginal_()
+{
+ auto elapsed = this->simulator_timer_.simulationTimeElapsed();
+ auto original_time_step = this->simulator_timer_.currentStepLength();
+ auto report_step = this->simulator_timer_.reportStepNum();
+ maybeUpdateTuning_(elapsed, original_time_step, report_step);
+ maybeModifySuggestedTimeStepAtBeginningOfReportStep_(original_time_step);
+
+ AdaptiveSimulatorTimer substep_timer{
+ this->simulator_timer_.startDateTime(),
+ original_time_step,
+ elapsed,
+ suggestedNextTimestep_(),
+ report_step,
+ maxTimeStep_()
+ };
+ SubStepIteration substepIteration{*this, substep_timer, original_time_step, /*final_step=*/true};
+ return substepIteration.run();
+}
+
+#ifdef RESERVOIR_COUPLING_ENABLED
+template <class TypeTag>
+template <class Solver>
+ReservoirCouplingMaster&
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+reservoirCouplingMaster_()
+{
+ return *adaptive_time_stepping_.reservoir_coupling_master_;
+}
+#endif
+
+#ifdef RESERVOIR_COUPLING_ENABLED
+template <class TypeTag>
+template <class Solver>
+ReservoirCouplingSlave&
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+reservoirCouplingSlave_()
+{
+ return *this->adaptive_time_stepping_.reservoir_coupling_slave_;
+}
+#endif
+
+#ifdef RESERVOIR_COUPLING_ENABLED
+// Description of the reservoir coupling master and slave substep loop
+// -------------------------------------------------------------------
+// The master and slave processes attempts to reach the end of the report step using a series of substeps
+// (also called timesteps). Each substep have an upper limit that is roughly determined by a combination
+// of the keywords TUNING (through the TSINIT, TSMAXZ values), NEXSTEP, WCYCLE, and the start of the
+// next report step (the last substep needs to coincide with this time). Note that NEXTSTEP can be
+// updated from an ACTIONX keyword. Although, this comment focuses on the maximum substep limit,
+// note that there is also a lower limit on the substep length. And the substep sizes will be adjusted
+// automatically (or retried) based on the convergence behavior of the solver and other criteria.
+//
+// When using reservoir coupling, the upper limit on each substep is further limited to: a) not overshoot
+// next report date of a slave reservoir, and b) to keep the flow rate of the slave groups within
+// certain limits. To determine this potential further limitation on the substep, the following procedure
+// is used at the beginning of each master substep:
+// - First, the slaves sends the master the date of their next report step
+// - The master receives the date of the next report step from the slaves
+// - If necessary, the master computes a modified substep length based on the received dates
+// TODO: explain how the substep is limited to keep the flow rate of the slave groups within certain
+// limits. Since this is not implemented yet, this part of the procedure is not explained here.
+//
+// Then, after the master has determined the substep length using the above procedure, it sends it
+// to the slaves. The slaves will now use the end of this substep as a fixed point (like a mini report
+// step), that they will try to reach using a single substep or multiple substeps. The end of the
+// substep received from the master will either conincide with the end of its current report step, or
+// it will end before (it cannot not end after since the master has already adjusted the step to not
+// exceed any report time in a slave). If the end of this substep is earlier in time than its next report
+// date, the slave will enter a loop and wait for the master to send it a new substep.
+// The loop is finished when the end of the received time step conincides with the end of its current
+// report step.
+
+template <class TypeTag>
+template <class Solver>
+SimulatorReport
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+runStepReservoirCouplingMaster_()
+{
+ bool substep_done = false;
+ int iteration = 0;
+ const double original_time_step = this->simulator_timer_.currentStepLength();
+ double current_time{this->simulator_timer_.simulationTimeElapsed()};
+ double step_end_time = current_time + original_time_step;
+ auto current_step_length = original_time_step;
+ SimulatorReport report;
+ while(!substep_done) {
+ reservoirCouplingMaster_().receiveNextReportDateFromSlaves();
+ if (iteration == 0) {
+ maybeUpdateTuning_(current_time, current_step_length, /*substep=*/0);
}
- catch (const ConvergenceMonitorFailure& e) {
- substepReport = solver.failureReport();
- causeOfFailure = "Convergence monitor failure";
+ current_step_length = reservoirCouplingMaster_().maybeChopSubStep(
+ current_step_length, current_time);
+ reservoirCouplingMaster_().sendNextTimeStepToSlaves(current_step_length);
+ if (iteration == 0) {
+ maybeModifySuggestedTimeStepAtBeginningOfReportStep_(current_step_length);
}
- catch (const LinearSolverProblem& e) {
- substepReport = solver.failureReport();
- causeOfFailure = "Linear solver convergence failure";
-
- logException_(e, solverVerbose_);
- // since linearIterations is < 0 this will restart the solver
+ AdaptiveSimulatorTimer substep_timer{
+ this->simulator_timer_.startDateTime(),
+ /*stepLength=*/current_step_length,
+ /*elapsedTime=*/current_time,
+ /*timeStepEstimate=*/suggestedNextTimestep_(),
+ /*reportStep=*/this->simulator_timer_.reportStepNum(),
+ maxTimeStep_()
+ };
+ bool final_step = ReservoirCoupling::Seconds::compare_gt_or_eq(
+ current_time + current_step_length, step_end_time
+ );
+ SubStepIteration substepIteration{*this, substep_timer, current_step_length, final_step};
+ auto sub_steps_report = substepIteration.run();
+ report += sub_steps_report;
+ current_time += current_step_length;
+ if (final_step) {
+ break;
}
- catch (const NumericalProblem& e) {
- substepReport = solver.failureReport();
- causeOfFailure = "Solver convergence failure - Numerical problem encountered";
+ iteration++;
+ }
+ return report;
+}
+#endif
- logException_(e, solverVerbose_);
- // since linearIterations is < 0 this will restart the solver
+#ifdef RESERVOIR_COUPLING_ENABLED
+template <class TypeTag>
+template <class Solver>
+SimulatorReport
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+runStepReservoirCouplingSlave_()
+{
+ bool substep_done = false;
+ int iteration = 0;
+ const double original_time_step = this->simulator_timer_.currentStepLength();
+ double current_time{this->simulator_timer_.simulationTimeElapsed()};
+ double step_end_time = current_time + original_time_step;
+ SimulatorReport report;
+ while(!substep_done) {
+ reservoirCouplingSlave_().sendNextReportDateToMasterProcess();
+ auto timestep = reservoirCouplingSlave_().receiveNextTimeStepFromMaster();
+ if (iteration == 0) {
+ maybeUpdateTuning_(current_time, original_time_step, /*substep=*/0);
+ maybeModifySuggestedTimeStepAtBeginningOfReportStep_(timestep);
}
- catch (const std::runtime_error& e) {
- substepReport = solver.failureReport();
-
- logException_(e, solverVerbose_);
- // also catch linear solver not converged
+ AdaptiveSimulatorTimer substep_timer{
+ this->simulator_timer_.startDateTime(),
+ /*step_length=*/timestep,
+ /*elapsed_time=*/current_time,
+ /*time_step_estimate=*/suggestedNextTimestep_(),
+ this->simulator_timer_.reportStepNum(),
+ maxTimeStep_()
+ };
+ bool final_step = ReservoirCoupling::Seconds::compare_gt_or_eq(
+ current_time + timestep, step_end_time
+ );
+ SubStepIteration substepIteration{*this, substep_timer, timestep, final_step};
+ auto sub_steps_report = substepIteration.run();
+ report += sub_steps_report;
+ current_time += timestep;
+ if (final_step) {
+ substep_done = true;
+ break;
}
- catch (const Dune::ISTLError& e) {
- substepReport = solver.failureReport();
+ iteration++;
+ }
+ return report;
+}
- logException_(e, solverVerbose_);
- // also catch errors in ISTL AMG that occur when time step is too large
- }
- catch (const Dune::MatrixBlockError& e) {
- substepReport = solver.failureReport();
+#endif
+
+template <class TypeTag>
+template <class Solver>
+double
+AdaptiveTimeStepping<TypeTag>::SubStepper<Solver>::
+suggestedNextTimestep_() const
+{
+ return this->adaptive_time_stepping_.suggestedNextStep();
+}
- logException_(e, solverVerbose_);
- // this can be thrown by ISTL's ILU0 in block mode, yet is not an ISTLError
+
+
+/************************************************
+ * Private class SubStepIteration public methods
+ ************************************************/
+
+template<class TypeTag>
+template<class Solver>
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+SubStepIteration(
+ SubStepper<Solver>& substepper,
+ AdaptiveSimulatorTimer& substep_timer,
+ const double original_time_step,
+ bool final_step
+)
+ : substepper_{substepper}
+ , substep_timer_{substep_timer}
+ , original_time_step_{original_time_step}
+ , final_step_{final_step}
+ , adaptive_time_stepping_{substepper.getAdaptiveTimerStepper()}
+{
+}
+
+template <class TypeTag>
+template <class Solver>
+SimulatorReport
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+run()
+{
+ auto& simulator = solver_().model().simulator();
+ auto& problem = simulator.problem();
+ // counter for solver restarts
+ int restarts = 0;
+ SimulatorReport report;
+
+ // sub step time loop
+ while (!this->substep_timer_.done()) {
+ maybeUpdateTuningAndTimeStep_();
+ const double dt = this->substep_timer_.currentStepLength();
+ if (timeStepVerbose_()) {
+ detail::logTimer(this->substep_timer_);
}
+ auto substep_report = runSubStep_();
+
//Pass substep to eclwriter for summary output
- simulator.problem().setSubStepReport(substepReport);
-
- report += substepReport;
-
- bool continue_on_uncoverged_solution = ignoreConvergenceFailure_ &&
- !substepReport.converged &&
- dt <= minTimeStep_;
-
- if (continue_on_uncoverged_solution && solverVerbose_) {
- const auto msg = fmt::format(
- "Solver failed to converge but timestep "
- "{} is smaller or equal to {}\n"
- "which is the minimum threshold given "
- "by option --solver-min-time-step\n",
- dt, minTimeStep_
- );
- OpmLog::problem(msg);
- }
+ problem.setSubStepReport(substep_report);
- if (substepReport.converged || continue_on_uncoverged_solution) {
-
- // advance by current dt
- ++substepTimer;
-
- // create object to compute the time error, simply forwards the call to the model
- SolutionTimeErrorSolverWrapper<Solver> relativeChange(solver);
-
- // compute new time step estimate
- const int iterations = useNewtonIteration_ ? substepReport.total_newton_iterations
- : substepReport.total_linear_iterations;
- double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
- substepTimer.simulationTimeElapsed());
-
- assert(dtEstimate > 0);
- // limit the growth of the timestep size by the growth factor
- dtEstimate = std::min(dtEstimate, double(maxGrowth_ * dt));
- assert(dtEstimate > 0);
- // further restrict time step size growth after convergence problems
- if (restarts > 0) {
- dtEstimate = std::min(growthFactor_ * dt, dtEstimate);
- // solver converged, reset restarts counter
- restarts = 0;
- }
+ report += substep_report;
- if (timestepVerbose_) {
- std::ostringstream ss;
- substepReport.reportStep(ss);
- OpmLog::info(ss.str());
- }
+ if (substep_report.converged || checkContinueOnUnconvergedSolution_(dt)) {
+ ++this->substep_timer_; // advance by current dt
- // write data if outputWriter was provided
- // if the time step is done we do not need
- // to write it as this will be done by the simulator
- // anyway.
- if (!substepTimer.done()) {
- time::StopWatch perfTimer;
- perfTimer.start();
+ const int iterations = getNumIterations_(substep_report);
+ auto dt_estimate = timeStepControlComputeEstimate_(
+ dt, iterations, this->substep_timer_.simulationTimeElapsed());
- problem.writeOutput(true);
+ assert(dt_estimate > 0);
+ dt_estimate = maybeRestrictTimeStepGrowth_(dt, dt_estimate, restarts);
+ restarts = 0; // solver converged, reset restarts counter
- report.success.output_write_time += perfTimer.secsSinceStart();
+ maybeReportSubStep_(substep_report);
+ if (this->final_step_ && this->substep_timer_.done()) {
+ // if the time step is done we do not need to write it as this will be done
+ // by the simulator anyway.
+ }
+ else {
+ report.success.output_write_time += writeOutput_();
}
// set new time step length
- substepTimer.provideTimeStepEstimate(dtEstimate);
-
- report.success.converged = substepTimer.done();
- substepTimer.setLastStepFailed(false);
+ setTimeStep_(dt_estimate);
+ report.success.converged = this->substep_timer_.done();
+ this->substep_timer_.setLastStepFailed(false);
}
else { // in case of no convergence
- substepTimer.setLastStepFailed(true);
-
- // If we have restarted (i.e. cut the timestep) too
- // many times, we have failed and throw an exception.
- if (restarts >= solverRestartMax_) {
- const auto msg = fmt::format(
- "Solver failed to converge after cutting timestep {} times.",
- restarts
- );
- if (solverVerbose_) {
- OpmLog::error(msg);
- }
- // Use throw directly to prevent file and line
- throw TimeSteppingBreakdown{msg};
+ this->substep_timer_.setLastStepFailed(true);
+ checkTimeStepMaxRestartLimit_(restarts);
+
+ const double new_time_step = restartFactor_() * dt;
+ checkTimeStepMinLimit_(new_time_step);
+ bool wells_shut = false;
+ if (new_time_step > minTimeStepBeforeClosingWells_()) {
+ chopTimeStep_(new_time_step);
+ } else {
+ wells_shut = chopTimeStepOrCloseFailingWells_(new_time_step);
+ }
+ if (wells_shut) {
+ setTimeStep_(dt); // retry the old timestep
}
+ else {
+ restarts++; // only increase if no wells were shut
+ }
+ }
+ problem.setNextTimeStepSize(this->substep_timer_.currentStepLength());
+ }
+ updateSuggestedNextStep_();
+ return report;
+}
- // The new, chopped timestep.
- const double newTimeStep = restartFactor_ * dt;
+/************************************************
+ * Private class SubStepIteration private methods
+ ************************************************/
- // If we have restarted (i.e. cut the timestep) too
- // much, we have failed and throw an exception.
- if (newTimeStep < minTimeStep_) {
- const auto msg = fmt::format(
- "Solver failed to converge after cutting timestep to {}\n"
- "which is the minimum threshold given by option --solver-min-time-step\n",
- minTimeStep_
- );
- if (solverVerbose_) {
- OpmLog::error(msg);
- }
- // Use throw directly to prevent file and line
- throw TimeSteppingBreakdown{msg};
- }
- // Define utility function for chopping timestep.
- auto chopTimestep = [&]() {
- substepTimer.provideTimeStepEstimate(newTimeStep);
- if (solverVerbose_) {
- const auto msg = fmt::format(
- "{}\nTimestep chopped to {} days\n",
- causeOfFailure,
- unit::convert::to(substepTimer.currentStepLength(), unit::day)
- );
- OpmLog::problem(msg);
- }
- ++restarts;
- };
+template<class TypeTag>
+template<class Solver>
+bool
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+checkContinueOnUnconvergedSolution_(double dt) const
+{
+ bool continue_on_uncoverged_solution = ignoreConvergenceFailure_() && dt <= minTimeStep_();
+ if (continue_on_uncoverged_solution && solverVerbose_()) {
+ // NOTE: This method is only called if the solver failed to converge.
+ const auto msg = fmt::format(
+ "Solver failed to converge but timestep {} is smaller or equal to {}\n"
+ "which is the minimum threshold given by option --solver-min-time-step\n",
+ dt, minTimeStep_()
+ );
+ OpmLog::problem(msg);
+ }
+ return continue_on_uncoverged_solution;
+}
- const double minimumChoppedTimestep = minTimeStepBeforeShuttingProblematicWells_;
- if (newTimeStep > minimumChoppedTimestep) {
- chopTimestep();
- } else {
- // We are below the threshold, and will check if there are any
- // wells we should close rather than chopping again.
- std::set<std::string> failing_wells = detail::consistentlyFailingWells(solver.model().stepReports());
- if (failing_wells.empty()) {
- // Found no wells to close, chop the timestep as above.
- chopTimestep();
- } else {
- // Close all consistently failing wells that are not under group control
- std::vector<std::string> shut_wells;
- for (const auto& well : failing_wells) {
- bool was_shut = solver.model().wellModel().forceShutWellByName(
- well, substepTimer.simulationTimeElapsed(), /*dont_shut_grup_wells =*/ true);
- if (was_shut) {
- shut_wells.push_back(well);
- }
- }
- // If no wells are closed we also try to shut wells under group control
- if (shut_wells.empty()) {
- for (const auto& well : failing_wells) {
- bool was_shut = solver.model().wellModel().forceShutWellByName(
- well, substepTimer.simulationTimeElapsed(), /*dont_shut_grup_wells =*/ false);
- if (was_shut) {
- shut_wells.push_back(well);
- }
- }
- }
- // If still no wells are closed we must fall back to chopping again
- if (shut_wells.empty()) {
- chopTimestep();
- } else {
- substepTimer.provideTimeStepEstimate(dt);
- if (solverVerbose_) {
- const std::string msg =
- fmt::format("\nProblematic well(s) were shut: {}"
- "(retrying timestep)\n",
- fmt::join(shut_wells, " "));
- OpmLog::problem(msg);
- }
- }
+template<class TypeTag>
+template<class Solver>
+void
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+checkTimeStepMaxRestartLimit_(const int restarts) const
+{
+ // If we have restarted (i.e. cut the timestep) too
+ // many times, we have failed and throw an exception.
+ if (restarts >= solverRestartMax_()) {
+ const auto msg = fmt::format(
+ "Solver failed to converge after cutting timestep {} times.", restarts
+ );
+ if (solverVerbose_()) {
+ OpmLog::error(msg);
+ }
+ // Use throw directly to prevent file and line
+ throw TimeSteppingBreakdown{msg};
+ }
+}
+
+template<class TypeTag>
+template<class Solver>
+void
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+checkTimeStepMinLimit_(const int new_time_step) const
+{
+ // If we have restarted (i.e. cut the timestep) too
+ // much, we have failed and throw an exception.
+ if (new_time_step < minTimeStep_()) {
+ const auto msg = fmt::format(
+ "Solver failed to converge after cutting timestep to {}\n"
+ "which is the minimum threshold given by option --solver-min-time-step\n",
+ minTimeStep_()
+ );
+ if (solverVerbose_()) {
+ OpmLog::error(msg);
+ }
+ // Use throw directly to prevent file and line
+ throw TimeSteppingBreakdown{msg};
+ }
+}
+
+template<class TypeTag>
+template<class Solver>
+void
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+chopTimeStep_(const double new_time_step)
+{
+ setTimeStep_(new_time_step);
+ if (solverVerbose_()) {
+ const auto msg = fmt::format("{}\nTimestep chopped to {} days\n",
+ this->cause_of_failure_,
+ unit::convert::to(this->substep_timer_.currentStepLength(), unit::day));
+ OpmLog::problem(msg);
+ }
+}
+
+template<class TypeTag>
+template<class Solver>
+bool
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+chopTimeStepOrCloseFailingWells_(const int new_time_step)
+{
+ bool wells_shut = false;
+ // We are below the threshold, and will check if there are any wells we should close
+ // rather than chopping again.
+ std::set<std::string> failing_wells = detail::consistentlyFailingWells(
+ solver_().model().stepReports());
+ if (failing_wells.empty()) {
+ // Found no wells to close, chop the timestep
+ chopTimeStep_(new_time_step);
+ } else {
+ // Close all consistently failing wells that are not under group control
+ std::vector<std::string> shut_wells;
+ for (const auto& well : failing_wells) {
+ bool was_shut = solver_().model().wellModel().forceShutWellByName(
+ well, this->substep_timer_.simulationTimeElapsed(), /*dont_shut_grup_wells =*/ true);
+ if (was_shut) {
+ shut_wells.push_back(well);
+ }
+ }
+ // If no wells are closed we also try to shut wells under group control
+ if (shut_wells.empty()) {
+ for (const auto& well : failing_wells) {
+ bool was_shut = solver_().model().wellModel().forceShutWellByName(
+ well, this->substep_timer_.simulationTimeElapsed(), /*dont_shut_grup_wells =*/ false);
+ if (was_shut) {
+ shut_wells.push_back(well);
}
}
}
- problem.setNextTimeStepSize(substepTimer.currentStepLength());
+ // If still no wells are closed we must fall back to chopping again
+ if (shut_wells.empty()) {
+ chopTimeStep_(new_time_step);
+ } else {
+ wells_shut = true;
+ if (solverVerbose_()) {
+ const std::string msg =
+ fmt::format("\nProblematic well(s) were shut: {}"
+ "(retrying timestep)\n",
+ fmt::join(shut_wells, " "));
+ OpmLog::problem(msg);
+ }
+ }
}
+ return wells_shut;
+}
- // store estimated time step for next reportStep
- suggestedNextTimestep_ = substepTimer.currentStepLength();
- if (timestepVerbose_) {
- std::ostringstream ss;
- substepTimer.report(ss);
- ss << "Suggested next step size = " << unit::convert::to(suggestedNextTimestep_, unit::day) << " (days)" << std::endl;
- OpmLog::debug(ss.str());
- }
+template<class TypeTag>
+template<class Solver>
+boost::posix_time::ptime
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+currentDateTime_() const
+{
+ return simulatorTimer_().currentDateTime();
+}
- if (! std::isfinite(suggestedNextTimestep_)) { // check for NaN
- suggestedNextTimestep_ = timestep;
+template<class TypeTag>
+template<class Solver>
+int
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+getNumIterations_(const SimulatorReportSingle &substep_report) const
+{
+ if (useNewtonIteration_()) {
+ return substep_report.total_newton_iterations;
+ }
+ else {
+ return substep_report.total_linear_iterations;
}
- return report;
}
template<class TypeTag>
-void AdaptiveTimeStepping<TypeTag>::
-updateTUNING(double max_next_tstep, const Tuning& tuning)
+template<class Solver>
+double
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+growthFactor_() const
{
- restartFactor_ = tuning.TSFCNV;
- growthFactor_ = tuning.TFDIFF;
- maxGrowth_ = tuning.TSFMAX;
- maxTimeStep_ = tuning.TSMAXZ;
- updateNEXTSTEP(max_next_tstep);
- timestepAfterEvent_ = tuning.TMAXWC;
+ return this->adaptive_time_stepping_.growth_factor_;
}
template<class TypeTag>
-void AdaptiveTimeStepping<TypeTag>::
-updateNEXTSTEP(double max_next_tstep)
+template<class Solver>
+bool
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+ignoreConvergenceFailure_() const
{
- // \Note Only update next suggested step if TSINIT was explicitly set in TUNING or NEXTSTEP is active.
- if (max_next_tstep > 0) {
- suggestedNextTimestep_ = max_next_tstep;
+ return adaptive_time_stepping_.ignore_convergence_failure_;
+}
+
+template<class TypeTag>
+template<class Solver>
+double
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+maxGrowth_() const
+{
+ return this->adaptive_time_stepping_.max_growth_;
+}
+
+template<class TypeTag>
+template<class Solver>
+void
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+maybeReportSubStep_(SimulatorReportSingle substep_report) const
+{
+ if (timeStepVerbose_()) {
+ std::ostringstream ss;
+ substep_report.reportStep(ss);
+ OpmLog::info(ss.str());
}
}
template<class TypeTag>
-template<class Serializer>
-void AdaptiveTimeStepping<TypeTag>::
-serializeOp(Serializer& serializer)
+template<class Solver>
+double
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+maybeRestrictTimeStepGrowth_(const double dt, double dt_estimate, const int restarts) const
{
- serializer(timeStepControlType_);
- switch (timeStepControlType_) {
- case TimeStepControlType::HardCodedTimeStep:
- allocAndSerialize<HardcodedTimeStepControl>(serializer);
- break;
- case TimeStepControlType::PIDAndIterationCount:
- allocAndSerialize<PIDAndIterationCountTimeStepControl>(serializer);
- break;
- case TimeStepControlType::SimpleIterationCount:
- allocAndSerialize<SimpleIterationCountTimeStepControl>(serializer);
- break;
- case TimeStepControlType::PID:
- allocAndSerialize<PIDTimeStepControl>(serializer);
- break;
+ // limit the growth of the timestep size by the growth factor
+ dt_estimate = std::min(dt_estimate, double(maxGrowth_() * dt));
+ assert(dt_estimate > 0);
+ // further restrict time step size growth after convergence problems
+ if (restarts > 0) {
+ dt_estimate = std::min(growthFactor_() * dt, dt_estimate);
}
- serializer(restartFactor_);
- serializer(growthFactor_);
- serializer(maxGrowth_);
- serializer(maxTimeStep_);
- serializer(minTimeStep_);
- serializer(ignoreConvergenceFailure_);
- serializer(solverRestartMax_);
- serializer(solverVerbose_);
- serializer(timestepVerbose_);
- serializer(suggestedNextTimestep_);
- serializer(fullTimestepInitially_);
- serializer(timestepAfterEvent_);
- serializer(useNewtonIteration_);
- serializer(minTimeStepBeforeShuttingProblematicWells_);
+
+ return dt_estimate;
}
+// The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicitBlackoil::runStep()
+// It has to be called for each substep since TUNING might have been changed for next sub step due
+// to ACTIONX (via NEXTSTEP) or WCYCLE keywords.
template<class TypeTag>
-AdaptiveTimeStepping<TypeTag>
-AdaptiveTimeStepping<TypeTag>::
-serializationTestObjectHardcoded()
+template<class Solver>
+void
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+maybeUpdateTuningAndTimeStep_()
{
- return serializationTestObject_<HardcodedTimeStepControl>();
+ // TODO: This function is currently only called if NEXTSTEP is activated from ACTIONX or
+ // if the WCYCLE keyword needs to modify the current timestep. So this method should rather
+ // be named maybeUpdateTimeStep_() or similar, since it should not update the tuning. However,
+ // the current definition of the maybeUpdateTuning_() callback is actually calling
+ // adaptiveTimeStepping_->updateTUNING(max_next_tstep, tuning) which is updating the tuning
+ // see SimulatorFullyImplicitBlackoil::runStep() for more details.
+ auto old_value = suggestedNextTimestep_();
+ if (this->substepper_.maybeUpdateTuning_(this->substep_timer_.simulationTimeElapsed(),
+ this->substep_timer_.currentStepLength(),
+ this->substep_timer_.currentStepNum()))
+ {
+ // Either NEXTSTEP and WCYCLE wants to change the current time step, but they cannot
+ // change the current time step directly. Instead, they change the suggested next time step
+ // by calling updateNEXTSTEP() via the maybeUpdateTuning() callback. We now need to update
+ // the current time step to the new suggested time step and reset the suggested time step
+ // to the old value.
+ setTimeStep_(suggestedNextTimestep_());
+ setSuggestedNextStep_(old_value);
+ }
}
template<class TypeTag>
-AdaptiveTimeStepping<TypeTag>
-AdaptiveTimeStepping<TypeTag>::
-serializationTestObjectPID()
+template<class Solver>
+double
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+minTimeStepBeforeClosingWells_() const
{
- return serializationTestObject_<PIDTimeStepControl>();
+ return this->adaptive_time_stepping_.min_time_step_before_shutting_problematic_wells_;
}
template<class TypeTag>
-AdaptiveTimeStepping<TypeTag>
-AdaptiveTimeStepping<TypeTag>::
-serializationTestObjectPIDIt()
+template<class Solver>
+double
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+minTimeStep_() const
{
- return serializationTestObject_<PIDAndIterationCountTimeStepControl>();
+ return this->adaptive_time_stepping_.min_time_step_;
}
template<class TypeTag>
-AdaptiveTimeStepping<TypeTag>
-AdaptiveTimeStepping<TypeTag>::
-serializationTestObjectSimple()
+template<class Solver>
+double
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+restartFactor_() const
{
- return serializationTestObject_<SimpleIterationCountTimeStepControl>();
+ return this->adaptive_time_stepping_.restart_factor_;
}
template<class TypeTag>
-bool
-AdaptiveTimeStepping<TypeTag>::
-operator==(const AdaptiveTimeStepping<TypeTag>& rhs) const
+template<class Solver>
+SimulatorReportSingle
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+runSubStep_()
{
- if (timeStepControlType_ != rhs.timeStepControlType_ ||
- (timeStepControl_ && !rhs.timeStepControl_) ||
- (!timeStepControl_ && rhs.timeStepControl_)) {
- return false;
+ SimulatorReportSingle substep_report;
+
+ auto handleFailure = [this, &substep_report]
+ (const std::string& failure_reason, const std::exception& e, bool log_exception = true)
+ {
+ substep_report = solver_().failureReport();
+ this->cause_of_failure_ = failure_reason;
+ if (log_exception && solverVerbose_()) {
+ std::string message;
+ message = "Caught Exception: ";
+ message += e.what();
+ OpmLog::debug(message);
+ }
+ };
+
+ try {
+ substep_report = solver_().step(this->substep_timer_);
+ if (solverVerbose_()) {
+ // report number of linear iterations
+ OpmLog::debug("Overall linear iterations used: "
+ + std::to_string(substep_report.total_linear_iterations));
+ }
}
-
- bool result = false;
- switch (timeStepControlType_) {
- case TimeStepControlType::HardCodedTimeStep:
- result = castAndComp<HardcodedTimeStepControl>(rhs);
- break;
- case TimeStepControlType::PIDAndIterationCount:
- result = castAndComp<PIDAndIterationCountTimeStepControl>(rhs);
- break;
- case TimeStepControlType::SimpleIterationCount:
- result = castAndComp<SimpleIterationCountTimeStepControl>(rhs);
- break;
- case TimeStepControlType::PID:
- result = castAndComp<PIDTimeStepControl>(rhs);
- break;
+ catch (const TooManyIterations& e) {
+ handleFailure("Solver convergence failure - Iteration limit reached", e);
+ }
+ catch (const ConvergenceMonitorFailure& e) {
+ handleFailure("Convergence monitor failure", e, /*log_exception=*/false);
+ }
+ catch (const LinearSolverProblem& e) {
+ handleFailure("Linear solver convergence failure", e);
+ }
+ catch (const NumericalProblem& e) {
+ handleFailure("Solver convergence failure - Numerical problem encountered", e);
+ }
+ catch (const std::runtime_error& e) {
+ handleFailure("Runtime error encountered", e);
+ }
+ catch (const Dune::ISTLError& e) {
+ handleFailure("ISTL error - Time step too large", e);
+ }
+ catch (const Dune::MatrixBlockError& e) {
+ handleFailure("Matrix block error", e);
}
- return result
- && this->restartFactor_ == rhs.restartFactor_
- && this->growthFactor_ == rhs.growthFactor_
- && this->maxGrowth_ == rhs.maxGrowth_
- && this->maxTimeStep_ == rhs.maxTimeStep_
- && this->minTimeStep_ == rhs.minTimeStep_
- && this->ignoreConvergenceFailure_ == rhs.ignoreConvergenceFailure_
- && this->solverRestartMax_== rhs.solverRestartMax_
- && this->solverVerbose_ == rhs.solverVerbose_
- && this->fullTimestepInitially_ == rhs.fullTimestepInitially_
- && this->timestepAfterEvent_ == rhs.timestepAfterEvent_
- && this->useNewtonIteration_ == rhs.useNewtonIteration_
- && this->minTimeStepBeforeShuttingProblematicWells_ ==
- rhs.minTimeStepBeforeShuttingProblematicWells_;
+ return substep_report;
}
template<class TypeTag>
-template<class Controller>
-AdaptiveTimeStepping<TypeTag>
-AdaptiveTimeStepping<TypeTag>::
-serializationTestObject_()
+template<class Solver>
+void
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+setTimeStep_(double dt_estimate)
{
- AdaptiveTimeStepping<TypeTag> result;
+ this->substep_timer_.provideTimeStepEstimate(dt_estimate);
+}
+
+template<class TypeTag>
+template<class Solver>
+Solver&
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+solver_() const
+{
+ return this->substepper_.solver_;
+}
- result.restartFactor_ = 1.0;
- result.growthFactor_ = 2.0;
- result.maxGrowth_ = 3.0;
- result.maxTimeStep_ = 4.0;
- result.minTimeStep_ = 5.0;
- result.ignoreConvergenceFailure_ = true;
- result.solverRestartMax_ = 6;
- result.solverVerbose_ = true;
- result.timestepVerbose_ = true;
- result.suggestedNextTimestep_ = 7.0;
- result.fullTimestepInitially_ = true;
- result.useNewtonIteration_ = true;
- result.minTimeStepBeforeShuttingProblematicWells_ = 9.0;
- result.timeStepControlType_ = Controller::Type;
- result.timeStepControl_ = std::make_unique<Controller>(Controller::serializationTestObject());
- return result;
+template<class TypeTag>
+template<class Solver>
+int
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+solverRestartMax_() const
+{
+ return this->adaptive_time_stepping_.solver_restart_max_;
}
template<class TypeTag>
-void AdaptiveTimeStepping<TypeTag>::
-init_(const UnitSystem& unitSystem)
+template<class Solver>
+void
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+setSuggestedNextStep_(double step)
{
- std::tie(timeStepControlType_,
- timeStepControl_,
- useNewtonIteration_) = detail::createController(unitSystem);
+ this->adaptive_time_stepping_.setSuggestedNextStep(step);
+}
- // make sure growth factor is something reasonable
- if (growthFactor_ < 1.0) {
- OPM_THROW(std::runtime_error,
- "Growth factor cannot be less than 1.");
+template <class TypeTag>
+template <class Solver>
+const SimulatorTimer&
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+simulatorTimer_() const
+{
+ return this->substepper_.simulator_timer_;
+}
+
+template <class TypeTag>
+template <class Solver>
+bool
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+solverVerbose_() const
+{
+ return this->adaptive_time_stepping_.solver_verbose_;
+}
+
+template<class TypeTag>
+template<class Solver>
+boost::posix_time::ptime
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+startDateTime_() const
+{
+ return simulatorTimer_().startDateTime();
+}
+
+template <class TypeTag>
+template <class Solver>
+double
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+suggestedNextTimestep_() const
+{
+ return this->adaptive_time_stepping_.suggestedNextStep();
+}
+
+template <class TypeTag>
+template <class Solver>
+double
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+timeStepControlComputeEstimate_(const double dt, const int iterations, double elapsed) const
+{
+ // create object to compute the time error, simply forwards the call to the model
+ SolutionTimeErrorSolverWrapper<Solver> relative_change{solver_()};
+ return this->adaptive_time_stepping_.time_step_control_->computeTimeStepSize(
+ dt, iterations, relative_change, elapsed);
+}
+
+template <class TypeTag>
+template <class Solver>
+bool
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+timeStepVerbose_() const
+{
+ return this->adaptive_time_stepping_.timestep_verbose_;
+}
+
+// The suggested time step is the stepsize that will be used as a first try for
+// the next sub step. It is updated at the end of each substep. It can also be
+// updated by the TUNING or NEXTSTEP keywords at the beginning of each report step or
+// at the beginning of each substep by the ACTIONX keyword (via NEXTSTEP), this is
+// done by the maybeUpdateTuning_() method which is called at the beginning of each substep
+// (and the begginning of each report step). Note that the WCYCLE keyword can also update the
+// suggested time step via the maybeUpdateTuning_() method.
+template <class TypeTag>
+template <class Solver>
+void
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+updateSuggestedNextStep_()
+{
+ auto suggested_next_step = this->substep_timer_.currentStepLength();
+ if (! std::isfinite(suggested_next_step)) { // check for NaN
+ suggested_next_step = this->original_time_step_;
}
+ if (timeStepVerbose_()) {
+ std::ostringstream ss;
+ this->substep_timer_.report(ss);
+ ss << "Suggested next step size = "
+ << unit::convert::to(suggested_next_step, unit::day) << " (days)" << std::endl;
+ OpmLog::debug(ss.str());
+ }
+ setSuggestedNextStep_(suggested_next_step);
+}
+
+template <class TypeTag>
+template <class Solver>
+bool
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+useNewtonIteration_() const
+{
+ return this->adaptive_time_stepping_.use_newton_iteration_;
+}
+
+template <class TypeTag>
+template <class Solver>
+double
+AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
+writeOutput_() const
+{
+ time::StopWatch perf_timer;
+ perf_timer.start();
+ auto& problem = solver_().model().simulator().problem();
+ problem.writeOutput(true);
+ return perf_timer.secsSinceStart();
+}
+
+/************************************************
+ * Private class SolutionTimeErrorSolverWrapper
+ * **********************************************/
+
+template<class TypeTag>
+template<class Solver>
+AdaptiveTimeStepping<TypeTag>::SolutionTimeErrorSolverWrapper<Solver>::SolutionTimeErrorSolverWrapper(
+ const Solver& solver
+)
+ : solver_{solver}
+{}
+
+template<class TypeTag>
+template<class Solver>
+double AdaptiveTimeStepping<TypeTag>::SolutionTimeErrorSolverWrapper<Solver>::relativeChange() const
+{
+ // returns: || u^n+1 - u^n || / || u^n+1 ||
+ return solver_.model().relativeChange();
}
} // namespace Opm
-#endif
+#endif // OPM_ADAPTIVE_TIME_STEPPING_IMPL_HPP
diff --git a/opm/simulators/timestepping/SimulatorTimerInterface.hpp b/opm/simulators/timestepping/SimulatorTimerInterface.hpp
index 6875db8711d..3948609bffe 100644
--- a/opm/simulators/timestepping/SimulatorTimerInterface.hpp
+++ b/opm/simulators/timestepping/SimulatorTimerInterface.hpp
@@ -40,14 +40,15 @@ namespace Opm
/// destructor
virtual ~SimulatorTimerInterface() {}
- /// Current step number. This is the number of timesteps that
- /// has been completed from the start of the run. The time
- /// after initialization but before the simulation has started
- /// is timestep number zero.
- virtual int currentStepNum() const = 0;
+ // -----------------------------------------------------------
+ // Pure virtual functions to be implemented by derived classes
+ // -----------------------------------------------------------
- /// Current report step number. This might differ from currentStepNum in case of sub stepping
- virtual int reportStepNum() const { return currentStepNum(); }
+ /// advance time by currentStepLength
+ virtual void advance() = 0 ;
+
+ /// return copy of current timer instance
+ virtual std::unique_ptr< SimulatorTimerInterface > clone () const = 0;
/// Current step length. This is the length of the step
/// the simulator will take in the next iteration.
@@ -55,6 +56,25 @@ namespace Opm
/// @note if done(), it is an error to call currentStepLength().
virtual double currentStepLength() const = 0;
+ /// Current step number. This is the number of timesteps that
+ /// has been completed from the start of the run. The time
+ /// after initialization but before the simulation has started
+ /// is timestep number zero.
+ virtual int currentStepNum() const = 0;
+
+ /// Return true if timer indicates that simulation of timer interval is finished
+ virtual bool done() const = 0;
+
+ /// Whether the current step is the first step.
+ virtual bool initialStep() const = 0;
+
+ /// Return true if last time step failed
+ virtual bool lastStepFailed() const = 0;
+
+ /// Time elapsed since the start of the simulation until the
+ /// beginning of the current time step [s].
+ virtual double simulationTimeElapsed() const = 0;
+
/// Previous step length. This is the length of the step that
/// was taken to arrive at this time.
///
@@ -63,30 +83,13 @@ namespace Opm
/// it is an error to call stepLengthTaken().
virtual double stepLengthTaken () const = 0;
- /// Previous report step length. This is the length of the step that
- /// was taken to arrive at this report step time.
- ///
- /// @note if no increments have been done (i.e. the timer is
- /// still in its constructed state and reportStepNum() == 0),
- /// it is an error to call stepLengthTaken().
- virtual double reportStepLengthTaken () const { return stepLengthTaken(); }
-
- /// Time elapsed since the start of the simulation until the
- /// beginning of the current time step [s].
- virtual double simulationTimeElapsed() const = 0;
-
- /// advance time by currentStepLength
- virtual void advance() = 0 ;
-
- /// Return true if timer indicates that simulation of timer interval is finished
- virtual bool done() const = 0;
-
- /// Whether the current step is the first step.
- virtual bool initialStep() const = 0;
-
/// Return start date of simulation
virtual boost::posix_time::ptime startDateTime() const = 0;
+ // -----------------------------------------------------------------
+ // Virtual functions (not pure) to allow for default implementations
+ // -----------------------------------------------------------------
+
/// Return the current time as a posix time object.
virtual boost::posix_time::ptime currentDateTime() const;
@@ -94,11 +97,17 @@ namespace Opm
/// 1970) until the current time step begins [s].
virtual time_t currentPosixTime() const;
- /// Return true if last time step failed
- virtual bool lastStepFailed() const = 0;
+ /// Previous report step length. This is the length of the step that
+ /// was taken to arrive at this report step time.
+ ///
+ /// @note if no increments have been done (i.e. the timer is
+ /// still in its constructed state and reportStepNum() == 0),
+ /// it is an error to call stepLengthTaken().
+ virtual double reportStepLengthTaken () const { return stepLengthTaken(); }
+
+ /// Current report step number. This might differ from currentStepNum in case of sub stepping
+ virtual int reportStepNum() const { return currentStepNum(); }
- /// return copy of current timer instance
- virtual std::unique_ptr< SimulatorTimerInterface > clone () const = 0;
};
diff --git a/opm/simulators/utils/readDeck.cpp b/opm/simulators/utils/readDeck.cpp
index bdb16c75a7c..aae46fcd094 100644
--- a/opm/simulators/utils/readDeck.cpp
+++ b/opm/simulators/utils/readDeck.cpp
@@ -54,7 +54,11 @@
#include <opm/input/eclipse/Schedule/Action/State.hpp>
#include <opm/input/eclipse/Schedule/ArrayDimChecker.hpp>
+#include <opm/input/eclipse/Schedule/ResCoup/ReservoirCouplingInfo.hpp>
+#include <opm/input/eclipse/Schedule/ResCoup/MasterGroup.hpp>
+#include <opm/input/eclipse/Schedule/ResCoup/Slaves.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
+
#include <opm/input/eclipse/Schedule/UDQ/UDQConfig.hpp>
#include <opm/input/eclipse/Schedule/UDQ/UDQState.hpp>
#include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>
@@ -151,7 +155,7 @@ namespace {
if (schedule == nullptr) {
schedule = std::make_shared<Opm::Schedule>
(deck, eclipseState, parseContext, errorGuard,
- std::move(python), lowActionParsingStrictness, /*slave_mode=*/false,
+ std::move(python), lowActionParsingStrictness, /*slaveMode=*/false,
keepKeywords, outputInterval, init_state);
}
@@ -180,14 +184,14 @@ namespace {
std::unique_ptr<Opm::UDQState>& udqState,
std::unique_ptr<Opm::Action::State>& actionState,
std::unique_ptr<Opm::WellTestState>& wtestState,
- Opm::ErrorGuard& errorGuard)
+ Opm::ErrorGuard& errorGuard,
+ const bool slaveMode)
{
if (schedule == nullptr) {
schedule = std::make_shared<Opm::Schedule>
(deck, eclipseState, parseContext,
- errorGuard, std::move(python), lowActionParsingStrictness, keepKeywords);
+ errorGuard, std::move(python), lowActionParsingStrictness, slaveMode, keepKeywords);
}
-
udqState = std::make_unique<Opm::UDQState>
((*schedule)[0].udq().params().undefinedValue());
@@ -233,6 +237,24 @@ namespace {
#endif
}
+ void inconsistentScheduleError(const std::string& message)
+ {
+ OPM_THROW(std::logic_error,
+ fmt::format("Inconsistent SCHEDULE section: {}", message));
+ }
+
+ void checkScheduleKeywordConsistency(const Opm::Schedule& schedule)
+ {
+ const auto& final_state = schedule.back();
+ const auto& rescoup = final_state.rescoup();
+ if (rescoup.slaveCount() > 0 && rescoup.masterGroupCount() == 0) {
+ inconsistentScheduleError("SLAVES keyword without GRUPMAST keyword");
+ }
+ if (rescoup.slaveCount() == 0 && rescoup.masterGroupCount() > 0) {
+ inconsistentScheduleError("GRUPMAST keyword without SLAVES keyword");
+ }
+ }
+
void readOnIORank(Opm::Parallel::Communication comm,
const std::string& deckFilename,
const Opm::ParseContext* parseContext,
@@ -249,7 +271,8 @@ namespace {
const bool lowActionParsingStrictness,
const bool keepKeywords,
const std::optional<int>& outputInterval,
- Opm::ErrorGuard& errorGuard)
+ Opm::ErrorGuard& errorGuard,
+ const bool slaveMode)
{
OPM_TIMEBLOCK(readDeck);
if (((schedule == nullptr) || (summaryConfig == nullptr)) &&
@@ -282,9 +305,10 @@ namespace {
lowActionParsingStrictness, keepKeywords,
std::move(python),
schedule, udqState, actionState, wtestState,
- errorGuard);
+ errorGuard, slaveMode);
}
+ checkScheduleKeywordConsistency(*schedule);
eclipseState->appendAqufluxSchedule(schedule->getAquiferFluxSchedule());
if (Opm::OpmLog::hasBackend("STDOUT_LOGGER")) {
@@ -480,7 +504,6 @@ Opm::setupLogging(Parallel::Communication& comm,
}
logFileStream << ".PRT";
debugFileStream << ".DBG";
-
FileOutputMode output;
{
static std::map<std::string, FileOutputMode> stringToOutputMode =
@@ -539,10 +562,10 @@ void Opm::readDeck(Opm::Parallel::Communication comm,
const bool initFromRestart,
const bool checkDeck,
const bool keepKeywords,
- const std::optional<int>& outputInterval)
+ const std::optional<int>& outputInterval,
+ const bool slaveMode)
{
auto errorGuard = std::make_unique<ErrorGuard>();
-
int parseSuccess = 1; // > 0 is success
std::string failureMessage;
@@ -573,7 +596,7 @@ void Opm::readDeck(Opm::Parallel::Communication comm,
eclipseState, schedule, udqState, actionState, wtestState,
summaryConfig, std::move(python), initFromRestart,
checkDeck, treatCriticalAsNonCritical, lowActionParsingStrictness,
- keepKeywords, outputInterval, *errorGuard);
+ keepKeywords, outputInterval, *errorGuard, slaveMode);
// Update schedule so that re-parsing after actions use same strictness
assert(schedule);
diff --git a/opm/simulators/utils/readDeck.hpp b/opm/simulators/utils/readDeck.hpp
index 87da15315fe..dc7692ad471 100644
--- a/opm/simulators/utils/readDeck.hpp
+++ b/opm/simulators/utils/readDeck.hpp
@@ -99,7 +99,8 @@ void readDeck(Parallel::Communication comm,
bool initFromRestart,
bool checkDeck,
bool keepKeywords,
- const std::optional<int>& outputInterval);
+ const std::optional<int>& outputInterval,
+ bool slaveMode);
void verifyValidCellGeometry(Parallel::Communication comm,
const EclipseState& eclipseState);