24#ifndef DUMUX_FVSATURATION2P_HH
25#define DUMUX_FVSATURATION2P_HH
74template<
class TypeTag>
84 dim = GridView::dimension, dimWorld = GridView::dimensionworld
104 using PrimaryVariables =
typename SolutionTypes::PrimaryVariables;
110 pw = Indices::pressureW,
111 pn = Indices::pressureNw,
112 pGlobal = Indices::pressureGlobal,
113 vw = Indices::velocityW,
114 vn = Indices::velocityNw,
115 vt = Indices::velocityTotal,
116 sw = Indices::saturationW,
117 sn = Indices::saturationNw
121 wPhaseIdx = Indices::wPhaseIdx,
122 nPhaseIdx = Indices::nPhaseIdx,
123 saturationIdx = Indices::saturationIdx,
124 satEqIdx = Indices::satEqIdx,
125 numPhases = getPropValue<TypeTag, Properties::NumPhases>()
130 using Element =
typename GridView::Traits::template Codim<0>::Entity;
131 using Intersection =
typename GridView::Intersection;
133 using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
134 using DimVector = Dune::FieldVector<Scalar, dim>;
139 return *capillaryFlux_;
144 return *capillaryFlux_;
149 return *gravityFlux_;
154 return *gravityFlux_;
169 void getFlux(Scalar&
update,
const Intersection& intersection, CellData& cellDataI);
175 void getSource(Scalar&
update,
const Element& element, CellData& cellDataI);
192 int size = problem_.gridView().size(0);
193 transportedQuantity.resize(size);
194 for (
int i = 0; i < size; i++)
196 switch (saturationType_)
200 transportedQuantity[i] = problem_.variables().cellData(i).saturation(wPhaseIdx);
205 transportedQuantity[i] = problem_.variables().cellData(i).saturation(nPhaseIdx);
220 int size = problem_.gridView().size(0);
221 for (
int i = 0; i < size; i++)
223 CellData& cellData = problem_.variables().cellData(i);
224 switch (saturationType_)
228 Scalar sat = transportedQuantity[i];
230 cellData.setSaturation(wPhaseIdx, sat);
231 cellData.setSaturation(nPhaseIdx, 1 - sat);
236 Scalar sat = transportedQuantity[i];
238 cellData.setSaturation(wPhaseIdx,1 -sat);
239 cellData.setSaturation(nPhaseIdx, sat);
252 template<
class DataEntry>
255 return (entry > -1e-6 && entry < 1.0 + 1e-6);
268 asImp_().updateSaturationSolution(updateVec);
279 asImp_().updateSaturationSolution(updateVec, dt);
289 Scalar dt = problem_.timeManager().timeStepSize();
290 int size = problem_.gridView().size(0);
291 for (
int i = 0; i < size; i++)
293 asImp_().updateSaturationSolution(i, updateVec[i][0], dt);
305 int size = problem_.gridView().size(0);
306 for (
int i = 0; i < size; i++)
308 asImp_().updateSaturationSolution(i, updateVec[i][0], dt);
323 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
325 switch (saturationType_)
329 Scalar sat = cellData.saturation(wPhaseIdx) + dt*
update;
331 cellData.setSaturation(wPhaseIdx, sat);
332 cellData.setSaturation(nPhaseIdx, 1 - sat);
337 Scalar sat = cellData.saturation(nPhaseIdx) + dt*
update;
339 cellData.setSaturation(wPhaseIdx,1 -sat);
340 cellData.setSaturation(nPhaseIdx, sat);
358 template<
class MultiWriter>
361 int size = problem_.gridView().size(0);
362 TransportSolutionType *
saturation = writer.allocateManagedBuffer(size);
363 TransportSolutionType *saturationSecond = 0;
365 if (vtkOutputLevel_ > 0)
367 saturationSecond = writer.allocateManagedBuffer(size);
370 for (
int i = 0; i < size; i++)
372 CellData& cellData = problem_.variables().cellData(i);
374 if (saturationType_ == sw)
376 (*saturation)[i] = cellData.saturation(wPhaseIdx);
377 if (vtkOutputLevel_ > 0)
379 (*saturationSecond)[i] = cellData.saturation(nPhaseIdx);
382 else if (saturationType_ == sn)
384 (*saturation)[i] = cellData.saturation(nPhaseIdx);
385 if (vtkOutputLevel_ > 0)
387 (*saturationSecond)[i] = cellData.saturation(wPhaseIdx);
392 if (saturationType_ == sw)
394 writer.attachCellData(*
saturation,
"wetting saturation");
395 if (vtkOutputLevel_ > 0)
397 writer.attachCellData(*saturationSecond,
"nonwetting saturation");
400 else if (saturationType_ == sn)
402 writer.attachCellData(*
saturation,
"nonwetting saturation");
403 if (vtkOutputLevel_ > 0)
405 writer.attachCellData(*saturationSecond,
"wetting saturation");
409 if (
velocity().calculateVelocityInTransport())
410 velocity().addOutputVtkFields(writer);
422 int eIdxGlobal = problem_.variables().index(element);
425 switch (saturationType_)
428 sat = problem_.variables().cellData(eIdxGlobal).saturation(wPhaseIdx);
431 sat = problem_.variables().cellData(eIdxGlobal).saturation(nPhaseIdx);
445 int eIdxGlobal = problem_.variables().index(element);
450 switch (saturationType_)
453 problem_.variables().cellData(eIdxGlobal).setSaturation(wPhaseIdx, sat);
454 problem_.variables().cellData(eIdxGlobal).setSaturation(nPhaseIdx, 1-sat);
457 problem_.variables().cellData(eIdxGlobal).setSaturation(nPhaseIdx, sat);
458 problem_.variables().cellData(eIdxGlobal).setSaturation(wPhaseIdx, 1-sat);
470 ParentType(problem), problem_(problem), threshold_(1e-6),
471 switchNormals_(
getParam<bool>(
"Impet.SwitchNormals"))
473 if (compressibility_ && velocityType_ == vt)
475 DUNE_THROW(Dune::NotImplemented,
476 "Total velocity - global pressure - model cannot be used with compressible fluids!");
478 if (saturationType_ != sw && saturationType_ != sn)
480 DUNE_THROW(Dune::NotImplemented,
"Saturation type not supported!");
482 if (pressureType_ != pw && pressureType_ != pn && pressureType_ != pGlobal)
484 DUNE_THROW(Dune::NotImplemented,
"Pressure type not supported!");
486 if (velocityType_ != vw && velocityType_ != vn && velocityType_ != vt)
488 DUNE_THROW(Dune::NotImplemented,
"Velocity type not supported!");
491 capillaryFlux_ = std::make_shared<CapillaryFlux>(problem);
492 gravityFlux_ = std::make_shared<GravityFlux>(problem);
493 velocity_ = std::make_shared<Velocity>(problem);
495 vtkOutputLevel_ = getParam<int>(
"Vtk.OutputLevel");
496 porosityThreshold_ = getParam<Scalar>(
"Impet.PorosityThreshold");
501 Implementation &asImp_()
502 {
return *
static_cast<Implementation *
>(
this); }
505 const Implementation &asImp_()
const
506 {
return *
static_cast<const Implementation *
>(
this); }
509 std::shared_ptr<Velocity> velocity_;
510 std::shared_ptr<CapillaryFlux> capillaryFlux_;
511 std::shared_ptr<GravityFlux> gravityFlux_;
514 Scalar porosityThreshold_;
516 static const bool compressibility_ = getPropValue<TypeTag, Properties::EnableCompressibility>();
517 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
518 static const int velocityType_ = getPropValue<TypeTag, Properties::VelocityFormulation>();
519 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
521 Scalar density_[numPhases];
522 Scalar viscosity_[numPhases];
524 const Scalar threshold_;
525 const bool switchNormals_;
539template<
class TypeTag>
542 auto elementI = intersection.inside();
543 auto elementJ = intersection.outside();
545 const CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(elementJ));
548 const GlobalPosition& globalPosI = elementI.geometry().center();
549 const GlobalPosition& globalPosJ = elementJ.geometry().center();
552 Scalar
volume = elementI.geometry().volume();
554 Scalar
porosity = max(problem_.spatialParams().porosity(elementI), porosityThreshold_);
556 if (compressibility_)
558 viscosity_[wPhaseIdx] = cellDataI.viscosity(wPhaseIdx);
559 viscosity_[nPhaseIdx] = cellDataI.viscosity(nPhaseIdx);
563 int isIndex = intersection.indexInInside();
565 GlobalPosition unitOuterNormal = intersection.centerUnitOuterNormal();
567 unitOuterNormal *= -1.0;
569 Scalar faceArea = intersection.geometry().volume();
571 if (velocity().calculateVelocityInTransport() && !cellDataI.fluxData().haveVelocity(isIndex))
572 velocity().calculateVelocity(intersection, cellDataI);
575 Scalar factorW = (cellDataI.fluxData().velocity(wPhaseIdx, isIndex) * unitOuterNormal) * faceArea;
576 Scalar factorNw = (cellDataI.fluxData().velocity(nPhaseIdx, isIndex) * unitOuterNormal) * faceArea;
577 Scalar factorTotal = factorW + factorNw;
580 GlobalPosition distVec = globalPosJ - globalPosI;
582 Scalar dist = distVec.two_norm();
584 bool takeNeighbor = (elementI.level() < elementJ.level());
587 (takeNeighbor) ? !cellDataJ.fluxData().isUpwindCell(wPhaseIdx, intersection.indexInOutside()) :
588 cellDataI.fluxData().isUpwindCell(wPhaseIdx, isIndex);
590 (takeNeighbor) ? !cellDataJ.fluxData().isUpwindCell(nPhaseIdx, intersection.indexInOutside()) :
591 cellDataI.fluxData().isUpwindCell(nPhaseIdx, isIndex);
599 lambdaW = cellDataI.mobility(wPhaseIdx);
600 if (compressibility_)
602 lambdaW /= cellDataI.density(wPhaseIdx);
607 lambdaW = cellDataJ.mobility(wPhaseIdx);
608 if (compressibility_)
610 lambdaW /= cellDataJ.density(wPhaseIdx);
616 lambdaNw = cellDataI.mobility(nPhaseIdx);
617 if (compressibility_)
619 lambdaNw /= cellDataI.density(nPhaseIdx);
624 lambdaNw = cellDataJ.mobility(nPhaseIdx);
625 if (compressibility_)
627 lambdaNw /= cellDataJ.density(nPhaseIdx);
631 switch (velocityType_)
636 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
637 viscosity_[nPhaseIdx], factorTotal, intersection);
640 Scalar satWI = cellDataI.saturation(wPhaseIdx);
641 Scalar satWJ = cellDataJ.saturation(wPhaseIdx);
643 Scalar pcI = cellDataI.capillaryPressure();
644 Scalar pcJ = cellDataJ.capillaryPressure();
647 GlobalPosition pcGradient = unitOuterNormal;
648 pcGradient *= (pcI - pcJ) / dist;
652 capillaryFlux().getFlux(flux, intersection, satWI, satWJ, pcGradient);
653 Scalar capillaryFlux = (flux * unitOuterNormal * faceArea);
656 gravityFlux().getFlux(flux, intersection, satWI, satWJ);
657 Scalar gravityFlux = (flux * unitOuterNormal * faceArea);
659 switch (saturationType_)
664 factorTotal *= lambdaW / (lambdaW + lambdaNw);
670 factorTotal *= lambdaNw / (lambdaW + lambdaNw);
676 factorTotal -= capillaryFlux;
677 factorTotal += gravityFlux;
680 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
681 viscosity_[nPhaseIdx], 10 * capillaryFlux, intersection);
682 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
683 viscosity_[nPhaseIdx], 10 * gravityFlux, intersection);
689 if (compressibility_)
691 factorW /= cellDataI.density(wPhaseIdx);
692 factorNw /= cellDataI.density(nPhaseIdx);
696 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
697 viscosity_[nPhaseIdx], factorW, intersection, wPhaseIdx);
698 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
699 viscosity_[nPhaseIdx], factorNw, intersection, nPhaseIdx);
705 switch (velocityType_)
711 switch (saturationType_)
732template<
class TypeTag>
735 auto elementI = intersection.inside();
738 const GlobalPosition& globalPosI = elementI.geometry().center();
741 const GlobalPosition& globalPosJ = intersection.geometry().center();
744 Scalar
volume = elementI.geometry().volume();
746 Scalar
porosity = max(problem_.spatialParams().porosity(elementI), porosityThreshold_);
748 if (compressibility_)
750 viscosity_[wPhaseIdx] = cellDataI.viscosity(wPhaseIdx);
751 viscosity_[nPhaseIdx] = cellDataI.viscosity(nPhaseIdx);
755 int isIndex = intersection.indexInInside();
757 GlobalPosition unitOuterNormal = intersection.centerUnitOuterNormal();
759 unitOuterNormal *= -1.0;
761 Scalar faceArea = intersection.geometry().volume();
763 if (velocity().calculateVelocityInTransport())
764 velocity().calculateVelocityOnBoundary(intersection, cellDataI);
767 Scalar factorW = (cellDataI.fluxData().velocity(wPhaseIdx, isIndex) * unitOuterNormal) * faceArea;
768 Scalar factorNw = (cellDataI.fluxData().velocity(nPhaseIdx, isIndex) * unitOuterNormal) * faceArea;
769 Scalar factorTotal = factorW + factorNw;
772 GlobalPosition distVec = globalPosJ - globalPosI;
775 Scalar dist = distVec.two_norm();
778 BoundaryTypes bcType;
779 problem_.boundaryTypes(bcType, intersection);
780 PrimaryVariables boundValues(0.0);
782 if (bcType.isDirichlet(satEqIdx))
784 problem_.dirichlet(boundValues, intersection);
786 Scalar satBound = boundValues[saturationIdx];
789 Scalar satWI = cellDataI.saturation(wPhaseIdx);
790 Scalar satWBound = 0;
791 switch (saturationType_)
795 satWBound = satBound;
800 satWBound = 1 - satBound;
805 const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(elementI.geometry().center());
806 const Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
812 if (cellDataI.fluxData().isUpwindCell(wPhaseIdx, isIndex))
814 lambdaW = cellDataI.mobility(wPhaseIdx);
815 if (compressibility_)
817 lambdaW /= cellDataI.density(wPhaseIdx);
822 if (compressibility_)
824 lambdaW = fluidMatrixInteraction.krw(satWBound)
829 lambdaW = fluidMatrixInteraction.krw(satWBound)
830 / viscosity_[wPhaseIdx];
834 if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, isIndex))
836 lambdaNw = cellDataI.mobility(nPhaseIdx);
837 if (compressibility_)
839 lambdaNw /= cellDataI.density(nPhaseIdx);
844 if (compressibility_)
846 lambdaNw = fluidMatrixInteraction.krn(satWBound)
851 lambdaNw = fluidMatrixInteraction.krn(satWBound)
852 / viscosity_[nPhaseIdx];
856 switch (velocityType_)
860 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
861 viscosity_[nPhaseIdx], factorTotal, intersection);
863 Scalar pcI = cellDataI.capillaryPressure();
866 GlobalPosition pcGradient = unitOuterNormal;
867 pcGradient *= (pcI - pcBound) / dist;
871 capillaryFlux().getFlux(flux, intersection, satWI, satWBound, pcGradient);
872 Scalar capillaryFlux = flux * unitOuterNormal * faceArea;
875 gravityFlux().getFlux(flux, intersection, satWI, satWBound);
876 Scalar gravityFlux = flux * unitOuterNormal * faceArea;
878 switch (saturationType_)
883 factorTotal *= lambdaW / (lambdaW + lambdaNw);
889 factorTotal *= lambdaNw / (lambdaW + lambdaNw);
891 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
892 viscosity_[nPhaseIdx], factorW, intersection, wPhaseIdx);
893 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
894 viscosity_[nPhaseIdx], factorNw, intersection, nPhaseIdx);
900 factorTotal -= capillaryFlux;
901 factorTotal += gravityFlux;
904 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
905 viscosity_[nPhaseIdx], 10 * capillaryFlux, intersection);
906 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
907 viscosity_[nPhaseIdx], 10 * gravityFlux, intersection);
913 if (compressibility_)
915 factorW /= cellDataI.density(wPhaseIdx);
916 factorNw /= cellDataI.density(nPhaseIdx);
920 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
921 viscosity_[nPhaseIdx], factorW, intersection, wPhaseIdx);
922 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
923 viscosity_[nPhaseIdx], factorNw, intersection, nPhaseIdx);
930 if (bcType.isNeumann(satEqIdx))
932 problem_.neumann(boundValues, intersection);
933 factorW = boundValues[wPhaseIdx];
934 factorNw = boundValues[nPhaseIdx];
936 factorNw *= faceArea;
939 Scalar lambdaW, lambdaNw;
941 lambdaW = cellDataI.mobility(wPhaseIdx);
942 lambdaNw = cellDataI.mobility(nPhaseIdx);
943 if (compressibility_)
945 lambdaW /= cellDataI.density(wPhaseIdx);
946 lambdaNw /= cellDataI.density(nPhaseIdx);
947 factorW /= cellDataI.density(wPhaseIdx);
948 factorNw /= cellDataI.density(nPhaseIdx);
952 factorW /= density_[wPhaseIdx];
953 factorNw /= density_[nPhaseIdx];
956 switch (velocityType_)
961 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
962 viscosity_[nPhaseIdx], factorW + factorNw, intersection);
967 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
968 viscosity_[nPhaseIdx], factorW, intersection, wPhaseIdx);
969 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
970 viscosity_[nPhaseIdx], factorNw, intersection, nPhaseIdx);
977 if (bcType.isOutflow(satEqIdx))
980 Scalar lambdaW = cellDataI.mobility(wPhaseIdx);
981 Scalar lambdaNw = cellDataI.mobility(nPhaseIdx);
982 if (compressibility_)
984 lambdaW /= cellDataI.density(wPhaseIdx);
985 lambdaNw /= cellDataI.density(nPhaseIdx);
988 if (velocityType_ == vt)
990 switch (saturationType_)
995 factorTotal *= lambdaW / (lambdaW + lambdaNw);
996 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
997 viscosity_[nPhaseIdx], factorTotal, intersection);
1003 factorTotal *= lambdaNw / (lambdaW + lambdaNw);
1004 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
1005 viscosity_[nPhaseIdx], factorTotal, intersection);
1012 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
1013 viscosity_[nPhaseIdx], factorW, intersection, wPhaseIdx);
1014 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
1015 viscosity_[nPhaseIdx], factorNw, intersection, nPhaseIdx);
1018 switch (velocityType_)
1024 switch (saturationType_)
1044template<
class TypeTag>
1048 Scalar
volume = element.geometry().volume();
1051 Scalar
porosity = max(problem_.spatialParams().porosity(element), porosityThreshold_);
1053 if (compressibility_)
1055 viscosity_[wPhaseIdx] = cellDataI.viscosity(wPhaseIdx);
1056 viscosity_[nPhaseIdx] = cellDataI.viscosity(nPhaseIdx);
1059 PrimaryVariables sourceVec(0.0);
1060 problem_.source(sourceVec, element);
1062 if (compressibility_)
1064 sourceVec[wPhaseIdx] /= cellDataI.density(wPhaseIdx);
1065 sourceVec[nPhaseIdx] /= cellDataI.density(nPhaseIdx);
1069 sourceVec[wPhaseIdx] /= density_[wPhaseIdx];
1070 sourceVec[nPhaseIdx] /= density_[nPhaseIdx];
1074 Scalar lambdaW = cellDataI.mobility(wPhaseIdx);
1075 Scalar lambdaNw = cellDataI.mobility(nPhaseIdx);
1076 if (compressibility_)
1078 lambdaW /= cellDataI.density(wPhaseIdx);
1079 lambdaNw /= cellDataI.density(nPhaseIdx);
1082 switch (saturationType_)
1086 if (sourceVec[wPhaseIdx] < 0 && cellDataI.saturation(wPhaseIdx) < threshold_)
1087 sourceVec[wPhaseIdx] = 0.0;
1089 update += sourceVec[wPhaseIdx] /
porosity;
1094 if (sourceVec[nPhaseIdx] < 0 && cellDataI.saturation(nPhaseIdx) < threshold_)
1095 sourceVec[nPhaseIdx] = 0.0;
1097 update += sourceVec[nPhaseIdx] /
porosity;
1102 switch (velocityType_)
1107 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx], viscosity_[nPhaseIdx],
1108 (sourceVec[wPhaseIdx] + sourceVec[nPhaseIdx]) * -1 *
volume, element);
1114 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
1115 viscosity_[nPhaseIdx], sourceVec[wPhaseIdx] * -1 *
volume, element,
1117 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
1118 viscosity_[nPhaseIdx], sourceVec[nPhaseIdx] * -1 *
volume, element,
1126template<
class TypeTag>
1129 ParentType::initialize();
1131 if (!compressibility_)
1133 const auto element = *problem_.gridView().template begin<0>();
1134 FluidState fluidState;
1135 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
1136 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
1137 fluidState.setTemperature(problem_.temperature(element));
1138 fluidState.setSaturation(wPhaseIdx, 1.);
1139 fluidState.setSaturation(nPhaseIdx, 0.);
1147 for (
const auto& element : elements(problem_.gridView()))
1149 PrimaryVariables initSol(0.0);
1150 problem_.initial(initSol, element);
1152 int eIdxGlobal = problem_.variables().index(element);
1154 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
1156 switch (saturationType_)
1160 cellData.setSaturation(wPhaseIdx, initSol[saturationIdx]);
1161 cellData.setSaturation(nPhaseIdx, 1 - initSol[saturationIdx]);
1166 cellData.setSaturation(wPhaseIdx,1 -initSol[saturationIdx]);
1167 cellData.setSaturation(nPhaseIdx, initSol[saturationIdx]);
1174 velocity_->initialize();
1175 capillaryFlux_->initialize();
1176 gravityFlux_->initialize();
1186template<
class TypeTag>
1190 for (
const auto& element : elements(problem_.gridView()))
1192 int eIdxGlobal = problem_.variables().index(element);
1194 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
1197 Scalar satW = cellData.saturation(wPhaseIdx);
1199 const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(element.geometry().center());
1200 const Scalar pc = fluidMatrixInteraction.pc(satW);
1202 cellData.setCapillaryPressure(pc);
1205 const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
1206 const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
1209 cellData.setMobility(wPhaseIdx, mobilityW);
1210 cellData.setMobility(nPhaseIdx, mobilityNw);
1213 cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNw));
1214 cellData.setFracFlowFunc(nPhaseIdx, mobilityNw / (mobilityW + mobilityNw));
Finite Volume discretization of a transport equation.
T getParam(Args &&... args)
A free function to get a parameter from the parameter tree singleton.
Definition: parameters.hh:348
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 saturation(int phaseIdx) noexcept
I/O name of saturation for multiphase systems.
Definition: name.hh:43
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
std::string porosity() noexcept
I/O name of porosity.
Definition: name.hh:139
Scalar volume(Shape shape, Scalar inscribedRadius)
Returns the volume of a given geometry based on the inscribed radius.
Definition: poreproperties.hh:73
The finite volume discretization of a saturation transport equation.
Definition: saturation.hh:76
const CapillaryFlux & capillaryFlux() const
Definition: saturation.hh:142
void getFlux(Scalar &update, const Intersection &intersection, CellData &cellDataI)
Function which calculates the flux update.
Definition: saturation.hh:540
Velocity & velocity()
Definition: saturation.hh:158
GravityFlux & gravityFlux()
Definition: saturation.hh:147
void updateTransportedQuantity(TransportSolutionType &updateVec)
Updates the primary transport variable.
Definition: saturation.hh:263
void updateSaturationSolution(int eIdxGlobal, Scalar update, Scalar dt)
Updates the saturation solution of a cell.
Definition: saturation.hh:321
void getFluxOnBoundary(Scalar &update, const Intersection &intersection, CellData &cellDataI)
Function which calculates the boundary flux update.
Definition: saturation.hh:733
void updateSaturationSolution(TransportSolutionType &updateVec)
Globally updates the saturation solution.
Definition: saturation.hh:287
CapillaryFlux & capillaryFlux()
Definition: saturation.hh:137
Velocity & velocity() const
Definition: saturation.hh:163
void updateSaturationSolution(TransportSolutionType &updateVec, Scalar dt)
Globally updates the saturation solution.
Definition: saturation.hh:303
void initialize()
Sets the initial solution.
Definition: saturation.hh:1127
void getTransportedQuantity(TransportSolutionType &transportedQuantity)
Writes the current values of the primary transport variable into the transportedQuantity-vector (come...
Definition: saturation.hh:190
const GravityFlux & gravityFlux() const
Definition: saturation.hh:152
void serializeEntity(std::ostream &outstream, const Element &element)
Function for serialization of the primary transport variable.
Definition: saturation.hh:420
void deserializeEntity(std::istream &instream, const Element &element)
Function for deserialization of the primary transport variable.
Definition: saturation.hh:443
void updateMaterialLaws()
Update the values of the material laws and constitutive relations.
Definition: saturation.hh:1187
void addOutputVtkFields(MultiWriter &writer)
Adds saturation output to the output file.
Definition: saturation.hh:359
FVSaturation2P(Problem &problem)
Constructs a FVSaturation2P object.
Definition: saturation.hh:469
void getSource(Scalar &update, const Element &element, CellData &cellDataI)
Function which calculates the source update.
Definition: saturation.hh:1045
void setTransportedQuantity(TransportSolutionType &transportedQuantity)
Writes the current values of the primary transport variable into the variable container.
Definition: saturation.hh:218
void updateTransportedQuantity(TransportSolutionType &updateVec, Scalar dt)
Updates the primary transport variable.
Definition: saturation.hh:277
bool inPhysicalRange(DataEntry &entry)
Check if saturation is in physical range.
Definition: saturation.hh:253
The finite volume discretization of a transport equation.
Definition: transport.hh:52
bool enableLocalTimeStepping()
Definition: transport.hh:202
void update(const Scalar t, Scalar &dt, TransportSolutionType &updateVec, bool impet=false)
Calculate the update vector.
Definition: transport.hh:270
void innerUpdate(TransportSolutionType &updateVec)
Definition: transport.hh:564
Specifies the properties for immiscible 2p transport.