Merge pull request #2698 from totto82/applyPVTWsalt

Pass the salt concentration to the water PVT

ebos/equil/initstateequil.hh

@@ -172,7 +172,8 @@ private:
     double
     density(const double press) const
     {
-        double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
+        double saltConcentration = 0.0; // TODO allow for non-zero initial salt concentration
+        double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, saltConcentration);
         rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
         return rho;
     }

opm/simulators/aquifers/AquiferCarterTracy.hpp

@@ -224,10 +224,11 @@ protected:
         // We use the temperature of the first cell connected to the aquifer
         // Here we copy the fluidstate of the first cell, so we do not accidentally mess up the reservoir fs
         fs_aquifer.assign(iq0.fluidState());
-        Eval temperature_aq, pa0_mean;
+        Eval temperature_aq, pa0_mean, saltConcentration_aq;
         temperature_aq = fs_aquifer.temperature(0);
+        saltConcentration_aq = fs_aquifer.saltConcentration();
         pa0_mean = this->pa0_;
-        Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean);
+        Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean, saltConcentration_aq);
         mu_w_ = mu_w_aquifer.value();
     }
 

opm/simulators/flow/MissingFeatures.cpp

@@ -124,6 +124,7 @@ namespace MissingFeatures {
             "AUTOCOAR",
             "AUTOREF",
             "BIGMODEL",
+            "BDENSITY",
             "BGGI",
             "BOGI",
             "BOUNDARY",
@@ -652,6 +653,7 @@ namespace MissingFeatures {
             "SALT",
             "SALTNODE",
             "SALTREST",
+            "SALTVD",
             "SCALELIM",
             "SCDATAB",
             "SCDETAB",

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1507,6 +1507,7 @@ namespace Opm
         // basically, it is a single value for all the segments
 
         EvalWell temperature;
+        EvalWell saltConcentration;
         // not sure how to handle the pvt region related to segment
         // for the current approach, we use the pvt region of the first perforated cell
         // although there are some text indicating using the pvt region of the lowest
@@ -1519,6 +1520,7 @@ namespace Opm
             const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
             const auto& fs = intQuants.fluidState();
             temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
+            saltConcentration = extendEval(fs.saltConcentration());
             pvt_region_index = fs.pvtRegionIndex();
         }
 
@@ -1547,9 +1549,9 @@ namespace Opm
             if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
                 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
                 b[waterCompIdx] =
-                    FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
+                    FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, saltConcentration);
                 visc[waterCompIdx] =
-                    FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure);
+                    FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure, saltConcentration);
             }
 
             EvalWell rv(0.0);
@@ -2740,6 +2742,7 @@ namespace Opm
     getSegmentSurfaceVolume(const Simulator& ebos_simulator, const int seg_idx) const
     {
         EvalWell temperature;
+        EvalWell saltConcentration;
         int pvt_region_index;
         {
             // using the pvt region of first perforated cell
@@ -2748,6 +2751,7 @@ namespace Opm
             const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
             const auto& fs = intQuants.fluidState();
             temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
+            saltConcentration = extendEval(fs.saltConcentration());
             pvt_region_index = fs.pvtRegionIndex();
         }
 
@@ -2762,7 +2766,7 @@ namespace Opm
         if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
             const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
             b[waterCompIdx] =
-                FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
+                FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, saltConcentration);
         }
 
         EvalWell rv(0.0);

opm/simulators/wells/RateConverter.hpp

@@ -448,6 +448,7 @@ namespace Opm {
                     ra.rs = 0.0;
                     ra.rv = 0.0;
                     ra.pv = 0.0;
+                    ra.saltConcentration = 0.0;
 
                 }
 
@@ -490,6 +491,7 @@ namespace Opm {
                     auto& rs  = ra.rs;
                     auto& rv  = ra.rv;
                     auto& pv  = ra.pv;
+                    auto& saltConcentration = ra.saltConcentration;
 
                     // sum p, rs, rv, and T.
                     double hydrocarbonPV = pv_cell*hydrocarbon;
@@ -499,6 +501,7 @@ namespace Opm {
                         rs += fs.Rs().value()*hydrocarbonPV;
                         rv += fs.Rv().value()*hydrocarbonPV;
                         T += fs.temperature(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
+                        saltConcentration += fs.saltConcentration().value()*hydrocarbonPV;
                     }
                 }
 
@@ -509,17 +512,20 @@ namespace Opm {
                       auto& rs  = ra.rs;
                       auto& rv  = ra.rv;
                       auto& pv  = ra.pv;
+                      auto& saltConcentration = ra.saltConcentration;
                       // communicate sums
                       p = comm.sum(p);
                       T = comm.sum(T);
                       rs = comm.sum(rs);
                       rv = comm.sum(rv);
                       pv = comm.sum(pv);
+                      saltConcentration = comm.sum(saltConcentration);
                       // compute average
                       p /= pv;
                       T /= pv;
                       rs /= pv;
                       rv /= pv;
+                      saltConcentration /=pv;
                 }
             }
 
@@ -565,6 +571,7 @@ namespace Opm {
                 const auto& ra = attr_.attributes(r);
                 const double p = ra.pressure;
                 const double T = ra.temperature;
+                const double saltConcentration = ra.saltConcentration;
 
                 const int   iw = Details::PhasePos::water(pu);
                 const int   io = Details::PhasePos::oil  (pu);
@@ -575,7 +582,7 @@ namespace Opm {
                 if (Details::PhaseUsed::water(pu)) {
                     // q[w]_r = q[w]_s / bw
 
-                    const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p);
+                    const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
 
                     coeff[iw] = 1.0 / bw;
                 }
@@ -645,6 +652,7 @@ namespace Opm {
                 const auto& ra = attr_.attributes(r);
                 const double p = ra.pressure;
                 const double T = ra.temperature;
+                const double saltConcentration = ra.saltConcentration;
 
                 const int   iw = Details::PhasePos::water(pu);
                 const int   io = Details::PhasePos::oil  (pu);
@@ -653,7 +661,7 @@ namespace Opm {
                 if (Details::PhaseUsed::water(pu)) {
                     // q[w]_r = q[w]_s / bw
 
-                    const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p);
+                    const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
 
                     voidage_rates[iw] = surface_rates[iw] / bw;
                 }
@@ -755,6 +763,7 @@ namespace Opm {
                     , rs(0.0)
                     , rv(0.0)
                     , pv(0.0)
+                    , saltConcentration(0.0)
                 {}
 
                 double pressure;
@@ -762,6 +771,7 @@ namespace Opm {
                 double rs;
                 double rv;
                 double pv;
+                double saltConcentration;
             };
 
             Details::RegionAttributes<RegionId, Attributes> attr_;

opm/simulators/wells/StandardWell_impl.hpp

@@ -1838,11 +1838,12 @@ namespace Opm
             const double p_above = perf == 0 ? well_state.bhp()[w] : well_state.perfPress()[first_perf_ + perf - 1];
             const double p_avg = (well_state.perfPress()[first_perf_ + perf] + p_above)/2;
             const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
+            const double saltConcentration = fs.saltConcentration().value();
 
             if (waterPresent) {
                 const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
                 b_perf[ waterCompIdx + perf * num_components_] =
-                FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
+                FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, saltConcentration);
             }
 
             if (gasPresent) {