Merge pull request #2197 from bska/load-aquifer-data

Load aquifer data

ebos/eclbaseaquifermodel.hh

@@ -27,8 +27,14 @@
 #ifndef EWOMS_ECL_BASE_AQUIFER_MODEL_HH
 #define EWOMS_ECL_BASE_AQUIFER_MODEL_HH
 
+#include <opm/output/data/Aquifer.hpp>
+
 #include <opm/models/utils/propertysystem.hh>
 
+#include <exception>
+#include <stdexcept>
+#include <vector>
+
 BEGIN_PROPERTIES
 
 NEW_PROP_TAG(Simulator);
@@ -64,6 +70,23 @@ public:
     void initialSolutionApplied()
     { }
 
+    /*!
+     * \brief Called if aquifers are being initialized from values retrieved
+     *        from a restart file.
+     *
+     * \param[in] aquiferSoln Set of aquifer-related initial values, mostly
+     *        pertaining to analytic aquifers.  Contains at minimum the
+     *        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)
+    {
+        throw std::logic_error {
+            "Initialization from restart data not supported "
+            "for base aquifer model"
+        };
+    }
+
     /*!
      * \brief This method is called when a new episode (report step) starts.
      */

ebos/eclproblem.hh

@@ -1748,6 +1748,9 @@ public:
     EclWellModel& wellModel()
     { return wellModel_; }
 
+    EclAquiferModel& mutableAquiferModel()
+    { return aquiferModel_; }
+
     // temporary solution to facilitate output of initial state from flow
     const InitialFluidState& initialFluidState(unsigned globalDofIdx) const
     { return initialFluidStates_[globalDofIdx]; }

ebos/eclwriter.hh

@@ -454,6 +454,9 @@ public:
 
             // initialize the well model from restart values
             simulator_.problem().wellModel().initFromRestartFile(restartValues);
+
+            if (!restartValues.aquifer.empty())
+                simulator_.problem().mutableAquiferModel().initFromRestart(restartValues.aquifer);
         }
     }
 

opm/simulators/aquifers/AquiferCarterTracy.hpp

@@ -23,6 +23,11 @@
 
 #include <opm/simulators/aquifers/AquiferInterface.hpp>
 
+#include <opm/output/data/Aquifer.hpp>
+
+#include <exception>
+#include <stdexcept>
+
 namespace Opm
 {
 
@@ -145,6 +150,14 @@ namespace Opm
                 }
             }
 
+            void assignRestartData(const data::AquiferData& /* xaq */) override
+            {
+                throw std::runtime_error {
+                    "Restart-based initialization not currently supported "
+                    "for Carter-Tracey analytic aquifers"
+                };
+            }
+
             inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td)
             {
                 // We use the opm-common numeric linear interpolator

opm/simulators/aquifers/AquiferFetkovich.hpp

@@ -23,6 +23,11 @@ along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 
 #include <opm/simulators/aquifers/AquiferInterface.hpp>
 
+#include <opm/output/data/Aquifer.hpp>
+
+#include <exception>
+#include <stdexcept>
+
 namespace Opm
 {
 
@@ -146,6 +151,19 @@ namespace Opm
       }
     }
 
+    void assignRestartData(const data::AquiferData& xaq) override
+    {
+      if (xaq.type != data::AquiferType::Fetkovich)
+      {
+        throw std::invalid_argument {
+          "Analytic aquifer data for unexpected aquifer type "
+          "passed to Fetkovich aquifer"
+        };
+      }
+
+      this->aquifer_pressure_ = xaq.pressure;
+    }
+
     inline Eval dpai(int idx)
     {
       const Eval dp = aquifer_pressure_ - Base::pressure_current_.at(idx)
@@ -176,6 +194,11 @@ namespace Opm
     inline void calculateAquiferCondition()
     {
       Base::rhow_.resize(Base::cell_idx_.size(),0.);
+
+      if (this->solution_set_from_restart_) {
+        return;
+      }
+
       if (!aqufetp_data_.p0)
       {
         Base::pa0_ = calculateReservoirEquilibrium();

opm/simulators/aquifers/AquiferInterface.hpp

@@ -27,6 +27,8 @@
 #include <opm/parser/eclipse/EclipseState/Aquancon.hpp>
 #include <opm/common/utility/numeric/linearInterpolation.hpp>
 
+#include <opm/output/data/Aquifer.hpp>
+
 #include <opm/material/common/MathToolbox.hpp>
 #include <opm/material/densead/Math.hpp>
 #include <opm/material/densead/Evaluation.hpp>
@@ -74,6 +76,26 @@ namespace Opm
     // Deconstructor
     virtual ~AquiferInterface() {}
 
+    void initFromRestart(const std::vector<data::AquiferData>& aquiferSoln)
+    {
+      auto xaqPos = std::find_if(aquiferSoln.begin(), aquiferSoln.end(),
+        [this](const data::AquiferData& xaq) -> bool
+      {
+        return xaq.aquiferID == this->connection_.aquiferID;
+      });
+
+      if (xaqPos == aquiferSoln.end()) {
+        // No restart value applies to this aquifer.  Nothing to do.
+        return;
+      }
+
+      this->assignRestartData(*xaqPos);
+
+      this->W_flux_ = xaqPos->volume;
+      this->pa0_    = xaqPos->initPressure;
+      this->solution_set_from_restart_ = true;
+    }
+
     void initialSolutionApplied()
     {
       initQuantities(connection_);
@@ -129,7 +151,10 @@ namespace Opm
     inline void initQuantities(const Aquancon::AquanconOutput& connection)
     {
       // We reset the cumulative flux at the start of any simulation, so, W_flux = 0
-      W_flux_ = 0.;
+      if (!this->solution_set_from_restart_)
+      {
+        W_flux_ = 0.;
+      }
 
       // We next get our connections to the aquifer and initialize these quantities using the initialize_connections function
       initializeConnections(connection);
@@ -207,8 +232,12 @@ namespace Opm
 
     Eval W_flux_;
 
+    bool solution_set_from_restart_{false};
+
     virtual void initializeConnections(const Aquancon::AquanconOutput& connection) =0;
 
+    virtual void assignRestartData(const data::AquiferData& xaq) = 0;
+
     virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0;
 
     virtual void calculateAquiferCondition() = 0;

opm/simulators/aquifers/BlackoilAquiferModel.hpp

@@ -29,11 +29,16 @@
 #include <opm/parser/eclipse/EclipseState/AquiferCT.hpp>
 #include <opm/parser/eclipse/EclipseState/Aquifetp.hpp>
 #include <opm/parser/eclipse/EclipseState/Aquancon.hpp>
-#include <opm/simulators/timestepping/SimulatorTimer.hpp>
+
+#include <opm/output/data/Aquifer.hpp>
+
 #include <opm/simulators/aquifers/AquiferCarterTracy.hpp>
 #include <opm/simulators/aquifers/AquiferFetkovich.hpp>
+
 #include <opm/material/densead/Math.hpp>
 
+#include <vector>
+
 namespace Opm {
 
         /// Class for handling the blackoil well model.
@@ -47,6 +52,8 @@ namespace Opm {
             explicit BlackoilAquiferModel(Simulator& simulator);
 
             void initialSolutionApplied();
+            void initFromRestart(const std::vector<data::AquiferData>& aquiferSoln);
+
             void beginEpisode();
             void beginTimeStep();
             void beginIteration();

opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp

@@ -29,6 +29,26 @@ namespace Opm {
     }
   }
 
+  template<typename TypeTag>
+  void
+  BlackoilAquiferModel<TypeTag>::initFromRestart(const std::vector<data::AquiferData>& aquiferSoln)
+  {
+    if(aquiferCarterTracyActive())
+    {
+      for (auto& aquifer : aquifers_CarterTracy)
+      {
+        aquifer.initFromRestart(aquiferSoln);
+      }
+    }
+    if(aquiferFetkovichActive())
+    {
+      for (auto& aquifer : aquifers_Fetkovich)
+      {
+        aquifer.initFromRestart(aquiferSoln);
+      }
+    }
+  }
+
   template<typename TypeTag>
   void
   BlackoilAquiferModel<TypeTag>::beginEpisode()