56 GetPropType<TypeTag, Properties::GridView>,
57 GetPropType<TypeTag, Properties::ElementMapper>,
58 GetPropType<TypeTag, Properties::Scalar>>
61 GetPropType<TypeTag, Properties::GridView>,
62 GetPropType<TypeTag, Properties::ElementMapper>,
63 GetPropType<TypeTag, Properties::Scalar>>;
64 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
65 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
66 using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
67 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
68 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
70 enum { numPhases = FluidSystem::numPhases };
74 :
BaseType(simulator.vanguard().cartesianIndexMapper(),
75 simulator.vanguard().gridView(),
76 simulator.model().elementMapper(),
77 simulator.vanguard().eclState())
78 , simulator_(simulator)
88 if (this->enableThresholdPressure_ && !this->thpresDefault_.empty()) {
89 this->computeDefaultThresholdPressures_();
90 this->applyExplicitThresholdPressures_();
96 void computeDefaultThresholdPressures_()
98 const auto& vanguard = simulator_.vanguard();
99 const auto& gridView = vanguard.gridView();
101 using Toolbox = MathToolbox<Evaluation>;
104 ElementContext elemCtx(simulator_);
105 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
106 elemCtx.updateAll(elem);
107 const auto& stencil = elemCtx.stencil(0);
109 for (
unsigned scvfIdx = 0; scvfIdx < stencil.numInteriorFaces(); ++ scvfIdx) {
110 const auto& face = stencil.interiorFace(scvfIdx);
112 unsigned i = face.interiorIndex();
113 unsigned j = face.exteriorIndex();
115 unsigned insideElemIdx = elemCtx.globalSpaceIndex(i, 0);
116 unsigned outsideElemIdx = elemCtx.globalSpaceIndex(j, 0);
118 unsigned equilRegionInside = this->elemEquilRegion_[insideElemIdx];
119 unsigned equilRegionOutside = this->elemEquilRegion_[outsideElemIdx];
121 if (equilRegionInside == equilRegionOutside)
126 const Evaluation& trans = simulator_.problem().transmissibility(elemCtx, i, j);
127 Scalar faceArea = face.area();
128 if (std::abs(faceArea*getValue(trans)) < 1e-18)
134 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, 0);
135 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
136 unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
137 const auto& up = elemCtx.intensiveQuantities(upIdx, 0);
139 if (up.mobility(phaseIdx) > 0.0) {
140 Scalar phaseVal = Toolbox::value(extQuants.pressureDifference(phaseIdx));
141 pth = std::max(pth, std::abs(phaseVal));
145 int offset1 = equilRegionInside*this->numEquilRegions_ + equilRegionOutside;
146 int offset2 = equilRegionOutside*this->numEquilRegions_ + equilRegionInside;
148 this->thpresDefault_[offset1] = std::max(this->thpresDefault_[offset1], pth);
149 this->thpresDefault_[offset2] = std::max(this->thpresDefault_[offset2], pth);
155 for (
unsigned i = 0; i < this->thpresDefault_.size(); ++i)
156 this->thpresDefault_[i] = gridView.comm().max(this->thpresDefault_[i]);
158 this->logPressures();
161 const Simulator& simulator_;