OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

changed: use updateProperty_ in updateCompositionChangeLimits_()

Browse Source
This commit is contained in:
Arne Morten Kvarving 2022-10-05 11:26:44 +02:00
parent 2d33932160
commit 26b8fe01a3
1 changed files with 58 additions and 91 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

149
ebos/eclproblem.hh
View File

@ -2108,103 +2108,70 @@ private:
{
// update the "last Rs" values for all elements, including the ones in the ghost
// and overlap regions
const auto& simulator = this->simulator();
int episodeIdx = this->episodeIndex();
std::array<bool,3> active{this->drsdtConvective_(episodeIdx),
this->drsdtActive_(episodeIdx),
this->drvdtActive_(episodeIdx)};
if (!active[0] && !active[1] && !active[2])
return;
OPM_BEGIN_PARALLEL_TRY_CATCH();
if (this->drsdtConvective_(episodeIdx)) {
// This implements the convective DRSDT as described in
// Sandve et al. "Convective dissolution in field scale CO2 storage simulations using the OPM Flow simulator"
// Submitted to TCCS 11, 2021
Scalar g = this->gravity_[dim - 1];
ElementContext elemCtx(simulator);
const auto& vanguard = simulator.vanguard();
auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
unsigned compressedDofIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const DimMatrix& perm = intrinsicPermeability(compressedDofIdx);
const Scalar permz = perm[dim - 1][dim - 1]; // The Z permeability
Scalar distZ = vanguard.cellThickness(compressedDofIdx);
const auto& iq = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = iq.fluidState();
Scalar t = getValue(fs.temperature(FluidSystem::oilPhaseIdx));
Scalar p = getValue(fs.pressure(FluidSystem::oilPhaseIdx));
Scalar so = getValue(fs.saturation(FluidSystem::oilPhaseIdx));
Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(),t,p);
Scalar saturatedInvB = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(),t,p);
Scalar rsZero = 0.0;
Scalar pureDensity = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(),t,p,rsZero) * FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex());
Scalar saturatedDensity = saturatedInvB * (FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex()) + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, fs.pvtRegionIndex()));
Scalar deltaDensity = saturatedDensity - pureDensity;
Scalar rs = getValue(fs.Rs());
Scalar visc = FluidSystem::oilPvt().viscosity(fs.pvtRegionIndex(),t,p,rs);
Scalar poro = getValue(iq.porosity());
// Note that for so = 0 this gives no limits (inf) for the dissolution rate
// Also we restrict the effect of convective mixing to positive density differences
// i.e. we only allow for fingers moving downward
this->convectiveDrs_[compressedDofIdx] = permz * rssat * max(0.0, deltaDensity) * g / ( so * visc * distZ * poro);
}
}
this->updateProperty_("EclProblem::updateCompositionChangeLimits_()) failed:",
[this,episodeIdx,active](unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
auto& simulator = this->simulator();
auto& vanguard = simulator.vanguard();
if (active[0]) {
// This implements the convective DRSDT as described in
// Sandve et al. "Convective dissolution in field scale CO2 storage simulations using the OPM Flow simulator"
// Submitted to TCCS 11, 2021
const Scalar g = this->gravity_[dim - 1];
const DimMatrix& perm = intrinsicPermeability(compressedDofIdx);
const Scalar permz = perm[dim - 1][dim - 1]; // The Z permeability
const Scalar distZ = vanguard.cellThickness(compressedDofIdx);
const auto& fs = iq.fluidState();
const Scalar t = getValue(fs.temperature(FluidSystem::oilPhaseIdx));
const Scalar p = getValue(fs.pressure(FluidSystem::oilPhaseIdx));
const Scalar so = getValue(fs.saturation(FluidSystem::oilPhaseIdx));
const Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(),t,p);
const Scalar saturatedInvB = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(),t,p);
const Scalar rsZero = 0.0;
const Scalar pureDensity = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(),t,p,rsZero) * FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex());
const Scalar saturatedDensity = saturatedInvB * (FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex()) + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, fs.pvtRegionIndex()));
const Scalar deltaDensity = saturatedDensity - pureDensity;
const Scalar rs = getValue(fs.Rs());
const Scalar visc = FluidSystem::oilPvt().viscosity(fs.pvtRegionIndex(),t,p,rs);
const Scalar poro = getValue(iq.porosity());
// Note that for so = 0 this gives no limits (inf) for the dissolution rate
// Also we restrict the effect of convective mixing to positive density differences
// i.e. we only allow for fingers moving downward
this->convectiveDrs_[compressedDofIdx] = permz * rssat * max(0.0, deltaDensity) * g / ( so * visc * distZ * poro);
}
if (this->drsdtActive_(episodeIdx)) {
ElementContext elemCtx(simulator);
const auto& vanguard = simulator.vanguard();
auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& elem = *elemIt;
if (active[1]) {
const auto& fs = iq.fluidState();
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
using FluidState = typename std::decay<decltype(fs)>::type;
unsigned compressedDofIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& iq = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = iq.fluidState();
int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
const auto& oilVaporizationControl = vanguard.schedule()[episodeIdx].oilvap();
if (oilVaporizationControl.getOption(pvtRegionIdx) || fs.saturation(gasPhaseIdx) > freeGasMinSaturation_)
this->lastRs_[compressedDofIdx] =
BlackOil::template getRs_<FluidSystem,
FluidState,
Scalar>(fs, iq.pvtRegionIndex());
else
this->lastRs_[compressedDofIdx] = std::numeric_limits<Scalar>::infinity();
}
using FluidState = typename std::decay<decltype(fs)>::type;
int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
const auto& oilVaporizationControl = simulator.vanguard().schedule()[episodeIdx].oilvap();
if (oilVaporizationControl.getOption(pvtRegionIdx) || fs.saturation(gasPhaseIdx) > freeGasMinSaturation_)
this->lastRs_[compressedDofIdx] =
BlackOil::template getRs_<FluidSystem,
FluidState,
Scalar>(fs, iq.pvtRegionIndex());
else
this->lastRs_[compressedDofIdx] = std::numeric_limits<Scalar>::infinity();
}
}
// update the "last Rv" values for all elements, including the ones in the ghost
// and overlap regions
if (this->drvdtActive_(episodeIdx)) {
ElementContext elemCtx(simulator);
const auto& vanguard = simulator.vanguard();
auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
unsigned compressedDofIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& iq = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = iq.fluidState();
using FluidState = typename std::decay<decltype(fs)>::type;
this->lastRv_[compressedDofIdx] =
BlackOil::template getRv_<FluidSystem,
FluidState,
Scalar>(fs, iq.pvtRegionIndex());
}
}
OPM_END_PARALLEL_TRY_CATCH("EclProblem::_updateCompositionLayers() failed: ", this->simulator().vanguard().grid().comm());
if (active[2]) {
const auto& fs = iq.fluidState();
using FluidState = typename std::decay<decltype(fs)>::type;
this->lastRv_[compressedDofIdx] =
BlackOil::template getRv_<FluidSystem,
FluidState,
Scalar>(fs, iq.pvtRegionIndex());
}
});
}
bool updateMaxOilSaturation_()