OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5620 from hakonhagland/timestepping2

Browse Source 
Reservoir coupling: Implement time stepping
This commit is contained in:
Markus Blatt 2025-01-21 14:38:24 +01:00 committed by GitHub
commit 990c3f0248
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
22 changed files with 2958 additions and 623 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
CMakeLists_files.cmake
View File

@ -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

2
opm/simulators/flow/FlowGenericVanguard.cpp
View File

@ -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();
}

12
opm/simulators/flow/FlowMain.hpp
View File

@ -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;
}

64
opm/simulators/flow/Main.cpp
View File

@ -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

8
opm/simulators/flow/Main.hpp
View File

@ -150,7 +150,8 @@ public:
~Main();
void setArgvArgc_(const std::string& filename);
void maybeSaveReservoirCouplingSlaveLogFilename_();
void maybeRedirectReservoirCouplingSlaveOutput_();
void initMPI();
int runDynamic()
@ -413,6 +414,7 @@ protected:
keepKeywords,
getNumThreads(),
Parameters::Get<Parameters::EclOutputInterval>(),
Parameters::Get<Parameters::Slave>(),
cmdline_params,
Opm::moduleVersion(),
Opm::compileTimestamp());
@ -697,6 +699,7 @@ private:
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 @@ private:
// 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

116
opm/simulators/flow/ReservoirCoupling.cpp Normal file
View File

@ -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

113
opm/simulators/flow/ReservoirCoupling.hpp Normal file
View File

@ -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

217
opm/simulators/flow/ReservoirCouplingMaster.cpp Normal file
View File

@ -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(
&timestep,
/*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

91
opm/simulators/flow/ReservoirCouplingMaster.hpp Normal file
View File

@ -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

276
opm/simulators/flow/ReservoirCouplingSlave.cpp Normal file
View File

@ -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(
&timestep,
/*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(&timestep, /*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

65
opm/simulators/flow/ReservoirCouplingSlave.hpp Normal file
View File

@ -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

335
opm/simulators/flow/ReservoirCouplingSpawnSlaves.cpp Normal file
View File

@ -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

72
opm/simulators/flow/ReservoirCouplingSpawnSlaves.hpp Normal file
View File

@ -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

4
opm/simulators/flow/SimulatorFullyImplicitBlackoil.cpp
View File

@ -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

105
opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp
View File

@ -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 @@ public:
/// \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 @@ public:
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 @@ public:
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 @@ public:
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 @@ protected:
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_;
};

31
opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp
View File

@ -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

15
opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp
View File

@ -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_;
};

337
opm/simulators/timestepping/AdaptiveTimeStepping.hpp
View File

@ -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;
SolutionTimeErrorSolverWrapper(const Solver& solver);
double relativeChange() const;
private:
const Solver& solver_;
};
//! \brief contructor taking parameter object
explicit AdaptiveTimeStepping(const UnitSystem& unitSystem,
const double max_next_tstep = -1.0,
const bool terminalOutput = true);
// 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;
//! \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);
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
);
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; }
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();
bool operator==(const AdaptiveTimeStepping<TypeTag>& rhs) const;
AdaptiveTimeStepping<TypeTag>& getAdaptiveTimerStepper();
SimulatorReport run();
friend class AdaptiveTimeStepping<TypeTag>::template SubStepIteration<Solver>;
private:
template<class Controller>
static AdaptiveTimeStepping<TypeTag> serializationTestObject_();
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;
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
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_;
};
}
template <class Solver>
class SubStepIteration {
public:
SubStepIteration(
SubStepper<Solver>& substepper,
AdaptiveSimulatorTimer& substep_timer,
const double original_time_step,
const bool final_step
);
SimulatorReport run();
private:
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_;
};
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

1584
opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp
View File

File diff suppressed because it is too large Load Diff

73
opm/simulators/timestepping/SimulatorTimerInterface.hpp
View File

@ -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;
};

43
opm/simulators/utils/readDeck.cpp
View File

@ -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);

3
opm/simulators/utils/readDeck.hpp
View File

@ -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);