22#ifndef DUMUX_FVMPFAO2DPRESSURE2P_HH
23#define DUMUX_FVMPFAO2DPRESSURE2P_HH
65template<
class TypeTag>
73 dim = GridView::dimension, dimWorld = GridView::dimensionworld
89 using PrimaryVariables =
typename SolutionTypes::PrimaryVariables;
90 using ScalarSolutionType =
typename SolutionTypes::ScalarSolution;
96 pw = Indices::pressureW,
97 pn = Indices::pressureNw,
98 sw = Indices::saturationW,
99 sn = Indices::saturationNw
103 wPhaseIdx = Indices::wPhaseIdx,
104 nPhaseIdx = Indices::nPhaseIdx,
105 pressureIdx = Indices::pressureIdx,
106 saturationIdx = Indices::saturationIdx,
107 pressEqIdx = Indices::pressureEqIdx,
108 satEqIdx = Indices::satEqIdx,
109 numPhases = getPropValue<TypeTag, Properties::NumPhases>()
116 dirichletDirichlet = 1,
117 dirichletNeumann = 2,
121 using Element =
typename GridView::Traits::template Codim<0>::Entity;
122 using IntersectionIterator =
typename GridView::IntersectionIterator;
123 using Intersection =
typename GridView::Intersection;
125 using LocalPosition = Dune::FieldVector<Scalar, dim>;
126 using GlobalPosition =
typename Element::Geometry::GlobalCoordinate;
127 using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
128 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
130 using DimVector = Dune::FieldVector<Scalar, dim>;
134 using GlobalInteractionVolumeVector = std::vector<InteractionVolume>;
135 using InnerBoundaryVolumeFaces = std::vector<Dune::FieldVector<bool, 2*dim> >;
138 Intersection getNextIntersection_(
const Element&,
const IntersectionIterator&);
142 void initializeMatrix();
144 void storeInteractionVolumeInfo();
166 storeInteractionVolumeInfo();
178 const auto element = *problem_.gridView().template begin<0>();
179 FluidState fluidState;
180 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
181 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
182 fluidState.setTemperature(problem_.temperature(element));
183 fluidState.setSaturation(wPhaseIdx, 1.);
184 fluidState.setSaturation(nPhaseIdx, 0.);
195 storeInteractionVolumeInfo();
214 timeStep_ = problem_.timeManager().timeStepSize();
216 int size = problem_.gridView().size(0);
217 for (
int i = 0; i < size; i++)
221 switch (saturationType_)
224 sat = problem_.variables().cellData(i).saturation(wPhaseIdx);
227 sat = problem_.variables().cellData(i).saturation(nPhaseIdx);
232 maxError_ = max(maxError_, (sat - 1.0) / timeStep_);
236 maxError_ = max(maxError_, (-sat) / timeStep_);
254 for (
const auto& element : elements(problem_.gridView()))
267 int eIdxGlobal = problem_.variables().index(element);
268 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
270 switch (pressureType_)
274 Scalar potW = this->
pressure()[eIdxGlobal];
276 Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
277 Scalar potPc = cellData.capillaryPressure()
278 + gravityDiff * (density_[nPhaseIdx] - density_[wPhaseIdx]);
280 cellData.setPotential(wPhaseIdx, potW);
281 cellData.setPotential(nPhaseIdx, potW + potPc);
283 Scalar pressW = potW - gravityDiff * density_[wPhaseIdx];
285 cellData.setPressure(wPhaseIdx, pressW);
286 cellData.setPressure(nPhaseIdx, pressW + cellData.capillaryPressure());
292 Scalar potNw = this->
pressure()[eIdxGlobal];
294 Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
295 Scalar potPc = cellData.capillaryPressure()
296 + gravityDiff * (density_[nPhaseIdx] - density_[wPhaseIdx]);
298 cellData.setPotential(nPhaseIdx, potNw);
299 cellData.setPotential(wPhaseIdx, potNw - potPc);
301 Scalar pressNw = potNw - gravityDiff * density_[nPhaseIdx];
303 cellData.setPressure(wPhaseIdx, pressNw - cellData.capillaryPressure());
304 cellData.setPressure(nPhaseIdx, pressNw);
310 cellData.fluxData().resetVelocity();
323 template<
class MultiWriter>
326 int size = problem_.gridView().size(0);
327 ScalarSolutionType *potential = writer.allocateManagedBuffer(size);
331 if (pressureType_ == pw)
333 writer.attachCellData(*potential,
"wetting potential");
336 if (pressureType_ == pn)
338 writer.attachCellData(*potential,
"nonwetting potential");
341 if (vtkOutputLevel_ > 0)
343 ScalarSolutionType *
pressure = writer.allocateManagedBuffer(size);
344 ScalarSolutionType *pressureSecond = writer.allocateManagedBuffer(size);
345 ScalarSolutionType *potentialSecond = writer.allocateManagedBuffer(size);
346 ScalarSolutionType *pc = writer.allocateManagedBuffer(size);
348 for (
const auto& element : elements(problem_.gridView()))
350 int idx = problem_.variables().index(element);
351 CellData& cellData = problem_.variables().cellData(idx);
353 (*pc)[idx] = cellData.capillaryPressure();
355 if (pressureType_ == pw)
357 (*pressure)[idx] = cellData.pressure(wPhaseIdx);
358 (*potentialSecond)[idx] = cellData.potential(nPhaseIdx);
359 (*pressureSecond)[idx] = cellData.pressure(nPhaseIdx);
362 if (pressureType_ == pn)
364 (*pressure)[idx] = cellData.pressure(nPhaseIdx);
365 (*potentialSecond)[idx] = cellData.potential(wPhaseIdx);
366 (*pressureSecond)[idx] = cellData.pressure(wPhaseIdx);
370 if (pressureType_ == pw)
372 writer.attachCellData(*
pressure,
"wetting pressure");
373 writer.attachCellData(*pressureSecond,
"nonwetting pressure");
374 writer.attachCellData(*potentialSecond,
"nonwetting potential");
377 if (pressureType_ == pn)
379 writer.attachCellData(*
pressure,
"nonwetting pressure");
380 writer.attachCellData(*pressureSecond,
"wetting pressure");
381 writer.attachCellData(*potentialSecond,
"wetting potential");
384 writer.attachCellData(*pc,
"capillary pressure");
395 ParentType(problem), problem_(problem), gravity_(problem.gravity()), maxError_(
398 if (pressureType_ != pw && pressureType_ != pn)
400 DUNE_THROW(Dune::NotImplemented,
"Pressure type not supported!");
402 if (saturationType_ != sw && saturationType_ != sn)
404 DUNE_THROW(Dune::NotImplemented,
"Saturation type not supported!");
406 if (getPropValue<TypeTag, Properties::EnableCompressibility>())
408 DUNE_THROW(Dune::NotImplemented,
"Compressibility not supported!");
412 DUNE_THROW(Dune::NotImplemented,
"Dimension not supported!");
415 ErrorTermFactor_ = getParam<Scalar>(
"Impet.ErrorTermFactor");
416 ErrorTermLowerBound_ = getParam<Scalar>(
"Impet.ErrorTermLowerBound");
417 ErrorTermUpperBound_ = getParam<Scalar>(
"Impet.ErrorTermUpperBound");
419 density_[wPhaseIdx] = 0.;
420 density_[nPhaseIdx] = 0.;
421 viscosity_[wPhaseIdx] = 0.;
422 viscosity_[nPhaseIdx] = 0.;
424 vtkOutputLevel_ = getParam<int>(
"Vtk.OutputLevel");
446 Scalar evaluateErrorTerm_(CellData& cellData)
451 switch (saturationType_)
454 sat = cellData.saturation(wPhaseIdx);
457 sat = cellData.saturation(nPhaseIdx);
461 Scalar error = (sat > 1.0) ? sat - 1.0 : 0.0;
469 Scalar errorAbs = abs(error);
471 if ((errorAbs * timeStep_ > 1e-6) && (errorAbs > ErrorTermLowerBound_ * maxError_)
472 && (!problem_.timeManager().willBeFinished()))
474 return ErrorTermFactor_ * error;
484 const GravityVector& gravity_;
488 Scalar ErrorTermFactor_;
489 Scalar ErrorTermLowerBound_;
490 Scalar ErrorTermUpperBound_;
492 Scalar density_[numPhases];
493 Scalar viscosity_[numPhases];
497 static constexpr Scalar threshold_ = 1e-15;
499 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
501 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
503 static const int velocityType_ = getPropValue<TypeTag, Properties::VelocityFormulation>();
507template<
class TypeTag>
508typename FvMpfaO2dPressure2p<TypeTag>::Intersection
510 const IntersectionIterator& isIt)
512 auto isItBegin = problem_.gridView().ibegin(element);
513 const auto isEndIt = problem_.gridView().iend(element);
515 auto tempIsIt = isIt;
516 auto nextIsIt = ++tempIsIt;
519 switch (getPropValue<TypeTag, Properties::GridImplementation>())
522 case GridTypeIndices::yaspGrid:
524 if (nextIsIt == isEndIt)
526 nextIsIt = isItBegin;
530 nextIsIt = ++tempIsIt;
532 if (nextIsIt == isEndIt)
534 auto tempIsItBegin = isItBegin;
535 nextIsIt = ++tempIsItBegin;
542 case GridTypeIndices::aluGrid:
543 case GridTypeIndices::ugGrid:
545 if (nextIsIt == isEndIt)
546 nextIsIt = isItBegin;
552 DUNE_THROW(Dune::NotImplemented,
"GridType can not be used with MPFAO implementation!");
561template<
class TypeTag>
562void FvMpfaO2dPressure2p<TypeTag>::initializeMatrix()
565 for (
const auto& element : elements(problem_.gridView()))
568 int eIdxGlobalI = problem_.variables().index(element);
574 const auto isEndIt = problem_.gridView().iend(element);
575 for (
auto isIt = problem_.gridView().ibegin(element); isIt != isEndIt; ++isIt)
577 const auto& intersection = *isIt;
578 auto nextIntersection = getNextIntersection_(element, isIt);
580 if (intersection.neighbor())
584 if (intersection.neighbor() && nextIntersection.neighbor())
586 for (
const auto& innerIntersection
587 : intersections(problem_.gridView(), intersection.outside()))
588 for (
const auto& innerNextIntersection
589 : intersections(problem_.gridView(), nextIntersection.outside()))
591 if (innerIntersection.neighbor() && innerNextIntersection.neighbor())
593 if (innerIntersection.outside() == innerNextIntersection.outside()
594 && innerIntersection.outside() != intersection.inside())
604 this->A_.setrowsize(eIdxGlobalI, rowSize);
609 this->A_.endrowsizes();
612 for (
const auto& element : elements(problem_.gridView()))
615 int eIdxGlobalI = problem_.variables().index(element);
618 this->A_.addindex(eIdxGlobalI, eIdxGlobalI);
621 const auto isEndIt = problem_.gridView().iend(element);
622 for (
auto isIt = problem_.gridView().ibegin(element); isIt != isEndIt; ++isIt)
624 const auto& intersection = *isIt;
625 auto nextIntersection = getNextIntersection_(element, isIt);
627 if (intersection.neighbor())
630 int eIdxGlobalJ = problem_.variables().index(intersection.outside());
634 this->A_.addindex(eIdxGlobalI, eIdxGlobalJ);
637 if (intersection.neighbor() && nextIntersection.neighbor())
639 for (
const auto& innerIntersection
640 : intersections(problem_.gridView(), intersection.outside()))
641 for (
const auto& innerNextIntersection
642 : intersections(problem_.gridView(), nextIntersection.outside()))
644 if (innerIntersection.neighbor() && innerNextIntersection.neighbor())
646 auto innerOutside = innerIntersection.outside();
648 if (innerOutside == innerNextIntersection.outside()
649 && innerOutside != intersection.inside())
651 int eIdxGlobalJ = problem_.variables().index(innerOutside);
653 this->A_.addindex(eIdxGlobalI, eIdxGlobalJ);
662 this->A_.endindices();
690template<
class TypeTag>
691void FvMpfaO2dPressure2p<TypeTag>::storeInteractionVolumeInfo()
703 BoundaryTypes bcType;
706 for (
const auto& element : elements(problem_.gridView()))
710 int eIdxGlobal1 = problem_.variables().index(element);
712 const GlobalPosition& globalPos1 =
element.geometry().center();
715 DimMatrix K1(problem_.spatialParams().intrinsicPermeability(element));
717 const auto isIt12End = problem_.gridView().iend(element);
718 for (
auto isIt12 = problem_.gridView().ibegin(element); isIt12 != isIt12End; ++isIt12)
720 const auto& intersection12 = *isIt12;
721 auto intersection14 = getNextIntersection_(element, isIt12);
723 int indexInInside12 = intersection12.indexInInside();
724 int indexInInside14 = intersection14.indexInInside();
729 const auto refElement = referenceElement(element);
731 GlobalPosition corner1234(0);
733 int globalVertIdx1234 = 0;
736 for (
int i = 0; i < intersection12.geometry().corners(); ++i)
738 bool finished =
false;
740 const GlobalPosition& isIt12corner = intersection12.geometry().corner(i);
742 int localVertIdx12corner = refElement.subEntity(indexInInside12, dim - 1, i, dim);
744 int globalVertIdx12corner = problem_.variables().index(
element.template subEntity<dim>(localVertIdx12corner));
746 for (
int j = 0; j < intersection14.geometry().corners(); ++j)
748 int localVertIdx14corner = refElement.subEntity(indexInInside14, dim - 1, j, dim);
750 int globalVertIdx14corner = problem_.variables().index(
element.template subEntity<dim>(localVertIdx14corner));
752 if (globalVertIdx12corner == globalVertIdx14corner)
754 corner1234 = isIt12corner;
756 globalVertIdx1234 = globalVertIdx12corner;
769 if (interactionVolumes_[globalVertIdx1234].isStored())
775 interactionVolumes_[globalVertIdx1234].setStored();
780 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element, 0);
781 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection12.indexInInside(), 0, 0);
782 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection14.indexInInside(), 0, 1);
785 const GlobalPosition& globalPosFace12 = intersection12.geometry().center();
788 Scalar faceVol12 = intersection12.geometry().volume() / 2.0;
791 DimVector unitOuterNormal12 = intersection12.centerUnitOuterNormal();
794 const GlobalPosition& globalPosFace41 = intersection14.geometry().center();
797 Scalar faceVol41 = intersection14.geometry().volume() / 2.0;
800 DimVector unitOuterNormal14 = intersection14.centerUnitOuterNormal();
804 R.mv(globalPos1 - globalPosFace12, nu14);
807 R.mv(globalPosFace41 - globalPos1, nu12);
809 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu12, K1, 0, 0);
810 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu14, K1, 0, 1);
811 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal12, 0, 0);
812 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal14, 0, 1);
813 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol12, 0, 0);
814 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 0, 1);
820 interactionVolumes_[globalVertIdx1234].setDF(abs(nu14 * Rnu12), 0);
823 if (intersection12.neighbor())
826 auto element2 = intersection12.outside();
828 int eIdxGlobal2 = problem_.variables().index(element2);
831 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element2, 1);
832 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection12.indexInOutside(), 1, 1);
834 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal12, 1, 1);
835 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol12, 1, 1);
838 const GlobalPosition& globalPos2 = element2.geometry().center();
841 DimMatrix K2(problem_.spatialParams().intrinsicPermeability(element2));
844 if (intersection14.neighbor())
848 auto element4 = intersection14.outside();
851 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element4, 3);
852 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection14.indexInOutside(), 3, 0);
854 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal14, 3, 0);
855 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 3, 0);
860 const GlobalPosition& globalPos4 = element4.geometry().center();
863 DimMatrix K4(problem_.spatialParams().intrinsicPermeability(element4));
866 GlobalPosition globalPos3(0);
868 GlobalPosition globalPosFace23(0);
869 GlobalPosition globalPosFace34(0);
871 for (
const auto& intersection2
872 : intersections(problem_.gridView(), element2))
874 bool finished =
false;
876 for (
const auto& intersection4
877 : intersections(problem_.gridView(), element4))
879 if (intersection2.neighbor() && intersection4.neighbor())
881 auto element32 = intersection2.outside();
882 auto element34 = intersection4.outside();
885 if (element32 == element34 && element32 != element)
888 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element32, 2);
890 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection2.indexInInside(), 1,
892 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection2.indexInOutside(), 2,
894 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection4.indexInInside(), 3,
896 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection4.indexInOutside(), 2,
900 globalPos3 = element32.geometry().center();
902 globalPosFace23 = intersection2.geometry().center();
903 globalPosFace34 = intersection4.geometry().center();
905 Scalar faceVol23 = intersection2.geometry().volume() / 2.0;
906 Scalar faceVol34 = intersection4.geometry().volume() / 2.0;
909 DimVector unitOuterNormal23 = intersection2.centerUnitOuterNormal();
911 DimVector unitOuterNormal43 = intersection4.centerUnitOuterNormal();
913 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal23, 1, 0);
914 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal23, 2, 1);
915 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal43, 2, 0);
916 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal43, 3, 1);
917 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol23, 1, 0);
918 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol23, 2, 1);
919 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol34, 2, 0);
920 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol34, 3, 1);
924 problem_.spatialParams().intrinsicPermeability(element32));
928 R.umv(globalPosFace12 - globalPos2, nu23);
931 R.umv(globalPosFace23 - globalPos2, nu21);
934 R.umv(globalPosFace34 - globalPos3, nu32);
937 R.umv(globalPos3 - globalPosFace23, nu34);
940 R.umv(globalPos4 - globalPosFace34, nu41);
943 R.umv(globalPos4 - globalPosFace41, nu43);
945 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu23, K2, 1, 0);
946 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu21, K2, 1, 1);
947 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu34, K3, 2, 0);
948 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu32, K3, 2, 1);
949 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu41, K4, 3, 0);
950 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu43, K4, 3, 1);
955 interactionVolumes_[globalVertIdx1234].setDF(abs(nu23 * Rnu21), 1);
959 interactionVolumes_[globalVertIdx1234].setDF(abs(nu32 * Rnu34), 2);
963 interactionVolumes_[globalVertIdx1234].setDF(abs(nu41 * Rnu43), 3);
981 problem_.boundaryTypes(bcType, intersection14);
982 PrimaryVariables boundValues(0.0);
984 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 3);
985 if (bcType.isNeumann(pressEqIdx))
987 problem_.neumann(boundValues, intersection14);
988 boundValues *= faceVol41;
989 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 3);
991 if (bcType.hasDirichlet())
993 problem_.dirichlet(boundValues, intersection14);
994 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 3);
1001 GlobalPosition globalPosFace23(0);
1004 Scalar faceVol23 = 0;
1007 DimVector unitOuterNormal23(0);
1009 bool finished =
false;
1011 for (
const auto& intersection2
1012 : intersections(problem_.gridView(), element2))
1014 if (intersection2.boundary())
1016 for (
int i = 0; i < intersection2.geometry().corners(); ++i)
1018 int localVertIdx2corner = refElement.subEntity(intersection2.indexInInside(), dim - 1, i,
1021 int globalVertIdx2corner = problem_.variables().index(element2.template subEntity<dim>(localVertIdx2corner));
1023 if (globalVertIdx2corner == globalVertIdx1234)
1025 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection2.indexInInside(), 1,
1028 globalPosFace23 = intersection2.geometry().center();
1030 faceVol23 = intersection2.geometry().volume() / 2.0;
1032 unitOuterNormal23 = intersection2.centerUnitOuterNormal();
1034 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal23, 1, 0);
1035 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol23, 1, 0);
1037 problem_.boundaryTypes(bcType, intersection2);
1040 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 1);
1041 if (bcType.isNeumann(pressEqIdx))
1043 problem_.neumann(boundValues, intersection2);
1044 boundValues *= faceVol23;
1045 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 1);
1047 if (bcType.hasDirichlet())
1049 problem_.dirichlet(boundValues, intersection2);
1050 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 1);
1053 interactionVolumes_[globalVertIdx1234].setOutsideFace(2);
1055 innerBoundaryVolumeFaces_[eIdxGlobal1][intersection12.indexInInside()] =
true;
1056 innerBoundaryVolumeFaces_[eIdxGlobal2][intersection12.indexInOutside()] =
true;
1060 R.umv(globalPosFace12 - globalPos2, nu23);
1063 R.umv(globalPosFace23 - globalPos2, nu21);
1065 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu23, K2, 1, 0);
1066 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu21, K2, 1, 1);
1071 interactionVolumes_[globalVertIdx1234].setDF(abs(nu23 * Rnu21), 1);
1087 Dune::NotImplemented,
1088 "fvmpfao2pfaboundpressure2p.hh, l. 997: boundary shape not available as interaction volume shape");
1096 problem_.boundaryTypes(bcType, *isIt12);
1097 PrimaryVariables boundValues(0.0);
1099 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 0);
1100 if (bcType.isNeumann(pressEqIdx))
1102 problem_.neumann(boundValues, *isIt12);
1103 boundValues *= faceVol12;
1104 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 0);
1106 if (bcType.hasDirichlet())
1108 problem_.dirichlet(boundValues, *isIt12);
1109 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 0);
1113 if (intersection14.boundary())
1115 problem_.boundaryTypes(bcType, intersection14);
1118 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 3);
1119 if (bcType.isNeumann(pressEqIdx))
1121 problem_.neumann(boundValues, intersection14);
1122 boundValues *= faceVol41;
1123 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 3);
1125 if (bcType.hasDirichlet())
1127 problem_.dirichlet(boundValues, intersection14);
1128 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 3);
1131 interactionVolumes_[globalVertIdx1234].setOutsideFace(1);
1132 interactionVolumes_[globalVertIdx1234].setOutsideFace(2);
1140 auto element4 = intersection14.outside();
1141 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection14.indexInOutside(), 3, 0);
1144 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element4, 3);
1146 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal14, 3, 0);
1147 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 3, 0);
1150 const GlobalPosition& globalPos4 = element4.geometry().center();
1152 int eIdxGlobal4 = problem_.variables().index(element4);
1154 bool finished =
false;
1157 for (
const auto& intersection4
1158 : intersections(problem_.gridView(), element4))
1160 if (intersection4.boundary())
1162 for (
int i = 0; i < intersection4.geometry().corners(); ++i)
1164 int localVertIdx4corner = refElement.subEntity(intersection4.indexInInside(), dim - 1, i,
1167 int globalVertIdx4corner = problem_.variables().index(element4.template subEntity<dim>(localVertIdx4corner));
1169 if (globalVertIdx4corner == globalVertIdx1234)
1171 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection4.indexInInside(), 3,
1174 const GlobalPosition& globalPosFace34 = intersection4.geometry().center();
1176 Scalar faceVol34 = intersection4.geometry().volume() / 2.0;
1178 DimVector unitOuterNormal43 = intersection4.centerUnitOuterNormal();
1180 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal43, 3, 1);
1181 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol34, 3, 1);
1183 problem_.boundaryTypes(bcType, intersection4);
1186 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 2);
1187 if (bcType.isNeumann(pressEqIdx))
1189 problem_.neumann(boundValues, intersection4);
1190 boundValues *= faceVol34;
1191 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 2);
1193 if (bcType.hasDirichlet())
1195 problem_.dirichlet(boundValues, intersection4);
1196 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 2);
1199 interactionVolumes_[globalVertIdx1234].setOutsideFace(1);
1201 innerBoundaryVolumeFaces_[eIdxGlobal1][intersection14.indexInInside()] =
true;
1202 innerBoundaryVolumeFaces_[eIdxGlobal4][intersection14.indexInOutside()] =
true;
1206 problem_.spatialParams().intrinsicPermeability(element4));
1210 R.umv(globalPos4 - globalPosFace34, nu41);
1213 R.umv(globalPos4 - globalPosFace41, nu43);
1215 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu41, K4, 3, 0);
1216 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu43, K4, 3, 1);
1221 interactionVolumes_[globalVertIdx1234].setDF(abs(nu41 * Rnu43), 3);
1237 Dune::NotImplemented,
1238 "fvmpfao2pfaboundpressure2p.hh, l. 1164: boundary shape not available as interaction volume shape");
1251template<
class TypeTag>
1252void FvMpfaO2dPressure2p<TypeTag>::assemble()
1259 for (
const auto& vertex : vertices(problem_.gridView()))
1261 int vIdxGlobal = problem_.variables().index(vertex);
1263 InteractionVolume& interactionVolume = interactionVolumes_[vIdxGlobal];
1265 if (interactionVolume.isInnerVolume())
1268 auto element1 = interactionVolume.getSubVolumeElement(0);
1269 auto element2 = interactionVolume.getSubVolumeElement(1);
1270 auto element3 = interactionVolume.getSubVolumeElement(2);
1271 auto element4 = interactionVolume.getSubVolumeElement(3);
1274 const GlobalPosition& globalPos1 = element1.geometry().center();
1275 const GlobalPosition& globalPos2 = element2.geometry().center();
1276 const GlobalPosition& globalPos3 = element3.geometry().center();
1277 const GlobalPosition& globalPos4 = element4.geometry().center();
1280 Scalar volume1 = element1.geometry().volume();
1281 Scalar volume2 = element2.geometry().volume();
1282 Scalar volume3 = element3.geometry().volume();
1283 Scalar volume4 = element4.geometry().volume();
1286 int eIdxGlobal1 = problem_.variables().index(element1);
1287 int eIdxGlobal2 = problem_.variables().index(element2);
1288 int eIdxGlobal3 = problem_.variables().index(element3);
1289 int eIdxGlobal4 = problem_.variables().index(element4);
1292 CellData& cellData1 = problem_.variables().cellData(eIdxGlobal1);
1293 CellData& cellData2 = problem_.variables().cellData(eIdxGlobal2);
1294 CellData& cellData3 = problem_.variables().cellData(eIdxGlobal3);
1295 CellData& cellData4 = problem_.variables().cellData(eIdxGlobal4);
1298 PrimaryVariables source(0.0);
1299 problem_.source(source, element1);
1300 this->f_[eIdxGlobal1] += volume1 / (4.0)
1301 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1302 problem_.source(source, element2);
1303 this->f_[eIdxGlobal2] += volume2 / (4.0)
1304 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1305 problem_.source(source, element3);
1306 this->f_[eIdxGlobal3] += volume3 / (4.0)
1307 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1308 problem_.source(source, element4);
1309 this->f_[eIdxGlobal4] += volume4 / (4.0)
1310 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1312 this->f_[eIdxGlobal1] += evaluateErrorTerm_(cellData1) * volume1 / (4.0);
1313 this->f_[eIdxGlobal2] += evaluateErrorTerm_(cellData2) * volume2 / (4.0);
1314 this->f_[eIdxGlobal3] += evaluateErrorTerm_(cellData3) * volume3 / (4.0);
1315 this->f_[eIdxGlobal4] += evaluateErrorTerm_(cellData4) * volume4 / (4.0);
1318 Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
1319 lambda1[nPhaseIdx] = cellData1.mobility(nPhaseIdx);
1322 Scalar lambdaTotal1 = lambda1[wPhaseIdx] + lambda1[nPhaseIdx];
1325 Dune::FieldVector<Scalar, numPhases> lambda2(cellData2.mobility(wPhaseIdx));
1326 lambda2[nPhaseIdx] = cellData2.mobility(nPhaseIdx);
1329 Scalar lambdaTotal2 = lambda2[wPhaseIdx] + lambda2[nPhaseIdx];
1332 Dune::FieldVector<Scalar, numPhases> lambda3(cellData3.mobility(wPhaseIdx));
1333 lambda3[nPhaseIdx] = cellData3.mobility(nPhaseIdx);
1336 Scalar lambdaTotal3 = lambda3[wPhaseIdx] + lambda3[nPhaseIdx];
1339 Dune::FieldVector<Scalar, numPhases> lambda4(cellData4.mobility(wPhaseIdx));
1340 lambda4[nPhaseIdx] = cellData4.mobility(nPhaseIdx);
1343 Scalar lambdaTotal4 = lambda4[wPhaseIdx] + lambda4[nPhaseIdx];
1345 Scalar gn12nu14 = interactionVolume.getNtkrkNu_df(lambdaTotal1, 0, 0, 1);
1346 Scalar gn12nu12 = interactionVolume.getNtkrkNu_df(lambdaTotal1, 0, 0, 0);
1347 Scalar gn14nu14 = interactionVolume.getNtkrkNu_df(lambdaTotal1, 0, 1, 1);
1348 Scalar gn14nu12 = interactionVolume.getNtkrkNu_df(lambdaTotal1, 0, 1, 0);
1349 Scalar gn12nu23 = interactionVolume.getNtkrkNu_df(lambdaTotal2, 1, 1, 0);
1350 Scalar gn12nu21 = interactionVolume.getNtkrkNu_df(lambdaTotal2, 1, 1, 1);
1351 Scalar gn23nu23 = interactionVolume.getNtkrkNu_df(lambdaTotal2, 1, 0, 0);
1352 Scalar gn23nu21 = interactionVolume.getNtkrkNu_df(lambdaTotal2, 1, 0, 1);
1353 Scalar gn43nu32 = interactionVolume.getNtkrkNu_df(lambdaTotal3, 2, 0, 1);
1354 Scalar gn43nu34 = interactionVolume.getNtkrkNu_df(lambdaTotal3, 2, 0, 0);
1355 Scalar gn23nu32 = interactionVolume.getNtkrkNu_df(lambdaTotal3, 2, 1, 1);
1356 Scalar gn23nu34 = interactionVolume.getNtkrkNu_df(lambdaTotal3, 2, 1, 0);
1357 Scalar gn43nu41 = interactionVolume.getNtkrkNu_df(lambdaTotal4, 3, 1, 0);
1358 Scalar gn43nu43 = interactionVolume.getNtkrkNu_df(lambdaTotal4, 3, 1, 1);
1359 Scalar gn14nu41 = interactionVolume.getNtkrkNu_df(lambdaTotal4, 3, 0, 0);
1360 Scalar gn14nu43 = interactionVolume.getNtkrkNu_df(lambdaTotal4, 3, 0, 1);
1363 Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim> C(0), F(0), A(0), B(0);
1366 C[0][0] = -gn12nu12;
1367 C[0][3] = -gn12nu14;
1369 C[1][1] = -gn23nu23;
1372 C[3][2] = -gn14nu43;
1375 F[0][0] = gn12nu12 + gn12nu14;
1376 F[1][1] = -gn23nu21 + gn23nu23;
1377 F[2][2] = -gn43nu34 - gn43nu32;
1378 F[3][3] = gn14nu43 - gn14nu41;
1380 A[0][0] = gn12nu12 + gn12nu21;
1381 A[0][1] = -gn12nu23;
1383 A[1][0] = -gn23nu21;
1384 A[1][1] = gn23nu23 + gn23nu32;
1386 A[2][1] = -gn43nu32;
1387 A[2][2] = -gn43nu34 - gn43nu43;
1389 A[3][0] = -gn14nu12;
1391 A[3][3] = -gn14nu41 - gn14nu14;
1393 B[0][0] = gn12nu12 + gn12nu14;
1394 B[0][1] = gn12nu21 - gn12nu23;
1395 B[1][1] = -gn23nu21 + gn23nu23;
1396 B[1][2] = gn23nu34 + gn23nu32;
1397 B[2][2] = -gn43nu34 - gn43nu32;
1398 B[2][3] = -gn43nu43 + gn43nu41;
1399 B[3][0] = -gn14nu12 - gn14nu14;
1400 B[3][3] = gn14nu43 - gn14nu41;
1404 F += C.rightmultiply(B.leftmultiply(A));
1405 Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim> T(F);
1410 this->A_[eIdxGlobal1][eIdxGlobal1] += T[0][0] + T[3][0];
1411 this->A_[eIdxGlobal1][eIdxGlobal2] += T[0][1] + T[3][1];
1412 this->A_[eIdxGlobal1][eIdxGlobal3] += T[0][2] + T[3][2];
1413 this->A_[eIdxGlobal1][eIdxGlobal4] += T[0][3] + T[3][3];
1415 this->A_[eIdxGlobal2][eIdxGlobal1] += -T[0][0] + T[1][0];
1416 this->A_[eIdxGlobal2][eIdxGlobal2] += -T[0][1] + T[1][1];
1417 this->A_[eIdxGlobal2][eIdxGlobal3] += -T[0][2] + T[1][2];
1418 this->A_[eIdxGlobal2][eIdxGlobal4] += -T[0][3] + T[1][3];
1420 this->A_[eIdxGlobal3][eIdxGlobal1] -= T[1][0] + T[2][0];
1421 this->A_[eIdxGlobal3][eIdxGlobal2] -= T[1][1] + T[2][1];
1422 this->A_[eIdxGlobal3][eIdxGlobal3] -= T[1][2] + T[2][2];
1423 this->A_[eIdxGlobal3][eIdxGlobal4] -= T[1][3] + T[2][3];
1425 this->A_[eIdxGlobal4][eIdxGlobal1] += T[2][0] - T[3][0];
1426 this->A_[eIdxGlobal4][eIdxGlobal2] += T[2][1] - T[3][1];
1427 this->A_[eIdxGlobal4][eIdxGlobal3] += T[2][2] - T[3][2];
1428 this->A_[eIdxGlobal4][eIdxGlobal4] += T[2][3] - T[3][3];
1430 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 0)])
1432 this->A_[eIdxGlobal1][eIdxGlobal1] += T[0][0];
1433 this->A_[eIdxGlobal1][eIdxGlobal2] += T[0][1];
1434 this->A_[eIdxGlobal1][eIdxGlobal3] += T[0][2];
1435 this->A_[eIdxGlobal1][eIdxGlobal4] += T[0][3];
1437 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 1)])
1439 this->A_[eIdxGlobal1][eIdxGlobal1] += T[3][0];
1440 this->A_[eIdxGlobal1][eIdxGlobal2] += T[3][1];
1441 this->A_[eIdxGlobal1][eIdxGlobal3] += T[3][2];
1442 this->A_[eIdxGlobal1][eIdxGlobal4] += T[3][3];
1445 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 0)])
1447 this->A_[eIdxGlobal2][eIdxGlobal1] += T[1][0];
1448 this->A_[eIdxGlobal2][eIdxGlobal2] += T[1][1];
1449 this->A_[eIdxGlobal2][eIdxGlobal3] += T[1][2];
1450 this->A_[eIdxGlobal2][eIdxGlobal4] += T[1][3];
1452 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 1)])
1454 this->A_[eIdxGlobal2][eIdxGlobal1] += -T[0][0];
1455 this->A_[eIdxGlobal2][eIdxGlobal2] += -T[0][1];
1456 this->A_[eIdxGlobal2][eIdxGlobal3] += -T[0][2];
1457 this->A_[eIdxGlobal2][eIdxGlobal4] += -T[0][3];
1459 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 0)])
1461 this->A_[eIdxGlobal3][eIdxGlobal1] -= T[2][0];
1462 this->A_[eIdxGlobal3][eIdxGlobal2] -= T[2][1];
1463 this->A_[eIdxGlobal3][eIdxGlobal3] -= T[2][2];
1464 this->A_[eIdxGlobal3][eIdxGlobal4] -= T[2][3];
1466 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 1)])
1468 this->A_[eIdxGlobal3][eIdxGlobal1] -= T[1][0];
1469 this->A_[eIdxGlobal3][eIdxGlobal2] -= T[1][1];
1470 this->A_[eIdxGlobal3][eIdxGlobal3] -= T[1][2];
1471 this->A_[eIdxGlobal3][eIdxGlobal4] -= T[1][3];
1473 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 0)])
1475 this->A_[eIdxGlobal4][eIdxGlobal1] += -T[3][0];
1476 this->A_[eIdxGlobal4][eIdxGlobal2] += -T[3][1];
1477 this->A_[eIdxGlobal4][eIdxGlobal3] += -T[3][2];
1478 this->A_[eIdxGlobal4][eIdxGlobal4] += -T[3][3];
1480 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 1)])
1482 this->A_[eIdxGlobal4][eIdxGlobal1] += T[2][0];
1483 this->A_[eIdxGlobal4][eIdxGlobal2] += T[2][1];
1484 this->A_[eIdxGlobal4][eIdxGlobal3] += T[2][2];
1485 this->A_[eIdxGlobal4][eIdxGlobal4] += T[2][3];
1490 Dune::FieldVector<Scalar, 2 * dim> pc(0);
1491 pc[0] = cellData1.capillaryPressure();
1492 pc[1] = cellData2.capillaryPressure();
1493 pc[2] = cellData3.capillaryPressure();
1494 pc[3] = cellData4.capillaryPressure();
1496 Dune::FieldVector<Scalar, 2 * dim> gravityDiff(0);
1500 gravityDiff[0] = (problem_.bBoxMax() - globalPos1) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1501 gravityDiff[1] = (problem_.bBoxMax() - globalPos2) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1502 gravityDiff[2] = (problem_.bBoxMax() - globalPos3) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1503 gravityDiff[3] = (problem_.bBoxMax() - globalPos4) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1507 if (pc[0] == 0 && pc[1] == 0 && pc[2] == 0 && pc[3] == 0)
1512 Dune::FieldVector<Scalar, 2 * dim> pcFlux(0);
1518 Scalar pcPotential12 = pcFlux[0];
1519 Scalar pcPotential14 = pcFlux[3];
1520 Scalar pcPotential32 = -pcFlux[1];
1521 Scalar pcPotential34 = -pcFlux[2];
1526 Dune::FieldVector<Scalar, numPhases> lambda12Upw(0.0);
1527 lambda12Upw[wPhaseIdx] = (pcPotential12 >= 0) ? lambda1[wPhaseIdx] : lambda2[wPhaseIdx];
1528 lambda12Upw[nPhaseIdx] = (pcPotential12 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
1531 Dune::FieldVector<Scalar, numPhases> lambda14Upw(0.0);
1532 lambda14Upw[wPhaseIdx] = (pcPotential14 >= 0) ? lambda1[wPhaseIdx] : lambda4[wPhaseIdx];
1533 lambda14Upw[nPhaseIdx] = (pcPotential14 >= 0) ? lambda1[nPhaseIdx] : lambda4[nPhaseIdx];
1536 Dune::FieldVector<Scalar, numPhases> lambda32Upw(0.0);
1537 lambda32Upw[wPhaseIdx] = (pcPotential32 >= 0) ? lambda3[wPhaseIdx] : lambda2[wPhaseIdx];
1538 lambda32Upw[nPhaseIdx] = (pcPotential32 >= 0) ? lambda3[nPhaseIdx] : lambda2[nPhaseIdx];
1541 Dune::FieldVector<Scalar, numPhases> lambda34Upw(0.0);
1542 lambda34Upw[wPhaseIdx] = (pcPotential34 >= 0) ? lambda3[wPhaseIdx] : lambda4[wPhaseIdx];
1543 lambda34Upw[nPhaseIdx] = (pcPotential34 >= 0) ? lambda3[nPhaseIdx] : lambda4[nPhaseIdx];
1545 for (
int i = 0; i < numPhases; i++)
1547 Scalar lambdaT12 = lambda12Upw[wPhaseIdx] + lambda12Upw[nPhaseIdx];
1548 Scalar lambdaT14 = lambda14Upw[wPhaseIdx] + lambda14Upw[nPhaseIdx];
1549 Scalar lambdaT32 = lambda32Upw[wPhaseIdx] + lambda32Upw[nPhaseIdx];
1550 Scalar lambdaT34 = lambda34Upw[wPhaseIdx] + lambda34Upw[nPhaseIdx];
1551 Scalar fracFlow12 = (lambdaT12 > threshold_) ? lambda12Upw[i] / (lambdaT12) : 0.0;
1552 Scalar fracFlow14 = (lambdaT14 > threshold_) ? lambda14Upw[i] / (lambdaT14) : 0.0;
1553 Scalar fracFlow32 = (lambdaT32 > threshold_) ? lambda32Upw[i] / (lambdaT32) : 0.0;
1554 Scalar fracFlow34 = (lambdaT34 > threshold_) ? lambda34Upw[i] / (lambdaT34) : 0.0;
1556 Dune::FieldVector<Scalar, 2 * dim> pcFluxReal(pcFlux);
1558 pcFluxReal[0] *= fracFlow12;
1559 pcFluxReal[1] *= fracFlow32;
1560 pcFluxReal[2] *= fracFlow34;
1561 pcFluxReal[3] *= fracFlow14;
1567 switch (pressureType_)
1574 this->f_[eIdxGlobal1] -= (pcFluxReal[0] + pcFluxReal[3]);
1575 this->f_[eIdxGlobal2] -= (pcFluxReal[1] - pcFluxReal[0]);
1576 this->f_[eIdxGlobal3] -= (-pcFluxReal[2] - pcFluxReal[1]);
1577 this->f_[eIdxGlobal4] -= (-pcFluxReal[3] + pcFluxReal[2]);
1579 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 0)])
1581 this->f_[eIdxGlobal1] -= pcFluxReal[0];
1583 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 1)])
1585 this->f_[eIdxGlobal1] -= pcFluxReal[3];
1587 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 0)])
1589 this->f_[eIdxGlobal2] -= pcFluxReal[1];
1591 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 1)])
1593 this->f_[eIdxGlobal2] += pcFluxReal[0];
1595 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 0)])
1597 this->f_[eIdxGlobal3] += pcFluxReal[2];
1599 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 1)])
1601 this->f_[eIdxGlobal3] += pcFluxReal[1];
1603 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 0)])
1605 this->f_[eIdxGlobal4] += pcFluxReal[3];
1607 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 1)])
1609 this->f_[eIdxGlobal4] -= pcFluxReal[2];
1619 this->f_[eIdxGlobal1] += (pcFluxReal[0] + pcFluxReal[1]);
1620 this->f_[eIdxGlobal2] += (pcFluxReal[1] - pcFluxReal[0]);
1621 this->f_[eIdxGlobal3] += (-pcFluxReal[2] - pcFluxReal[1]);
1622 this->f_[eIdxGlobal4] += (-pcFluxReal[3] + pcFluxReal[2]);
1624 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 0)])
1626 this->f_[eIdxGlobal1] += pcFluxReal[0];
1628 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 1)])
1630 this->f_[eIdxGlobal1] += pcFluxReal[3];
1632 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 0)])
1634 this->f_[eIdxGlobal2] += pcFluxReal[1];
1636 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 1)])
1638 this->f_[eIdxGlobal2] -= pcFluxReal[0];
1640 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 0)])
1642 this->f_[eIdxGlobal3] -= pcFluxReal[2];
1644 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 1)])
1646 this->f_[eIdxGlobal3] -= pcFluxReal[1];
1648 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 0)])
1650 this->f_[eIdxGlobal4] -= pcFluxReal[3];
1652 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 1)])
1654 this->f_[eIdxGlobal4] += pcFluxReal[2];
1666 for (
int elemIdx = 0; elemIdx < 2 * dim; elemIdx++)
1668 bool isOutside =
false;
1669 for (
int fIdx = 0; fIdx < dim; fIdx++)
1671 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
1672 if (interactionVolume.isOutsideFace(intVolFaceIdx))
1683 auto element = interactionVolume.getSubVolumeElement(elemIdx);
1686 const GlobalPosition& globalPos =
element.geometry().center();
1692 int eIdxGlobal = problem_.variables().index(element);
1695 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
1698 DimMatrix
permeability(problem_.spatialParams().intrinsicPermeability(element));
1701 PrimaryVariables source(0);
1702 problem_.source(source, element);
1703 this->f_[eIdxGlobal] +=
volume / (4.0)
1704 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1706 this->f_[eIdxGlobal] += evaluateErrorTerm_(cellData) *
volume / (4.0);
1709 Dune::FieldVector<Scalar, numPhases> lambda(cellData.mobility(wPhaseIdx));
1710 lambda[nPhaseIdx] = cellData.mobility(nPhaseIdx);
1712 Scalar pc = cellData.capillaryPressure();
1714 Scalar gravityDiff = (problem_.bBoxMax() - globalPos) * gravity_
1715 * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1719 for (
int fIdx = 0; fIdx < dim; fIdx++)
1721 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
1723 if (interactionVolume.isBoundaryFace(intVolFaceIdx))
1726 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(pressEqIdx))
1728 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
1730 const auto refElement = referenceElement(element);
1732 const LocalPosition& localPos = refElement.position(boundaryFaceIdx, 1);
1734 const GlobalPosition& globalPosFace =
element.geometry().global(localPos);
1736 DimVector distVec(globalPosFace - globalPos);
1737 Scalar dist = distVec.two_norm();
1738 DimVector unitDistVec(distVec);
1739 unitDistVec /= dist;
1741 Scalar faceArea = interactionVolume.getFaceArea(elemIdx, fIdx);
1744 Scalar satWBound = cellData.saturation(wPhaseIdx);
1746 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(satEqIdx))
1748 Scalar satBound = interactionVolume.getDirichletValues(intVolFaceIdx)[saturationIdx];
1749 switch (saturationType_)
1753 satWBound = satBound;
1758 satWBound = 1 - satBound;
1765 const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(
element.geometry().center());
1766 Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
1768 Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
1769 * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1771 pcBound += gravityDiffBound;
1773 Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
1774 lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
1775 lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
1776 lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
1778 Scalar potentialBound = interactionVolume.getDirichletValues(intVolFaceIdx)[pressureIdx];
1779 Scalar gdeltaZ = (problem_.bBoxMax()-globalPosFace) * gravity_;
1782 Scalar potentialDiffW = 0;
1783 Scalar potentialDiffNw = 0;
1784 switch (pressureType_)
1788 potentialBound += density_[wPhaseIdx]*gdeltaZ;
1789 potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBound) / dist;
1790 potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBound - pcBound) / dist;
1795 potentialBound += density_[nPhaseIdx]*gdeltaZ;
1796 potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBound + pcBound) / dist;
1797 potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBound) / dist;
1802 Scalar lambdaTotal = (potentialDiffW >= 0.) ? lambda[wPhaseIdx] : lambdaBound[wPhaseIdx];
1803 lambdaTotal += (potentialDiffNw >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
1805 DimVector permTimesNormal(0);
1809 Scalar entry = lambdaTotal * (unitDistVec * permTimesNormal) / dist * faceArea;
1814 switch (pressureType_)
1819 DimVector pcGradient = unitDistVec;
1820 pcGradient *= (pc - pcBound) / dist;
1823 pcFlux = 0.5 * (lambda[nPhaseIdx] + lambdaBound[nPhaseIdx])
1824 * (permTimesNormal * pcGradient) * faceArea;
1831 DimVector pcGradient = unitDistVec;
1832 pcGradient *= (pc - pcBound) / dist;
1835 pcFlux = 0.5 * (lambda[wPhaseIdx] + lambdaBound[wPhaseIdx])
1836 * (permTimesNormal * pcGradient) * faceArea;
1844 this->A_[eIdxGlobal][eIdxGlobal] += entry;
1845 this->f_[eIdxGlobal] += entry * potentialBound;
1847 if (pc == 0 && pcBound == 0)
1852 for (
int i = 0; i < numPhases; i++)
1854 switch (pressureType_)
1861 this->f_[eIdxGlobal] -= pcFlux;
1870 this->f_[eIdxGlobal] += pcFlux;
1878 else if (interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressEqIdx))
1880 Scalar J = interactionVolume.getNeumannValues(intVolFaceIdx)[wPhaseIdx]
1881 / density_[wPhaseIdx];
1882 J += interactionVolume.getNeumannValues(intVolFaceIdx)[nPhaseIdx] / density_[nPhaseIdx];
1883 this->f_[eIdxGlobal] -= J;
1887 std::cout <<
"interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressEqIdx)"
1888 << interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressEqIdx) <<
"\n";
1889 DUNE_THROW(Dune::NotImplemented,
1890 "No valid boundary condition type defined for pressure equation!");
1901 if (problem_.gridView().comm().size() > 1)
1904 for (
const auto& element : elements(problem_.gridView()))
1906 if (
element.partitionType() == Dune::InteriorEntity)
1910 int eIdxGlobalI = problem_.variables().index(element);
1912 this->A_[eIdxGlobalI] = 0.0;
1913 this->A_[eIdxGlobalI][eIdxGlobalI] = 1.0;
1914 this->f_[eIdxGlobalI] = this->
pressure()[eIdxGlobalI];
1927template<
class TypeTag>
1931 for (
const auto& element : elements(problem_.gridView()))
1933 int eIdxGlobal = problem_.variables().index(element);
1935 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
1937 const Scalar satW = cellData.saturation(wPhaseIdx);
1939 const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(element.geometry().center());
1940 const Scalar pc = fluidMatrixInteraction.pc(satW);
1942 cellData.setCapillaryPressure(pc);
1945 const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
1946 const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
1949 cellData.setMobility(wPhaseIdx, mobilityW);
1950 cellData.setMobility(nPhaseIdx, mobilityNw);
1953 cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNw));
1954 cellData.setFracFlowFunc(nPhaseIdx, mobilityNw / (mobilityW + mobilityNw));
Class including the information of an interaction volume of a MPFA O-method that does not change with...
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 pressure(int phaseIdx) noexcept
I/O name of pressure for multiphase systems.
Definition: name.hh:34
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Scalar volume(Shape shape, Scalar inscribedRadius)
Returns the volume of a given geometry based on the inscribed radius.
Definition: poreproperties.hh:73
Finite volume MPFA O-method discretization of a two-phase flow pressure equation of the sequential IM...
Definition: omethod/2dpressure.hh:67
void updateMaterialLaws()
Constitutive functions are initialized and stored in the variables object.
Definition: omethod/2dpressure.hh:1928
GlobalInteractionVolumeVector interactionVolumes_
Definition: omethod/2dpressure.hh:480
void storePressureSolution(const Element &element)
Stores the pressure solution of a cell.
Definition: omethod/2dpressure.hh:265
InnerBoundaryVolumeFaces innerBoundaryVolumeFaces_
Definition: omethod/2dpressure.hh:481
void updateInteractionVolumeInfo()
Updates interaction volumes.
Definition: omethod/2dpressure.hh:158
void addOutputVtkFields(MultiWriter &writer)
Adds pressure output to the output file.
Definition: omethod/2dpressure.hh:324
void update()
Pressure update.
Definition: omethod/2dpressure.hh:210
void initialize()
Initializes the pressure model.
Definition: omethod/2dpressure.hh:174
void storePressureSolution()
Globally stores the pressure solution.
Definition: omethod/2dpressure.hh:251
FvMpfaO2dPressure2p(Problem &problem)
Constructs a FvMpfaO2dPressure2p object.
Definition: omethod/2dpressure.hh:394
Class including the information of an interaction volume of a MPFA O-method that does not change with...
Definition: ointeractionvolume.hh:38
The finite volume base class for the solution of a pressure equation.
Definition: sequential/cellcentered/pressure.hh:49
void initialize()
Initialize pressure model.
Definition: sequential/cellcentered/pressure.hh:213
PressureSolution & pressure()
Returns the vector containing the pressure solution.
Definition: sequential/cellcentered/pressure.hh:120
void solve()
Solves the global system of equations to get the spatial distribution of the pressure.
Definition: sequential/cellcentered/pressure.hh:527
Specifies the properties for immiscible 2p diffusion/pressure models.
Properties for a MPFA method.
Finite Volume Diffusion Model.