24#ifndef DUMUX_FVTRANSPORT2P2C_HH
25#define DUMUX_FVTRANSPORT2P2C_HH
28#include <unordered_map>
30#include <dune/grid/common/gridenums.hh>
31#include <dune/common/float_cmp.hh>
58template<
class TypeTag>
80 dim = GridView::dimension, dimWorld = GridView::dimensionworld
84 pw = Indices::pressureW,
85 pn = Indices::pressureN
89 wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
90 wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx,
91 contiWEqIdx=Indices::contiWEqIdx, contiNEqIdx=Indices::contiNEqIdx,
92 NumPhases = getPropValue<TypeTag, Properties::NumPhases>(),
93 NumComponents = getPropValue<TypeTag, Properties::NumComponents>()
96 using Element =
typename GridView::Traits::template Codim<0>::Entity;
97 using Intersection =
typename GridView::Intersection;
99 using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
100 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
101 using PhaseVector = Dune::FieldVector<Scalar, NumPhases>;
102 using ComponentVector = Dune::FieldVector<Scalar, NumComponents>;
124 void innerUpdate(TransportSolutionType& updateVec);
127 virtual void update(
const Scalar t, Scalar& dt, TransportSolutionType& updateVec,
bool impet =
false);
137 const Intersection& intersection, CellData& cellDataI);
141 const Intersection& intersection,
const CellData& cellDataI);
144 const Intersection& intersection,
145 FluidState& BCfluidState,
146 PhaseVector& pressBound);
156 int transportedQuantities = getPropValue<TypeTag, Properties::NumEq>() - 1;
158 for (
int eqNumber = 0; eqNumber < transportedQuantities; eqNumber++)
169 template<
class MultiWriter>
172 if(
problem().vtkOutputLevel()>3)
175 int size =
problem_.gridView().size(0);
176 ScalarSolutionType *totalC1PV = writer.allocateManagedBuffer(size);
177 ScalarSolutionType *totalC2PV = writer.allocateManagedBuffer(size);
180 writer.attachCellData(*totalC1PV,
"total Concentration w-Comp - vector");
181 writer.attachCellData(*totalC2PV,
"total Concentration n-Comp - vector");
188 int eIdxGlobal =
problem().variables().index(element);
195 int eIdxGlobal =
problem().variables().index(element);
196 CellData& cellData =
problem().variables().cellData(eIdxGlobal);
214 transportedQuantity.resize((getPropValue<TypeTag, Properties::NumEq>() - 1));
215 for(
int compIdx = 0; compIdx < (getPropValue<TypeTag, Properties::NumEq>() - 1); compIdx++)
216 transportedQuantity[compIdx].resize(
problem_.gridView().size(0));
244 template<
class DataEntry>
247 int numComp = getPropValue<TypeTag, Properties::NumEq>() - 1;
248 for(
int compIdx = 0; compIdx < numComp; compIdx++)
250 if (entry[compIdx] < -1.0e-6)
282 Scalar cFLFactor = getParam<Scalar>(
"Impet.CFLFactor");
284 subCFLFactor_ = min(getParam<Scalar>(
"Impet.SubCFLFactor"), cFLFactor);
285 verbosity_ = getParam<int>(
"TimeManager.SubTimestepVerbosity");
290 std::cout<<
"max CFL-Number of "<<cFLFactor<<
", max Sub-CFL-Number of "
291 <<
subCFLFactor_<<
": Enable local time-stepping!" << std::endl;
305 static const int pressureType = getPropValue<TypeTag, Properties::PressureFormulation>();
335 Implementation &asImp_()
336 {
return *
static_cast<Implementation *
>(
this); }
339 const Implementation &asImp_()
const
340 {
return *
static_cast<const Implementation *
>(
this); }
361template<
class TypeTag>
363 TransportSolutionType& updateVec,
bool impet)
368 unsigned int size = problem_.gridView().size(0);
369 if (localTimeStepping_)
371 if (this->timeStepData_.size() != size)
372 this->timeStepData_.resize(size);
379 updateVec.resize(getPropValue<TypeTag, Properties::NumComponents>());
380 updateVec[wCompIdx].resize(problem_.gridView().size(0));
381 updateVec[nCompIdx].resize(problem_.gridView().size(0));
382 updateVec[wCompIdx] = 0;
383 updateVec[nCompIdx] = 0;
386 int restrictingCell = -1;
388 ComponentVector entries(0.);
391 for (
const auto& element : elements(problem().gridView()))
394 int eIdxGlobalI = problem().variables().index(element);
395 CellData& cellDataI = problem().variables().cellData(eIdxGlobalI);
398 double sumfactorin = 0;
399 double sumfactorout = 0;
402 for (
const auto& intersection : intersections(problem().gridView(), element))
404 int indexInInside = intersection.indexInInside();
407 if (intersection.neighbor())
408 asImp_().getFlux(entries, timestepFlux, intersection, cellDataI);
411 if (intersection.boundary())
412 asImp_().getFluxOnBoundary(entries, timestepFlux, intersection, cellDataI);
415 if (localTimeStepping_)
418 if (localData.
faceTargetDt[indexInInside] < accumulatedDt_ + dtThreshold_)
420 localData.
faceFluxes[indexInInside] = entries;
426 updateVec[wCompIdx][eIdxGlobalI] += entries[wCompIdx];
427 updateVec[nCompIdx][eIdxGlobalI] += entries[nCompIdx];
431 sumfactorin += timestepFlux[0];
432 sumfactorout += timestepFlux[1];
436 if (localTimeStepping_)
439 for (
int i=0; i < 2*dim; i++)
441 updateVec[wCompIdx][eIdxGlobalI] += localData.
faceFluxes[i][wCompIdx];
442 updateVec[nCompIdx][eIdxGlobalI] += localData.
faceFluxes[i][nCompIdx];
447 PrimaryVariables q(NAN);
448 problem().source(q, element);
449 updateVec[wCompIdx][eIdxGlobalI] += q[contiWEqIdx];
450 updateVec[nCompIdx][eIdxGlobalI] += q[contiNEqIdx];
454 sumfactorin = max(sumfactorin,sumfactorout)
455 / problem().spatialParams().porosity(element);
458 if (localTimeStepping_)
460 timeStepData_[eIdxGlobalI].dt = 1./sumfactorin;
461 if ( 1./sumfactorin < dt)
464 restrictingCell= eIdxGlobalI;
469 if ( 1./sumfactorin < dt)
472 restrictingCell= eIdxGlobalI;
480 using ElementMapper =
typename SolutionTypes::ElementMapper;
482 for (
int i = 0; i < updateVec.size(); i++)
484 DataHandle dataHandle(problem_.variables().elementMapper(), updateVec[i]);
485 problem_.gridView().template communicate<DataHandle>(dataHandle,
486 Dune::InteriorBorder_All_Interface,
487 Dune::ForwardCommunication);
490 if (localTimeStepping_)
494 TimeDataHandle timeDataHandle(problem_.elementMapper(), timeStepData_);
495 problem_.gridView().template communicate<TimeDataHandle>(timeDataHandle,
496 Dune::InteriorBorder_All_Interface,
497 Dune::ForwardCommunication);
500 dt = problem_.gridView().comm().min(dt);
505 Dune::dinfo <<
"Timestep restricted by CellIdx " << restrictingCell <<
" leads to dt = "
506 <<dt * getParam<Scalar>(
"Impet.CFLFactor")<< std::endl;
507 if (averagedFaces_ != 0)
508 Dune::dinfo <<
" Averageing done for " << averagedFaces_ <<
" faces. "<< std::endl;
517template<
class TypeTag>
520 if (this->enableLocalTimeStepping())
521 this->innerUpdate(updateVector);
523 updateConcentrations(updateVector, problem().timeManager().timeStepSize());
533template<
class TypeTag>
536 updateConcentrations(updateVector, dt);
546template<
class TypeTag>
550 for (
int i = 0; i< problem().gridView().size(0); i++)
552 CellData& cellDataI = problem().variables().cellData(i);
553 for(
int compIdx = 0; compIdx < getPropValue<TypeTag, Properties::NumComponents>(); compIdx++)
555 totalConcentration_[compIdx][i] += (updateVector[compIdx][i]*=dt);
556 cellDataI.setMassConcentration(compIdx, totalConcentration_[compIdx][i]);
577template<
class TypeTag>
579 Dune::FieldVector<Scalar, 2>& timestepFlux,
580 const Intersection& intersection,
586 auto elementI = intersection.inside();
587 int eIdxGlobalI = problem().variables().index(elementI);
590 const GlobalPosition globalPos = elementI.geometry().center();
591 const GlobalPosition& gravity_ = problem().gravity();
593 Scalar
volume = elementI.geometry().volume();
596 Scalar pressI = problem().pressureModel().pressure(eIdxGlobalI);
597 Scalar pcI = cellDataI.capillaryPressure();
598 DimMatrix K_I(problem().spatialParams().intrinsicPermeability(elementI));
600 const auto fluidMatrixInteraction = problem().spatialParams().fluidMatrixInteractionAtPos(elementI.geometry().center());
601 PhaseVector SmobI(0.);
603 SmobI[wPhaseIdx] = max((cellDataI.saturation(wPhaseIdx)
604 - fluidMatrixInteraction.pcSwCurve().effToAbsParams().swr())
606 SmobI[nPhaseIdx] = max((cellDataI.saturation(nPhaseIdx)
607 - fluidMatrixInteraction.pcSwCurve().effToAbsParams().snr())
610 Scalar densityWI (0.), densityNWI(0.);
611 densityWI= cellDataI.density(wPhaseIdx);
612 densityNWI = cellDataI.density(nPhaseIdx);
615 GlobalPosition unitOuterNormal = intersection.centerUnitOuterNormal();
617 unitOuterNormal *= -1.0;
618 Scalar faceArea = intersection.geometry().volume();
625 Dune::FieldVector<Scalar, 2> factor (0.);
626 Dune::FieldVector<Scalar, 2> updFactor (0.);
628 PhaseVector potential(0.);
631 auto neighbor = intersection.outside();
632 int eIdxGlobalJ = problem().variables().index(neighbor);
633 CellData& cellDataJ = problem().variables().cellData(eIdxGlobalJ);
636 const GlobalPosition& globalPosNeighbor = neighbor.geometry().center();
639 GlobalPosition distVec = globalPosNeighbor - globalPos;
641 Scalar dist = distVec.two_norm();
643 GlobalPosition unitDistVec(distVec);
647 Scalar densityWJ (0.), densityNWJ(0.);
648 densityWJ = cellDataJ.density(wPhaseIdx);
649 densityNWJ = cellDataJ.density(nPhaseIdx);
652 double densityW_mean = (densityWI + densityWJ) * 0.5;
653 double densityNW_mean = (densityNWI + densityNWJ) * 0.5;
655 double pressJ = problem().pressureModel().pressure(eIdxGlobalJ);
656 Scalar pcJ = cellDataJ.capillaryPressure();
659 DimMatrix meanK_(0.);
662 problem().spatialParams().intrinsicPermeability(neighbor));
663 Dune::FieldVector<Scalar,dim> K(0);
664 meanK_.umv(unitDistVec,K);
667 switch (pressureType)
671 potential[wPhaseIdx] = (pressI - pressJ) / (dist);
672 potential[nPhaseIdx] = (pressI - pressJ + pcI - pcJ) / (dist);
677 potential[wPhaseIdx] = (pressI - pressJ - pcI + pcJ) / (dist);
678 potential[nPhaseIdx] = (pressI - pressJ) / (dist);
683 potential[nPhaseIdx] += (gravity_ * unitDistVec) * densityNW_mean;
684 potential[wPhaseIdx] += (gravity_ * unitDistVec) * densityW_mean;
686 potential[wPhaseIdx] *= fabs(K * unitOuterNormal);
687 potential[nPhaseIdx] *= fabs(K * unitOuterNormal);
690 Dune::FieldVector<bool, NumPhases> doUpwinding(
true);
691 PhaseVector lambda(0.);
692 for(
int phaseIdx = 0; phaseIdx < NumPhases; phaseIdx++)
695 if(phaseIdx == wPhaseIdx)
696 contiEqIdx = contiWEqIdx;
698 contiEqIdx = contiNEqIdx;
700 if(!impet_ || restrictFluxInTransport_==0)
702 if (potential[phaseIdx] > 0.)
704 lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
705 cellDataI.setUpwindCell(intersection.indexInInside(), contiEqIdx,
true);
708 else if(potential[phaseIdx] < 0.)
710 lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
711 cellDataI.setUpwindCell(intersection.indexInInside(), contiEqIdx,
false);
715 doUpwinding[phaseIdx] =
false;
716 cellDataI.setUpwindCell(intersection.indexInInside(), contiEqIdx,
false);
717 cellDataJ.setUpwindCell(intersection.indexInOutside(), contiEqIdx,
false);
722 bool wasUpwindCell = cellDataI.isUpwindCell(intersection.indexInInside(), contiEqIdx);
724 if (potential[phaseIdx] > 0. && wasUpwindCell)
725 lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
726 else if (potential[phaseIdx] < 0. && !wasUpwindCell)
727 lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
729 else if(restrictFluxInTransport_ == 2)
730 doUpwinding[phaseIdx] =
false;
734 if (potential[phaseIdx] > 0. && (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellDataJ.mobility(phaseIdx), 0.0, 1.0e-30)
735 || (cellDataI.wasRefined() && cellDataJ.wasRefined() && elementI.father() == neighbor.father())))
736 lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
737 else if (potential[phaseIdx] < 0. && (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellDataI.mobility(phaseIdx), 0.0, 1.0e-30)
738 || (cellDataI.wasRefined() && cellDataJ.wasRefined() && elementI.father() == neighbor.father())))
739 lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
741 doUpwinding[phaseIdx] =
false;
747 if(!doUpwinding[phaseIdx])
750 if(cellDataI.mobility(phaseIdx)+cellDataJ.mobility(phaseIdx)==0.)
752 potential[phaseIdx] = 0;
757 fluxEntries[wCompIdx] -= potential[phaseIdx] * faceArea /
volume
758 *
harmonicMean(cellDataI.massFraction(phaseIdx, wCompIdx) * cellDataI.mobility(phaseIdx) * cellDataI.density(phaseIdx),
759 cellDataJ.massFraction(phaseIdx, wCompIdx) * cellDataJ.mobility(phaseIdx) * cellDataJ.density(phaseIdx));
760 fluxEntries[nCompIdx] -= potential[phaseIdx] * faceArea /
volume
761 *
harmonicMean(cellDataI.massFraction(phaseIdx, nCompIdx) * cellDataI.mobility(phaseIdx) * cellDataI.density(phaseIdx),
762 cellDataJ.massFraction(phaseIdx, nCompIdx) * cellDataJ.mobility(phaseIdx) * cellDataJ.density(phaseIdx));
766 timestepFlux[0] += max(0.,
767 - potential[phaseIdx] * faceArea /
volume
768 *
harmonicMean(cellDataI.mobility(phaseIdx),cellDataJ.mobility(phaseIdx)));
770 timestepFlux[1] += max(0.,
771 potential[phaseIdx] * faceArea /
volume
772 *
harmonicMean(cellDataI.mobility(phaseIdx),cellDataJ.mobility(phaseIdx))/SmobI[phaseIdx]);
775 if(!(cellDataI.wasRefined() && cellDataJ.wasRefined() && elementI.father() == neighbor.father())
776 && eIdxGlobalI > eIdxGlobalJ)
779 #if DUNE_MINIMAL_DEBUG_LEVEL < 3
781 if(eIdxGlobalI > eIdxGlobalJ)
782 Dune::dinfo <<
"harmonicMean flux of phase" << phaseIdx <<
" used from cell" << eIdxGlobalI<<
" into " << eIdxGlobalJ
783 <<
" ; TE upwind I = "<< cellDataI.isUpwindCell(intersection.indexInInside(), contiEqIdx)
784 <<
" but pot = "<< potential[phaseIdx] << std::endl;
789 potential[phaseIdx] = 0;
795 double velocityJIw = max((-lambda[wPhaseIdx] * potential[wPhaseIdx]) * faceArea /
volume, 0.0);
796 double velocityIJw = max(( lambda[wPhaseIdx] * potential[wPhaseIdx]) * faceArea /
volume, 0.0);
797 double velocityJIn = max((-lambda[nPhaseIdx] * potential[nPhaseIdx]) * faceArea /
volume, 0.0);
798 double velocityIJn = max(( lambda[nPhaseIdx] * potential[nPhaseIdx]) * faceArea /
volume, 0.0);
801 timestepFlux[0] += velocityJIw + velocityJIn;
803 double foutw = velocityIJw/SmobI[wPhaseIdx];
804 double foutn = velocityIJn/SmobI[nPhaseIdx];
807 if (isnan(foutw) || isinf(foutw) || foutw < 0) foutw = 0;
808 if (isnan(foutn) || isinf(foutn) || foutn < 0) foutn = 0;
809 timestepFlux[1] += foutw + foutn;
811 fluxEntries[wCompIdx] +=
812 velocityJIw * cellDataJ.massFraction(wPhaseIdx, wCompIdx) * densityWJ
813 - velocityIJw * cellDataI.massFraction(wPhaseIdx, wCompIdx) * densityWI
814 + velocityJIn * cellDataJ.massFraction(nPhaseIdx, wCompIdx) * densityNWJ
815 - velocityIJn * cellDataI.massFraction(nPhaseIdx, wCompIdx) * densityNWI;
816 fluxEntries[nCompIdx] +=
817 velocityJIw * cellDataJ.massFraction(wPhaseIdx, nCompIdx) * densityWJ
818 - velocityIJw * cellDataI.massFraction(wPhaseIdx, nCompIdx) * densityWI
819 + velocityJIn * cellDataJ.massFraction(nPhaseIdx, nCompIdx) * densityNWJ
820 - velocityIJn * cellDataI.massFraction(nPhaseIdx, nCompIdx) * densityNWI;
883template<
class TypeTag>
885 Dune::FieldVector<Scalar, 2>& timestepFlux,
886 const Intersection& intersection,
887 const CellData& cellDataI)
891 auto elementI = intersection.inside();
892 int eIdxGlobalI = problem().variables().index(elementI);
895 const GlobalPosition globalPos = elementI.geometry().center();
898 Scalar
volume = elementI.geometry().volume();
899 const GlobalPosition& gravity_ = problem().gravity();
901 Scalar pressI = problem().pressureModel().pressure(eIdxGlobalI);
902 Scalar pcI = cellDataI.capillaryPressure();
903 DimMatrix K_I(problem().spatialParams().intrinsicPermeability(elementI));
905 if(regulateBoundaryPermeability)
907 int axis = intersection.indexInInside() / 2;
908 if(K_I[axis][axis] < minimalBoundaryPermeability)
909 K_I[axis][axis] = minimalBoundaryPermeability;
912 const auto fluidMatrixInteraction = problem().spatialParams().fluidMatrixInteractionAtPos(elementI.geometry().center());
913 Scalar SwmobI = max((cellDataI.saturation(wPhaseIdx)
914 - fluidMatrixInteraction.pcSwCurve().effToAbsParams().swr())
916 Scalar SnmobI = max((cellDataI.saturation(nPhaseIdx)
917 - fluidMatrixInteraction.pcSwCurve().effToAbsParams().snr())
920 Scalar densityWI (0.), densityNWI(0.);
921 densityWI= cellDataI.density(wPhaseIdx);
922 densityNWI = cellDataI.density(nPhaseIdx);
925 GlobalPosition unitOuterNormal = intersection.centerUnitOuterNormal();
927 unitOuterNormal *= -1.0;
928 Scalar faceArea = intersection.geometry().volume();
931 Dune::FieldVector<Scalar, 2> factor (0.);
932 Dune::FieldVector<Scalar, 2> updFactor (0.);
934 PhaseVector potential(0.);
936 const GlobalPosition& globalPosFace = intersection.geometry().center();
939 GlobalPosition distVec = globalPosFace - globalPos;
941 Scalar dist = distVec.two_norm();
943 GlobalPosition unitDistVec(distVec);
947 FluidState BCfluidState;
950 BoundaryTypes bcTypes;
951 problem().boundaryTypes(bcTypes, intersection);
954 if (bcTypes.isDirichlet(contiWEqIdx))
957 PhaseVector pressBound(0.);
961 this->evalBoundary(globalPosFace,
967 Scalar densityWBound = BCfluidState.density(wPhaseIdx);
968 Scalar densityNWBound = BCfluidState.density(nPhaseIdx);
971 if(getPropValue<TypeTag, Properties::EnableCapillarity>())
972 pcBound = (BCfluidState.pressure(nPhaseIdx) - BCfluidState.pressure(wPhaseIdx));
974 double densityW_mean = (densityWI + densityWBound) / 2;
975 double densityNW_mean = (densityNWI + densityNWBound) / 2;
978 Dune::FieldVector<Scalar,dim> K(0);
979 K_I.umv(unitDistVec,K);
982 switch (pressureType)
986 potential[wPhaseIdx] = (pressI - pressBound[wPhaseIdx]) / dist;
987 potential[nPhaseIdx] = (pressI + pcI - pressBound[wPhaseIdx] - pcBound)
993 potential[wPhaseIdx] = (pressI - pcI - pressBound[nPhaseIdx] + pcBound)
995 potential[nPhaseIdx] = (pressI - pressBound[nPhaseIdx]) / dist;
999 potential[wPhaseIdx] += (gravity_ * unitDistVec) * densityW_mean;
1000 potential[nPhaseIdx] += (gravity_ * unitDistVec) * densityNW_mean;
1002 potential[wPhaseIdx] *= fabs(K * unitOuterNormal);
1003 potential[nPhaseIdx] *= fabs(K * unitOuterNormal);
1005 const auto fluidMatrixInteraction = problem().spatialParams().fluidMatrixInteractionAtPos(elementI.geometry().center());
1008 PhaseVector lambda(0.);
1009 if (potential[wPhaseIdx] >= 0.)
1010 lambda[wPhaseIdx] = cellDataI.mobility(wPhaseIdx);
1013 if(getPropValue<TypeTag, Properties::BoundaryMobility>()==Indices::satDependent)
1014 lambda[wPhaseIdx] = BCfluidState.saturation(wPhaseIdx) / viscosityWBound;
1016 lambda[wPhaseIdx] = fluidMatrixInteraction.krw(BCfluidState.saturation(wPhaseIdx)) / viscosityWBound;
1018 if (potential[nPhaseIdx] >= 0.)
1019 lambda[nPhaseIdx] = cellDataI.mobility(nPhaseIdx);
1022 if(getPropValue<TypeTag, Properties::BoundaryMobility>()==Indices::satDependent)
1023 lambda[nPhaseIdx] = BCfluidState.saturation(nPhaseIdx) / viscosityNWBound;
1025 lambda[nPhaseIdx] = fluidMatrixInteraction.krn(BCfluidState.saturation(wPhaseIdx)) / viscosityNWBound;
1029 double velocityJIw = max((-lambda[wPhaseIdx] * potential[wPhaseIdx]) * faceArea /
volume, 0.0);
1030 double velocityIJw = max(( lambda[wPhaseIdx] * potential[wPhaseIdx]) * faceArea /
volume, 0.0);
1031 double velocityJIn = max((-lambda[nPhaseIdx] * potential[nPhaseIdx]) * faceArea /
volume, 0.0);
1032 double velocityIJn = max(( lambda[nPhaseIdx] * potential[nPhaseIdx]) * faceArea /
volume, 0.0);
1035 timestepFlux[0] = velocityJIw + velocityJIn;
1037 double foutw = velocityIJw/SwmobI;
1038 double foutn = velocityIJn/SnmobI;
1041 if (isnan(foutw) || isinf(foutw) || foutw < 0) foutw = 0;
1042 if (isnan(foutn) || isinf(foutn) || foutn < 0) foutn = 0;
1043 timestepFlux[1] = foutw + foutn;
1045 fluxEntries[wCompIdx] =
1046 + velocityJIw * BCfluidState.massFraction(wPhaseIdx, wCompIdx) * densityWBound
1047 - velocityIJw * cellDataI.massFraction(wPhaseIdx, wCompIdx) * densityWI
1048 + velocityJIn * BCfluidState.massFraction(nPhaseIdx, wCompIdx) * densityNWBound
1049 - velocityIJn * cellDataI.massFraction(nPhaseIdx, wCompIdx) * densityNWI ;
1050 fluxEntries[nCompIdx] =
1051 velocityJIw * BCfluidState.massFraction(wPhaseIdx, nCompIdx) * densityWBound
1052 - velocityIJw * cellDataI.massFraction(wPhaseIdx, nCompIdx) * densityWI
1053 + velocityJIn * BCfluidState.massFraction(nPhaseIdx, nCompIdx) * densityNWBound
1054 - velocityIJn * cellDataI.massFraction(nPhaseIdx, nCompIdx) * densityNWI ;
1056 else if (bcTypes.isNeumann(contiWEqIdx))
1059 PrimaryVariables J(NAN);
1060 problem().neumann(J, intersection);
1061 fluxEntries[wCompIdx] = - J[contiWEqIdx] * faceArea /
volume;
1062 fluxEntries[nCompIdx] = - J[contiNEqIdx] * faceArea /
volume;
1065 #define cflIgnoresNeumann
1066 #ifdef cflIgnoresNeumann
1067 timestepFlux[0] = 0;
1068 timestepFlux[1] = 0;
1070 double inflow = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW;
1073 timestepFlux[0] = updFactor[wCompIdx] / densityW
1074 + updFactor[nCompIdx] / densityNW;
1075 timestepFlux[1] = -(updFactor[wCompIdx] / densityW /SwmobI
1076 + updFactor[nCompIdx] / densityNW / SnmobI);
1080 timestepFlux[0] = -(updFactor[wCompIdx] / densityW
1081 + updFactor[nCompIdx] / densityNW);
1082 timestepFlux[1] = updFactor[wCompIdx] / densityW /SwmobI
1083 + updFactor[nCompIdx] / densityNW / SnmobI;
1105template<
class TypeTag>
1107 const Intersection& intersection,
1108 FluidState &BCfluidState,
1109 PhaseVector &pressBound)
1114 auto element = intersection.inside();
1116 PrimaryVariables primaryVariablesOnBoundary(0.);
1117 problem().dirichlet(primaryVariablesOnBoundary, intersection);
1119 const auto fluidMatrixInteraction = problem().spatialParams().fluidMatrixInteractionAtPos(element.geometry().center());
1122 typename Indices::BoundaryFormulation bcType;
1123 problem().boundaryFormulation(bcType, intersection);
1126 Scalar satBound = primaryVariablesOnBoundary[contiWEqIdx];
1128 if(getPropValue<TypeTag, Properties::EnableCapillarity>())
1130 Scalar pcBound = fluidMatrixInteraction.pc(satBound);
1131 switch (pressureType)
1135 pressBound[wPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
1136 pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx] + pcBound;
1141 pressBound[wPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx] - pcBound;
1142 pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
1148 pressBound[wPhaseIdx] = pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
1151 problem().temperatureAtPos(globalPosFace));
1153 else if (bcType == Indices::concentration)
1156 pressBound[wPhaseIdx] = pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
1157 Scalar Z0Bound = primaryVariablesOnBoundary[contiWEqIdx];
1159 problem().temperatureAtPos(globalPosFace));
1161 if(getPropValue<TypeTag, Properties::EnableCapillarity>())
1163 Scalar pcBound = fluidMatrixInteraction.pc(BCfluidState.saturation(wPhaseIdx));
1166 for(
int iter=0; iter < maxiter; iter++)
1169 switch (pressureType)
1173 pressBound[wPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
1174 pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx]
1180 pressBound[wPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx]
1182 pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
1188 Scalar oldPc = pcBound;
1191 problem().temperatureAtPos(globalPosFace));
1192 pcBound = fluidMatrixInteraction.pc(BCfluidState.saturation(wPhaseIdx));
1196 if (abs(oldPc - pcBound) < 10.0)
1202 DUNE_THROW(Dune::NotImplemented,
"Boundary Formulation neither Concentration nor Saturation??");
1205template<
class TypeTag>
1208 dt = std::numeric_limits<Scalar>::max();
1211 for (
const auto& element : elements(problem_.gridView()))
1214 if (element.partitionType() != Dune::InteriorEntity)
1221 int eIdxGlobalI = problem_.variables().index(element);
1226 using FaceDt = std::unordered_map<int, Scalar>;
1230 for (
const auto& intersection : intersections(problem_.gridView(), element))
1232 int indexInInside = intersection.indexInInside();
1234 if (intersection.neighbor())
1236 auto neighbor = intersection.outside();
1237 int eIdxGlobalJ = problem_.variables().index(neighbor);
1239 int levelI = element.level();
1240 int levelJ = neighbor.level();
1243 if (eIdxGlobalI < eIdxGlobalJ && levelI <= levelJ)
1247 int indexInOutside = intersection.indexInOutside();
1249 if (localDataI.
faceTargetDt[indexInInside] < accumulatedDt_ + dtThreshold_
1250 || localDataJ.
faceTargetDt[indexInOutside] < accumulatedDt_ + dtThreshold_)
1252 Scalar timeStep = min(localDataI.
dt, localDataJ.
dt);
1254 if (levelI < levelJ)
1256 typename FaceDt::iterator it = faceDt.find(indexInInside);
1257 if (it != faceDt.end())
1259 it->second = min(it->second, timeStep);
1263 faceDt.insert(std::make_pair(indexInInside, timeStep));
1268 localDataI.
faceTargetDt[indexInInside] += subCFLFactor_ * timeStep;
1269 localDataJ.
faceTargetDt[indexInOutside] += subCFLFactor_ * timeStep;
1272 dt = min(dt, timeStep);
1276 else if (intersection.boundary())
1278 if (localDataI.
faceTargetDt[indexInInside] < accumulatedDt_ + dtThreshold_)
1280 localDataI.
faceTargetDt[indexInInside] += subCFLFactor_ * localDataI.
dt;
1281 dt = min(dt, subCFLFactor_ * localDataI.
dt);
1285 if (faceDt.size() > 0)
1287 typename FaceDt::iterator it = faceDt.begin();
1288 for (;it != faceDt.end();++it)
1290 localDataI.
faceTargetDt[it->first] += subCFLFactor_ * it->second;
1293 for (
const auto& intersection : intersections(problem_.gridView(), element))
1295 if (intersection.neighbor())
1297 int indexInInside = intersection.indexInInside();
1299 it = faceDt.find(indexInInside);
1300 if (it != faceDt.end())
1302 auto neighbor = intersection.outside();
1303 int eIdxGlobalJ = problem_.variables().index(neighbor);
1307 int indexInOutside = intersection.indexInOutside();
1309 localDataJ.
faceTargetDt[indexInOutside] += subCFLFactor_ * it->second;
1319 using ElementMapper =
typename SolutionTypes::ElementMapper;
1322 TimeDataHandle timeDataHandle(problem_.elementMapper(), timeStepData_);
1323 problem_.gridView().template communicate<TimeDataHandle>(timeDataHandle,
1324 Dune::InteriorBorder_All_Interface,
1325 Dune::ForwardCommunication);
1327 dt = problem_.gridView().comm().min(dt);
1331template<
class TypeTag>
1334 if (localTimeStepping_)
1336 Scalar realDt = problem_.timeManager().timeStepSize();
1338 Scalar subDt = realDt;
1340 updatedTargetDt_(subDt);
1342 Scalar accumulatedDtOld = accumulatedDt_;
1343 accumulatedDt_ += subDt;
1345 Scalar t = problem_.timeManager().time();
1347 if (accumulatedDt_ < realDt)
1352 Scalar dtCorrection = min(0.0, realDt - accumulatedDt_);
1353 subDt += dtCorrection;
1356 std::cout<<
" Sub-time-step size: "<<subDt<< std::endl;
1358 bool stopTimeStep =
false;
1359 int size = problem_.gridView().size(0);
1360 for (
int i = 0; i < size; i++)
1363 int transportedQuantities = getPropValue<TypeTag, Properties::NumEq>() - 1;
1364 for (
int eqNumber = 0; eqNumber < transportedQuantities; eqNumber++)
1366 newVal[eqNumber] = totalConcentration_[eqNumber][i];
1367 newVal[eqNumber] += updateVec[eqNumber][i] * subDt;
1369 if (!asImp_().inPhysicalRange(newVal))
1371 stopTimeStep =
true;
1380 rank = problem_.gridView().comm().rank();
1382 rank = problem_.gridView().comm().max(rank);
1383 problem_.gridView().comm().broadcast(&stopTimeStep,1,rank);
1387 if (stopTimeStep && accumulatedDtOld > dtThreshold_)
1389 problem_.timeManager().setTimeStepSize(accumulatedDtOld);
1394 asImp_().updateTransportedQuantity(updateVec, subDt);
1398 if (accumulatedDt_ >= realDt)
1403 problem_.pressureModel().updateMaterialLaws();
1404 problem_.model().updateTransport(t, subDt, updateVec);
1406 updatedTargetDt_(subDt);
1408 accumulatedDtOld = accumulatedDt_;
1409 accumulatedDt_ += subDt;
1414 asImp_().updateTransportedQuantity(updateVec, realDt);
1417 resetTimeStepData_();
Define some often used mathematical functions.
Determines the pressures and saturations of all fluid phases given the total mass of all components.
Contains a class to exchange entries of a vector.
constexpr Scalar harmonicMean(Scalar x, Scalar y, Scalar wx=1.0, Scalar wy=1.0) noexcept
Calculate the (weighted) harmonic mean of two scalar values.
Definition: math.hh:69
void harmonicMeanMatrix(Dune::FieldMatrix< Scalar, m, n > &K, const Dune::FieldMatrix< Scalar, m, n > &Ki, const Dune::FieldMatrix< Scalar, m, n > &Kj)
Calculate the harmonic mean of a fixed-size matrix.
Definition: math.hh:109
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
Scalar volume(Shape shape, Scalar inscribedRadius)
Returns the volume of a given geometry based on the inscribed radius.
Definition: poreproperties.hh:73
Flash calculation routines for compositional sequential models.
Definition: compositionalflash.hh:51
static void saturationFlash2p2c(FluidState &fluidState, const Scalar saturation, const PhaseVector &phasePressure, const Scalar temperature)
Definition: compositionalflash.hh:219
static void concentrationFlash2p2c(FluidState &fluidState, const Scalar Z0, const PhaseVector &phasePressure, const Scalar temperature)
Definition: compositionalflash.hh:91
A data handle class to exchange entries of a vector.
Definition: vectorcommdatahandle.hh:78
Compositional transport step in a Finite Volume discretization.
Definition: fvtransport.hh:60
virtual void update(const Scalar t, Scalar &dt, TransportSolutionType &updateVec, bool impet=false)
Calculate the update vector and determine timestep size This method calculates the update vector of ...
Definition: fvtransport.hh:362
void evalBoundary(GlobalPosition globalPosFace, const Intersection &intersection, FluidState &BCfluidState, PhaseVector &pressBound)
Evaluate the boundary conditions.
Definition: fvtransport.hh:1106
void deserializeEntity(std::istream &instream, const Element &element)
Function needed for restart option of the transport model: Read in.
Definition: fvtransport.hh:193
bool regulateBoundaryPermeability
Enables regulation of permeability in the direction of a Dirichlet Boundary Condition.
Definition: fvtransport.hh:309
void updateConcentrations(TransportSolutionType &updateVector, Scalar dt)
Updates the concentrations once an update is calculated.
Definition: fvtransport.hh:547
bool localTimeStepping_
Definition: fvtransport.hh:329
Scalar subCFLFactor_
Definition: fvtransport.hh:328
int restrictFluxInTransport_
Restriction of flux on new pressure field if direction reverses from the pressure equation.
Definition: fvtransport.hh:307
FVTransport2P2C(Problem &problem)
Constructs a FVTransport2P2C object Currently, the compositional transport scheme can not be applied ...
Definition: fvtransport.hh:272
const Scalar dtThreshold_
Threshold for sub-time-stepping routine.
Definition: fvtransport.hh:316
void initialize()
Set the initial values before the first pressure equation This method is called before first pressure...
Definition: fvtransport.hh:153
void addOutputVtkFields(MultiWriter &writer)
Write transport variables into the output files.
Definition: fvtransport.hh:170
void getSource(Scalar &update, const Element &element, CellData &cellDataI)
Definition: fvtransport.hh:235
Scalar & totalConcentration(int compIdx, int eIdxGlobal)
Return the the total concentration stored in the transport vector To get real cell values,...
Definition: fvtransport.hh:230
Problem & problem_
Definition: fvtransport.hh:300
void innerUpdate(TransportSolutionType &updateVec)
Definition: fvtransport.hh:1332
int averagedFaces_
number of faces were flux was restricted
Definition: fvtransport.hh:302
void updatedTargetDt_(Scalar &dt)
Definition: fvtransport.hh:1206
void resetTimeStepData_()
Definition: fvtransport.hh:322
Dune::FieldVector< Scalar, 2 > EntryType
Definition: fvtransport.hh:107
bool impet_
indicating if we are in an estimate (false) or real impet (true) step.
Definition: fvtransport.hh:301
void serializeEntity(std::ostream &outstream, const Element &element)
Function needed for restart option of the transport model: Write out.
Definition: fvtransport.hh:186
void getTransportedQuantity(TransportSolutionType &transportedQuantity)
Return the vector of the transported quantity For an immiscible IMPES scheme, this is the saturation....
Definition: fvtransport.hh:211
void getFlux(ComponentVector &fluxEntries, EntryType ×tepFlux, const Intersection &intersection, CellData &cellDataI)
Get flux at an interface between two cells The flux through is calculated according to the underlyin...
Definition: fvtransport.hh:578
Problem & problem()
Acess function for the current problem.
Definition: fvtransport.hh:121
std::vector< LocalTimesteppingData > timeStepData_
Stores data for sub-time-stepping.
Definition: fvtransport.hh:318
int verbosity_
Definition: fvtransport.hh:330
bool inPhysicalRange(DataEntry &entry)
Function to control the abort of the transport-sub-time-stepping depending on a physical parameter ra...
Definition: fvtransport.hh:245
TransportSolutionType totalConcentration_
private vector of transported primary variables
Definition: fvtransport.hh:299
Scalar accumulatedDt_
Current time-interval in sub-time-stepping routine.
Definition: fvtransport.hh:314
virtual ~FVTransport2P2C()
Definition: fvtransport.hh:294
void updateTransportedQuantity(TransportSolutionType &updateVector)
Updates the transported quantity once an update is calculated.
Definition: fvtransport.hh:518
void getFluxOnBoundary(ComponentVector &fluxEntries, EntryType ×tepFlux, const Intersection &intersection, const CellData &cellDataI)
Get flux on Boundary.
Definition: fvtransport.hh:884
static const int pressureType
gives kind of pressure used ( , , )
Definition: fvtransport.hh:305
bool switchNormals
Definition: fvtransport.hh:312
bool enableLocalTimeStepping()
Function to check if local time stepping is activated.
Definition: fvtransport.hh:259
Dune::FieldVector< Scalar, 2 > TimeStepFluxType
Definition: fvtransport.hh:108
Scalar minimalBoundaryPermeability
Minimal limit for the boundary permeability.
Definition: fvtransport.hh:311
Definition: fvtransport.hh:112
Dune::FieldVector< EntryType, 2 *dim > faceFluxes
Definition: fvtransport.hh:113
Dune::FieldVector< Scalar, 2 *dim > faceTargetDt
Definition: fvtransport.hh:114
LocalTimesteppingData()
Definition: fvtransport.hh:116
Scalar dt
Definition: fvtransport.hh:115
Defines the properties required for the sequential 2p2c models.