28#ifndef OPM_FLOW_THRESHOLD_PRESSURE_HPP
29#define OPM_FLOW_THRESHOLD_PRESSURE_HPP
31#include <opm/material/densead/Evaluation.hpp>
32#include <opm/material/densead/Math.hpp>
54template <
class TypeTag>
56 GetPropType<TypeTag, Properties::GridView>,
57 GetPropType<TypeTag, Properties::ElementMapper>,
58 GetPropType<TypeTag, Properties::Scalar>>
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)
89 this->computeDefaultThresholdPressures_();
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);
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;
161 const Simulator& simulator_;
This class calculates the threshold pressure for grid faces according to the Eclipse Reference Manual...
Definition: FlowThresholdPressure.hpp:59
FlowThresholdPressure(const Simulator &simulator)
Definition: FlowThresholdPressure.hpp:73
void finishInit()
Actually compute the threshold pressures over a face as a pre-compute step.
Definition: FlowThresholdPressure.hpp:85
Definition: GenericThresholdPressure.hpp:43
bool enableThresholdPressure_
Definition: GenericThresholdPressure.hpp:113
void logPressures()
Definition: GenericThresholdPressure_impl.hpp:252
void finishInit()
Actually compute the threshold pressures over a face as a pre-compute step.
Definition: GenericThresholdPressure_impl.hpp:101
std::vector< unsigned short > elemEquilRegion_
Definition: GenericThresholdPressure.hpp:107
unsigned numEquilRegions_
Definition: GenericThresholdPressure.hpp:106
std::vector< GetPropType< TypeTag, Properties::Scalar > > thpresDefault_
Definition: GenericThresholdPressure.hpp:104
void applyExplicitThresholdPressures_()
Definition: GenericThresholdPressure_impl.hpp:150
Declare the properties used by the infrastructure code of the finite volume discretizations.
Defines the common properties required by the porous medium multi-phase models.
Definition: blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:235
The Opm property system, traits with inheritance.