Prepare to Save/Restore Richer Set of Aquifer Values

This commit adds logic to communicate more dynamic aquifer values
between the simulation and I/O layers.  In particular, we ensure
that we allocate the 'aquFet' and 'aquCT' substructures of the
dynamic aquifer data as appropriate and that we collect this
information on the I/O rank as part of the restart output process.
We furthermore make the 'ParallelRestart' facility aware of dynamic
aquifer data in preparation of loading these values from the restart
file.

ebos/eclbaseaquifermodel.hh

@@ -72,7 +72,7 @@ public:
      *        aquifer pressure and the base run's total produced liquid
      *        volume from the model's aquifers.
      */
-    void initFromRestart(const std::vector<data::AquiferData>& aquiferSoln OPM_UNUSED)
+    void initFromRestart(const data::Aquifers& aquiferSoln OPM_UNUSED)
     {
         throw std::logic_error {
             "Initialization from restart data not supported "

ebos/eclgenericwriter.cc

@@ -392,6 +392,7 @@ doWriteOutput(const int                     reportStepNum,
               data::Solution&&              localCellData,
               data::Wells&&                 localWellData,
               data::GroupAndNetworkValues&& localGroupAndNetworkData,
+              data::Aquifers&&              localAquiferData,
               const Action::State& actionState,
               const UDQState& udqState,
               const SummaryState& summaryState,
@@ -410,7 +411,10 @@ doWriteOutput(const int                     reportStepNum,
                    : std::move(localWellData),
 
         isParallel ? this->collectToIORank_.globalGroupAndNetworkData()
-                   : std::move(localGroupAndNetworkData)
+                   : std::move(localGroupAndNetworkData),
+
+        isParallel ? this->collectToIORank_.globalAquiferData()
+                   : std::move(localAquiferData)
     };
 
     if (eclState_.getSimulationConfig().useThresholdPressure()) {

ebos/eclgenericwriter.hh

@@ -82,6 +82,7 @@ protected:
                        data::Solution&&              localCellData,
                        data::Wells&&                 localWellData,
                        data::GroupAndNetworkValues&& localGroupAndNetworkData,
+                       data::Aquifers&&              localAquiferData,
                        const Action::State& actionState,
                        const UDQState& udqState,
                        const SummaryState& summaryState,

ebos/eclproblem.hh

@@ -1885,6 +1885,9 @@ public:
     EclWellModel& wellModel()
     { return wellModel_; }
 
+    const EclAquiferModel& aquiferModel() const
+    { return aquiferModel_; }
+
     EclAquiferModel& mutableAquiferModel()
     { return aquiferModel_; }
 

ebos/eclwriter.hh

@@ -141,7 +141,6 @@ public:
     /*!
      * \brief collect and pass data and pass it to eclIO writer
      */
-
     void evalSummaryState(bool isSubStep)
     {
         const int reportStepNum = simulator_.episodeIndex() + 1;
@@ -176,7 +175,7 @@ public:
             .groupAndNetworkData(reportStepNum, simulator_.vanguard().schedule());
 
 
-        const auto localAquiferData = simulator_.problem().mutableAquiferModel().aquiferData();
+        const auto localAquiferData = simulator_.problem().aquiferModel().aquiferData();
 
         this->prepareLocalCellData(isSubStep, reportStepNum);
 
@@ -221,7 +220,6 @@ public:
                           eclOutputModule_->initialInplace());
     }
 
-
     void writeOutput(bool isSubStep)
     {
         const int reportStepNum = simulator_.episodeIndex() + 1;
@@ -234,6 +232,8 @@ public:
         auto localGroupAndNetworkData = simulator_.problem().wellModel()
             .groupAndNetworkData(reportStepNum, simulator_.vanguard().schedule());
 
+        auto localAquiferData = simulator_.problem().aquiferModel().aquiferData();
+
         data::Solution localCellData = {};
         if (! isSubStep) {
             this->eclOutputModule_->assignToSolution(localCellData);
@@ -248,7 +248,7 @@ public:
                                            eclOutputModule_->getWBPData(),
                                            localWellData,
                                            localGroupAndNetworkData,
-                                           {});
+                                           localAquiferData);
         }
 
         if (this->collectToIORank_.isIORank()) {
@@ -258,6 +258,7 @@ public:
                                 std::move(localCellData),
                                 std::move(localWellData),
                                 std::move(localGroupAndNetworkData),
+                                std::move(localAquiferData),
                                 this->actionState(),
                                 this->udqState(),
                                 this->summaryState(),

opm/simulators/aquifers/AquiferFetkovich.hpp

@@ -83,8 +83,9 @@ public:
         data.volume = this->W_flux_.value();
         data.initPressure = this->pa0_;
         data.type = data::AquiferType::Fetkovich;
-        // Not handling std::shared_ptr<FetkovichData> aquFet for now,
+
-        // because we do not need it yet
+        data.aquFet = std::make_shared<data::FetkovichData>();
+
         return data;
     }
 

opm/simulators/aquifers/AquiferInterface.hpp

@@ -88,19 +88,16 @@ public:
     {
     }
 
-    void initFromRestart(const std::vector<data::AquiferData>& aquiferSoln)
+    void initFromRestart(const data::Aquifers& aquiferSoln)
     {
-        auto xaqPos
+        auto xaqPos = aquiferSoln.find(this->aquiferID());
-            = std::find_if(aquiferSoln.begin(), aquiferSoln.end(), [this](const data::AquiferData& xaq) -> bool {
-                   return xaq.aquiferID == this->aquiferID();
-              });
-
         if (xaqPos == aquiferSoln.end())
             return;
 
-        this->assignRestartData(*xaqPos);
+        this->assignRestartData(xaqPos->second);
-        this->W_flux_ = xaqPos->volume;
+
-        this->pa0_ = xaqPos->initPressure;
+        this->W_flux_ = xaqPos->second.volume;
+        this->pa0_ = xaqPos->second.initPressure;
         this->solution_set_from_restart_ = true;
     }
 

opm/simulators/aquifers/AquiferNumerical.hpp

@@ -72,7 +72,7 @@ public:
         }
     }
 
-    void initFromRestart([[maybe_unused]]const std::vector<data::AquiferData>& aquiferSoln)
+    void initFromRestart([[maybe_unused]]const data::Aquifers& aquiferSoln)
     {
         // NOT handling Restart for now
     }

opm/simulators/aquifers/BlackoilAquiferModel.hpp

@@ -87,7 +87,7 @@ public:
     explicit BlackoilAquiferModel(Simulator& simulator);
 
     void initialSolutionApplied();
-    void initFromRestart(const std::vector<data::AquiferData>& aquiferSoln);
+    void initFromRestart(const data::Aquifers& aquiferSoln);
 
     void beginEpisode();
     void beginTimeStep();

opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp

@@ -59,19 +59,21 @@ BlackoilAquiferModel<TypeTag>::initialSolutionApplied()
 
 template <typename TypeTag>
 void
-BlackoilAquiferModel<TypeTag>::initFromRestart(const std::vector<data::AquiferData>& aquiferSoln)
+BlackoilAquiferModel<TypeTag>::initFromRestart(const data::Aquifers& aquiferSoln)
 {
-    if (aquiferCarterTracyActive()) {
+    if (this->aquiferCarterTracyActive()) {
-        for (auto& aquifer : aquifers_CarterTracy) {
+        for (auto& aquifer : this->aquifers_CarterTracy) {
             aquifer.initFromRestart(aquiferSoln);
         }
     }
-    if (aquiferFetkovichActive()) {
+
-        for (auto& aquifer : aquifers_Fetkovich) {
+    if (this->aquiferFetkovichActive()) {
+        for (auto& aquifer : this->aquifers_Fetkovich) {
             aquifer.initFromRestart(aquiferSoln);
         }
     }
-    if (aquiferNumericalActive()) {
+
+    if (this->aquiferNumericalActive()) {
         for (auto& aquifer : this->aquifers_numerical) {
             aquifer.initFromRestart(aquiferSoln);
         }

opm/simulators/utils/ParallelRestart.cpp

@@ -22,9 +22,12 @@
 #endif
 
 #include "ParallelRestart.hpp"
-#include <ctime>
+#include <cassert>
 #include <cstring>
+#include <ctime>
+#include <memory>
 #include <dune/common/parallel/mpitraits.hh>
+#include <opm/output/data/Aquifer.hpp>
 #include <opm/output/data/Cells.hpp>
 #include <opm/output/data/Groups.hpp>
 #include <opm/output/data/GuideRateValue.hpp>
@@ -214,13 +217,37 @@ std::size_t packSize(const std::array<T,N>& data, Dune::MPIHelper::MPICommunicat
     return N*packSize(data[0], comm);
 }
 
+HANDLE_AS_POD(data::CarterTracyData)
 HANDLE_AS_POD(data::Connection)
 HANDLE_AS_POD(data::CurrentControl)
+HANDLE_AS_POD(data::FetkovichData)
 HANDLE_AS_POD(data::GroupConstraints)
 HANDLE_AS_POD(data::NodeData)
 HANDLE_AS_POD(data::Rates)
 HANDLE_AS_POD(data::Segment)
 
+std::size_t packSize(const data::AquiferData& data, Dune::MPIHelper::MPICommunicator comm)
+{
+    const auto type = 0ull;
+
+    const auto base = packSize(data.aquiferID, comm)
+        + packSize(data.pressure, comm)
+        + packSize(data.fluxRate, comm)
+        + packSize(data.volume, comm)
+        + packSize(data.initPressure, comm)
+        + packSize(data.datumDepth, comm)
+        + packSize(type, comm);
+
+    if (data.aquFet != nullptr) {
+        return base + packSize(*data.aquFet, comm);
+    }
+    else if (data.aquCT != nullptr) {
+        return base + packSize(*data.aquCT, comm);
+    }
+
+    return base;
+}
+
 std::size_t packSize(const data::GuideRateValue&, Dune::MPIHelper::MPICommunicator comm)
 {
     const auto nItem = static_cast<std::size_t>(data::GuideRateValue::Item::NumItems);
@@ -289,6 +316,7 @@ std::size_t packSize(const RestartValue& data, Dune::MPIHelper::MPICommunicator
     return packSize(data.solution, comm)
         +  packSize(data.wells, comm)
         +  packSize(data.grp_nwrk, comm)
+        +  packSize(data.aquifer, comm)
         +  packSize(data.extra, comm);
 }
 
@@ -485,6 +513,29 @@ void pack(const std::unordered_map<T1,T2,H,P,A>& data, std::vector<char>& buffer
     }
 }
 
+void pack(const data::AquiferData& data, std::vector<char>& buffer, int& position,
+          Dune::MPIHelper::MPICommunicator comm)
+{
+    const auto type =
+          (data.aquFet != nullptr)*(1ull << 0)
+        + (data.aquCT  != nullptr)*(1ull << 1);
+
+    pack(data.aquiferID, buffer, position, comm);
+    pack(data.pressure, buffer, position, comm);
+    pack(data.fluxRate, buffer, position, comm);
+    pack(data.volume, buffer, position, comm);
+    pack(data.initPressure, buffer, position, comm);
+    pack(data.datumDepth, buffer, position, comm);
+    pack(type, buffer, position, comm);
+
+    if (data.aquFet != nullptr) {
+        pack(*data.aquFet, buffer, position, comm);
+    }
+    else if (data.aquCT != nullptr) {
+        pack(*data.aquCT, buffer, position, comm);
+    }
+}
+
 void pack(const data::GuideRateValue& data, std::vector<char>& buffer, int& position,
           Dune::MPIHelper::MPICommunicator comm)
 {
@@ -582,6 +633,7 @@ void pack(const RestartValue& data, std::vector<char>& buffer, int& position,
     pack(data.solution, buffer, position, comm);
     pack(data.wells, buffer, position, comm);
     pack(data.grp_nwrk, buffer, position, comm);
+    pack(data.aquifer, buffer, position, comm);
     pack(data.extra, buffer, position, comm);
 }
 
@@ -813,6 +865,31 @@ void unpack(data::Well& data, std::vector<char>& buffer, int& position,
     unpack(data.guide_rates, buffer, position, comm);
 }
 
+void unpack(data::AquiferData& data, std::vector<char>& buffer, int& position,
+            Dune::MPIHelper::MPICommunicator comm)
+{
+    auto type = 0ull;
+
+    unpack(data.aquiferID, buffer, position, comm);
+    unpack(data.pressure, buffer, position, comm);
+    unpack(data.fluxRate, buffer, position, comm);
+    unpack(data.volume, buffer, position, comm);
+    unpack(data.initPressure, buffer, position, comm);
+    unpack(data.datumDepth, buffer, position, comm);
+    unpack(type, buffer, position, comm);
+
+    if (type == 1ull) {
+        data.type = data::AquiferType::Fetkovich;
+        data.aquFet = std::make_shared<data::FetkovichData>();
+        unpack(*data.aquFet, buffer, position, comm);
+    }
+    else if (type == 2ull) {
+        data.type = data::AquiferType::CarterTracy;
+        data.aquCT = std::make_shared<data::CarterTracyData>();
+        unpack(*data.aquCT, buffer, position, comm);
+    }
+}
+
 void unpack(data::GuideRateValue& data, std::vector<char>& buffer, int& position,
             Dune::MPIHelper::MPICommunicator comm)
 {
@@ -893,6 +970,7 @@ void unpack(RestartValue& data, std::vector<char>& buffer, int& position,
     unpack(data.solution, buffer, position, comm);
     unpack(data.wells, buffer, position, comm);
     unpack(data.grp_nwrk, buffer, position, comm);
+    unpack(data.aquifer, buffer, position, comm);
     unpack(data.extra, buffer, position, comm);
 }
 
@@ -964,6 +1042,7 @@ INSTANTIATE_PACK(std::map<std::string,std::map<std::pair<int,int>,int>>)
 INSTANTIATE_PACK(std::map<std::string,int>)
 INSTANTIATE_PACK(std::map<std::string,double>)
 INSTANTIATE_PACK(std::map<int,int>)
+INSTANTIATE_PACK(std::map<int,data::AquiferData>)
 INSTANTIATE_PACK(std::unordered_map<std::string,size_t>)
 INSTANTIATE_PACK(std::unordered_map<std::string,std::string>)
 INSTANTIATE_PACK(std::unordered_set<std::string>)
@@ -979,10 +1058,7 @@ RestartValue loadParallelRestart(const EclipseIO* eclIO, Action::State& actionSt
                                  Dune::CollectiveCommunication<Dune::MPIHelper::MPICommunicator> comm)
 {
 #if HAVE_MPI
-    data::Solution sol;
+    RestartValue restartValues{};
-    data::Wells wells;
-    data::GroupAndNetworkValues grp_nwrk;
-    RestartValue restartValues(sol, wells, grp_nwrk);
 
     if (eclIO)
     {

opm/simulators/utils/ParallelRestart.hpp

@@ -49,9 +49,12 @@ class RestartValue;
 
 namespace data
 {
+struct AquiferData;
+struct CarterTracyData;
 struct CellData;
 struct Connection;
 struct CurrentControl;
+struct FetkovichData;
 class GroupAndNetworkValues;
 struct GroupConstraints;
 struct GroupData;
@@ -331,9 +334,12 @@ void unpack(char* str, std::size_t length, std::vector<char>& buffer, int& posit
   void unpack(T& data, std::vector<char>& buffer, int& position, \
               Dune::MPIHelper::MPICommunicator comm);
 
+ADD_PACK_PROTOTYPES(data::AquiferData)
+ADD_PACK_PROTOTYPES(data::CarterTracyData)
 ADD_PACK_PROTOTYPES(data::CellData)
 ADD_PACK_PROTOTYPES(data::Connection)
 ADD_PACK_PROTOTYPES(data::CurrentControl)
+ADD_PACK_PROTOTYPES(data::FetkovichData)
 ADD_PACK_PROTOTYPES(data::Rates)
 ADD_PACK_PROTOTYPES(data::Segment)
 ADD_PACK_PROTOTYPES(data::Solution)

tests/test_ParallelRestart.cpp

@@ -24,6 +24,7 @@
 
 #include <boost/test/unit_test.hpp>
 
+#include <memory>
 #include <tuple>
 #include <utility>
 
@@ -116,6 +117,7 @@
 #include <opm/parser/eclipse/EclipseState/Tables/TableContainer.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/TableSchema.hpp>
+#include <opm/output/data/Aquifer.hpp>
 #include <opm/output/eclipse/RestartValue.hpp>
 #include <opm/simulators/utils/ParallelRestart.hpp>
 #include <ebos/eclmpiserializer.hh>
@@ -260,6 +262,39 @@ Opm::data::NodeData getNodeData()
         123.457
     };
 }
+
+Opm::data::AquiferData getFetkovichAquifer(const int aquiferID = 1)
+{
+    auto aquifer = Opm::data::AquiferData {
+        aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5, Opm::data::AquiferType::Fetkovich
+    };
+
+    aquifer.aquFet = std::make_shared<Opm::data::FetkovichData>();
+
+    aquifer.aquFet->initVolume = 1.23;
+    aquifer.aquFet->prodIndex = 45.67;
+    aquifer.aquFet->timeConstant = 890.123;
+
+    return aquifer;
+}
+
+Opm::data::AquiferData getCarterTracyAquifer(const int aquiferID = 5)
+{
+    auto aquifer = Opm::data::AquiferData {
+        aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5, Opm::data::AquiferType::CarterTracy
+    };
+
+    aquifer.aquCT = std::make_shared<Opm::data::CarterTracyData>();
+
+    aquifer.aquCT->timeConstant = 987.65;
+    aquifer.aquCT->influxConstant = 43.21;
+    aquifer.aquCT->waterDensity = 1014.5;
+    aquifer.aquCT->waterViscosity = 0.00318;
+    aquifer.aquCT->dimensionless_time = 42.0;
+    aquifer.aquCT->dimensionless_pressure = 2.34;
+
+    return aquifer;
+}
 }
 
 
@@ -313,6 +348,35 @@ BOOST_AUTO_TEST_CASE(Rates)
     DO_CHECKS(data::Rates)
 }
 
+BOOST_AUTO_TEST_CASE(dataFetkovichData)
+{
+    const auto val1 = getFetkovichAquifer();
+    const auto val2 = PackUnpack(val1);
+
+    DO_CHECKS(data::FetkovichData)
+}
+
+BOOST_AUTO_TEST_CASE(dataCarterTracyData)
+{
+    const auto val1 = getCarterTracyAquifer();
+    const auto val2 = PackUnpack(val1);
+
+    DO_CHECKS(data::CarterTracyData)
+}
+
+BOOST_AUTO_TEST_CASE(dataAquifers)
+{
+    const auto val1 = Opm::data::Aquifers {
+        { 1, getFetkovichAquifer(1) },
+        { 4, getFetkovichAquifer(4) },
+        { 5, getCarterTracyAquifer(5) },
+    };
+
+    const auto val2 = PackUnpack(val1);
+
+    DO_CHECKS(data::Aquifers)
+}
+
 BOOST_AUTO_TEST_CASE(dataGuideRateValue)
 {
     using Item = Opm::data::GuideRateValue::Item;
@@ -436,8 +500,13 @@ BOOST_AUTO_TEST_CASE(RestartValue)
         }
     };
 
+    auto aquifers = Opm::data::Aquifers {
+        { 2, getCarterTracyAquifer(2) },
+        {11, getFetkovichAquifer(11) },
+    };
+
     const auto val1 = Opm::RestartValue {
-        getSolution(), std::move(wells1), std::move(grp_nwrk_1)
+        getSolution(), std::move(wells1), std::move(grp_nwrk_1), std::move(aquifers)
     };
     const auto val2 = PackUnpack(val1);