diff --git a/ebos/eclequilinitializer.hh b/ebos/eclequilinitializer.hh index 7513445a7..460f5c269 100644 --- a/ebos/eclequilinitializer.hh +++ b/ebos/eclequilinitializer.hh @@ -148,6 +148,10 @@ public: fluidState.setDensity(phaseIdx, rho); } + + // set salt concentration + if (enableBrine) + fluidState.setSaltConcentration(initialState.saltConcentration()[elemIdx]); } } diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index b0df9fdaf..f1b6663e8 100644 --- a/ebos/eclproblem.hh +++ b/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(); } diff --git a/ebos/equil/equilibrationhelpers.hh b/ebos/equil/equilibrationhelpers.hh index ae3410294..f5027ecbc 100644 --- a/ebos/equil/equilibrationhelpers.hh +++ b/ebos/equil/equilibrationhelpers.hh @@ -612,6 +612,8 @@ private: */ class EquilReg { + using TabulatedFunction = Opm::Tabulated1DFunction; + public: /** * Constructor. @@ -624,10 +626,12 @@ public: EquilReg(const Opm::EquilRecord& rec, std::shared_ptr rs, std::shared_ptr 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 rs_; /**< RS calculator */ std::shared_ptr rv_; /**< RV calculator */ + const TabulatedFunction& saltVdTable_; const int pvtIdx_; }; diff --git a/ebos/equil/initstateequil.hh b/ebos/equil/initstateequil.hh index 45ada5c09..0a30ecb0c 100644 --- a/ebos/equil/initstateequil.hh +++ b/ebos/equil/initstateequil.hh @@ -46,6 +46,7 @@ #include #include #include +#include #include #include @@ -148,31 +149,37 @@ namespace PhasePressODE { template class Water { +using TabulatedFunction = Opm::Tabulated1DFunction; 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 + // The initializing algorithm can give depths outside the range due to numerical noise i.e. we extrapolate + double saltConcentration = saltVdTable_.eval(depth, /*extrapolate=*/true); double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, saltConcentration); rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); return rho; @@ -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(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(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 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,42 @@ private: this->temperature_ = eclState.fieldProps().get_double("TEMPI"); } + template + void updateInitialSaltConcentration_(const Opm::EclipseState& eclState, const RMap& reg, const Grid& grid) + { + const int numEquilReg = rsFunc_.size(); + saltVdTable_.resize(numEquilReg); + const auto& tables = eclState.getTableManager(); + const Opm::TableContainer& saltvdTables = tables.getSaltvdTables(); + + // If no saltvd table is given, we create a trivial table for the density calculations + if (saltvdTables.empty()) { + std::vector x = {0.0,1.0}; + std::vector y = {0.0,0.0}; + for (auto& table : this->saltVdTable_) { + table.setXYContainers(x, y); + } + } else { + for (size_t i = 0; i < saltvdTables.size(); ++i) { + const Opm::SaltvdTable& saltvdTable = saltvdTables.getTable(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 > rsFunc_; std::vector< std::shared_ptr > rvFunc_; + using TabulatedFunction = Opm::Tabulated1DFunction; + std::vector saltVdTable_; std::vector regionPvtIdx_; Vec temperature_; + Vec saltConcentration_; PVec pp_; PVec sat_; Vec rs_; @@ -1768,7 +1812,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 diff --git a/opm/simulators/flow/MissingFeatures.cpp b/opm/simulators/flow/MissingFeatures.cpp index b02f53ece..9a19a538f 100644 --- a/opm/simulators/flow/MissingFeatures.cpp +++ b/opm/simulators/flow/MissingFeatures.cpp @@ -653,7 +653,6 @@ namespace MissingFeatures { "SALT", "SALTNODE", "SALTREST", - "SALTVD", "SCALELIM", "SCDATAB", "SCDETAB", diff --git a/tests/test_equil.cc b/tests/test_equil.cc index a278ff64c..42e36f1a4 100644 --- a/tests/test_equil.cc +++ b/tests/test_equil.cc @@ -182,6 +182,11 @@ void test_PhasePressure() using TypeTag = TTAG(TestEquilTypeTag); using FluidSystem = GET_PROP_TYPE(TypeTag, FluidSystem); + using TabulatedFunction = Opm::Tabulated1DFunction; + std::vector x = {0.0,100.0}; + std::vector y = {0.0,0.0}; + TabulatedFunction trivialSaltVdTable = TabulatedFunction(2, x, y); + auto simulator = initSimulator("equil_base.DATA"); initDefaultFluidSystem(); @@ -189,6 +194,7 @@ void test_PhasePressure() record, std::make_shared(), std::make_shared(), + trivialSaltVdTable, 0 }; @@ -235,26 +241,35 @@ void test_CellSubset() const Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ), mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) }; + using TabulatedFunction = Opm::Tabulated1DFunction; + std::vector x = {0.0,100.0}; + std::vector y = {0.0,0.0}; + TabulatedFunction trivialSaltVdTable = TabulatedFunction(2, x, y); + const Opm::EQUIL::EquilReg region[] = { Opm::EQUIL::EquilReg(record[0], std::make_shared(), std::make_shared(), + trivialSaltVdTable, 0) , Opm::EQUIL::EquilReg(record[0], std::make_shared(), std::make_shared(), + trivialSaltVdTable, 0) , Opm::EQUIL::EquilReg(record[1], std::make_shared(), std::make_shared(), + trivialSaltVdTable, 0) , Opm::EQUIL::EquilReg(record[1], std::make_shared(), std::make_shared(), + trivialSaltVdTable, 0) }; @@ -330,26 +345,35 @@ void test_RegMapping() auto simulator = initSimulator("equil_base.DATA"); initDefaultFluidSystem(); + using TabulatedFunction = Opm::Tabulated1DFunction; + std::vector x = {0.0,100.0}; + std::vector y = {0.0,0.0}; + TabulatedFunction trivialSaltVdTable = TabulatedFunction(2, x, y); + const Opm::EQUIL::EquilReg region[] = { Opm::EQUIL::EquilReg(record[0], std::make_shared(), std::make_shared(), + trivialSaltVdTable, 0) , Opm::EQUIL::EquilReg(record[0], std::make_shared(), std::make_shared(), + trivialSaltVdTable, 0) , Opm::EQUIL::EquilReg(record[1], std::make_shared(), std::make_shared(), + trivialSaltVdTable, 0) , Opm::EQUIL::EquilReg(record[1], std::make_shared(), std::make_shared(), + trivialSaltVdTable, 0) };