22#ifndef DUMUX_MPFAL2DPRESSUREVELOCITY2P_HH
23#define DUMUX_MPFAL2DPRESSUREVELOCITY2P_HH
25#include <dune/common/float_cmp.hh>
53 dim = GridView::dimension, dimWorld = GridView::dimensionworld
61 using PrimaryVariables =
typename SolutionTypes::PrimaryVariables;
70 using IndexSet =
typename GridView::IndexSet;
71 using Intersection =
typename GridView::Intersection;
77 pw = Indices::pressureW,
78 pn = Indices::pressureNw,
79 sw = Indices::saturationW,
80 sn = Indices::saturationNw,
81 wPhaseIdx = Indices::wPhaseIdx,
82 nPhaseIdx = Indices::nPhaseIdx,
83 pressureIdx = Indices::pressureIdx,
84 saturationIdx = Indices::saturationIdx,
85 pressEqIdx = Indices::pressureEqIdx,
86 satEqIdx = Indices::satEqIdx,
87 numPhases = getPropValue<TypeTag, Properties::NumPhases>()
90 using Element =
typename GridView::template Codim<0>::Entity;
92 using GlobalPosition =
typename Element::Geometry::GlobalCoordinate;
93 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
102 ParentType(problem), problem_(problem), velocity_(problem)
104 density_[wPhaseIdx] = 0.;
105 density_[nPhaseIdx] = 0.;
106 viscosity_[wPhaseIdx] = 0.;
107 viscosity_[nPhaseIdx] = 0.;
109 calcVelocityInTransport_ = getParam<bool>(
"MPFA.CalcVelocityInTransport");
137 const auto element = *problem_.gridView().template begin<0>();
138 FluidState fluidState;
139 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
140 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
141 fluidState.setTemperature(problem_.temperature(element));
142 fluidState.setSaturation(wPhaseIdx, 1.);
143 fluidState.setSaturation(nPhaseIdx, 0.);
150 velocity_.initialize();
176 return calcVelocityInTransport_;
189 template<
class MultiWriter>
193 velocity_.addOutputVtkFields(writer);
200 Scalar density_[numPhases];
201 Scalar viscosity_[numPhases];
202 bool calcVelocityInTransport_;
205 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
207 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
217template<
class TypeTag>
221 for (
const auto& vertex : vertices(problem_.gridView()))
223 int vIdxGlobal = problem_.variables().index(vertex);
235 int eIdxGlobal1 = problem_.variables().index(element1);
236 int eIdxGlobal2 = problem_.variables().index(element2);
237 int eIdxGlobal3 = problem_.variables().index(element3);
238 int eIdxGlobal4 = problem_.variables().index(element4);
241 CellData& cellData1 = problem_.variables().cellData(eIdxGlobal1);
242 CellData& cellData2 = problem_.variables().cellData(eIdxGlobal2);
243 CellData& cellData3 = problem_.variables().cellData(eIdxGlobal3);
244 CellData& cellData4 = problem_.variables().cellData(eIdxGlobal4);
246 velocity_.calculateInnerInteractionVolumeVelocity(interactionVolume, cellData1, cellData2,
247 cellData3, cellData4, this->innerBoundaryVolumeFaces_);
253 for (
int elemIdx = 0; elemIdx < 2 * dim; elemIdx++)
255 bool isOutside =
false;
256 for (
int fIdx = 0; fIdx < dim; fIdx++)
270 int eIdxGlobal = problem_.variables().index(interactionVolume.
getSubVolumeElement(elemIdx));
272 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
274 velocity_.calculateBoundaryInteractionVolumeVelocity(interactionVolume, cellData, elemIdx);
293template<
class TypeTag>
296 int numVertices = intersection.geometry().corners();
298 auto elementI = intersection.inside();
299 auto elementJ = intersection.outside();
301 int eIdxGlobalI = problem_.variables().index(elementI);
302 int eIdxGlobalJ = problem_.variables().index(elementJ);
304 CellData& cellDataJ = problem_.variables().cellData(eIdxGlobalJ);
306 const auto refElement = referenceElement(elementI);
308 int indexInInside = intersection.indexInInside();
309 int indexInOutside = intersection.indexInOutside();
311 Dune::FieldVector<CellData, 4> cellDataTemp;
313 for (
int vIdx = 0; vIdx < numVertices; vIdx++)
315 int localVertIdx = refElement.subEntity(indexInInside, dim - 1, vIdx, dim);
317 int vIdxGlobal = problem_.variables().index(elementI.template subEntity<dim>(localVertIdx));
331 eIdxGlobal[0] = problem_.variables().index(element1);
332 eIdxGlobal[1] = problem_.variables().index(element2);
333 eIdxGlobal[2] = problem_.variables().index(element3);
334 eIdxGlobal[3] = problem_.variables().index(element4);
337 cellDataTemp[0] = problem_.variables().cellData(eIdxGlobal[0]);
338 cellDataTemp[1] = problem_.variables().cellData(eIdxGlobal[1]);
339 cellDataTemp[2] = problem_.variables().cellData(eIdxGlobal[2]);
340 cellDataTemp[3] = problem_.variables().cellData(eIdxGlobal[3]);
342 velocity_.calculateInnerInteractionVolumeVelocity(interactionVolume, cellDataTemp[0], cellDataTemp[1],
343 cellDataTemp[2], cellDataTemp[3], this->innerBoundaryVolumeFaces_);
345 for (
int i = 0; i < 4; i++)
347 if (eIdxGlobal[i] == eIdxGlobalI)
349 cellData.fluxData().setVelocity(wPhaseIdx, indexInInside,
350 cellDataTemp[i].fluxData().velocity(wPhaseIdx, indexInInside));
351 cellData.fluxData().setVelocity(nPhaseIdx, indexInInside,
352 cellDataTemp[i].fluxData().velocity(nPhaseIdx, indexInInside));
353 cellData.fluxData().setUpwindPotential(wPhaseIdx, indexInInside,
354 cellDataTemp[i].fluxData().upwindPotential(wPhaseIdx, indexInInside));
355 cellData.fluxData().setUpwindPotential(nPhaseIdx, indexInInside,
356 cellDataTemp[i].fluxData().upwindPotential(nPhaseIdx, indexInInside));
358 else if (eIdxGlobal[i] == eIdxGlobalJ)
360 cellDataJ.fluxData().setVelocity(wPhaseIdx, indexInOutside,
361 cellDataTemp[i].fluxData().velocity(wPhaseIdx, indexInOutside));
362 cellDataJ.fluxData().setVelocity(nPhaseIdx, indexInOutside,
363 cellDataTemp[i].fluxData().velocity(nPhaseIdx, indexInOutside));
364 cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, indexInOutside,
365 cellDataTemp[i].fluxData().upwindPotential(wPhaseIdx, indexInOutside));
366 cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, indexInOutside,
367 cellDataTemp[i].fluxData().upwindPotential(nPhaseIdx, indexInOutside));
372 cellData.fluxData().setVelocityMarker(indexInInside);
373 cellDataJ.fluxData().setVelocityMarker(indexInOutside);
384template<
class TypeTag>
387 auto element = intersection.inside();
390 int isIndex = intersection.indexInInside();
393 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
395 BoundaryTypes bcType;
397 problem_.boundaryTypes(bcType, intersection);
398 PrimaryVariables boundValues(0.0);
400 if (bcType.isDirichlet(pressEqIdx))
402 problem_.dirichlet(boundValues, intersection);
405 const GlobalPosition& globalPosI = element.geometry().center();
408 const GlobalPosition& globalPosJ = intersection.geometry().center();
411 Scalar lambdaWI = cellData.mobility(wPhaseIdx);
412 Scalar lambdaNwI = cellData.mobility(nPhaseIdx);
415 Scalar pcI = cellData.capillaryPressure();
418 GlobalPosition distVec = globalPosJ - globalPosI;
421 Scalar dist = distVec.two_norm();
425 DimMatrix meanPermeability(0);
427 problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(element));
434 if (bcType.isDirichlet(satEqIdx))
436 switch (saturationType_)
440 satW = boundValues[saturationIdx];
445 satW = 1 - boundValues[saturationIdx];
452 satW = cellData.saturation(wPhaseIdx);
455 const Scalar pressBound = boundValues[pressureIdx];
457 const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(element.geometry().center());
458 const Scalar pcBound = fluidMatrixInteraction.pc(satW);
461 Scalar pressWBound = 0;
462 Scalar pressNwBound = 0;
463 if (pressureType_ == pw)
465 pressWBound = pressBound;
466 pressNwBound = pressBound + pcBound;
468 else if (pressureType_ == pn)
470 pressWBound = pressBound - pcBound;
471 pressNwBound = pressBound;
474 const Scalar lambdaWBound = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
475 const Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
477 Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
478 Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
481 potentialDiffW = (cellData.pressure(wPhaseIdx) - pressWBound);
482 potentialDiffNw = (cellData.pressure(nPhaseIdx) - pressNwBound);
484 potentialDiffW += density_[wPhaseIdx] * (distVec * problem_.gravity());
485 potentialDiffNw += density_[nPhaseIdx] * (distVec * problem_.gravity());
488 cellData.fluxData().setUpwindPotential(wPhaseIdx, isIndex, potentialDiffW);
489 cellData.fluxData().setUpwindPotential(nPhaseIdx, isIndex, potentialDiffNw);
492 Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWBound;
493 lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
494 Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwBound;
495 lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
501 Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
502 Dune::FieldVector<Scalar, dimWorld> velocityNw(unitOuterNormal);
506 Scalar areaScaling = (unitOuterNormal * distVec);
508 Scalar gravityTermW = (problem_.gravity() * distVec) * density_[wPhaseIdx] * areaScaling;
509 Scalar gravityTermNw = (problem_.gravity() * distVec) * density_[nPhaseIdx] * areaScaling;
512 switch (pressureType_)
516 velocityW *= lambdaW * scalarPerm * ((cellData.pressure(wPhaseIdx) - pressBound) / dist + gravityTermW);
517 velocityNw *= lambdaNw * scalarPerm * ((cellData.pressure(wPhaseIdx) - pressBound) / dist + gravityTermNw)
518 + 0.5 * (lambdaNwI + lambdaNwBound) * scalarPerm * (pcI - pcBound) / dist;
523 velocityW *= lambdaW * scalarPerm * ((cellData.pressure(nPhaseIdx) - pressBound) / dist + gravityTermW)
524 - 0.5 * (lambdaWI + lambdaWBound) * scalarPerm * (pcI - pcBound) / dist;
525 velocityNw *= lambdaNw * scalarPerm * ((cellData.pressure(nPhaseIdx) - pressBound) / dist + gravityTermNw);
531 cellData.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW);
532 cellData.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNw);
533 cellData.fluxData().setVelocityMarker(isIndex);
537 else if (bcType.isNeumann(pressEqIdx))
539 problem_.neumann(boundValues, intersection);
541 Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
542 Dune::FieldVector<Scalar, dimWorld> velocityNw(unitOuterNormal);
544 velocityW *= boundValues[wPhaseIdx];
545 velocityNw *= boundValues[nPhaseIdx];
547 velocityW /= density_[wPhaseIdx];
548 velocityNw /= density_[nPhaseIdx];
551 cellData.fluxData().setUpwindPotential(wPhaseIdx, isIndex, boundValues[wPhaseIdx]);
552 cellData.fluxData().setUpwindPotential(nPhaseIdx, isIndex, boundValues[nPhaseIdx]);
554 cellData.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW);
555 cellData.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNw);
556 cellData.fluxData().setVelocityMarker(isIndex);
560 DUNE_THROW(Dune::NotImplemented,
"No valid boundary condition type defined for pressure equation!");
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property
Definition: propertysystem.hh:141
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:150
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
Finite volume MPFA L-method discretization of a two-phase flow pressure equation of the sequential IM...
Definition: lmethod/2dpressure.hh:72
void update()
Pressure update.
Definition: lmethod/2dpressure.hh:288
void initialize()
Initializes the pressure model.
Definition: lmethod/2dpressure.hh:191
FVMPFALInteractionVolume< TypeTag > InteractionVolume
Type of the interaction volume objects.
Definition: lmethod/2dpressure.hh:148
void storePressureSolution()
Globally stores the pressure solution.
Definition: lmethod/2dpressure.hh:225
void addOutputVtkFields(MultiWriter &writer)
Adds pressure output to the output file.
Definition: lmethod/2dpressure.hh:337
void updateMaterialLaws()
constitutive functions are initialized and stored in the variables object
Definition: lmethod/2dpressure.hh:1828
Class for the calculation of velocities from the pressure solution of an IMPES scheme using a MPFA L-...
Definition: lmethod/2dpressurevelocity.hh:46
void updateVelocity()
Function for updating the velocity field if iterations are necessary in the transport solution.
Definition: lmethod/2dpressurevelocity.hh:120
void calculateVelocityOnBoundary(const Intersection &intersection, CellData &cellData)
Calculates the velocity at a boundary.
Definition: lmethod/2dpressurevelocity.hh:385
FvMpfaL2dPressureVelocity2p(Problem &problem)
Constructs a FvMpfaL2dPressureVelocity2p object.
Definition: lmethod/2dpressurevelocity.hh:101
void initialize()
Initializes pressure and velocity.
Definition: lmethod/2dpressurevelocity.hh:135
void addOutputVtkFields(MultiWriter &writer)
Adds velocity output to the output file.
Definition: lmethod/2dpressurevelocity.hh:190
bool calculateVelocityInTransport()
Indicates if velocity is reconstructed in the pressure step or in the transport step.
Definition: lmethod/2dpressurevelocity.hh:174
void calculateVelocity()
Calculates the velocities at all cell-cell interfaces.
Definition: lmethod/2dpressurevelocity.hh:218
void update()
Pressure and velocity update.
Definition: lmethod/2dpressurevelocity.hh:162
Class for calculating 2-d velocities from cell-wise constant pressure values.
Definition: lmethod/2dvelocity.hh:58
Class including the information of an interaction volume of a MPFA L-method that does not change with...
Definition: linteractionvolume.hh:41
bool isInnerVolume()
Returns true if the interaction volume is completely inside the model domain.
Definition: linteractionvolume.hh:316
Element getSubVolumeElement(int subVolumeIdx)
Get an element of the interaction volume.
Definition: linteractionvolume.hh:281
int getFaceIndexFromSubVolume(int subVolumeIdx, int subVolumeFaceIdx)
Map from local interaction volume numbering on element to numbering on interaction volume.
Definition: linteractionvolume.hh:270
bool isOutsideFace(int subVolumeFaceIdx)
Returns true if an interaction volume flux face is outside the model domain.
Definition: linteractionvolume.hh:301
The finite volume base class for the solution of a pressure equation.
Definition: sequential/cellcentered/pressure.hh:49
Finite volume MPFA O-method discretization of a two-phase pressure equation of the sequential IMPES m...
Velocity calculation using a 2-d MPFA O-method.