Merge pull request #4519 from bska/unify-analytic-aquifer-construction

Refactor Construction of Analytic Aquifer Objects

opm/simulators/aquifers/AquiferConstantFlux.hpp

@@ -20,41 +20,50 @@
 #ifndef OPM_AQUIFERCONSTANTFLUX_HPP
 #define OPM_AQUIFERCONSTANTFLUX_HPP
 
-#include <opm/material/common/MathToolbox.hpp>
-#include <opm/material/densead/Evaluation.hpp>
-
 #include <opm/simulators/aquifers/AquiferInterface.hpp>
 
 #include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
 #include <opm/input/eclipse/EclipseState/Aquifer/AquiferFlux.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
+
+#include <opm/common/ErrorMacros.hpp>
+
+#include <opm/material/common/MathToolbox.hpp>
+#include <opm/material/densead/Evaluation.hpp>
+
+#include <algorithm>
+#include <numeric>
+#include <stdexcept>
 
 namespace Opm {
+
 template<typename TypeTag>
-class AquiferConstantFlux : public AquiferInterface<TypeTag> {
+class AquiferConstantFlux : public AquiferInterface<TypeTag>
+{
 public:
     using RateVector = GetPropType<TypeTag, Properties::RateVector>;
     using Simulator = GetPropType<TypeTag, Properties::Simulator>;
     using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
+
     static constexpr int numEq = BlackoilIndices::numEq;
     using Eval = DenseAd::Evaluation<double, /*size=*/numEq>;
 
-    AquiferConstantFlux(const SingleAquiferFlux& aquifer,
+    AquiferConstantFlux(const std::vector<Aquancon::AquancCell>& connections,
-                        const std::vector<Aquancon::AquancCell>& connections,
+                        const Simulator&                         ebos_simulator,
-                        const Simulator& ebos_simulator)
+                        const SingleAquiferFlux&                 aquifer)
         : AquiferInterface<TypeTag>(aquifer.id, ebos_simulator)
-         , connections_(connections)
+        , connections_             (connections)
-         , aquifer_data_(aquifer)
+        , aquifer_data_            (aquifer)
+        , connection_flux_         (connections_.size(), Eval{0})
     {
-        // init_cumulative_flux is the flux volume from previoius running
         this->initializeConnections();
-        connection_flux_.resize(this->connections_.size(), {0});
     }
 
     static AquiferConstantFlux serializationTestObject(const Simulator& ebos_simulator)
     {
-        AquiferConstantFlux<TypeTag> result({}, {}, ebos_simulator);
+        AquiferConstantFlux<TypeTag> result({}, ebos_simulator, {});
         result.cumulative_flux_ = 1.0;
 
         return result;
@@ -62,7 +71,8 @@ public:
 
     virtual ~AquiferConstantFlux() = default;
 
-    void updateAquifer(const SingleAquiferFlux& aquifer) {
+    void updateAquifer(const SingleAquiferFlux& aquifer)
+    {
         aquifer_data_ = aquifer;
     }
 
@@ -72,31 +82,32 @@ public:
     void initialSolutionApplied() override {
     }
 
-    void beginTimeStep() override {
-    }
 
-    void endTimeStep() override {
+    void beginTimeStep() override
-        this->flux_rate_ = 0.;
+    {}
-        for (const auto& q : this->connection_flux_) {
-            this->flux_rate_ += Opm::getValue(q);
-        }
 
-        this->cumulative_flux_ += this->flux_rate_ * this->ebos_simulator_.timeStepSize();
+    void endTimeStep() override
+    {
+        this->flux_rate_ = this->totalFluxRate();
+        this->cumulative_flux_ +=
+            this->flux_rate_ * this->ebos_simulator_.timeStepSize();
     }
 
     data::AquiferData aquiferData() const override
     {
         data::AquiferData data;
+
         data.aquiferID = this->aquifer_data_.id;
-        // pressure for constant flux aquifer is 0
+
-        data.pressure = 0.;
+        // Pressure for constant flux aquifer is 0
-        data.fluxRate = 0.;
+        data.pressure = 0.0;
-        for (const auto& q : this->connection_flux_) {
+        data.fluxRate = this->totalFluxRate();
-            data.fluxRate += q.value();
+
-        }
         data.volume = this->cumulative_flux_;
+
         // not totally sure whether initPressure matters
-        data.initPressure = 0.;
+        data.initPressure = 0.0;
+
         return data;
     }
 
@@ -104,19 +115,17 @@ public:
                      const unsigned cellIdx,
                      const unsigned timeIdx) override
     {
-        const auto& model = this->ebos_simulator_.model();
-
         const int idx = this->cellToConnectionIdx_[cellIdx];
-        if (idx < 0)
+        if (idx < 0) {
             return;
-
-        const auto* intQuantsPtr = model.cachedIntensiveQuantities(cellIdx, timeIdx);
-        if (intQuantsPtr == nullptr) {
-            throw std::logic_error("Invalid intensive quantities cache detected in AquiferAnalytical::addToSource()");
         }
 
-        const double fw = this->aquifer_data_.flux;
+        const auto& model = this->ebos_simulator_.model();
+
+        const auto fw = this->aquifer_data_.flux;
+
         this->connection_flux_[idx] = fw * this->connections_[idx].effective_facearea;
+
         rates[BlackoilIndices::conti0EqIdx + compIdx_()]
                 += this->connection_flux_[idx] / model.dofTotalVolume(cellIdx);
     }
@@ -134,17 +143,19 @@ public:
 
 private:
     const std::vector<Aquancon::AquancCell>& connections_;
+
     SingleAquiferFlux aquifer_data_;
-    std::vector<int> cellToConnectionIdx_;
+    std::vector<Eval> connection_flux_{};
-    std::vector<Eval> connection_flux_;
+    std::vector<int> cellToConnectionIdx_{};
-    double flux_rate_ {};
+    double flux_rate_{};
-    double cumulative_flux_ = 0.;
+    double cumulative_flux_{};
 
     void initializeConnections() {
         this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
         for (std::size_t idx = 0; idx < this->connections_.size(); ++idx) {
             const auto global_index = this->connections_[idx].global_index;
-            const int cell_index = this->ebos_simulator_.vanguard().compressedIndexForInterior(global_index);
+            const int cell_index = this->ebos_simulator_.vanguard()
+                .compressedIndexForInterior(global_index);
 
             if (cell_index < 0) {
                 continue;
@@ -152,8 +163,10 @@ private:
 
             this->cellToConnectionIdx_[cell_index] = idx;
         }
-        // TODO: at the moment, we are using the effective_facearea from the parser. Should we update the facearea here if
+
-        //  the grid changed during the preprocessing?
+        // TODO: At the moment, we are using the effective_facearea from the
+        // parser.  Should we update the facearea here if the grid changed
+        // during the preprocessing?
     }
 
     // TODO: this is a function from AquiferAnalytical
@@ -164,7 +177,18 @@ private:
 
         return FluidSystem::waterCompIdx;
     }
+
+    double totalFluxRate() const
+    {
+        return std::accumulate(this->connection_flux_.begin(),
+                               this->connection_flux_.end(), 0.0,
+                               [](const double rate, const auto& q)
+                               {
+                                   return rate + getValue(q);
+                               });
+    }
 };
-}
+
+} // namespace Opm
 
 #endif //OPM_AQUIFERCONSTANTFLUX_HPP

opm/simulators/aquifers/BlackoilAquiferModel.hpp

@@ -48,6 +48,7 @@
 
 #include <vector>
 #include <type_traits>
+#include <string_view>
 
 namespace Opm
 {
@@ -125,8 +126,20 @@ protected:
     // TODO: possibly better to use unorder_map here for aquifers
     std::vector<std::unique_ptr<AquiferInterface<TypeTag>>> aquifers;
 
-    // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy
+    // This initialization function is used to connect the parser objects
+    // with the ones needed by AquiferCarterTracy
     void init();
+
+private:
+    void createDynamicAquifers(const int episode_index);
+
+    void initializeStaticAquifers();
+
+    template <typename AquiferType, typename AquiferData>
+    std::unique_ptr<AquiferType>
+    createAnalyticAquiferPointer(const AquiferData& aqData,
+                                 const int          aquiferID,
+                                 std::string_view   aqType) const;
 };
 
 

opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp

@@ -25,6 +25,7 @@
 #include <algorithm>
 #include <memory>
 #include <stdexcept>
+#include <string_view>
 
 namespace Opm
 {
@@ -60,33 +61,12 @@ template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::beginEpisode()
 {
-    // probably function name beginReportStep() is more appropriate.
+    // Probably function name beginReportStep() is more appropriate.
-    // basically, we want to update the aquifer related information from SCHEDULE setup in this section
+    //
-    // it is the beginning of a report step
+    // Basically, we want to update the aquifer related information from
-    const auto& connections = this->simulator_.vanguard().eclState().aquifer().connections();
+    // SCHEDULE setup in this section it is the beginning of a report step
-    const int report_step  = this->simulator_.episodeIndex();
+
-    const auto& aqufluxs = this->simulator_.vanguard().schedule()[report_step].aqufluxs;
+    this->createDynamicAquifers(this->simulator_.episodeIndex());
-    for (const auto& elem : aqufluxs) {
-        const int id = elem.first;
-        auto find = std::find_if(begin(this->aquifers), end(this->aquifers), [id](auto& v){ return v->aquiferID() == id; });
-        if (find == this->aquifers.end()) {
-            // the aquifer id does not exist in BlackoilAquiferModel yet
-            const auto& aquinfo = elem.second;
-            auto aqf = std::make_unique<AquiferConstantFlux<TypeTag>>(aquinfo, connections.getConnections(aquinfo.id), this->simulator_);
-            this->aquifers.push_back(std::move(aqf));
-        } else {
-            const auto& aquinfo = elem.second;
-            auto aqu = dynamic_cast<AquiferConstantFlux<TypeTag>*> (find->get());
-            if (!aqu) {
-                // if the aquifers can return types easily, we might be able to give a better message with type information
-                const auto msg = fmt::format("Aquifer {} is updated with constant flux aquifer keyword AQUFLUX at report step {},"
-                                                       " while it might be specified to be a different type of aquifer before this. We do not support the conversion between"
-                                                       " different types of aquifer.\n", id, report_step);
-                OPM_THROW(std::runtime_error, msg);
-            }
-            aqu->updateAquifer(aquinfo);
-        }
-    }
 }
 
 template <typename TypeTag>
@@ -167,62 +147,12 @@ BlackoilAquiferModel<TypeTag>::deserialize(Restarter& /* res */)
 
 // Initialize the aquifers in the deck
 template <typename TypeTag>
-void
+void BlackoilAquiferModel<TypeTag>::init()
-BlackoilAquiferModel<TypeTag>::init()
 {
-    const auto& aquifer = this->simulator_.vanguard().eclState().aquifer();
+    if (this->simulator_.vanguard().eclState().aquifer().active()) {
-
+        this->initializeStaticAquifers();
-    if (!aquifer.active()) {
-        return;
     }
 
-    const auto& connections = aquifer.connections();
-    for (const auto& aq : aquifer.ct()) {
-        if (!connections.hasAquiferConnections(aq.aquiferID)) {
-            auto msg = fmt::format("No valid connections for Carter-Tracy aquifer {}, aquifer {} will be ignored.",
-                                   aq.aquiferID, aq.aquiferID);
-            OpmLog::warning(msg);
-            continue;
-        }
-        auto aqf = std::make_unique<AquiferCarterTracy<TypeTag>>(connections.getConnections(aq.aquiferID),
-                                                                 this->simulator_, aq);
-        aquifers.push_back(std::move(aqf));
-    }
-
-    for (const auto& aq : aquifer.fetp()) {
-        if (!connections.hasAquiferConnections(aq.aquiferID)) {
-            auto msg = fmt::format("No valid connections for Fetkovich aquifer {}, aquifer {} will be ignored.",
-                                   aq.aquiferID, aq.aquiferID);
-            OpmLog::warning(msg);
-            continue;
-        }
-        auto aqf = std::make_unique<AquiferFetkovich<TypeTag>>(connections.getConnections(aq.aquiferID),
-                                                               this->simulator_, aq);
-        aquifers.push_back(std::move(aqf));
-    }
-
-    for (const auto& [id, aq] : aquifer.aquflux()) {
-        // make sure not dummy constant flux aquifers
-        if ( !aq.active ) continue;
-
-        if (!connections.hasAquiferConnections(id)) {
-            auto msg = fmt::format("No valid connections for constant flux aquifer {}, aquifer {} will be ignored.",
-                                   id, id);
-            OpmLog::warning(msg);
-            continue;
-        }
-
-        auto aqf = std::make_unique<AquiferConstantFlux<TypeTag>>(aq, connections.getConnections(id), this->simulator_);
-        this->aquifers.push_back(std::move(aqf));
-    }
-
-    if (aquifer.hasNumericalAquifer()) {
-        const auto& num_aquifers = aquifer.numericalAquifers().aquifers();
-        for ([[maybe_unused]]const auto& [id, aqu] : num_aquifers) {
-            auto aqf = std::make_unique<AquiferNumerical<TypeTag>>(aqu, this->simulator_);
-            aquifers.push_back(std::move(aqf));
-        }
-    }
 }
 
 template<typename TypeTag>
@@ -259,4 +189,119 @@ serializeOp(Serializer& serializer)
     }
 }
 
+template <typename TypeTag>
+void BlackoilAquiferModel<TypeTag>::initializeStaticAquifers()
+{
+    const auto& aquifer =
+        this->simulator_.vanguard().eclState().aquifer();
+
+    for (const auto& aquCT : aquifer.ct()) {
+        auto aquCTPtr = this->template createAnalyticAquiferPointer
+            <AquiferCarterTracy<TypeTag>>(aquCT, aquCT.aquiferID, "Carter-Tracy");
+
+        if (aquCTPtr != nullptr) {
+            this->aquifers.push_back(std::move(aquCTPtr));
+        }
+    }
+
+    for (const auto& aquFetp : aquifer.fetp()) {
+        auto aquFetpPtr = this->template createAnalyticAquiferPointer
+            <AquiferFetkovich<TypeTag>>(aquFetp, aquFetp.aquiferID, "Fetkovich");
+
+        if (aquFetpPtr != nullptr) {
+            this->aquifers.push_back(std::move(aquFetpPtr));
+        }
+    }
+
+    for (const auto& [id, aquFlux] : aquifer.aquflux()) {
+        // Make sure not dummy constant flux aquifers
+        if (! aquFlux.active) { continue; }
+
+        auto aquFluxPtr = this->template createAnalyticAquiferPointer
+            <AquiferConstantFlux<TypeTag>>(aquFlux, id, "Constant Flux");
+
+        if (aquFluxPtr != nullptr) {
+            this->aquifers.push_back(std::move(aquFluxPtr));
+        }
+    }
+
+    if (aquifer.hasNumericalAquifer()) {
+        for (const auto& aquNum : aquifer.numericalAquifers().aquifers()) {
+            auto aquNumPtr = std::make_unique<AquiferNumerical<TypeTag>>
+                (aquNum.second, this->simulator_);
+
+            this->aquifers.push_back(std::move(aquNumPtr));
+        }
+    }
+}
+
+template <typename TypeTag>
+template <typename AquiferType, typename AquiferData>
+std::unique_ptr<AquiferType>
+BlackoilAquiferModel<TypeTag>::
+createAnalyticAquiferPointer(const AquiferData& aqData,
+                             const int          aquiferID,
+                             std::string_view   aqType) const
+{
+    const auto& connections =
+        this->simulator_.vanguard().eclState().aquifer().connections();
+
+    if (! connections.hasAquiferConnections(aquiferID)) {
+        const auto msg = fmt::format("No valid connections for {} aquifer {}.  "
+                                     "Aquifer {} will be ignored.",
+                                     aqType, aquiferID, aquiferID);
+        OpmLog::warning(msg);
+
+        return {};
+    }
+
+    return std::make_unique<AquiferType>
+        (connections.getConnections(aquiferID), this->simulator_, aqData);
+}
+
+template <typename TypeTag>
+void BlackoilAquiferModel<TypeTag>::createDynamicAquifers(const int episode_index)
+{
+    const auto& sched = this->simulator_.vanguard().schedule()[episode_index];
+
+    for (const auto& [id, aquFlux] : sched.aqufluxs) {
+        auto aquPos =
+            std::find_if(std::begin(this->aquifers),
+                         std::end(this->aquifers),
+                [id](const auto& aquPtr)
+            {
+                return aquPtr->aquiferID() == id;
+            });
+
+        if (aquPos == std::end(this->aquifers)) {
+            // An aquifer with this 'id' does not yet exist in
+            // the collection managed by this object.  Create it.
+            auto aquFluxPtr = this->template createAnalyticAquiferPointer
+                <AquiferConstantFlux<TypeTag>>(aquFlux, id, "Constant Flux");
+
+            if (aquFluxPtr != nullptr) {
+                this->aquifers.push_back(std::move(aquFluxPtr));
+            }
+        }
+        else {
+            auto aquFluxPtr = dynamic_cast<AquiferConstantFlux<TypeTag>*>(aquPos->get());
+            if (aquFluxPtr == nullptr) {
+                // If the aquifers can return types easily, we might be able
+                // to give a better message with type information.
+                const auto msg =
+                    fmt::format("Aquifer {} is updated with constant flux "
+                                "aquifer keyword AQUFLUX at report step {}, "
+                                "while it might be specified to be a "
+                                "different type of aquifer before this. "
+                                "We do not support the conversion between "
+                                "different types of aquifer.\n", id, episode_index);
+
+                OPM_THROW(std::runtime_error, msg);
+            }
+
+            aquFluxPtr->updateAquifer(aquFlux);
+        }
+    }
+}
+
 } // namespace Opm

tests/test_RestartSerialization.cpp

@@ -491,7 +491,7 @@ BOOST_AUTO_TEST_CASE(AquiferConstantFlux)
     Opm::Serializer ser(packer);
     ser.pack(data_out);
     const size_t pos1 = ser.position();
-    decltype(data_out) data_in({}, {}, sim);
+    decltype(data_out) data_in({}, sim, {});
     ser.unpack(data_in);
     const size_t pos2 = ser.position();
     BOOST_CHECK_MESSAGE(pos1 == pos2, "Packed size differ from unpack size for AquiferConstantFlux");