Refactor Phase Saturation Derivation Procedure

This commit introduces a new helper class,

    Opm::EQUIL::Details::PhaseSaturations<>

that subsumes the responsibility of the existing helper function

    Opm::EQUIL::phaseSaturations<>()

and generalises that functionality to arbitrary depth points within
single cells.  This is in preparation of adding support for the N<0
case of the initial fluid in place procedure defined in the EQUIL
keyword.  The class consumes an already equlibrated pressure table
for the pertinent equilibration region, calculates capillary
pressure values and inverts Pc curves to derive saturation values.
If the capillary pressure curves are constant within a cell, then a
simple depth consideration with respect to the implied sharp phase
interface is used to derive saturation values.  We also preserve
existing support for SWATINIT-type initialisation of the water
saturation field.

Switch InitialStateComputer<>::calcPressSatRsRv() over to using the
pressure and saturation helper classes instead of the original
helper functions since this provides additional control.  Also
remove those helper functions to reduce risk of confusion over which
method to use.  Update the unit tests accordingly.
This commit is contained in:
Bård Skaflestad 2020-04-09 16:06:29 +02:00
parent 6243e62b69
commit fd2d8536eb
2 changed files with 780 additions and 418 deletions

File diff suppressed because it is too large Load Diff

View File

@ -168,56 +168,75 @@ static Opm::EquilRecord mkEquilRecord( double datd, double datp,
return Opm::EquilRecord( datd, datp, zwoc, pcow_woc, zgoc, pcgo_goc, true, true, 0); return Opm::EquilRecord( datd, datp, zwoc, pcow_woc, zgoc, pcgo_goc, true, true, 0);
} }
template <typename Simulator>
double centerDepth(const Simulator& sim, const std::size_t cell)
{
return Opm::UgGridHelpers::cellCenterDepth(sim.vanguard().grid(), cell);
}
namespace { namespace {
void test_PhasePressure() void test_PhasePressure()
{ {
typedef std::vector<double> PVal; const auto record = mkEquilRecord( 0, 1e5, 5, 0, 0, 0 );
typedef std::vector<PVal> PPress;
auto record = mkEquilRecord( 0, 1e5, 5, 0, 0, 0 ); using TypeTag = TTAG(TestEquilTypeTag);
using FluidSystem = GET_PROP_TYPE(TypeTag, FluidSystem);
typedef TTAG(TestEquilTypeTag) TypeTag;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
auto simulator = initSimulator<TypeTag>("equil_base.DATA"); auto simulator = initSimulator<TypeTag>("equil_base.DATA");
initDefaultFluidSystem<TypeTag>(); initDefaultFluidSystem<TypeTag>();
Opm::EQUIL::EquilReg const auto region = Opm::EQUIL::EquilReg {
region(record, record,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0); 0
};
std::vector<int> cells(simulator->vanguard().grid().size(0)); auto numCells = 0;
std::iota(cells.begin(), cells.end(), 0); auto vspan = std::array<double, 2>{};
{
auto cells = std::vector<int>(simulator->vanguard().grid().size(0));
std::iota(cells.begin(), cells.end(), 0);
const double grav = 10; Opm::EQUIL::Details::verticalExtent(simulator->vanguard().grid(),
const PPress ppress = Opm::EQUIL::phasePressures<FluidSystem>(simulator->vanguard().grid(), region, cells, grav); cells, numCells, vspan);
}
const int first = 0, last = simulator->vanguard().grid().size(0) - 1; const auto grav = 10.0;
const double reltol = 1.0e-8; auto ptable = Opm::EQUIL::Details::PressureTable<
CHECK_CLOSE(ppress[FluidSystem::waterPhaseIdx][first] , 90e3 , reltol); FluidSystem, Opm::EQUIL::EquilReg
CHECK_CLOSE(ppress[FluidSystem::waterPhaseIdx][last ] , 180e3 , reltol); >{ grav };
CHECK_CLOSE(ppress[FluidSystem::oilPhaseIdx][first] , 103.5e3 , reltol);
CHECK_CLOSE(ppress[FluidSystem::oilPhaseIdx][last ] , 166.5e3 , reltol); ptable.equilibrate(region, vspan);
const auto reltol = 1.0e-8;
const auto first = centerDepth(*simulator, 0);
const auto last = centerDepth(*simulator, numCells - 1);
CHECK_CLOSE(ptable.water(first), 90e3 , reltol);
CHECK_CLOSE(ptable.water(last) , 180e3 , reltol);
CHECK_CLOSE(ptable.oil (first), 103.5e3, reltol);
CHECK_CLOSE(ptable.oil (last) , 166.5e3, reltol);
} }
void test_CellSubset() void test_CellSubset()
{ {
typedef std::vector<double> PVal; using PVal = std::vector<double>;
typedef std::vector<PVal> PPress; using PPress = std::vector<PVal>;
using TypeTag = TTAG(TestEquilTypeTag);
using FluidSystem = GET_PROP_TYPE(TypeTag, FluidSystem);
typedef TTAG(TestEquilTypeTag) TypeTag;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
auto simulator = initSimulator<TypeTag>("equil_base.DATA"); auto simulator = initSimulator<TypeTag>("equil_base.DATA");
const auto& eclipseState = simulator->vanguard().eclState(); const auto& eclipseState = simulator->vanguard().eclState();
Opm::GridManager gm(eclipseState.getInputGrid()); Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid()); const UnstructuredGrid& grid = *(gm.c_grid());
initDefaultFluidSystem<TypeTag>(); initDefaultFluidSystem<TypeTag>();
Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ), const Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) }; mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
Opm::EQUIL::EquilReg region[] = const Opm::EQUIL::EquilReg region[] =
{ {
Opm::EQUIL::EquilReg(record[0], Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
@ -260,25 +279,33 @@ void test_CellSubset()
cells[ix].push_back(c); cells[ix].push_back(c);
} }
PPress ppress(2, PVal(simulator->vanguard().grid().size(0), 0)); auto numCells = 0;
for (std::vector< std::vector<int> >::const_iterator auto vspan = std::array<double, 2>{};
r = cells.begin(), e = cells.end();
r != e; ++r)
{ {
const int rno = int(r - cells.begin()); auto vspancells = std::vector<int>(simulator->vanguard().grid().size(0));
const double grav = 10; std::iota(vspancells.begin(), vspancells.end(), 0);
const PPress p =
Opm::EQUIL::phasePressures<FluidSystem>(simulator->vanguard().grid(), region[rno], *r, grav); Opm::EQUIL::Details::verticalExtent(simulator->vanguard().grid(),
vspancells, numCells, vspan);
}
const auto grav = 10.0;
auto ptable = Opm::EQUIL::Details::PressureTable<
FluidSystem, Opm::EQUIL::EquilReg
>{ grav };
auto ppress = PPress(2, PVal(numCells, 0.0));
for (auto r = cells.begin(), e = cells.end(); r != e; ++r) {
const int rno = int(r - cells.begin());
ptable.equilibrate(region[rno], vspan);
PVal::size_type i = 0; PVal::size_type i = 0;
for (std::vector<int>::const_iterator for (auto c = r->begin(), ce = r->end(); c != ce; ++c, ++i) {
c = r->begin(), ce = r->end(); const auto depth = centerDepth(*simulator, *c);
c != ce; ++c, ++i)
{
assert (i < p[0].size());
ppress[0][*c] = p[0][i]; ppress[0][*c] = ptable.water(depth);
ppress[1][*c] = p[1][i]; ppress[1][*c] = ptable.oil (depth);
} }
} }
@ -292,18 +319,20 @@ void test_CellSubset()
void test_RegMapping() void test_RegMapping()
{ {
typedef std::vector<double> PVal; const Opm::EquilRecord record[] = {
typedef std::vector<PVal> PPress; mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ),
};
Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ), using PVal = std::vector<double>;
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) }; using PPress = std::vector<PVal>;
using TypeTag = TTAG(TestEquilTypeTag);
using FluidSystem = GET_PROP_TYPE(TypeTag, FluidSystem);
typedef TTAG(TestEquilTypeTag) TypeTag;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
auto simulator = initSimulator<TypeTag>("equil_base.DATA"); auto simulator = initSimulator<TypeTag>("equil_base.DATA");
initDefaultFluidSystem<TypeTag>(); initDefaultFluidSystem<TypeTag>();
Opm::EQUIL::EquilReg region[] = const Opm::EQUIL::EquilReg region[] =
{ {
Opm::EQUIL::EquilReg(record[0], Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
@ -326,6 +355,21 @@ void test_RegMapping()
0) 0)
}; };
auto numCells = 0;
auto vspan = std::array<double, 2>{};
{
auto cells = std::vector<int>(simulator->vanguard().grid().size(0));
std::iota(cells.begin(), cells.end(), 0);
Opm::EQUIL::Details::verticalExtent(simulator->vanguard().grid(),
cells, numCells, vspan);
}
const auto grav = 10.0;
auto ptable = Opm::EQUIL::Details::PressureTable<
FluidSystem, Opm::EQUIL::EquilReg
>{ grav };
std::vector<int> eqlnum(simulator->vanguard().grid().size(0)); std::vector<int> eqlnum(simulator->vanguard().grid().size(0));
// [ 0 1; 2 3] // [ 0 1; 2 3]
{ {
@ -347,23 +391,18 @@ void test_RegMapping()
} }
} }
Opm::RegionMapping<> eqlmap(eqlnum); const Opm::RegionMapping<> eqlmap(eqlnum);
PPress ppress(2, PVal(simulator->vanguard().grid().size(0), 0)); auto ppress = PPress(2, PVal(numCells, 0.0));
for (const auto& r : eqlmap.activeRegions()) { for (const auto& r : eqlmap.activeRegions()) {
const auto& rng = eqlmap.cells(r); ptable.equilibrate(region[r], vspan);
const int rno = r;
const double grav = 10;
const PPress p =
Opm::EQUIL::phasePressures<FluidSystem>(simulator->vanguard().grid(), region[rno], rng, grav);
PVal::size_type i = 0; PVal::size_type i = 0;
for (const auto& c : rng) { for (const auto& c : eqlmap.cells(r)) {
assert (i < p[0].size()); const auto depth = centerDepth(*simulator, c);
ppress[0][c] = p[0][i]; ppress[0][c] = ptable.water(depth);
ppress[1][c] = p[1][i]; ppress[1][c] = ptable.oil (depth);
++i; ++i;
} }