Reservoir coupling: Implement time stepping

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

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

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,15 +368,23 @@ 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;
         }
 
         // 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;
         }
 

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

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

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