Merge pull request #2699 from totto82/pvtwsaltSaltvd

Implement SALTVD
This commit is contained in:
Bård Skaflestad 2020-08-18 12:32:26 +02:00 committed by GitHub
commit af9a9fa512
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 89 additions and 7 deletions

View File

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

View File

@ -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();
}

View File

@ -612,6 +612,8 @@ private:
*/
class EquilReg
{
using TabulatedFunction = 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_;
};

View File

@ -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,37 @@ namespace PhasePressODE {
template <class FluidSystem>
class Water
{
using TabulatedFunction = 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
// 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<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,42 @@ private:
this->temperature_ = eclState.fieldProps().get_double("TEMPI");
}
template <class RMap>
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<double> x = {0.0,1.0};
std::vector<double> 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<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 = Opm::Tabulated1DFunction<double>;
std::vector<TabulatedFunction> saltVdTable_;
std::vector<int> 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

View File

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

View File

@ -182,6 +182,11 @@ void test_PhasePressure()
using TypeTag = TTAG(TestEquilTypeTag);
using FluidSystem = GET_PROP_TYPE(TypeTag, FluidSystem);
using TabulatedFunction = Opm::Tabulated1DFunction<double>;
std::vector<double> x = {0.0,100.0};
std::vector<double> y = {0.0,0.0};
TabulatedFunction trivialSaltVdTable = TabulatedFunction(2, x, y);
auto simulator = initSimulator<TypeTag>("equil_base.DATA");
initDefaultFluidSystem<TypeTag>();
@ -189,6 +194,7 @@ void test_PhasePressure()
record,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
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<double>;
std::vector<double> x = {0.0,100.0};
std::vector<double> y = {0.0,0.0};
TabulatedFunction trivialSaltVdTable = TabulatedFunction(2, x, y);
const Opm::EQUIL::EquilReg region[] =
{
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
};
@ -330,26 +345,35 @@ void test_RegMapping()
auto simulator = initSimulator<TypeTag>("equil_base.DATA");
initDefaultFluidSystem<TypeTag>();
using TabulatedFunction = Opm::Tabulated1DFunction<double>;
std::vector<double> x = {0.0,100.0};
std::vector<double> y = {0.0,0.0};
TabulatedFunction trivialSaltVdTable = TabulatedFunction(2, x, y);
const Opm::EQUIL::EquilReg region[] =
{
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
};