implement Saltvd

ebos/eclequilinitializer.hh

@@ -148,6 +148,10 @@ public:
                 fluidState.setDensity(phaseIdx, rho);
 
             }
+
+            // set salt concentration
+            if (enableBrine)
+                fluidState.setSaltConcentration(initialState.saltConcentration()[elemIdx]);
         }
     }
 

ebos/eclproblem.hh

@@ -411,6 +411,7 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
     enum { enableExperiments = GET_PROP_VALUE(TypeTag, EnableExperiments) };
     enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
     enum { enablePolymer = GET_PROP_VALUE(TypeTag, EnablePolymer) };
+    enum { enableBrine = GET_PROP_VALUE(TypeTag, EnableBrine) };
     enum { enablePolymerMolarWeight = GET_PROP_VALUE(TypeTag, EnablePolymerMW) };
     enum { enableFoam = GET_PROP_VALUE(TypeTag, EnableFoam) };
     enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) };
@@ -1567,6 +1568,9 @@ public:
         if (enablePolymerMolarWeight)
             values[Indices::polymerMoleWeightIdx]= polymerMoleWeight_[globalDofIdx];
 
+        if (enableBrine)
+            values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration();
+
         values.checkDefined();
     }
 

ebos/equil/equilibrationhelpers.hh

@@ -612,6 +612,8 @@ private:
  */
 class EquilReg
 {
+    using TabulatedFunction = typename Opm::Tabulated1DFunction<double>;
+
 public:
     /**
      * Constructor.
@@ -624,10 +626,12 @@ public:
     EquilReg(const Opm::EquilRecord& rec,
              std::shared_ptr<Miscibility::RsFunction> rs,
              std::shared_ptr<Miscibility::RsFunction> rv,
+             const TabulatedFunction& saltVdTable,
              const int pvtIdx)
         : rec_    (rec)
         , rs_     (rs)
         , rv_     (rv)
+        , saltVdTable_ (saltVdTable)
         , pvtIdx_ (pvtIdx)
     {}
 
@@ -701,6 +705,8 @@ public:
     const CalcEvaporation&
     evaporationCalculator() const { return *this->rv_; }
 
+    const TabulatedFunction&
+    saltVdTable() const { return saltVdTable_;}
     /**
      * Retrieve pvtIdx of the region.
      */
@@ -710,6 +716,7 @@ private:
     Opm::EquilRecord rec_;     /**< Equilibration data */
     std::shared_ptr<Miscibility::RsFunction> rs_;      /**< RS calculator */
     std::shared_ptr<Miscibility::RsFunction> rv_;      /**< RV calculator */
+    const TabulatedFunction& saltVdTable_;
     const int pvtIdx_;
 };
 

ebos/equil/initstateequil.hh

@@ -46,6 +46,7 @@
 #include <opm/parser/eclipse/EclipseState/Tables/RvvdTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/PbvdTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/PdvdTable.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/SaltvdTable.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
 #include <opm/common/data/SimulationDataContainer.hpp>
 
@@ -148,31 +149,36 @@ namespace PhasePressODE {
 template <class FluidSystem>
 class Water
 {
+using TabulatedFunction = typename Opm::Tabulated1DFunction<double>;
 public:
     Water(const double temp,
+          const TabulatedFunction& saltVdTable,
           const int pvtRegionIdx,
           const double normGrav)
         : temp_(temp)
+        , saltVdTable_(saltVdTable)
         , pvtRegionIdx_(pvtRegionIdx)
         , g_(normGrav)
     {}
 
     double
-    operator()(const double /* depth */,
+    operator()(const double depth,
                const double press) const
     {
-        return this->density(press) * g_;
+        return this->density(depth, press) * g_;
     }
 
 private:
     const double temp_;
+    const TabulatedFunction& saltVdTable_;
     const int pvtRegionIdx_;
     const double g_;
 
     double
-    density(const double press) const
+    density(const double depth,
+            const double press) const
     {
-        double saltConcentration = 0.0; // TODO allow for non-zero initial salt concentration
+        double saltConcentration = saltVdTable_.eval(depth);
         double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, saltConcentration);
         rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
         return rho;
@@ -523,6 +529,7 @@ private:
     double gravity_;
     int    nsample_;
     double temperature_{ 273.15 + 20 };
+    double saltConcentration_{0.0};
 
     std::unique_ptr<OPress> oil_{};
     std::unique_ptr<GPress> gas_{};
@@ -749,7 +756,7 @@ makeWatPressure(const typename WPress::InitCond& ic,
                 const VSpan&                     span)
 {
     const auto drho = WatPressODE {
-        this->temperature_, reg.pvtIdx(), this->gravity_
+        this->temperature_, reg.saltVdTable(), reg.pvtIdx(), this->gravity_
     };
 
     this->wat_ = std::make_unique<WPress>(drho, ic, this->nsample_, span);
@@ -1576,6 +1583,7 @@ public:
                          const double grav = Opm::unit::gravity,
                          const bool applySwatInit = true)
         : temperature_(grid.size(/*codim=*/0)),
+          saltConcentration_(grid.size(/*codim=*/0)),
           pp_(FluidSystem::numPhases,
               std::vector<double>(grid.size(/*codim=*/0))),
           sat_(FluidSystem::numPhases,
@@ -1699,6 +1707,9 @@ public:
         // EXTRACT the initial temperature
         updateInitialTemperature_(eclipseState);
 
+        // EXTRACT the initial salt concentration
+        updateInitialSaltConcentration_(eclipseState, eqlmap, grid);
+
         // Compute pressures, saturations, rs and rv factors.
         calcPressSatRsRv(eqlmap, rec, materialLawManager, grid, grav);
 
@@ -1710,6 +1721,7 @@ public:
     typedef std::vector<Vec>    PVec; // One per phase.
 
     const Vec& temperature() const { return temperature_; }
+    const Vec& saltConcentration() const { return saltConcentration_; }
     const PVec& press() const { return pp_; }
     const PVec& saturation() const { return sat_; }
     const Vec& rs() const { return rs_; }
@@ -1721,10 +1733,31 @@ private:
         this->temperature_ = eclState.fieldProps().get_double("TEMPI");
     }
 
+    template <class RMap>
+    void updateInitialSaltConcentration_(const Opm::EclipseState& eclState, const RMap& reg, const Grid& grid)
+    {
+        const auto& tables = eclState.getTableManager();
+        const Opm::TableContainer& saltvdTables = tables.getSaltvdTables();
+        saltVdTable_.resize(saltvdTables.size());
+        for (size_t i = 0; i < saltvdTables.size(); ++i) {
+            const Opm::SaltvdTable& saltvdTable = saltvdTables.getTable<Opm::SaltvdTable>(i);
+            saltVdTable_[i].setXYContainers(saltvdTable.getDepthColumn(), saltvdTable.getSaltColumn());
+
+            const auto& cells = reg.cells(i);
+            for (const auto& cell : cells) {
+                const double depth = UgGridHelpers::cellCenterDepth(grid, cell);
+                this->saltConcentration_[cell] = saltVdTable_[i].eval(depth);
+            }
+        }
+    }
+
     std::vector< std::shared_ptr<Miscibility::RsFunction> > rsFunc_;
     std::vector< std::shared_ptr<Miscibility::RsFunction> > rvFunc_;
+    using TabulatedFunction = typename Opm::Tabulated1DFunction<double>;
+    std::vector<TabulatedFunction> saltVdTable_;
     std::vector<int> regionPvtIdx_;
     Vec temperature_;
+    Vec saltConcentration_;
     PVec pp_;
     PVec sat_;
     Vec rs_;
@@ -1768,7 +1801,7 @@ private:
             Details::verticalExtent(grid, cells, vspan);
 
             const auto eqreg = EquilReg {
-                rec[r], this->rsFunc_[r], this->rvFunc_[r], this->regionPvtIdx_[r]
+                rec[r], this->rsFunc_[r], this->rvFunc_[r], this->saltVdTable_[r], this->regionPvtIdx_[r]
             };
 
             // Ensure gas/oil and oil/water contacts are within the span for the

opm/simulators/flow/MissingFeatures.cpp

@@ -653,7 +653,6 @@ namespace MissingFeatures {
             "SALT",
             "SALTNODE",
             "SALTREST",
-            "SALTVD",
             "SCALELIM",
             "SCDATAB",
             "SCDETAB",