108 using Implementation = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
110 using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
111 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
112 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
113 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
114 using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
115 using GridView = GetPropType<TypeTag, Properties::GridView>;
116 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
118 enum { dimWorld = GridView::dimensionworld };
119 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
120 enum { numPhases = FluidSystem::numPhases };
121 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
122 enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
123 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
125 using Toolbox = MathToolbox<Evaluation>;
126 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
127 using EvalDimVector = Dune::FieldVector<Evaluation, dimWorld>;
128 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
136 throw std::invalid_argument(
"The ECL transmissibility module does not provide an explicit intrinsic permeability");
147 throw std::invalid_argument(
"The ECL transmissibility module does not provide explicit potential gradients");
157 {
return pressureDifference_[phaseIdx]; }
167 throw std::invalid_argument(
"The ECL transmissibility module does not provide explicit filter velocities");
180 {
return volumeFlux_[phaseIdx]; }
192 assert(phaseIdx < numPhases);
194 return upIdx_[phaseIdx];
206 assert(phaseIdx < numPhases);
208 return dnIdx_[phaseIdx];
211 void updateSolvent(
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx)
212 { asImp_().updateVolumeFluxTrans(elemCtx, scvfIdx, timeIdx); }
214 void updatePolymer(
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx)
215 { asImp_().updateShearMultipliers(elemCtx, scvfIdx, timeIdx); }
219 static void volumeAndPhasePressureDifferences(std::array<short, numPhases>& upIdx,
220 std::array<short, numPhases>& dnIdx,
222 Evaluation (&pressureDifferences)[numPhases],
223 const ElementContext& elemCtx,
227 const auto& problem = elemCtx.problem();
228 const auto& stencil = elemCtx.stencil(timeIdx);
229 const auto& scvf = stencil.interiorFace(scvfIdx);
230 unsigned interiorDofIdx = scvf.interiorIndex();
231 unsigned exteriorDofIdx = scvf.exteriorIndex();
233 assert(interiorDofIdx != exteriorDofIdx);
235 unsigned I = stencil.globalSpaceIndex(interiorDofIdx);
236 unsigned J = stencil.globalSpaceIndex(exteriorDofIdx);
237 Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
238 Scalar faceArea = scvf.area();
239 Scalar thpres = problem.thresholdPressure(I, J);
244 Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
246 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
247 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
254 Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
255 Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
259 Scalar distZ = zIn - zEx;
261 Scalar Vin = elemCtx.dofVolume(interiorDofIdx, 0);
262 Scalar Vex = elemCtx.dofVolume(exteriorDofIdx, 0);
264 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
265 if (!FluidSystem::phaseIsActive(phaseIdx))
267 calculatePhasePressureDiff_(upIdx[phaseIdx],
269 pressureDifferences[phaseIdx],
281 if (pressureDifferences[phaseIdx] == 0) {
286 const bool upwindIsInterior = (
static_cast<unsigned>(upIdx[phaseIdx]) == interiorDofIdx);
287 const IntensiveQuantities& up = upwindIsInterior ? intQuantsIn : intQuantsEx;
289 const Evaluation transMult = (intQuantsIn.rockCompTransMultiplier() + Toolbox::value(intQuantsEx.rockCompTransMultiplier()))/2;
291 const auto& materialLawManager = problem.materialLawManager();
292 FaceDir::DirEnum facedir = FaceDir::DirEnum::Unknown;
293 if (materialLawManager->hasDirectionalRelperms()) {
294 facedir = scvf.faceDirFromDirId();
296 if (upwindIsInterior)
298 pressureDifferences[phaseIdx]*up.mobility(phaseIdx, facedir)*transMult*(-trans/faceArea);
301 pressureDifferences[phaseIdx]*
302 (Toolbox::value(up.mobility(phaseIdx, facedir))*transMult*(-trans/faceArea));
306 template<
class EvalType>
307 static void calculatePhasePressureDiff_(
short& upIdx,
310 const IntensiveQuantities& intQuantsIn,
311 const IntensiveQuantities& intQuantsEx,
312 const unsigned phaseIdx,
313 const unsigned interiorDofIdx,
314 const unsigned exteriorDofIdx,
317 const unsigned globalIndexIn,
318 const unsigned globalIndexEx,
326 if (intQuantsIn.mobility(phaseIdx) <= 0.0 &&
327 intQuantsEx.mobility(phaseIdx) <= 0.0)
329 upIdx = interiorDofIdx;
330 dnIdx = exteriorDofIdx;
337 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
338 Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
339 Evaluation rhoAvg = (rhoIn + rhoEx)/2;
341 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
342 Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(phaseIdx));
345 pressureExterior += Toolbox::value(rhoAvg)*(distZg);
347 pressureExterior += rhoAvg*(distZg);
356 upIdx = exteriorDofIdx;
357 dnIdx = interiorDofIdx;
360 upIdx = interiorDofIdx;
361 dnIdx = exteriorDofIdx;
367 upIdx = interiorDofIdx;
368 dnIdx = exteriorDofIdx;
370 else if (Vin < Vex) {
371 upIdx = exteriorDofIdx;
372 dnIdx = interiorDofIdx;
378 if (globalIndexIn < globalIndexEx) {
379 upIdx = interiorDofIdx;
380 dnIdx = exteriorDofIdx;
383 upIdx = exteriorDofIdx;
384 dnIdx = interiorDofIdx;
412 Valgrind::SetUndefined(*
this);
414 volumeAndPhasePressureDifferences(upIdx_ , dnIdx_, volumeFlux_, pressureDifference_, elemCtx, scvfIdx, timeIdx);
420 template <
class Flu
idState>
424 const FluidState& exFluidState)
426 const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(scvfIdx);
427 const Scalar faceArea = scvf.area();
428 const Scalar zEx = scvf.integrationPos()[dimWorld - 1];
429 const auto& problem = elemCtx.problem();
430 const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(0, timeIdx);
431 const auto& intQuantsIn = elemCtx.intensiveQuantities(0, timeIdx);
443 pressureDifference_);
448 if constexpr (enableSolvent) {
449 if (upIdx_[gasPhaseIdx] == 0) {
450 const Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, scvfIdx);
451 const Scalar transModified = trans * Toolbox::value(intQuantsIn.rockCompTransMultiplier());
452 const auto solventFlux = pressureDifference_[gasPhaseIdx] * intQuantsIn.mobility(gasPhaseIdx) * (-transModified/faceArea);
453 asImp_().setSolventVolumeFlux(solventFlux);
455 asImp_().setSolventVolumeFlux(0.0);
464 template <
class Problem,
class Flu
idState,
class EvaluationContainer>
466 const unsigned globalSpaceIdx,
467 const IntensiveQuantities& intQuantsIn,
468 const unsigned bfIdx,
469 const double faceArea,
471 const FluidState& exFluidState,
472 std::array<short, numPhases>& upIdx,
473 std::array<short, numPhases>& dnIdx,
478 bool enableBoundaryMassFlux = problem.nonTrivialBoundaryConditions();
479 if (!enableBoundaryMassFlux)
482 Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, bfIdx);
487 Scalar g = problem.gravity()[dimWorld - 1];
494 Scalar zIn = problem.dofCenterDepth(globalSpaceIdx);
498 Scalar distZ = zIn - zEx;
500 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
501 if (!FluidSystem::phaseIsActive(phaseIdx))
506 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
507 const auto& rhoEx = exFluidState.density(phaseIdx);
508 Evaluation rhoAvg = (rhoIn + rhoEx)/2;
510 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
511 Evaluation pressureExterior = exFluidState.pressure(phaseIdx);
512 pressureExterior += rhoAvg*(distZ*g);
520 const unsigned interiorDofIdx = 0;
522 upIdx[phaseIdx] = -1;
523 dnIdx[phaseIdx] = interiorDofIdx;
526 upIdx[phaseIdx] = interiorDofIdx;
527 dnIdx[phaseIdx] = -1;
530 Evaluation transModified = trans;
532 if (upIdx[phaseIdx] == interiorDofIdx) {
536 const auto& up = intQuantsIn;
539 const Scalar transMult = Toolbox::value(up.rockCompTransMultiplier());
540 transModified *= transMult;
548 const auto& matParams = problem.materialLawParams(globalSpaceIdx);
549 std::array<typename FluidState::Scalar,numPhases> kr;
550 MaterialLaw::relativePermeabilities(kr, matParams, exFluidState);
552 const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx);
567 void calculateBoundaryFluxes_(
const ElementContext&,
unsigned,
unsigned)
571 Implementation& asImp_()
572 {
return *
static_cast<Implementation*
>(
this); }
574 const Implementation& asImp_()
const
575 {
return *
static_cast<const Implementation*
>(
this); }
578 Evaluation volumeFlux_[numPhases];
582 Evaluation pressureDifference_[numPhases];
585 std::array<short, numPhases> upIdx_;
586 std::array<short, numPhases> dnIdx_;
static void calculateBoundaryGradients_(const Problem &problem, const unsigned globalSpaceIdx, const IntensiveQuantities &intQuantsIn, const unsigned bfIdx, const double faceArea, const double zEx, const FluidState &exFluidState, std::array< short, numPhases > &upIdx, std::array< short, numPhases > &dnIdx, EvaluationContainer &volumeFlux, EvaluationContainer &pressureDifference)
Update the required gradients for boundary faces.
Definition NewTranFluxModule.hpp:465