28#ifndef EWOMS_CONVECTIVEMIXING_MODULE_HH
29#define EWOMS_CONVECTIVEMIXING_MODULE_HH
31#include <dune/common/fvector.hh>
33#include <opm/input/eclipse/Schedule/OilVaporizationProperties.hpp>
34#include <opm/input/eclipse/Schedule/Schedule.hpp>
36#include <opm/material/common/MathToolbox.hpp>
37#include <opm/material/common/Valgrind.hpp>
63template <
class TypeTag,
bool enableConvectiveMixing>
70template <
class TypeTag>
79 struct ConvectiveMixingModuleParam
83 static void beginEpisode(
const EclipseState&,
86 ConvectiveMixingModuleParam&)
90 template <
class Context>
95 const IntensiveQuantities&,
96 const IntensiveQuantities&,
98 const ConvectiveMixingModuleParam&)
101 template <
class Context>
113 const IntensiveQuantities&,
114 const IntensiveQuantities&,
120 const ConvectiveMixingModuleParam&)
124template <
class TypeTag>
134 using Toolbox = MathToolbox<Evaluation>;
136 enum { conti0EqIdx = Indices::conti0EqIdx };
137 enum { dimWorld = GridView::dimensionworld };
138 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
139 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
140 static constexpr bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
141 static constexpr unsigned contiEnergyEqIdx = Indices::contiEnergyEqIdx;
144 struct ConvectiveMixingModuleParam
152 static void beginEpisode(
const EclipseState& eclState,
153 const Schedule& schedule,
154 const int episodeIdx,
155 ConvectiveMixingModuleParam& info)
158 std::size_t numRegions = eclState.runspec().tabdims().getNumPVTTables();
159 const auto& control = schedule[episodeIdx].oilvap();
160 if (info.active_.empty()) {
161 info.active_.resize(numRegions);
162 info.Xhi_.resize(numRegions);
163 info.Psi_.resize(numRegions);
165 for (std::size_t i = 0; i < numRegions; ++i ) {
166 info.active_[i] = control.drsdtConvective(i);
167 if (control.drsdtConvective(i)) {
168 info.Xhi_[i] = control.getMaxDRSDT(i);
169 info.Psi_[i] = control.getPsi(i);
176 const IntensiveQuantities& intQuantsIn,
177 const IntensiveQuantities& intQuantsEx,
178 const unsigned phaseIdx,
179 const ConvectiveMixingModuleParam& info) {
181 if (info.active_.empty()) {
184 if (!info.active_[intQuantsIn.pvtRegionIndex()] || !info.active_[intQuantsEx.pvtRegionIndex()]) {
188 if (phaseIdx == FluidSystem::gasPhaseIdx) {
192 const auto& liquidPhaseIdx =
193 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
194 ? FluidSystem::waterPhaseIdx
195 : FluidSystem::oilPhaseIdx;
198 const auto& t_in = intQuantsIn.fluidState().temperature(liquidPhaseIdx);
199 const auto& p_in = intQuantsIn.fluidState().pressure(liquidPhaseIdx);
200 const auto& salt_in = intQuantsIn.fluidState().saltConcentration();
202 const auto& bLiquidIn =
203 FluidSystem::phaseIsActive(waterPhaseIdx)
204 ? FluidSystem::waterPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(),
205 t_in, p_in, Evaluation(0.0), salt_in)
206 : FluidSystem::oilPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(),
207 t_in, p_in, Evaluation(0.0));
209 const auto& refDensityLiquidIn =
210 FluidSystem::phaseIsActive(waterPhaseIdx)
211 ? FluidSystem::waterPvt().waterReferenceDensity(intQuantsIn.pvtRegionIndex())
212 : FluidSystem::oilPvt().oilReferenceDensity(intQuantsIn.pvtRegionIndex());
214 const auto& rho_in = bLiquidIn * refDensityLiquidIn;
216 const auto t_ex = Toolbox::value(intQuantsEx.fluidState().temperature(liquidPhaseIdx));
217 const auto p_ex = Toolbox::value(intQuantsEx.fluidState().pressure(liquidPhaseIdx));
218 const auto salt_ex = Toolbox::value(intQuantsEx.fluidState().saltConcentration());
220 const auto bLiquidEx =
221 FluidSystem::phaseIsActive(waterPhaseIdx)
222 ? FluidSystem::waterPvt().inverseFormationVolumeFactor(intQuantsEx.pvtRegionIndex(),
223 t_ex, p_ex, Scalar{0.0}, salt_ex)
224 : FluidSystem::oilPvt().inverseFormationVolumeFactor(intQuantsEx.pvtRegionIndex(),
225 t_ex, p_ex, Scalar{0.0});
227 const auto& refDensityLiquidEx =
228 FluidSystem::phaseIsActive(waterPhaseIdx)
229 ? FluidSystem::waterPvt().waterReferenceDensity(intQuantsEx.pvtRegionIndex())
230 : FluidSystem::oilPvt().oilReferenceDensity(intQuantsEx.pvtRegionIndex());
232 const auto rho_ex = bLiquidEx * refDensityLiquidEx;
234 rhoAvg = (rho_in + rho_ex) / 2;
237 template <
class Context>
239 const Context& elemCtx,
244 const auto& problem = elemCtx.problem();
245 const auto& stencil = elemCtx.stencil(timeIdx);
246 const auto& scvf = stencil.interiorFace(scvfIdx);
248 unsigned interiorDofIdx = scvf.interiorIndex();
249 unsigned exteriorDofIdx = scvf.exteriorIndex();
250 assert(interiorDofIdx != exteriorDofIdx);
252 const auto& globalIndexIn = stencil.globalSpaceIndex(interiorDofIdx);
253 const auto& globalIndexEx = stencil.globalSpaceIndex(exteriorDofIdx);
254 Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
255 Scalar faceArea = scvf.area();
256 const Scalar g = problem.gravity()[dimWorld - 1];
257 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
258 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
259 const Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
260 const Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
261 const Scalar distZ = zIn - zEx;
262 addConvectiveMixingFlux(flux,
270 problem.moduleParams().convectiveMixingModuleParam);
278 const IntensiveQuantities& intQuantsIn,
279 const IntensiveQuantities& intQuantsEx,
280 const unsigned globalIndexIn,
281 const unsigned globalIndexEx,
284 const Scalar faceArea,
285 const ConvectiveMixingModuleParam& info)
287 if (info.active_.empty()) {
291 if (!info.active_[intQuantsIn.pvtRegionIndex()] || !info.active_[intQuantsEx.pvtRegionIndex()]) {
295 const auto& liquidPhaseIdx =
296 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
297 ? FluidSystem::waterPhaseIdx
298 : FluidSystem::oilPhaseIdx;
301 const auto& t_in = intQuantsIn.fluidState().temperature(liquidPhaseIdx);
302 const auto& p_in = intQuantsIn.fluidState().pressure(liquidPhaseIdx);
303 const auto& rssat_in = intQuantsIn.saturatedDissolutionFactor();
304 const auto& salt_in = intQuantsIn.fluidState().saltSaturation();
306 const auto bLiquidSatIn =
307 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
308 ? FluidSystem::waterPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(),
309 t_in, p_in, rssat_in, salt_in)
310 : FluidSystem::oilPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(),
311 t_in, p_in, rssat_in);
313 const auto& densityLiquidIn =
314 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
315 ? FluidSystem::waterPvt().waterReferenceDensity(intQuantsIn.pvtRegionIndex())
316 : FluidSystem::oilPvt().oilReferenceDensity(intQuantsIn.pvtRegionIndex());
318 const auto rho_in = Opm::getValue(intQuantsIn.fluidState().invB(liquidPhaseIdx)) * densityLiquidIn;
319 const auto rho_sat_in = bLiquidSatIn *
321 rssat_in * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx,
322 intQuantsIn.pvtRegionIndex()));
325 const auto t_ex = Opm::getValue(intQuantsEx.fluidState().temperature(liquidPhaseIdx));
326 const auto p_ex = Opm::getValue(intQuantsEx.fluidState().pressure(liquidPhaseIdx));
327 const auto rssat_ex = Opm::getValue(intQuantsEx.saturatedDissolutionFactor());
328 const auto salt_ex = Opm::getValue(intQuantsEx.fluidState().saltSaturation());
329 const auto bLiquidSatEx =
330 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
331 ? FluidSystem::waterPvt().inverseFormationVolumeFactor(intQuantsEx.pvtRegionIndex(),
332 t_ex, p_ex, rssat_ex, salt_ex)
333 : FluidSystem::oilPvt().inverseFormationVolumeFactor(intQuantsEx.pvtRegionIndex(),
334 t_ex, p_ex, rssat_ex);
336 const auto& densityLiquidEx =
337 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
338 ? FluidSystem::waterPvt().waterReferenceDensity(intQuantsEx.pvtRegionIndex())
339 : FluidSystem::oilPvt().oilReferenceDensity(intQuantsEx.pvtRegionIndex());
341 const auto rho_ex = Opm::getValue(intQuantsEx.fluidState().invB(liquidPhaseIdx)) * densityLiquidEx;
342 const auto rho_sat_ex = bLiquidSatEx *
344 rssat_ex * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx,
345 intQuantsEx.pvtRegionIndex()));
347 const auto delta_rho = (rho_sat_ex + rho_sat_in - rho_in - rho_ex) / 2;
348 const auto pressure_difference_convective_mixing = delta_rho * distZg;
351 if (Opm::abs(pressure_difference_convective_mixing) > 1e-12) {
353 short interiorDofIdx = 0;
354 constexpr short exteriorDofIdx = 1;
357 if (pressure_difference_convective_mixing > 0) {
358 upIdx = exteriorDofIdx;
361 const auto& up = (upIdx == interiorDofIdx) ? intQuantsIn : intQuantsEx;
362 const auto& rssat_up = (upIdx == interiorDofIdx) ? rssat_in : rssat_ex;
363 unsigned globalUpIndex = (upIdx == interiorDofIdx) ? globalIndexIn : globalIndexEx;
365 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
366 ? up.fluidState().Rsw()
367 : up.fluidState().Rs();
369 const Evaluation& transMult = up.rockCompTransMultiplier();
370 const auto& invB = up.fluidState().invB(liquidPhaseIdx);
371 const auto& visc = up.fluidState().viscosity(liquidPhaseIdx);
374 const Evaluation RsupRestricted = Opm::min(Rsup, rssat_up*info.Psi_[up.pvtRegionIndex()]);
376 const auto convectiveFlux = -trans * transMult * info.Xhi_[up.pvtRegionIndex()] * invB *
377 pressure_difference_convective_mixing * RsupRestricted / (visc * faceArea);
378 unsigned activeGasCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx);
379 if (globalUpIndex == globalIndexIn) {
380 flux[conti0EqIdx + activeGasCompIdx] += convectiveFlux;
383 flux[conti0EqIdx + activeGasCompIdx] += Opm::getValue(convectiveFlux);
386 if constexpr (enableEnergy) {
387 const auto& h = up.fluidState().enthalpy(liquidPhaseIdx) *
388 FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, up.pvtRegionIndex());
389 if (globalUpIndex == globalIndexIn) {
390 flux[contiEnergyEqIdx] += convectiveFlux * h;
393 flux[contiEnergyEqIdx] += Opm::getValue(h) * Opm::getValue(convectiveFlux);
400template <
class TypeTag,
bool enableConvectiveMixingV>
410template <
class TypeTag>
424 const auto liquidPhaseIdx =
425 FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
426 ? FluidSystem::waterPhaseIdx
427 : FluidSystem::oilPhaseIdx;
429 const Evaluation SoMax = 0.0;
430 saturatedDissolutionFactor_ = FluidSystem::saturatedDissolutionFactor(asImp_().fluidState(),
432 asImp_().pvtRegionIndex(),
437 {
return saturatedDissolutionFactor_; }
441 {
return *
static_cast<Implementation*
>(
this); }
446template <
class TypeTag>
void updateSaturatedDissolutionFactor_()
Compute the intensive quantities needed to handle convective dissolution.
Definition: blackoilconvectivemixingmodule.hh:422
Implementation & asImp_()
Definition: blackoilconvectivemixingmodule.hh:440
const Evaluation & saturatedDissolutionFactor() const
Definition: blackoilconvectivemixingmodule.hh:436
Evaluation saturatedDissolutionFactor_
Definition: blackoilconvectivemixingmodule.hh:443
Provides the volumetric quantities required for the equations needed by the convective mixing (DRSDTC...
Definition: blackoilconvectivemixingmodule.hh:401
static void addConvectiveMixingFlux(RateVector &, const Context &, unsigned, unsigned)
Definition: blackoilconvectivemixingmodule.hh:102
static bool active(const Context &)
Definition: blackoilconvectivemixingmodule.hh:91
static void modifyAvgDensity(Evaluation &, const IntensiveQuantities &, const IntensiveQuantities &, const unsigned int, const ConvectiveMixingModuleParam &)
Definition: blackoilconvectivemixingmodule.hh:94
static void addConvectiveMixingFlux(RateVector &, const IntensiveQuantities &, const IntensiveQuantities &, const unsigned, const unsigned, const Scalar, const Scalar, const Scalar, const ConvectiveMixingModuleParam &)
Adds the convective mixing mass flux flux to the flux vector over a flux integration point.
Definition: blackoilconvectivemixingmodule.hh:112
static void addConvectiveMixingFlux(RateVector &flux, const IntensiveQuantities &intQuantsIn, const IntensiveQuantities &intQuantsEx, const unsigned globalIndexIn, const unsigned globalIndexEx, const Scalar distZg, const Scalar trans, const Scalar faceArea, const ConvectiveMixingModuleParam &info)
Adds the convective mixing mass flux flux to the flux vector over a flux integration point.
Definition: blackoilconvectivemixingmodule.hh:277
static void addConvectiveMixingFlux(RateVector &flux, const Context &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilconvectivemixingmodule.hh:238
static void modifyAvgDensity(Evaluation &rhoAvg, const IntensiveQuantities &intQuantsIn, const IntensiveQuantities &intQuantsEx, const unsigned phaseIdx, const ConvectiveMixingModuleParam &info)
Definition: blackoilconvectivemixingmodule.hh:175
Definition: blackoilconvectivemixingmodule.hh:64
Declare the properties used by the infrastructure code of the finite volume discretizations.
Defines the common properties required by the porous medium multi-phase models.
Definition: blackoilbioeffectsmodules.hh:43
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:233
std::vector< Scalar > Xhi_
Definition: blackoilconvectivemixingmodule.hh:147
std::vector< Scalar > Psi_
Definition: blackoilconvectivemixingmodule.hh:148
std::vector< bool > active_
Definition: blackoilconvectivemixingmodule.hh:146