22#ifndef DUMUX_MPFAL2DPRESSUREVELOCITIES2P_ADAPTIVE_HH
23#define DUMUX_MPFAL2DPRESSUREVELOCITIES2P_ADAPTIVE_HH
25#include <dune/common/float_cmp.hh>
55 dim = GridView::dimension, dimWorld = GridView::dimensionworld
63 using PrimaryVariables =
typename SolutionTypes::PrimaryVariables;
72 using Grid =
typename GridView::Grid;
73 using IndexSet =
typename GridView::IndexSet;
74 using Intersection =
typename GridView::Intersection;
76 using Geometry =
typename Grid::template Codim<0>::Entity::Geometry;
77 using JacobianTransposed =
typename Geometry::JacobianTransposed ;
83 pw = Indices::pressureW,
84 pn = Indices::pressureNw,
85 sw = Indices::saturationW,
86 sn = Indices::saturationNw
90 wPhaseIdx = Indices::wPhaseIdx,
91 nPhaseIdx = Indices::nPhaseIdx,
92 pressureIdx = Indices::pressureIdx,
93 saturationIdx = Indices::saturationIdx,
94 pressEqIdx = Indices::pressureEqIdx,
95 satEqIdx = Indices::satEqIdx,
96 numPhases = getPropValue<TypeTag, Properties::NumPhases>()
99 using Element =
typename GridView::template Codim<0>::Entity;
101 using GlobalPosition =
typename Element::Geometry::GlobalCoordinate;
102 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
103 using DimVector = Dune::FieldVector<Scalar, dim>;
112 ParentType(problem), problem_(problem), velocity_(problem)
114 density_[wPhaseIdx] = 0.;
115 density_[nPhaseIdx] = 0.;
116 viscosity_[wPhaseIdx] = 0.;
117 viscosity_[nPhaseIdx] = 0.;
119 calcVelocityInTransport_ = getParam<bool>(
"MPFA.CalcVelocityInTransport");
147 const auto element = *problem_.gridView().template begin<0>();
148 FluidState fluidState;
149 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
150 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
151 fluidState.setTemperature(problem_.temperature(element));
152 fluidState.setSaturation(wPhaseIdx, 1.);
153 fluidState.setSaturation(nPhaseIdx, 0.);
160 velocity_.initialize();
186 return calcVelocityInTransport_;
199 template<
class MultiWriter>
203 velocity_.addOutputVtkFields(writer);
210 Scalar density_[numPhases];
211 Scalar viscosity_[numPhases];
212 bool calcVelocityInTransport_;
215 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
217 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
226template<
class TypeTag>
230 for (
const auto& vertex : vertices(problem_.gridView()))
232 int vIdxGlobal = problem_.variables().index(vertex);
246 int eIdxGlobal1 = problem_.variables().index(element1);
247 int eIdxGlobal2 = problem_.variables().index(element2);
248 int eIdxGlobal3 = problem_.variables().index(element3);
249 int eIdxGlobal4 = problem_.variables().index(element4);
252 CellData& cellData1 = problem_.variables().cellData(eIdxGlobal1);
253 CellData& cellData2 = problem_.variables().cellData(eIdxGlobal2);
254 CellData& cellData3 = problem_.variables().cellData(eIdxGlobal3);
255 CellData& cellData4 = problem_.variables().cellData(eIdxGlobal4);
257 velocity_.calculateInnerInteractionVolumeVelocity(interactionVolume, cellData1, cellData2, cellData3,
258 cellData4, this->innerBoundaryVolumeFaces_);
268 int eIdxGlobal1 = problem_.variables().index(element1);
269 int eIdxGlobal2 = problem_.variables().index(element2);
270 int eIdxGlobal4 = problem_.variables().index(element4);
273 CellData& cellData1 = problem_.variables().cellData(eIdxGlobal1);
274 CellData& cellData2 = problem_.variables().cellData(eIdxGlobal2);
275 CellData& cellData4 = problem_.variables().cellData(eIdxGlobal4);
277 velocity_.calculateHangingNodeInteractionVolumeVelocity(interactionVolume, cellData1, cellData2,
278 cellData4, this->innerBoundaryVolumeFaces_);
283 DUNE_THROW(Dune::NotImplemented,
"Unknown interactionvolume type!");
290 for (
int elemIdx = 0; elemIdx < 2 * dim; elemIdx++)
292 bool isOutside =
false;
293 for (
int fIdx = 0; fIdx < dim; fIdx++)
307 int eIdxGlobal = problem_.variables().index(interactionVolume.
getSubVolumeElement(elemIdx));
309 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
311 velocity_.calculateBoundaryInteractionVolumeVelocity(interactionVolume, cellData, elemIdx);
328template<
class TypeTag>
331 int numVertices = intersection.geometry().corners();
333 auto elementI = intersection.inside();
334 auto elementJ = intersection.outside();
336 int levelI = elementI.level();
337 int levelJ = elementJ.level();
339 int eIdxGlobalI = problem_.variables().index(elementI);
340 int eIdxGlobalJ = problem_.variables().index(elementJ);
342 CellData& cellDataJ = problem_.variables().cellData(eIdxGlobalJ);
344 const auto refElement = referenceElement(elementI);
346 int indexInInside = intersection.indexInInside();
347 int indexInOutside = intersection.indexInOutside();
349 int fIdx = indexInInside;
352 fIdx = indexInOutside;
354 std::vector<CellData> cellDataTemp(0);
356 if (levelI != levelJ)
359 cellData.fluxData().setVelocity(wPhaseIdx, indexInInside, vel);
360 cellData.fluxData().setVelocity(nPhaseIdx, indexInInside, vel);
361 cellData.fluxData().setUpwindPotential(wPhaseIdx, indexInInside, 0);
362 cellData.fluxData().setUpwindPotential(nPhaseIdx, indexInInside, 0);
364 cellDataJ.fluxData().setVelocity(wPhaseIdx, indexInOutside, vel);
365 cellDataJ.fluxData().setVelocity(nPhaseIdx, indexInOutside, vel);
366 cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, indexInOutside, 0);
367 cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, indexInOutside, 0);
370 for (
int vIdx = 0; vIdx < numVertices; vIdx++)
372 int localVertIdx = refElement.subEntity(fIdx, dim - 1, vIdx, dim);
375 if (levelI >= levelJ)
377 vIdxGlobal = problem_.variables().index(elementI.template subEntity<dim>(localVertIdx));
381 vIdxGlobal = problem_.variables().index(elementJ.template subEntity<dim>(localVertIdx));
389 std::vector<int> eIdxGlobal(0);
398 eIdxGlobal.resize(4);
400 eIdxGlobal[0] = problem_.variables().index(element1);
401 eIdxGlobal[1] = problem_.variables().index(element2);
402 eIdxGlobal[2] = problem_.variables().index(element3);
403 eIdxGlobal[3] = problem_.variables().index(element4);
406 cellDataTemp.resize(4);
408 cellDataTemp[0] = problem_.variables().cellData(eIdxGlobal[0]);
409 cellDataTemp[1] = problem_.variables().cellData(eIdxGlobal[1]);
410 cellDataTemp[2] = problem_.variables().cellData(eIdxGlobal[2]);
411 cellDataTemp[3] = problem_.variables().cellData(eIdxGlobal[3]);
413 velocity_.calculateInnerInteractionVolumeVelocity(interactionVolume, cellDataTemp[0], cellDataTemp[1],
414 cellDataTemp[2], cellDataTemp[3], this->innerBoundaryVolumeFaces_);
422 eIdxGlobal.resize(3);
424 eIdxGlobal[0] = problem_.variables().index(element1);
425 eIdxGlobal[1] = problem_.variables().index(element2);
426 eIdxGlobal[2] = problem_.variables().index(element4);
429 cellDataTemp.resize(3);
431 cellDataTemp[0] = problem_.variables().cellData(eIdxGlobal[0]);
432 cellDataTemp[1] = problem_.variables().cellData(eIdxGlobal[1]);
433 cellDataTemp[2] = problem_.variables().cellData(eIdxGlobal[2]);
436 velocity_.calculateHangingNodeInteractionVolumeVelocity(interactionVolume, cellDataTemp[0], cellDataTemp[1],
437 cellDataTemp[2], this->innerBoundaryVolumeFaces_);
441 DUNE_THROW(Dune::NotImplemented,
"Unknown interactionvolume type!");
444 int size = cellDataTemp.size();
445 for (
int i = 0; i < size; i++)
447 if (eIdxGlobal[i] == eIdxGlobalI)
449 if (levelI >= levelJ)
451 cellData.fluxData().setVelocity(wPhaseIdx, indexInInside,
452 cellDataTemp[i].fluxData().velocity(wPhaseIdx, indexInInside));
453 cellData.fluxData().setVelocity(nPhaseIdx, indexInInside,
454 cellDataTemp[i].fluxData().velocity(nPhaseIdx, indexInInside));
455 cellData.fluxData().setUpwindPotential(wPhaseIdx, indexInInside,
456 cellDataTemp[i].fluxData().upwindPotential(wPhaseIdx, indexInInside));
457 cellData.fluxData().setUpwindPotential(nPhaseIdx, indexInInside,
458 cellDataTemp[i].fluxData().upwindPotential(nPhaseIdx, indexInInside));
462 cellDataJ.fluxData().setVelocity(wPhaseIdx, indexInOutside,
463 cellDataTemp[i].fluxData().velocity(wPhaseIdx, indexInInside));
464 cellDataJ.fluxData().setVelocity(nPhaseIdx, indexInOutside,
465 cellDataTemp[i].fluxData().velocity(nPhaseIdx, indexInInside));
466 cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, indexInOutside,
467 cellDataTemp[i].fluxData().upwindPotential(wPhaseIdx, indexInInside));
468 cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, indexInOutside,
469 cellDataTemp[i].fluxData().upwindPotential(nPhaseIdx, indexInInside));
475 else if (eIdxGlobal[i] == eIdxGlobalJ)
477 if (levelJ >= levelI)
479 cellDataJ.fluxData().setVelocity(wPhaseIdx, indexInOutside,
480 cellDataTemp[i].fluxData().velocity(wPhaseIdx, indexInOutside));
481 cellDataJ.fluxData().setVelocity(nPhaseIdx, indexInOutside,
482 cellDataTemp[i].fluxData().velocity(nPhaseIdx, indexInOutside));
483 cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, indexInOutside,
484 cellDataTemp[i].fluxData().upwindPotential(wPhaseIdx, indexInOutside));
485 cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, indexInOutside,
486 cellDataTemp[i].fluxData().upwindPotential(nPhaseIdx, indexInOutside));
490 cellData.fluxData().setVelocity(wPhaseIdx, indexInInside,
491 cellDataTemp[i].fluxData().velocity(wPhaseIdx, indexInOutside));
492 cellData.fluxData().setVelocity(nPhaseIdx, indexInInside,
493 cellDataTemp[i].fluxData().velocity(nPhaseIdx, indexInOutside));
494 cellData.fluxData().setUpwindPotential(wPhaseIdx, indexInInside,
495 cellDataTemp[i].fluxData().upwindPotential(wPhaseIdx, indexInOutside));
496 cellData.fluxData().setUpwindPotential(nPhaseIdx, indexInInside,
497 cellDataTemp[i].fluxData().upwindPotential(nPhaseIdx, indexInOutside));
505 if (levelI == levelJ)
507 cellData.fluxData().setVelocityMarker(indexInInside);
508 cellDataJ.fluxData().setVelocityMarker(indexInOutside);
520template<
class TypeTag>
523 auto element = intersection.inside();
526 int isIndex = intersection.indexInInside();
529 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
531 BoundaryTypes bcType;
533 problem_.boundaryTypes(bcType, intersection);
534 PrimaryVariables boundValues(0.0);
536 if (bcType.isDirichlet(pressEqIdx))
538 problem_.dirichlet(boundValues, intersection);
541 const GlobalPosition& globalPosI = element.geometry().center();
544 const GlobalPosition& globalPosJ = intersection.geometry().center();
547 Scalar lambdaWI = cellData.mobility(wPhaseIdx);
548 Scalar lambdaNwI = cellData.mobility(nPhaseIdx);
551 Scalar pcI = cellData.capillaryPressure();
554 GlobalPosition distVec = globalPosJ - globalPosI;
557 Scalar dist = distVec.two_norm();
561 DimMatrix meanPermeability(0);
563 problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(element));
570 if (bcType.isDirichlet(satEqIdx))
572 switch (saturationType_)
576 satW = boundValues[saturationIdx];
581 satW = 1 - boundValues[saturationIdx];
588 satW = cellData.saturation(wPhaseIdx);
591 const Scalar pressBound = boundValues[pressureIdx];
593 const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(element.geometry().center());
594 const Scalar pcBound = fluidMatrixInteraction.pc(satW);
597 Scalar pressWBound = 0;
598 Scalar pressNwBound = 0;
599 if (pressureType_ == pw)
601 pressWBound = pressBound;
602 pressNwBound = pressBound + pcBound;
604 else if (pressureType_ == pn)
606 pressWBound = pressBound - pcBound;
607 pressNwBound = pressBound;
610 const Scalar lambdaWBound = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
611 const Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
613 Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
614 Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
617 potentialDiffW = (cellData.pressure(wPhaseIdx) - pressWBound);
618 potentialDiffNw = (cellData.pressure(nPhaseIdx) - pressNwBound);
620 potentialDiffW += density_[wPhaseIdx] * (distVec * problem_.gravity());
621 potentialDiffNw += density_[nPhaseIdx] * (distVec * problem_.gravity());
624 cellData.fluxData().setUpwindPotential(wPhaseIdx, isIndex, potentialDiffW);
625 cellData.fluxData().setUpwindPotential(nPhaseIdx, isIndex, potentialDiffNw);
628 Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWBound;
629 lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
630 Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwBound;
631 lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
637 Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
638 Dune::FieldVector<Scalar, dimWorld> velocityNw(unitOuterNormal);
642 Scalar areaScaling = (unitOuterNormal * distVec);
645 Scalar gravityTermW = (problem_.gravity() * distVec) * density_[wPhaseIdx] * areaScaling;
646 Scalar gravityTermNw = (problem_.gravity() * distVec) * density_[nPhaseIdx] * areaScaling;
649 switch (pressureType_)
653 velocityW *= lambdaW * scalarPerm * ((cellData.pressure(wPhaseIdx) - pressBound) / dist + gravityTermW);
654 velocityNw *= lambdaNw * scalarPerm * ((cellData.pressure(wPhaseIdx) - pressBound) / dist + gravityTermNw)
655 + 0.5 * (lambdaNwI + lambdaNwBound) * scalarPerm * (pcI - pcBound) / dist;
660 velocityW *= lambdaW * scalarPerm * ((cellData.pressure(nPhaseIdx) - pressBound) / dist + gravityTermW)
661 - 0.5 * (lambdaWI + lambdaWBound) * scalarPerm * (pcI - pcBound) / dist;
662 velocityNw *= lambdaNw * scalarPerm * ((cellData.pressure(nPhaseIdx) - pressBound) / dist + gravityTermNw);
668 cellData.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW);
669 cellData.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNw);
670 cellData.fluxData().setVelocityMarker(isIndex);
674 else if (bcType.isNeumann(pressEqIdx))
676 problem_.neumann(boundValues, intersection);
678 Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
679 Dune::FieldVector<Scalar, dimWorld> velocityNw(unitOuterNormal);
681 velocityW *= boundValues[wPhaseIdx];
682 velocityNw *= boundValues[nPhaseIdx];
684 velocityW /= density_[wPhaseIdx];
685 velocityNw /= density_[nPhaseIdx];
688 cellData.fluxData().setUpwindPotential(wPhaseIdx, isIndex, boundValues[wPhaseIdx]);
689 cellData.fluxData().setUpwindPotential(nPhaseIdx, isIndex, boundValues[nPhaseIdx]);
691 cellData.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW);
692 cellData.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNw);
693 cellData.fluxData().setVelocityMarker(isIndex);
697 DUNE_THROW(Dune::NotImplemented,
"No valid boundary condition type defined for pressure equation!");
Velocity calculation on adaptive grids using a 2-d MPFA L-method.
Grid adaptive finite volume MPFA L-method discretization of a two-phase pressure equation of the sequ...
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
Grid adaptive finite volume MPFA L-method discretization of a two-phase flow pressure equation of the...
Definition: 2dpressureadaptive.hh:75
FVMPFALInteractionVolume< TypeTag > InteractionVolume
Type of the interaction volume objects.
Definition: 2dpressureadaptive.hh:152
void addOutputVtkFields(MultiWriter &writer)
Adds pressure output to the output file.
Definition: 2dpressureadaptive.hh:374
void storePressureSolution()
Globally stores the pressure solution.
Definition: 2dpressureadaptive.hh:231
void initialize()
Initializes the pressure model.
Definition: 2dpressureadaptive.hh:200
void updateMaterialLaws()
constitutive functions are initialized and stored in the variables object
Definition: 2dpressureadaptive.hh:2546
void update()
Pressure update.
Definition: 2dpressureadaptive.hh:298
Class for the calculation of velocities from the pressure solution of an IMPES scheme using a grid ad...
Definition: 2dpressurevelocityadaptive.hh:48
void addOutputVtkFields(MultiWriter &writer)
Adds velocity output to the output file.
Definition: 2dpressurevelocityadaptive.hh:200
void updateVelocity()
Function for updating the velocity field if iterations are necessary in the transport solution.
Definition: 2dpressurevelocityadaptive.hh:130
bool calculateVelocityInTransport()
Indicates if velocity is reconstructed in the pressure step or in the transport step.
Definition: 2dpressurevelocityadaptive.hh:184
FvMpfaL2dPressureVelocity2pAdaptive(Problem &problem)
Constructs a FvMpfaL2dPressureVelocity2pAdaptive object.
Definition: 2dpressurevelocityadaptive.hh:111
void initialize()
Initializes pressure and velocity.
Definition: 2dpressurevelocityadaptive.hh:145
void calculateVelocityOnBoundary(const Intersection &intersection, CellData &cellData)
Calculates the velocity at a boundary.
Definition: 2dpressurevelocityadaptive.hh:521
void calculateVelocity()
Calculates the velocities at all cell-cell interfaces.
Definition: 2dpressurevelocityadaptive.hh:227
void update()
Pressure and velocity update.
Definition: 2dpressurevelocityadaptive.hh:172
Class for calculating 2-d velocities from cell-wise constant pressure values.
Definition: 2dvelocityadaptive.hh:54
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
int getElementNumber()
Get number of stored elements.
Definition: linteractionvolume.hh:246
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