203 Valgrind::SetUndefined(*
this);
206 static const bool isEcfv = std::is_same<Discretization, EcfvDiscretization<TypeTag> >::value;
208 static_assert(isEcfv);
210 const auto& stencil = elemCtx.stencil(timeIdx);
211 const auto& scvf = stencil.interiorFace(scvfIdx);
213 interiorDofIdx_ = scvf.interiorIndex();
214 exteriorDofIdx_ = scvf.exteriorIndex();
215 assert(interiorDofIdx_ != exteriorDofIdx_);
217 unsigned I = stencil.globalSpaceIndex(interiorDofIdx_);
218 unsigned J = stencil.globalSpaceIndex(exteriorDofIdx_);
220 Scalar trans = transmissibility_(elemCtx, scvfIdx, timeIdx);
225 Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
227 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx_, timeIdx);
228 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx_, timeIdx);
230 Scalar zIn = dofCenterDepth_(elemCtx, interiorDofIdx_, timeIdx);
231 Scalar zEx = dofCenterDepth_(elemCtx, exteriorDofIdx_, timeIdx);
234 Scalar distZ = zIn - zEx;
236 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
237 if (!FluidSystem::phaseIsActive(phaseIdx))
242 if (intQuantsIn.mobility(phaseIdx) <= 0.0 &&
243 intQuantsEx.mobility(phaseIdx) <= 0.0)
245 upIdx_[phaseIdx] = interiorDofIdx_;
246 dnIdx_[phaseIdx] = exteriorDofIdx_;
247 pressureDifference_[phaseIdx] = 0.0;
248 volumeFlux_[phaseIdx] = 0.0;
254 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
255 Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
256 Evaluation rhoAvg = (rhoIn + rhoEx)/2;
258 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
259 Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(phaseIdx));
261 pressureExterior += rhoAvg*(distZ*g);
263 pressureDifference_[phaseIdx] = pressureExterior - pressureInterior;
269 if (pressureDifference_[phaseIdx] > 0.0) {
270 upIdx_[phaseIdx] = exteriorDofIdx_;
271 dnIdx_[phaseIdx] = interiorDofIdx_;
273 else if (pressureDifference_[phaseIdx] < 0.0) {
274 upIdx_[phaseIdx] = interiorDofIdx_;
275 dnIdx_[phaseIdx] = exteriorDofIdx_;
280 Scalar Vin = elemCtx.dofVolume(interiorDofIdx_, 0);
281 Scalar Vex = elemCtx.dofVolume(exteriorDofIdx_, 0);
283 upIdx_[phaseIdx] = interiorDofIdx_;
284 dnIdx_[phaseIdx] = exteriorDofIdx_;
286 else if (Vin < Vex) {
287 upIdx_[phaseIdx] = exteriorDofIdx_;
288 dnIdx_[phaseIdx] = interiorDofIdx_;
295 upIdx_[phaseIdx] = interiorDofIdx_;
296 dnIdx_[phaseIdx] = exteriorDofIdx_;
299 upIdx_[phaseIdx] = exteriorDofIdx_;
300 dnIdx_[phaseIdx] = interiorDofIdx_;
309 const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx);
311 if (upstreamIdx == interiorDofIdx_)
312 volumeFlux_[phaseIdx] =
313 pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-trans);
315 volumeFlux_[phaseIdx] =
316 pressureDifference_[phaseIdx]*(Toolbox::value(up.mobility(phaseIdx))*(-trans));
327 const FluidState& exFluidState)
329 const auto& stencil = elemCtx.stencil(timeIdx);
330 const auto& scvf = stencil.boundaryFace(scvfIdx);
332 interiorDofIdx_ = scvf.interiorIndex();
334 Scalar trans = transmissibilityBoundary_(elemCtx, scvfIdx, timeIdx);
339 Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
341 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx_, timeIdx);
348 Scalar zIn = dofCenterDepth_(elemCtx, interiorDofIdx_, timeIdx);
349 Scalar zEx = scvf.integrationPos()[dimWorld - 1];
353 Scalar distZ = zIn - zEx;
355 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
356 if (!FluidSystem::phaseIsActive(phaseIdx))
361 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
362 const auto& rhoEx = exFluidState.density(phaseIdx);
363 Evaluation rhoAvg = (rhoIn + rhoEx)/2;
365 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
366 Evaluation pressureExterior = exFluidState.pressure(phaseIdx);
367 pressureExterior += rhoAvg*(distZ*g);
369 pressureDifference_[phaseIdx] = pressureExterior - pressureInterior;
375 if (pressureDifference_[phaseIdx] > 0.0) {
376 upIdx_[phaseIdx] = -1;
377 dnIdx_[phaseIdx] = interiorDofIdx_;
380 upIdx_[phaseIdx] = interiorDofIdx_;
381 dnIdx_[phaseIdx] = -1;
385 if (upstreamIdx == interiorDofIdx_) {
390 const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx);
392 volumeFlux_[phaseIdx] =
393 pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-trans);
399 const auto& matParams =
400 elemCtx.problem().materialLawParams(elemCtx,
403 std::array<typename FluidState::Scalar,numPhases> kr;
404 MaterialLaw::relativePermeabilities(kr, matParams, exFluidState);
406 const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx);
407 volumeFlux_[phaseIdx] =
408 pressureDifference_[phaseIdx]*mob*(-trans);