Chase Type Specific Aquifer Data API Change

This commit switches to using the new 'typeData' interface for
representing type-specific aquifer data items.  In particular we use
the new 'typeData.is<>()' and 'typeData.get<>()' member functions to
query and access the data that is specific to each aquifer type
(e.g., Carter-Tracy or numerical).

While here, also expand the reported data items for numerical
aquifers to one initial pressure value for each aquifer cell.  This
is needed for restart purposes.

ebos/collecttoiorank.cc

@@ -37,6 +37,7 @@
 #include <ebos/femcpgridcompat.hh>
 #endif
 
+#include <algorithm>
 #include <cassert>
 #include <stdexcept>
 #include <string>
@@ -627,12 +628,20 @@ public:
 
             this->unpackCommon(data, aq);
 
-            if (data.aquFet != nullptr) {
+            if (auto const* aquFet = data.typeData.get<data::AquiferType::Fetkovich>();
-                this->unpackAquFet(std::move(data.aquFet), aq);
+                aquFet != nullptr)
+            {
+                this->unpackAquFet(*aquFet, aq);
             }
-
+            else if (auto const* aquCT = data.typeData.get<data::AquiferType::CarterTracy>();
-            if (data.aquCT != nullptr) {
+                     aquCT != nullptr)
-                this->unpackCarterTracy(std::move(data.aquCT), aq);
+            {
+                this->unpackCarterTracy(*aquCT, aq);
+            }
+            else if (auto const* aquNum = data.typeData.get<data::AquiferType::Numerical>();
+                     aquNum != nullptr)
+            {
+                this->unpackNumericAquifer(*aquNum, aq);
             }
         }
     }
@@ -648,14 +657,15 @@ private:
         aq.volume += data.volume;
     }
 
-    void unpackAquFet(std::shared_ptr<data::FetkovichData> aquFet, data::AquiferData& aq)
+    void unpackAquFet(const data::FetkovichData& aquFet, data::AquiferData& aq)
     {
-        if (aq.aquFet == nullptr) {
+        if (! aq.typeData.is<data::AquiferType::Fetkovich>()) {
-            aq.aquFet = std::move(aquFet);
+            auto* tData = aq.typeData.create<data::AquiferType::Fetkovich>();
+            *tData = aquFet;
         }
         else {
-            const auto& src = *aquFet;
+            const auto& src = aquFet;
-            auto&       dst = *aq.aquFet;
+            auto&       dst = *aq.typeData.getMutable<data::AquiferType::Fetkovich>();
 
             dst.initVolume   = std::max(dst.initVolume  , src.initVolume);
             dst.prodIndex    = std::max(dst.prodIndex   , src.prodIndex);
@@ -663,19 +673,46 @@ private:
         }
     }
 
-    void unpackCarterTracy(std::shared_ptr<data::CarterTracyData> aquCT, data::AquiferData& aq)
+    void unpackCarterTracy(const data::CarterTracyData& aquCT, data::AquiferData& aq)
     {
-        if (aq.aquCT == nullptr) {
+        if (! aq.typeData.is<data::AquiferType::CarterTracy>()) {
-            aq.aquCT = std::move(aquCT);
+            auto* tData = aq.typeData.create<data::AquiferType::CarterTracy>();
+            *tData = aquCT;
         }
         else {
-            const auto& src = *aquCT;
+            const auto& src = aquCT;
-            auto&       dst = *aq.aquCT;
+            auto&       dst = *aq.typeData.getMutable<data::AquiferType::CarterTracy>();
+
+            dst.timeConstant   = std::max(dst.timeConstant  , src.timeConstant);
+            dst.influxConstant = std::max(dst.influxConstant, src.influxConstant);
+            dst.waterDensity   = std::max(dst.waterDensity  , src.waterDensity);
+            dst.waterViscosity = std::max(dst.waterViscosity, src.waterViscosity);
 
             dst.dimensionless_time     = std::max(dst.dimensionless_time    , src.dimensionless_time);
             dst.dimensionless_pressure = std::max(dst.dimensionless_pressure, src.dimensionless_pressure);
         }
     }
+
+    void unpackNumericAquifer(const data::NumericAquiferData& aquNum, data::AquiferData& aq)
+    {
+        if (! aq.typeData.is<data::AquiferType::Numerical>()) {
+            auto* tData = aq.typeData.create<data::AquiferType::Numerical>();
+            *tData = aquNum;
+        }
+        else {
+            const auto& src = aquNum;
+            auto&       dst = *aq.typeData.getMutable<data::AquiferType::Numerical>();
+
+            std::transform(src.initPressure.begin(),
+                           src.initPressure.end(),
+                           dst.initPressure.begin(),
+                           dst.initPressure.begin(),
+                           [](const double p0_1, const double p0_2)
+                           {
+                               return std::max(p0_1, p0_2);
+                           });
+        }
+    }
 };
 
 

opm/simulators/aquifers/AquiferCarterTracy.hpp

@@ -81,15 +81,14 @@ public:
         }
         data.volume = this->W_flux_.value();
         data.initPressure = this->pa0_;
-        data.type = data::AquiferType::CarterTracy;
 
-        data.aquCT = std::make_shared<data::CarterTracyData>();
+        auto* aquCT = data.typeData.template create<data::AquiferType::CarterTracy>();
-        data.aquCT->timeConstant = this->aquct_data_.timeConstant();
+        aquCT->timeConstant = this->aquct_data_.timeConstant();
-        data.aquCT->influxConstant = this->aquct_data_.influxConstant();
+        aquCT->influxConstant = this->aquct_data_.influxConstant();
-        data.aquCT->waterDensity = this->aquct_data_.waterDensity();
+        aquCT->waterDensity = this->aquct_data_.waterDensity();
-        data.aquCT->waterViscosity = this->aquct_data_.waterViscosity();
+        aquCT->waterViscosity = this->aquct_data_.waterViscosity();
-        data.aquCT->dimensionless_time = this->dimensionless_time_;
+        aquCT->dimensionless_time = this->dimensionless_time_;
-        data.aquCT->dimensionless_pressure = this->dimensionless_pressure_;
+        aquCT->dimensionless_pressure = this->dimensionless_pressure_;
 
         return data;
     }

opm/simulators/aquifers/AquiferFetkovich.hpp

@@ -27,6 +27,7 @@ along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 
 #include <exception>
 #include <stdexcept>
+#include <utility>
 
 namespace Opm
 {
@@ -82,12 +83,11 @@ public:
                                         });
         data.volume = this->W_flux_.value();
         data.initPressure = this->pa0_;
-        data.type = data::AquiferType::Fetkovich;
 
-        data.aquFet = std::make_shared<data::FetkovichData>();
+        auto* aquFet = data.typeData.template create<data::AquiferType::Fetkovich>();
-        data.aquFet->initVolume = this->aqufetp_data_.initial_watvolume;
+        aquFet->initVolume = this->aqufetp_data_.initial_watvolume;
-        data.aquFet->prodIndex = this->aqufetp_data_.prod_index;
+        aquFet->prodIndex = this->aqufetp_data_.prod_index;
-        data.aquFet->timeConstant = this->aqufetp_data_.timeConstant();
+        aquFet->timeConstant = this->aqufetp_data_.timeConstant();
 
         return data;
     }
@@ -99,7 +99,7 @@ protected:
 
     void assignRestartData(const data::AquiferData& xaq) override
     {
-        if (xaq.type != data::AquiferType::Fetkovich) {
+        if (! xaq.typeData.is<data::AquiferType::Fetkovich>()) {
             throw std::invalid_argument {
                 "Analytic aquifer data for unexpected aquifer "
                 "type passed to Fetkovich aquifer"

opm/simulators/aquifers/AquiferNumerical.hpp

@@ -24,6 +24,8 @@
 #include <opm/output/data/Aquifer.hpp>
 #include <opm/parser/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
 
+#include <utility>
+
 namespace Opm
 {
 template <typename TypeTag>
@@ -57,17 +59,17 @@ public:
         , flux_rate_(0.)
         , cumulative_flux_(0.)
         , global_cell_(global_cell)
+        , init_pressure_(aquifer.numCells(), 0.0)
     {
         this->cell_to_aquifer_cell_idx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
 
         for (size_t idx = 0; idx < aquifer.numCells(); ++idx) {
-            const auto& cell = *(aquifer.getCellPrt(idx));
+            const auto* cell = aquifer.getCellPrt(idx);
-            const int global_idx = cell.global_index;
+
-            const auto search = cartesian_to_compressed.find(global_idx);
             // Due to parallelisation, the cell might not exist in the current process
+            auto search = cartesian_to_compressed.find(cell->global_index);
             if (search != cartesian_to_compressed.end()) {
-                const int cell_idx = cartesian_to_compressed.at(global_idx);
+                this->cell_to_aquifer_cell_idx_[search->second] = idx;
-                this->cell_to_aquifer_cell_idx_[cell_idx] = idx;
             }
         }
     }
@@ -87,19 +89,20 @@ public:
     data::AquiferData aquiferData() const
     {
         data::AquiferData data;
-        data.aquiferID = this->id_;
+        data.aquiferID = this->aquiferID();
-        data.initPressure = this->init_pressure_;
         data.pressure = this->pressure_;
         data.fluxRate = this->flux_rate_;
         data.volume = this->cumulative_flux_;
-        data.type = data::AquiferType::Numerical;
+
+        auto* aquNum = data.typeData.template create<data::AquiferType::Numerical>();
+        aquNum->initPressure = this->init_pressure_;
+
         return data;
     }
 
     void initialSolutionApplied()
     {
-        this->init_pressure_ = this->calculateAquiferPressure();
+        this->pressure_ = this->calculateAquiferPressure(this->init_pressure_);
-        this->pressure_ = this->init_pressure_;
         this->flux_rate_ = 0.;
         this->cumulative_flux_ = 0.;
     }
@@ -115,13 +118,19 @@ private:
     double flux_rate_; // aquifer influx rate
     double cumulative_flux_; // cumulative aquifer influx
     const int* global_cell_; // mapping to global index
-    double init_pressure_;
+    std::vector<double> init_pressure_{};
     double pressure_; // aquifer pressure
 
     // TODO: maybe unordered_map can also do the work to save memory?
     std::vector<int> cell_to_aquifer_cell_idx_;
 
     double calculateAquiferPressure() const
+    {
+        auto capture = std::vector<double>(this->init_pressure_.size(), 0.0);
+        return this->calculateAquiferPressure(capture);
+    }
+
+    double calculateAquiferPressure(std::vector<double>& cell_pressure) const
     {
         double sum_pressure_watervolume = 0.;
         double sum_watervolume = 0.;
@@ -158,11 +167,17 @@ private:
             const double water_volume = volume * porosity * water_saturation;
             sum_pressure_watervolume += water_volume * water_pressure_reservoir;
             sum_watervolume += water_volume;
+
+            cell_pressure[idx] = water_pressure_reservoir;
         }
 
         const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
         comm.sum(&sum_pressure_watervolume, 1);
         comm.sum(&sum_watervolume, 1);
+
+        // Ensure all processes have same notion of the aquifer cells' pressure values.
+        comm.sum(cell_pressure.data(), cell_pressure.size());
+
         return sum_pressure_watervolume / sum_watervolume;
     }
 

opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp

@@ -83,14 +83,12 @@ BlackoilAquiferModel<TypeTag>::initFromRestart(const data::Aquifers& aquiferSoln
 template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::beginEpisode()
-{
+{}
-}
 
 template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::beginIteration()
-{
+{}
-}
 
 template <typename TypeTag>
 void
@@ -131,8 +129,7 @@ BlackoilAquiferModel<TypeTag>::addToSource(RateVector& rates,
 template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::endIteration()
-{
+{}
-}
 
 template <typename TypeTag>
 void
@@ -155,11 +152,11 @@ BlackoilAquiferModel<TypeTag>::endTimeStep()
         }
     }
 }
+
 template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::endEpisode()
-{
+{}
-}
 
 template <typename TypeTag>
 template <class Restarter>
@@ -215,12 +212,14 @@ BlackoilAquiferModel<TypeTag>::init()
         }
     }
 }
+
 template <typename TypeTag>
 bool
 BlackoilAquiferModel<TypeTag>::aquiferCarterTracyActive() const
 {
     return !aquifers_CarterTracy.empty();
 }
+
 template <typename TypeTag>
 bool
 BlackoilAquiferModel<TypeTag>::aquiferFetkovichActive() const
@@ -236,7 +235,8 @@ BlackoilAquiferModel<TypeTag>::aquiferNumericalActive() const
 }
 
 template<typename TypeTag>
-data::Aquifers BlackoilAquiferModel<TypeTag>::aquiferData() const {
+data::Aquifers BlackoilAquiferModel<TypeTag>::aquiferData() const
+{
     data::Aquifers data;
     if (this->aquiferCarterTracyActive()) {
         for (const auto& aqu : this->aquifers_CarterTracy) {

opm/simulators/utils/ParallelRestart.cpp

@@ -226,6 +226,11 @@ HANDLE_AS_POD(data::NodeData)
 HANDLE_AS_POD(data::Rates)
 HANDLE_AS_POD(data::Segment)
 
+std::size_t packSize(const data::NumericAquiferData& data, Dune::MPIHelper::MPICommunicator comm)
+{
+    return packSize(data.initPressure, comm);
+}
+
 std::size_t packSize(const data::AquiferData& data, Dune::MPIHelper::MPICommunicator comm)
 {
     const auto type = 0ull;
@@ -238,11 +243,20 @@ std::size_t packSize(const data::AquiferData& data, Dune::MPIHelper::MPICommunic
         + packSize(data.datumDepth, comm)
         + packSize(type, comm);
 
-    if (data.aquFet != nullptr) {
+    if (auto const* aquFet = data.typeData.get<data::AquiferType::Fetkovich>();
-        return base + packSize(*data.aquFet, comm);
+        aquFet != nullptr)
+    {
+        return base + packSize(*aquFet, comm);
     }
-    else if (data.aquCT != nullptr) {
+    else if (auto const* aquCT = data.typeData.get<data::AquiferType::CarterTracy>();
-        return base + packSize(*data.aquCT, comm);
+             aquCT != nullptr)
+    {
+        return base + packSize(*aquCT, comm);
+    }
+    else if (auto const* aquNum = data.typeData.get<data::AquiferType::Numerical>();
+             aquNum != nullptr)
+    {
+        return base + packSize(*aquNum, comm);
     }
 
     return base;
@@ -513,12 +527,19 @@ void pack(const std::unordered_map<T1,T2,H,P,A>& data, std::vector<char>& buffer
     }
 }
 
+void pack(const data::NumericAquiferData& data, std::vector<char>& buffer, int& position,
+          Dune::MPIHelper::MPICommunicator comm)
+{
+    pack(data.initPressure, buffer, position, comm);
+}
+
 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.typeData.is<data::AquiferType::Fetkovich>()   * (1ull << 0))
-        + (data.aquCT  != nullptr)*(1ull << 1);
+        + (data.typeData.is<data::AquiferType::CarterTracy>() * (1ull << 1))
+        + (data.typeData.is<data::AquiferType::Numerical>()   * (1ull << 2));
 
     pack(data.aquiferID, buffer, position, comm);
     pack(data.pressure, buffer, position, comm);
@@ -528,11 +549,20 @@ void pack(const data::AquiferData& data, std::vector<char>& buffer, int& positio
     pack(data.datumDepth, buffer, position, comm);
     pack(type, buffer, position, comm);
 
-    if (data.aquFet != nullptr) {
+    if (auto const* aquFet = data.typeData.get<data::AquiferType::Fetkovich>();
-        pack(*data.aquFet, buffer, position, comm);
+        aquFet != nullptr)
+    {
+        pack(*aquFet, buffer, position, comm);
     }
-    else if (data.aquCT != nullptr) {
+    else if (auto const* aquCT = data.typeData.get<data::AquiferType::CarterTracy>();
-        pack(*data.aquCT, buffer, position, comm);
+             aquCT != nullptr)
+    {
+        pack(*aquCT, buffer, position, comm);
+    }
+    else if (auto const* aquNum = data.typeData.get<data::AquiferType::Numerical>();
+             aquNum != nullptr)
+    {
+        pack(*aquNum, buffer, position, comm);
     }
 }
 
@@ -865,6 +895,12 @@ void unpack(data::Well& data, std::vector<char>& buffer, int& position,
     unpack(data.guide_rates, buffer, position, comm);
 }
 
+void unpack(data::NumericAquiferData& data, std::vector<char>& buffer, int& position,
+            Dune::MPIHelper::MPICommunicator comm)
+{
+    unpack(data.initPressure, buffer, position, comm);
+}
+
 void unpack(data::AquiferData& data, std::vector<char>& buffer, int& position,
             Dune::MPIHelper::MPICommunicator comm)
 {
@@ -878,15 +914,17 @@ void unpack(data::AquiferData& data, std::vector<char>& buffer, int& position,
     unpack(data.datumDepth, buffer, position, comm);
     unpack(type, buffer, position, comm);
 
-    if (type == 1ull) {
+    if (type == (1ull << 0)) {
-        data.type = data::AquiferType::Fetkovich;
+        auto* aquFet = data.typeData.create<data::AquiferType::Fetkovich>();
-        data.aquFet = std::make_shared<data::FetkovichData>();
+        unpack(*aquFet, buffer, position, comm);
-        unpack(*data.aquFet, buffer, position, comm);
     }
-    else if (type == 2ull) {
+    else if (type == (1ull << 1)) {
-        data.type = data::AquiferType::CarterTracy;
+        auto* aquCT = data.typeData.create<data::AquiferType::CarterTracy>();
-        data.aquCT = std::make_shared<data::CarterTracyData>();
+        unpack(*aquCT, buffer, position, comm);
-        unpack(*data.aquCT, buffer, position, comm);
+    }
+    else if (type == (1ull << 2)) {
+        auto* aquNum = data.typeData.create<data::AquiferType::Numerical>();
+        unpack(*aquNum, buffer, position, comm);
     }
 }
 

opm/simulators/utils/ParallelRestart.hpp

@@ -61,6 +61,7 @@ struct GroupData;
 struct GroupGuideRates;
 class GuideRateValue;
 struct NodeData;
+struct NumericAquiferData;
 class Rates;
 struct Segment;
 class Solution;
@@ -349,6 +350,7 @@ ADD_PACK_PROTOTYPES(data::GroupGuideRates)
 ADD_PACK_PROTOTYPES(data::GroupData)
 ADD_PACK_PROTOTYPES(data::NodeData)
 ADD_PACK_PROTOTYPES(data::GroupAndNetworkValues)
+ADD_PACK_PROTOTYPES(data::NumericAquiferData)
 ADD_PACK_PROTOTYPES(data::Well)
 ADD_PACK_PROTOTYPES(data::Wells)
 ADD_PACK_PROTOTYPES(RestartKey)

tests/test_ParallelRestart.cpp

@@ -265,14 +265,14 @@ Opm::data::NodeData getNodeData()
 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
+        aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5
     };
 
-    aquifer.aquFet = std::make_shared<Opm::data::FetkovichData>();
+    auto* aquFet = aquifer.typeData.create<Opm::data::AquiferType::Fetkovich>();
 
-    aquifer.aquFet->initVolume = 1.23;
+    aquFet->initVolume = 1.23;
-    aquifer.aquFet->prodIndex = 45.67;
+    aquFet->prodIndex = 45.67;
-    aquifer.aquFet->timeConstant = 890.123;
+    aquFet->timeConstant = 890.123;
 
     return aquifer;
 }
@@ -280,17 +280,33 @@ Opm::data::AquiferData getFetkovichAquifer(const int aquiferID = 1)
 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
+        aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5
     };
 
-    aquifer.aquCT = std::make_shared<Opm::data::CarterTracyData>();
+    auto* aquCT = aquifer.typeData.create<Opm::data::AquiferType::CarterTracy>();
 
-    aquifer.aquCT->timeConstant = 987.65;
+    aquCT->timeConstant = 987.65;
-    aquifer.aquCT->influxConstant = 43.21;
+    aquCT->influxConstant = 43.21;
-    aquifer.aquCT->waterDensity = 1014.5;
+    aquCT->waterDensity = 1014.5;
-    aquifer.aquCT->waterViscosity = 0.00318;
+    aquCT->waterViscosity = 0.00318;
-    aquifer.aquCT->dimensionless_time = 42.0;
+    aquCT->dimensionless_time = 42.0;
-    aquifer.aquCT->dimensionless_pressure = 2.34;
+    aquCT->dimensionless_pressure = 2.34;
+
+    return aquifer;
+}
+
+Opm::data::AquiferData getNumericalAquifer(const int aquiferID = 2)
+{
+    auto aquifer = Opm::data::AquiferData {
+        aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5
+    };
+
+    auto* aquNum = aquifer.typeData.create<Opm::data::AquiferType::Numerical>();
+
+    aquNum->initPressure.push_back(1.234);
+    aquNum->initPressure.push_back(2.345);
+    aquNum->initPressure.push_back(3.4);
+    aquNum->initPressure.push_back(9.876);
 
     return aquifer;
 }
@@ -363,12 +379,21 @@ BOOST_AUTO_TEST_CASE(dataCarterTracyData)
     DO_CHECKS(data::CarterTracyData)
 }
 
+BOOST_AUTO_TEST_CASE(dataNumericAquiferData)
+{
+    const auto val1 = getNumericalAquifer();
+    const auto val2 = PackUnpack(val1);
+
+    DO_CHECKS(data::NumericAquiferData)
+}
+
 BOOST_AUTO_TEST_CASE(dataAquifers)
 {
     const auto val1 = Opm::data::Aquifers {
         { 1, getFetkovichAquifer(1) },
         { 4, getFetkovichAquifer(4) },
         { 5, getCarterTracyAquifer(5) },
+        { 9, getNumericalAquifer(9) },
     };
 
     const auto val2 = PackUnpack(val1);