added: (re)introduce AquiferInterface

now a base class for both analytical and numerical aquifers

opm/simulators/aquifers/AquiferAnalytical.hpp

@@ -22,6 +22,8 @@
 #ifndef OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
 #define OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
 
+#include <opm/simulators/aquifers/AquiferInterface.hpp>
+
 #include <opm/common/utility/numeric/linearInterpolation.hpp>
 #include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
 
@@ -45,7 +47,7 @@
 namespace Opm
 {
 template <typename TypeTag>
-class AquiferAnalytical
+class AquiferAnalytical : public AquiferInterface<TypeTag>
 {
 public:
     using Simulator = GetPropType<TypeTag, Properties::Simulator>;
@@ -81,9 +83,8 @@ public:
     AquiferAnalytical(int aqID,
                      const std::vector<Aquancon::AquancCell>& connections,
                      const Simulator& ebosSimulator)
-        : aquiferID_(aqID)
+        : AquiferInterface<TypeTag>(aqID, ebosSimulator)
         , connections_(connections)
-        , ebos_simulator_(ebosSimulator)
     {
     }
 
@@ -92,7 +93,7 @@ public:
     {
     }
 
-    void initFromRestart(const data::Aquifers& aquiferSoln)
+    void initFromRestart(const data::Aquifers& aquiferSoln) override
     {
         auto xaqPos = aquiferSoln.find(this->aquiferID());
         if (xaqPos == aquiferSoln.end())
@@ -105,16 +106,16 @@ public:
         this->solution_set_from_restart_ = true;
     }
 
-    void initialSolutionApplied()
+    void initialSolutionApplied() override
     {
         initQuantities();
     }
 
-    void beginTimeStep()
+    void beginTimeStep() override
     {
-        ElementContext elemCtx(ebos_simulator_);
+        ElementContext elemCtx(this->ebos_simulator_);
-        auto elemIt = ebos_simulator_.gridView().template begin<0>();
+        auto elemIt = this->ebos_simulator_.gridView().template begin<0>();
-        const auto& elemEndIt = ebos_simulator_.gridView().template end<0>();
+        const auto& elemEndIt = this->ebos_simulator_.gridView().template end<0>();
         OPM_BEGIN_PARALLEL_TRY_CATCH();
 
         for (; elemIt != elemEndIt; ++elemIt) {
@@ -129,27 +130,18 @@ public:
 
             elemCtx.updateIntensiveQuantities(0);
             const auto& iq = elemCtx.intensiveQuantities(0, 0);
-            pressure_previous_[idx] = getValue(iq.fluidState().pressure(phaseIdx_()));
+            pressure_previous_[idx] = getValue(iq.fluidState().pressure(this->phaseIdx_()));
         }
 
-        OPM_END_PARALLEL_TRY_CATCH("AquiferAnalytical::beginTimeStep() failed: ", ebos_simulator_.vanguard().grid().comm());
+        OPM_END_PARALLEL_TRY_CATCH("AquiferAnalytical::beginTimeStep() failed: ",
-    }
+                                   this->ebos_simulator_.vanguard().grid().comm());
-
-    template <class Context>
-    void addToSource(RateVector& rates,
-                     const Context& context,
-                     const unsigned spaceIdx,
-                     const unsigned timeIdx)
-    {
-        const unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
-        addToSource(rates, cellIdx, timeIdx);
     }
 
     void addToSource(RateVector& rates,
                      const unsigned cellIdx,
-                     const unsigned timeIdx)
+                     const unsigned timeIdx) override
     {
-        const auto& model = ebos_simulator_.model();
+        const auto& model = this->ebos_simulator_.model();
 
         const int idx = this->cellToConnectionIdx_[cellIdx];
         if (idx < 0)
@@ -162,46 +154,39 @@ public:
 
         // This is the pressure at td + dt
         this->updateCellPressure(this->pressure_current_, idx, *intQuantsPtr);
-        this->calculateInflowRate(idx, ebos_simulator_);
+        this->calculateInflowRate(idx, this->ebos_simulator_);
 
         rates[BlackoilIndices::conti0EqIdx + compIdx_()]
             += this->Qai_[idx] / model.dofTotalVolume(cellIdx);
     }
-    std::size_t size() const {
+
+    std::size_t size() const
+    {
         return this->connections_.size();
     }
 
-    int aquiferID() const { return this->aquiferID_; }
-
 protected:
-    inline Scalar gravity_() const
+    virtual void assignRestartData(const data::AquiferData& xaq) = 0;
+    virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0;
+    virtual void calculateAquiferCondition() = 0;
+    virtual void calculateAquiferConstants() = 0;
+    virtual Scalar aquiferDepth() const = 0;
+
+    Scalar gravity_() const
     {
-        return ebos_simulator_.problem().gravity()[2];
+        return this->ebos_simulator_.problem().gravity()[2];
     }
 
-    inline bool co2store_() const
+    int compIdx_() const
     {
-        return ebos_simulator_.vanguard().eclState().runspec().co2Storage();
+        if (this->co2store_())
-    }
-
-    inline int phaseIdx_() const
-    {
-        if(co2store_())
-            return FluidSystem::oilPhaseIdx;
-
-        return FluidSystem::waterPhaseIdx;
-    }
-
-    inline int compIdx_() const
-    {
-        if(co2store_())
             return FluidSystem::oilCompIdx;
 
         return FluidSystem::waterCompIdx;
     }
 
 
-    inline void initQuantities()
+    void initQuantities()
     {
         // We reset the cumulative flux at the start of any simulation, so, W_flux = 0
         if (!this->solution_set_from_restart_) {
@@ -219,45 +204,22 @@ protected:
         Qai_.resize(this->connections_.size(), Scalar{0});
     }
 
-    inline void
+    void updateCellPressure(std::vector<Eval>& pressure_water,
-    updateCellPressure(std::vector<Eval>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
+                            const int idx,
+                            const IntensiveQuantities& intQuants)
     {
         const auto& fs = intQuants.fluidState();
-        pressure_water.at(idx) = fs.pressure(phaseIdx_());
+        pressure_water.at(idx) = fs.pressure(this->phaseIdx_());
     }
 
-    inline void
+    void updateCellPressure(std::vector<Scalar>& pressure_water,
-    updateCellPressure(std::vector<Scalar>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
+                            const int idx,
+                            const IntensiveQuantities& intQuants)
     {
         const auto& fs = intQuants.fluidState();
-        pressure_water.at(idx) = fs.pressure(phaseIdx_()).value();
+        pressure_water.at(idx) = fs.pressure(this->phaseIdx_()).value();
     }
 
-    virtual void endTimeStep() = 0;
-
-    const int aquiferID_{};
-    const std::vector<Aquancon::AquancCell> connections_;
-    const Simulator& ebos_simulator_;
-
-    // Grid variables
-    std::vector<Scalar> faceArea_connected_;
-    std::vector<int> cellToConnectionIdx_;
-
-    // Quantities at each grid id
-    std::vector<Scalar> cell_depth_;
-    std::vector<Scalar> pressure_previous_;
-    std::vector<Eval> pressure_current_;
-    std::vector<Eval> Qai_;
-    std::vector<Scalar> alphai_;
-
-    Scalar Tc_{}; // Time constant
-    Scalar pa0_{}; // initial aquifer pressure
-    Scalar rhow_{};
-
-    Eval W_flux_;
-
-    bool solution_set_from_restart_ {false};
-
     void initializeConnections()
     {
         this->cell_depth_.resize(this->size(), this->aquiferDepth());
@@ -381,18 +343,8 @@ protected:
         this->W_flux_ *= this_area / tot_area;
     }
 
-    virtual void assignRestartData(const data::AquiferData& xaq) = 0;
-
-    virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0;
-
-    virtual void calculateAquiferCondition() = 0;
-
-    virtual void calculateAquiferConstants() = 0;
-
-    virtual Scalar aquiferDepth() const = 0;
-
     // This function is for calculating the aquifer properties from equilibrium state with the reservoir
-    virtual Scalar calculateReservoirEquilibrium()
+    Scalar calculateReservoirEquilibrium()
     {
         // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
         std::vector<Scalar> pw_aquifer;
@@ -415,8 +367,8 @@ protected:
             const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
             const auto& fs = iq0.fluidState();
 
-            water_pressure_reservoir = fs.pressure(phaseIdx_()).value();
+            water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
-            const auto water_density = fs.density(phaseIdx_());
+            const auto water_density = fs.density(this->phaseIdx_());
 
             const auto gdz =
                 this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
@@ -426,7 +378,7 @@ protected:
         }
 
         // We take the average of the calculated equilibrium pressures.
-        const auto& comm = ebos_simulator_.vanguard().grid().comm();
+        const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
 
         Scalar vals[2];
         vals[0] = std::accumulate(this->alphai_.begin(), this->alphai_.end(), Scalar{0});
@@ -437,7 +389,26 @@ protected:
         return vals[1] / vals[0];
     }
 
-    // This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer
+    const std::vector<Aquancon::AquancCell> connections_;
+
+    // Grid variables
+    std::vector<Scalar> faceArea_connected_;
+    std::vector<int> cellToConnectionIdx_;
+
+    // Quantities at each grid id
+    std::vector<Scalar> cell_depth_;
+    std::vector<Scalar> pressure_previous_;
+    std::vector<Eval> pressure_current_;
+    std::vector<Eval> Qai_;
+    std::vector<Scalar> alphai_;
+
+    Scalar Tc_{}; // Time constant
+    Scalar pa0_{}; // initial aquifer pressure
+    Scalar rhow_{};
+
+    Eval W_flux_;
+
+    bool solution_set_from_restart_ {false};
 };
 
 } // namespace Opm

opm/simulators/aquifers/AquiferInterface.hpp

@@ -0,0 +1,94 @@
+/*
+  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2017 Statoil ASA.
+  Copyright 2017 IRIS
+
+  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_AQUIFERINTERFACE_HEADER_INCLUDED
+#define OPM_AQUIFERINTERFACE_HEADER_INCLUDED
+
+#include <opm/output/data/Aquifer.hpp>
+
+namespace Opm
+{
+
+template <typename TypeTag>
+class AquiferInterface
+{
+public:
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
+    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
+
+    // Constructor
+    AquiferInterface(int aqID,
+                     const Simulator& ebosSimulator)
+        : aquiferID_(aqID)
+        , ebos_simulator_(ebosSimulator)
+    {
+    }
+
+    // Destructor
+    virtual ~AquiferInterface() = default;
+
+    virtual void initFromRestart(const data::Aquifers& aquiferSoln) = 0;
+
+    virtual void initialSolutionApplied() = 0;
+
+    virtual void beginTimeStep() = 0;
+    virtual void endTimeStep() = 0;
+
+    virtual data::AquiferData aquiferData() const = 0;
+
+    template <class Context>
+    void addToSource(RateVector& rates,
+                     const Context& context,
+                     const unsigned spaceIdx,
+                     const unsigned timeIdx)
+    {
+        const unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
+        addToSource(rates, cellIdx, timeIdx);
+    }
+
+    virtual void addToSource(RateVector& rates,
+                             const unsigned cellIdx,
+                             const unsigned timeIdx) = 0;
+
+    int aquiferID() const { return this->aquiferID_; }
+
+protected:
+    bool co2store_() const
+    {
+        return ebos_simulator_.vanguard().eclState().runspec().co2Storage();
+    }
+
+    int phaseIdx_() const
+    {
+        if (co2store_())
+            return FluidSystem::oilPhaseIdx;
+
+        return FluidSystem::waterPhaseIdx;
+    }
+
+    const int aquiferID_{};
+    const Simulator& ebos_simulator_;
+};
+
+} // namespace Opm
+
+#endif

opm/simulators/aquifers/AquiferNumerical.hpp

@@ -25,6 +25,7 @@
 
 #include <opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
 
+#include <opm/simulators/aquifers/AquiferInterface.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 
 #include <algorithm>
@@ -37,16 +38,17 @@
 namespace Opm
 {
 template <typename TypeTag>
-class AquiferNumerical
+class AquiferNumerical : public AquiferInterface<TypeTag>
 {
 public:
-    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
-    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
-    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
-
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
     using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
+    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
 
     enum { dimWorld = GridView::dimensionworld };
     enum { numPhases = FluidSystem::numPhases };
@@ -55,15 +57,14 @@ public:
     using Eval =  DenseAd::Evaluation<double, numEq>;
     using Toolbox = MathToolbox<Eval>;
 
-    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using typename AquiferInterface<TypeTag>::RateVector;
-    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
+
     // Constructor
     AquiferNumerical(const SingleNumericalAquifer& aquifer,
                      const std::unordered_map<int, int>& cartesian_to_compressed,
                      const Simulator& ebos_simulator,
                      const int* global_cell)
-        : id_             (aquifer.id())
+        : AquiferInterface<TypeTag>(aquifer.id(), ebos_simulator)
-        , ebos_simulator_ (ebos_simulator)
         , flux_rate_      (0.0)
         , cumulative_flux_(0.0)
         , global_cell_    (global_cell)
@@ -88,7 +89,7 @@ public:
         }
     }
 
-    void initFromRestart(const data::Aquifers& aquiferSoln)
+    void initFromRestart(const data::Aquifers& aquiferSoln) override
     {
         auto xaqPos = aquiferSoln.find(this->aquiferID());
         if (xaqPos == aquiferSoln.end())
@@ -107,14 +108,17 @@ public:
         this->solution_set_from_restart_ = true;
     }
 
-    void endTimeStep()
+    void beginTimeStep() override {}
+    void addToSource(RateVector&, const unsigned, const unsigned) override {}
+
+    void endTimeStep() override
     {
         this->pressure_ = this->calculateAquiferPressure();
         this->flux_rate_ = this->calculateAquiferFluxRate();
         this->cumulative_flux_ += this->flux_rate_ * this->ebos_simulator_.timeStepSize();
     }
 
-    data::AquiferData aquiferData() const
+    data::AquiferData aquiferData() const override
     {
         data::AquiferData data;
         data.aquiferID = this->aquiferID();
@@ -128,7 +132,7 @@ public:
         return data;
     }
 
-    void initialSolutionApplied()
+    void initialSolutionApplied() override
     {
         if (this->solution_set_from_restart_) {
             return;
@@ -139,38 +143,7 @@ public:
         this->cumulative_flux_ = 0.;
     }
 
-    int aquiferID() const
-    {
-        return static_cast<int>(this->id_);
-    }
-
 private:
-    const std::size_t id_;
-    const Simulator& ebos_simulator_;
-    double flux_rate_; // aquifer influx rate
-    double cumulative_flux_; // cumulative aquifer influx
-    const int* global_cell_; // mapping to global index
-    std::vector<double> init_pressure_{};
-    double pressure_; // aquifer pressure
-    bool solution_set_from_restart_ {false};
-    bool connects_to_reservoir_ {false};
-
-    // TODO: maybe unordered_map can also do the work to save memory?
-    std::vector<int> cell_to_aquifer_cell_idx_;
-
-    inline bool co2store_() const
-    {
-        return ebos_simulator_.vanguard().eclState().runspec().co2Storage();
-    }
-
-    inline int phaseIdx_() const
-    {
-        if(co2store_())
-            return FluidSystem::oilPhaseIdx;
-
-        return FluidSystem::waterPhaseIdx;
-    }
-
     void checkConnectsToReservoir()
     {
         ElementContext elem_ctx(this->ebos_simulator_);
@@ -325,6 +298,19 @@ private:
 
         return aquifer_flux;
     }
+
+    double flux_rate_; // aquifer influx rate
+    double cumulative_flux_; // cumulative aquifer influx
+    const int* global_cell_; // mapping to global index
+    std::vector<double> init_pressure_{};
+    double pressure_; // aquifer pressure
+    bool solution_set_from_restart_ {false};
+    bool connects_to_reservoir_ {false};
+
+    // TODO: maybe unordered_map can also do the work to save memory?
+    std::vector<int> cell_to_aquifer_cell_idx_;
 };
+
 } // namespace Opm
+
 #endif