24#ifndef DUMUX_FVVELOCITY2P_ADAPTIVE_HH
25#define DUMUX_FVVELOCITY2P_ADAPTIVE_HH
27#include <dune/common/float_cmp.hh>
38template<
class TypeTag>
48 using Indices =
typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
50 using FluidSystem =
typename GET_PROP_TYPE(TypeTag, FluidSystem);
51 using FluidState =
typename GET_PROP_TYPE(TypeTag, FluidState);
55 using Intersection =
typename GridView::Intersection;
59 dim = GridView::dimension, dimWorld = GridView::dimensionworld
62 using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
66 pw = Indices::pressureW,
67 pn = Indices::pressureNw,
68 pGlobal = Indices::pressureGlobal,
69 vw = Indices::velocityW,
70 vn = Indices::velocityNw,
71 vt = Indices::velocityTotal
78 using Element =
typename GridView::template Codim<0>::Entity;
80 using GlobalPosition =
typename Element::Geometry::GlobalCoordinate;
81 using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
82 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
90 :
ParentType(problem), problem_(problem), gravity_(problem.gravity())
94 DUNE_THROW(Dune::NotImplemented,
"Total velocity - global pressure - model cannot be used with compressible fluids!");
96 if (velocityType_ != vw && velocityType_ != vn && velocityType_ != vt)
98 DUNE_THROW(Dune::NotImplemented,
"Velocity type not supported!");
101 density_[wPhaseIdx] = 0.;
102 density_[nPhaseIdx] = 0.;
103 viscosity_[wPhaseIdx] = 0.;
104 viscosity_[nPhaseIdx] = 0.;
112 if (!compressibility_)
114 const auto element = *problem_.gridView().template begin<0>();
115 FluidState fluidState;
116 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
117 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
118 fluidState.setTemperature(problem_.temperature(element));
119 fluidState.setSaturation(wPhaseIdx, 1.);
120 fluidState.setSaturation(nPhaseIdx, 0.);
144 const GravityVector& gravity_;
145 Scalar density_[numPhases];
146 Scalar viscosity_[numPhases];
152 static const int pressureType_ =
GET_PROP_VALUE(TypeTag, PressureFormulation);
164template<
class TypeTag>
167 auto elementI = intersection.inside();
168 auto elementJ = intersection.outside();
170 if (elementI.level() == elementJ.level())
172 ParentType::calculateVelocity(intersection, cellData);
174 else if (elementJ.level() == elementI.level() + 1 && dim == 2)
176 CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(elementJ));
179 const GlobalPosition& globalPosI = elementI.geometry().center();
180 const GlobalPosition& globalPosJ = elementJ.geometry().center();
183 Scalar lambdaWI = cellData.mobility(wPhaseIdx);
184 Scalar lambdaNwI = cellData.mobility(nPhaseIdx);
185 Scalar fractionalWI = cellData.fracFlowFunc(wPhaseIdx);
186 Scalar fractionalNwI = cellData.fracFlowFunc(nPhaseIdx);
187 Scalar lambdaWJ = cellDataJ.mobility(wPhaseIdx);
188 Scalar lambdaNwJ = cellDataJ.mobility(nPhaseIdx);
189 Scalar fractionalWJ = cellDataJ.fracFlowFunc(wPhaseIdx);
190 Scalar fractionalNwJ = cellDataJ.fracFlowFunc(nPhaseIdx);
193 Scalar pcI = cellData.capillaryPressure();
194 Scalar pcJ = cellDataJ.capillaryPressure();
197 int isIndexI = intersection.indexInInside();
199 Scalar faceArea = intersection.geometry().volume();
200 Scalar faceAreaSum = faceArea;
203 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
207 int isIndexJ = intersection.indexInOutside();
210 auto elementK = intersection.outside();
217 for (
const auto& intersectionI :
intersections(problem_.gridView(), elementI))
219 if (intersectionI.neighbor())
221 auto neighbor2 = intersectionI.outside();
225 if (neighbor2 != elementJ && intersectionI.indexInInside() == isIndexI)
227 globalIdxK = problem_.variables().index(neighbor2);
228 elementK = neighbor2;
229 faceAreaSum += intersectionI.geometry().volume();
236 CellData& cellDataK = problem_.variables().cellData(globalIdxK);
239 const GlobalPosition& globalPosInterface = intersection.geometry().center();
241 Dune::FieldVector < Scalar, dimWorld > distVec = globalPosInterface - globalPosI;
242 Scalar lI = distVec * unitOuterNormal;
243 distVec = globalPosJ - globalPosInterface;
244 Scalar lJ = distVec * unitOuterNormal;
247 DimMatrix permeabilityI(0);
248 DimMatrix permeabilityJ(0);
249 DimMatrix permeabilityK(0);
251 problem_.spatialParams().meanK(permeabilityI,
252 problem_.spatialParams().intrinsicPermeability(elementI));
253 problem_.spatialParams().meanK(permeabilityJ,
254 problem_.spatialParams().intrinsicPermeability(elementJ));
255 problem_.spatialParams().meanK(permeabilityK,
256 problem_.spatialParams().intrinsicPermeability(elementK));
259 Scalar kI, kJ, kK, ng, kMean;
260 Dune::FieldVector < Scalar, dim > permI(0);
261 Dune::FieldVector < Scalar, dim > permJ(0);
262 Dune::FieldVector < Scalar, dim > permK(0);
264 permeabilityI.mv(unitOuterNormal, permI);
265 permeabilityJ.mv(unitOuterNormal, permJ);
266 permeabilityK.mv(unitOuterNormal, permK);
269 kI = unitOuterNormal * permI;
270 kJ = unitOuterNormal * permJ;
271 kK = unitOuterNormal * permK;
272 kMean = kI * kJ * kK * l / (kJ * kK * lI + kI * (kJ + kK) / 2 * lJ);
274 ng = gravity_ * unitOuterNormal;
276 Scalar fractionalWK = cellDataK.fracFlowFunc(wPhaseIdx);
277 Scalar fractionalNwK = cellDataK.fracFlowFunc(nPhaseIdx);
279 Scalar pcK = cellDataK.capillaryPressure();
280 Scalar pcJK = (pcJ + pcK) / 2;
284 Scalar potentialDiffWIJ = cellData.fluxData().upwindPotential(wPhaseIdx, isIndexI);
285 Scalar potentialDiffNwIJ = cellData.fluxData().upwindPotential(nPhaseIdx, isIndexI);
286 Scalar potentialDiffWIK = potentialDiffWIJ;
287 Scalar potentialDiffNwIK = potentialDiffNwIJ;
291 Scalar rhoMeanWIJ = density_[wPhaseIdx];
292 Scalar rhoMeanNwIJ = density_[nPhaseIdx];
293 Scalar rhoMeanWIK = density_[wPhaseIdx];
294 Scalar rhoMeanNwIK = density_[nPhaseIdx];
296 if (compressibility_)
298 rhoMeanWIJ = (lJ * cellData.density(wPhaseIdx) + lI * cellDataJ.density(wPhaseIdx)) / l;
299 rhoMeanNwIJ = (lJ * cellData.density(nPhaseIdx) + lI * cellDataJ.density(nPhaseIdx)) / l;
300 rhoMeanWIK = (lJ * cellData.density(wPhaseIdx) + lI * cellDataK.density(wPhaseIdx)) / l;
301 rhoMeanNwIK = (lJ * cellData.density(nPhaseIdx) + lI * cellDataK.density(nPhaseIdx)) / l;
304 Scalar fMeanWIJ = (lJ * fractionalWI + lI * fractionalWJ) / l;
305 Scalar fMeanNwIJ = (lJ * fractionalNwI + lI * fractionalNwJ) / l;
306 Scalar fMeanWIK = (lJ * fractionalWI + lI * fractionalWK) / l;
307 Scalar fMeanNwIK = (lJ * fractionalNwI + lI * fractionalNwK) / l;
309 Scalar densityWIJ = density_[wPhaseIdx];
310 Scalar densityNwIJ = density_[nPhaseIdx];
311 Scalar densityWIK = density_[wPhaseIdx];
312 Scalar densityNwIK = density_[nPhaseIdx];
314 if (compressibility_)
317 densityWIJ = (potentialDiffWIJ > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
318 densityNwIJ = (potentialDiffNwIJ > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
319 densityWIK = (potentialDiffWIK > 0.) ? cellData.density(wPhaseIdx) : cellDataK.density(wPhaseIdx);
320 densityNwIK = (potentialDiffNwIK > 0.) ? cellData.density(nPhaseIdx) : cellDataK.density(nPhaseIdx);
322 densityWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIJ, 0.0, 1.0e-30)) ? rhoMeanWIJ : densityWIJ;
323 densityNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIJ, 0.0, 1.0e-30)) ? rhoMeanNwIJ : densityNwIJ;
324 densityWIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIK, 0.0, 1.0e-30)) ? rhoMeanWIK : densityWIK;
325 densityNwIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIK, 0.0, 1.0e-30)) ? rhoMeanNwIK : densityNwIK;
328 Scalar fractionalWIJ = (potentialDiffWIJ > 0.) ? fractionalWI : fractionalWJ;
329 Scalar fractionalNwIJ = (potentialDiffNwIJ > 0.) ? fractionalNwI : fractionalNwJ;
330 Scalar fractionalWIK = (potentialDiffWIK > 0.) ? fractionalWI : fractionalWK;
331 Scalar fractionalNwIK = (potentialDiffNwIK > 0.) ? fractionalNwI : fractionalNwK;
333 fractionalWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIJ, 0.0, 1.0e-30)) ? fMeanWIJ : fractionalWIJ;
334 fractionalNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIJ, 0.0, 1.0e-30)) ? fMeanNwIJ : fractionalNwIJ;
335 fractionalWIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIK, 0.0, 1.0e-30)) ? fMeanWIK : fractionalWIK;
336 fractionalNwIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIK, 0.0, 1.0e-30)) ? fMeanNwIK : fractionalNwIK;
338 switch (pressureType_)
342 Scalar pressJK = (cellDataJ.globalPressure() + cellDataK.globalPressure()) / 2;
344 potentialDiffWIJ = (cellData.globalPressure() - fMeanNwIJ * pcI - (pressJK - fMeanNwIJ * pcJK)) / l
345 + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng;
346 potentialDiffNwIJ = (cellData.globalPressure() + fMeanWIJ * pcI - (pressJK + fMeanWIJ * pcJK)) / l
347 + (densityNwIJ - lJ / l * (kI + kK) / kI * (densityNwIK - densityNwIJ) / 2) * ng;
348 potentialDiffWIK = (cellData.globalPressure() - fMeanNwIK * pcI - (pressJK - fMeanNwIK * pcJK)) / l
349 + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng;
350 potentialDiffNwIK = (cellData.globalPressure() + fMeanWIK * pcI - (pressJK + fMeanWIK * pcJK)) / l
351 + (densityNwIK - lJ / l * (kI + kK) / kI * (densityNwIJ - densityNwIK) / 2) * ng;
356 potentialDiffWIJ = (cellData.pressure(wPhaseIdx)
357 - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l
358 + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng;
359 potentialDiffNwIJ = (cellData.pressure(nPhaseIdx)
360 - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l
361 + (densityNwIJ - lJ / l * (kI + kK) / kI * (densityNwIK - densityNwIJ) / 2) * ng;
362 potentialDiffWIK = (cellData.pressure(wPhaseIdx)
363 - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l
364 + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng;
365 potentialDiffNwIK = (cellData.pressure(nPhaseIdx)
366 - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l
367 + (densityNwIK - lJ / l * (kI + kK) / kI * (densityNwIJ - densityNwIK) / 2) * ng;
375 cellData.fluxData().addUpwindPotential(wPhaseIdx, isIndexI, potentialDiffWIJ);
376 cellData.fluxData().addUpwindPotential(nPhaseIdx, isIndexI, potentialDiffNwIJ);
377 cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, isIndexJ, -potentialDiffWIJ);
378 cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, isIndexJ, -potentialDiffNwIJ);
381 Scalar lambdaWIJ = (potentialDiffWIJ > 0.) ? lambdaWI : lambdaWJ;
382 lambdaWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIJ, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
383 Scalar lambdaNwIJ = (potentialDiffNwIJ > 0.) ? lambdaNwI : lambdaNwJ;
384 lambdaNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIJ, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNwIJ;
386 if (compressibility_)
388 densityWIJ = (potentialDiffWIJ > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
389 densityNwIJ = (potentialDiffNwIJ > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
390 densityWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIJ, 0.0, 1.0e-30)) ? rhoMeanWIJ : densityWIJ;
391 densityNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIJ, 0.0, 1.0e-30)) ? rhoMeanNwIJ : densityNwIJ;
392 densityWIK = (potentialDiffWIK > 0.) ? cellData.density(wPhaseIdx) : cellDataK.density(wPhaseIdx);
393 densityNwIK = (potentialDiffNwIK > 0.) ? cellData.density(nPhaseIdx) : cellDataK.density(nPhaseIdx);
394 densityWIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIK, 0.0, 1.0e-30)) ? rhoMeanWIK : densityWIK;
395 densityNwIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIK, 0.0, 1.0e-30)) ? rhoMeanNwIK : densityNwIK;
399 Dune::FieldVector < Scalar, dimWorld > velocityW(unitOuterNormal);
400 Dune::FieldVector < Scalar, dimWorld > velocityNw(unitOuterNormal);
401 Dune::FieldVector < Scalar, dimWorld > gravityTermW(unitOuterNormal);
402 Dune::FieldVector < Scalar, dimWorld > gravityTermNw(unitOuterNormal);
404 gravityTermW *= lambdaWIJ * kMean * ng;
405 gravityTermW *= densityWIJ - (lJ / l) * (kI + kK) / kI * (densityWIK - densityWIJ) / 2;
406 gravityTermNw *= lambdaNwIJ * kMean * ng;
407 gravityTermNw *= densityNwIJ - (lJ / l) * (kI + kK) / kI * (densityNwIK - densityNwIJ) / 2;
409 switch (pressureType_)
413 velocityW *= lambdaWIJ * kMean * (cellData.pressure(wPhaseIdx) -
414 (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx)) / 2.0) / l;
415 velocityNw *= lambdaNwIJ * kMean * (cellData.pressure(nPhaseIdx) -
416 (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)) / 2.0) / l;
418 velocityW += gravityTermW;
419 velocityNw += gravityTermNw;
424 velocityNw *= lambdaNwIJ * kMean * (cellData.pressure(nPhaseIdx) -
425 (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)) / 2.0) / l;
426 velocityW *= lambdaWIJ * kMean * (cellData.pressure(wPhaseIdx) -
427 (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx)) / 2.0) / l;
429 velocityW += gravityTermW;
430 velocityNw += gravityTermNw;
435 Scalar pressJK = (cellDataJ.globalPressure() + cellDataK.globalPressure()) / 2;
437 velocityW *= lambdaWIJ * kMean * (cellData.globalPressure() - pressJK) / l;
438 velocityW += gravityTermW;
439 velocityW += gravityTermNw;
445 cellDataJ.fluxData().setVelocity(wPhaseIdx, isIndexJ, velocityW);
446 cellDataJ.fluxData().setVelocity(nPhaseIdx, isIndexJ, velocityNw);
447 cellDataJ.fluxData().setVelocityMarker(isIndexJ);
450 velocityW *= faceArea/faceAreaSum;
451 velocityNw *= faceArea/faceAreaSum;
452 cellData.fluxData().addVelocity(wPhaseIdx, isIndexI, velocityW);
453 cellData.fluxData().addVelocity(nPhaseIdx, isIndexI, velocityNw);
455 else if (elementI.level() > elementJ.level() && dim == 3)
457 int globalIdxJ = problem_.variables().index(elementJ);
459 CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
462 const GlobalPosition& globalPosI = elementI.geometry().center();
463 const GlobalPosition& globalPosJ = elementJ.geometry().center();
466 Scalar lambdaWI = cellData.mobility(wPhaseIdx);
467 Scalar lambdaNwI = cellData.mobility(nPhaseIdx);
468 Scalar lambdaWJ = cellDataJ.mobility(wPhaseIdx);
469 Scalar lambdaNwJ = cellDataJ.mobility(nPhaseIdx);
472 Scalar pcI = cellData.capillaryPressure();
473 Scalar pcJ = cellDataJ.capillaryPressure();
476 int isIndexI = intersection.indexInInside();
477 int isIndexJ = intersection.indexInOutside();
480 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
483 GlobalPosition distVec = globalPosJ - globalPosI;
486 Scalar dist = distVec.two_norm();
489 DimMatrix meanPermeability(0);
491 problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(elementI),
492 problem_.spatialParams().intrinsicPermeability(elementJ));
498 Scalar potentialDiffW = cellData.potential(wPhaseIdx) - cellDataJ.potential(wPhaseIdx);
499 Scalar potentialDiffNw = cellData.potential(nPhaseIdx) - cellDataJ.potential(nPhaseIdx);
501 if (compressibility_)
503 density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
504 density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
506 density_[wPhaseIdx] =
507 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
509 density_[nPhaseIdx] =
510 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
513 potentialDiffW = (cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx));
514 potentialDiffNw = (cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx));
517 potentialDiffW += density_[wPhaseIdx] * (distVec * gravity_);
518 potentialDiffNw += density_[nPhaseIdx] * (distVec * gravity_);
522 cellData.fluxData().setUpwindPotential(wPhaseIdx, isIndexI, potentialDiffW);
523 cellData.fluxData().setUpwindPotential(nPhaseIdx, isIndexI, potentialDiffNw);
525 cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, isIndexJ, -potentialDiffW);
526 cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, isIndexJ, -potentialDiffNw);
529 Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWJ;
530 lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30 )) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
531 Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwJ;
532 lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNw;
534 if (compressibility_)
536 density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
537 density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
539 density_[wPhaseIdx] =
540 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
542 density_[nPhaseIdx] =
543 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
550 Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
551 Dune::FieldVector<Scalar, dimWorld> velocityNw(unitOuterNormal);
555 Scalar areaScaling = (unitOuterNormal * distVec);
558 Scalar gravityTermW = (gravity_ * distVec) * density_[wPhaseIdx] * areaScaling;
559 Scalar gravityTermNw = (gravity_ * distVec) * density_[nPhaseIdx] * areaScaling;
562 switch (pressureType_)
566 velocityW *= lambdaW * scalarPerm
567 * ((cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + gravityTermW);
568 velocityNw *= lambdaNw * scalarPerm
569 * ((cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + gravityTermNw)
570 + 0.5 * (lambdaNwI + lambdaNwJ) * scalarPerm * (pcI - pcJ) / dist;
575 velocityW *= lambdaW * scalarPerm
576 * ((cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist + gravityTermW)
577 - 0.5 * (lambdaWI + lambdaWJ) * scalarPerm * (pcI - pcJ) / dist;
578 velocityNw *= lambdaNw * scalarPerm
579 * ((cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist + gravityTermNw);
584 velocityW *= (lambdaW + lambdaNw) * scalarPerm * (cellData.globalPressure() - cellDataJ.globalPressure()) / dist
585 + scalarPerm * (lambdaW * gravityTermW + lambdaNw * gravityTermNw);
592 cellData.fluxData().setVelocity(wPhaseIdx, isIndexI, velocityW);
593 cellData.fluxData().setVelocity(nPhaseIdx, isIndexI, velocityNw);
594 cellData.fluxData().setVelocityMarker(isIndexI);
597 Scalar weightingFactor = pow(0.5, (dim - 1)*(elementI.level() - elementJ.level()));
599 velocityW *= weightingFactor;
600 velocityNw *= weightingFactor;
602 cellDataJ.fluxData().addVelocity(wPhaseIdx, isIndexJ, velocityW);
603 cellDataJ.fluxData().addVelocity(nPhaseIdx, isIndexJ, velocityNw);
604 cellDataJ.fluxData().setVelocityMarker(isIndexJ);
#define GET_PROP_VALUE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:282
#define GET_PROP_TYPE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:283
Dune::IteratorRange< typename MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper >::Intersections::const_iterator > intersections(const MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper > &glue)
Range generator to iterate with range-based for loops over all intersections as follows: for (const a...
Definition: glue.hh:62
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Property tag EnableCompressibility
Returns whether compressibility is allowed.
Definition: porousmediumflow/2p/sequential/properties.hh:55
Property tag SaturationFormulation
The formulation of the saturation model.
Definition: porousmediumflow/2p/sequential/properties.hh:53
Property tag NumPhases
Number of phases in the system.
Definition: porousmediumflow/sequential/properties.hh:69
Property tag VelocityFormulation
The type of velocity reconstructed for the transport model.
Definition: porousmediumflow/2p/sequential/properties.hh:54
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
std::string permeability() noexcept
I/O name of permeability.
Definition: name.hh:143
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Determines the velocity from a finite volume solution of the pressure equation of a sequential model ...
Definition: 2p/sequential/diffusion/cellcentered/velocity.hh:61
void initialize()
For initialization.
Definition: 2p/sequential/diffusion/cellcentered/velocity.hh:143
Determines the velocity from a finite volume solution of the pressure equation of a sequential model ...
Definition: velocityadaptive.hh:40
void initialize()
For initialization.
Definition: velocityadaptive.hh:108
bool calculateVelocityInTransport()
Indicates if velocity is reconstructed in the pressure step or in the transport step.
Definition: velocityadaptive.hh:137
FVVelocity2PAdaptive(Problem &problem)
Constructs a FVVelocity2PAdaptive object.
Definition: velocityadaptive.hh:89
void calculateVelocity(const Intersection &intersection, CellData &cellData)
Calculates the velocity at a cell-cell interface.
Definition: velocityadaptive.hh:165
Velocity field from a finite volume solution of a pressure equation.