119 enum { dimWorld = GridView::dimensionworld };
120 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
121 enum { numPhases = FluidSystem::numPhases };
129 using Toolbox = MathToolbox<Evaluation>;
130 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
131 using EvalDimVector = Dune::FieldVector<Evaluation, dimWorld>;
132 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
142 throw std::invalid_argument(
"The ECL transmissibility module does not provide an explicit intrinsic permeability");
153 throw std::invalid_argument(
"The ECL transmissibility module does not provide explicit potential gradients");
163 {
return pressureDifference_[phaseIdx]; }
173 throw std::invalid_argument(
"The ECL transmissibility module does not provide explicit filter velocities");
186 {
return volumeFlux_[phaseIdx]; }
198 assert(phaseIdx < numPhases);
200 return upIdx_[phaseIdx];
212 assert(phaseIdx < numPhases);
214 return dnIdx_[phaseIdx];
217 void updateSolvent(
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx)
218 { asImp_().updateVolumeFluxTrans(elemCtx, scvfIdx, timeIdx); }
220 void updatePolymer(
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx)
221 { asImp_().updateShearMultipliers(elemCtx, scvfIdx, timeIdx); }
225 static void volumeAndPhasePressureDifferences(std::array<short, numPhases>& upIdx,
226 std::array<short, numPhases>& dnIdx,
228 Evaluation (&pressureDifferences)[numPhases],
229 const ElementContext& elemCtx,
233 const auto& problem = elemCtx.problem();
234 const auto& stencil = elemCtx.stencil(timeIdx);
235 const auto& scvf = stencil.interiorFace(scvfIdx);
236 unsigned interiorDofIdx = scvf.interiorIndex();
237 unsigned exteriorDofIdx = scvf.exteriorIndex();
239 assert(interiorDofIdx != exteriorDofIdx);
241 unsigned I = stencil.globalSpaceIndex(interiorDofIdx);
242 unsigned J = stencil.globalSpaceIndex(exteriorDofIdx);
243 Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
244 Scalar faceArea = scvf.area();
245 Scalar thpres = problem.thresholdPressure(I, J);
250 Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
252 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
253 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
260 Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
261 Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
265 Scalar distZ = zIn - zEx;
267 Scalar Vin = elemCtx.dofVolume(interiorDofIdx, 0);
268 Scalar Vex = elemCtx.dofVolume(exteriorDofIdx, 0);
270 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
271 if (!FluidSystem::phaseIsActive(phaseIdx))
273 calculatePhasePressureDiff_(upIdx[phaseIdx],
275 pressureDifferences[phaseIdx],
287 problem.moduleParams());
289 const bool upwindIsInterior = (
static_cast<unsigned>(upIdx[phaseIdx]) == interiorDofIdx);
290 const IntensiveQuantities& up = upwindIsInterior ? intQuantsIn : intQuantsEx;
292 const Evaluation transMult = (intQuantsIn.rockCompTransMultiplier() + Toolbox::value(intQuantsEx.rockCompTransMultiplier()))/2;
294 const auto& materialLawManager = problem.materialLawManager();
295 FaceDir::DirEnum facedir = FaceDir::DirEnum::Unknown;
296 if (materialLawManager->hasDirectionalRelperms()) {
297 facedir = scvf.faceDirFromDirId();
299 if (upwindIsInterior)
301 pressureDifferences[phaseIdx]*up.mobility(phaseIdx, facedir)*transMult*(-trans/faceArea);
304 pressureDifferences[phaseIdx]*
305 (Toolbox::value(up.mobility(phaseIdx, facedir))*transMult*(-trans/faceArea));
309 template<
class EvalType>
310 static void calculatePhasePressureDiff_(
short& upIdx,
313 const IntensiveQuantities& intQuantsIn,
314 const IntensiveQuantities& intQuantsEx,
315 const unsigned phaseIdx,
316 const unsigned interiorDofIdx,
317 const unsigned exteriorDofIdx,
320 const unsigned globalIndexIn,
321 const unsigned globalIndexEx,
324 const ModuleParams& moduleParams)
329 if (intQuantsIn.mobility(phaseIdx) <= 0.0 &&
330 intQuantsEx.mobility(phaseIdx) <= 0.0)
332 upIdx = interiorDofIdx;
333 dnIdx = exteriorDofIdx;
340 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
341 Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
342 Evaluation rhoAvg = (rhoIn + rhoEx)/2;
344 if constexpr(enableConvectiveMixing) {
345 ConvectiveMixingModule::modifyAvgDensity(rhoAvg, intQuantsIn, intQuantsEx, phaseIdx, moduleParams.convectiveMixingModuleParam);
348 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
349 Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(phaseIdx));
352 pressureExterior += Toolbox::value(rhoAvg)*(distZg);
354 pressureExterior += rhoAvg*(distZg);
363 upIdx = exteriorDofIdx;
364 dnIdx = interiorDofIdx;
367 upIdx = interiorDofIdx;
368 dnIdx = exteriorDofIdx;
374 upIdx = interiorDofIdx;
375 dnIdx = exteriorDofIdx;
377 else if (Vin < Vex) {
378 upIdx = exteriorDofIdx;
379 dnIdx = interiorDofIdx;
385 if (globalIndexIn < globalIndexEx) {
386 upIdx = interiorDofIdx;
387 dnIdx = exteriorDofIdx;
390 upIdx = exteriorDofIdx;
391 dnIdx = interiorDofIdx;
419 Valgrind::SetUndefined(*
this);
421 volumeAndPhasePressureDifferences(upIdx_ , dnIdx_, volumeFlux_, pressureDifference_, elemCtx, scvfIdx, timeIdx);
427 template <
class Flu
idState>
431 const FluidState& exFluidState)
433 const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(scvfIdx);
434 const Scalar faceArea = scvf.area();
435 const Scalar zEx = scvf.integrationPos()[dimWorld - 1];
436 const auto& problem = elemCtx.problem();
437 const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(0, timeIdx);
438 const auto& intQuantsIn = elemCtx.intensiveQuantities(0, timeIdx);
450 pressureDifference_);
455 if constexpr (enableSolvent) {
456 if (upIdx_[gasPhaseIdx] == 0) {
457 const Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, scvfIdx);
458 const Scalar transModified = trans * Toolbox::value(intQuantsIn.rockCompTransMultiplier());
459 const auto solventFlux = pressureDifference_[gasPhaseIdx] * intQuantsIn.mobility(gasPhaseIdx) * (-transModified/faceArea);
460 asImp_().setSolventVolumeFlux(solventFlux);
462 asImp_().setSolventVolumeFlux(0.0);
471 template <
class Problem,
class Flu
idState,
class EvaluationContainer>
473 const unsigned globalSpaceIdx,
474 const IntensiveQuantities& intQuantsIn,
475 const unsigned bfIdx,
476 const double faceArea,
478 const FluidState& exFluidState,
479 std::array<short, numPhases>& upIdx,
480 std::array<short, numPhases>& dnIdx,
485 bool enableBoundaryMassFlux = problem.nonTrivialBoundaryConditions();
486 if (!enableBoundaryMassFlux)
489 Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, bfIdx);
494 Scalar g = problem.gravity()[dimWorld - 1];
501 Scalar zIn = problem.dofCenterDepth(globalSpaceIdx);
505 Scalar distZ = zIn - zEx;
507 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
508 if (!FluidSystem::phaseIsActive(phaseIdx))
513 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
514 const auto& rhoEx = exFluidState.density(phaseIdx);
515 Evaluation rhoAvg = (rhoIn + rhoEx)/2;
517 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
518 Evaluation pressureExterior = exFluidState.pressure(phaseIdx);
519 pressureExterior += rhoAvg*(distZ*g);
527 const unsigned interiorDofIdx = 0;
529 upIdx[phaseIdx] = -1;
530 dnIdx[phaseIdx] = interiorDofIdx;
533 upIdx[phaseIdx] = interiorDofIdx;
534 dnIdx[phaseIdx] = -1;
537 Evaluation transModified = trans;
539 if (upIdx[phaseIdx] == interiorDofIdx) {
543 const auto& up = intQuantsIn;
546 const Scalar transMult = Toolbox::value(up.rockCompTransMultiplier());
547 transModified *= transMult;
555 const auto& matParams = problem.materialLawParams(globalSpaceIdx);
556 std::array<typename FluidState::Scalar,numPhases> kr;
557 MaterialLaw::relativePermeabilities(kr, matParams, exFluidState);
559 const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx);
574 void calculateBoundaryFluxes_(
const ElementContext&,
unsigned,
unsigned)
578 Implementation& asImp_()
579 {
return *
static_cast<Implementation*
>(
this); }
581 const Implementation& asImp_()
const
582 {
return *
static_cast<const Implementation*
>(
this); }
585 Evaluation volumeFlux_[numPhases];
589 Evaluation pressureDifference_[numPhases];
592 std::array<short, numPhases> upIdx_;
593 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:472