86 return std::max(val1, val2);
91 auto equilRegion1Idx = elemEquilRegion_[elem1Idx];
92 auto equilRegion2Idx = elemEquilRegion_[elem2Idx];
94 if (equilRegion1Idx == equilRegion2Idx) {
98 return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx];
101template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
105 const auto& simConfig = eclState_.getSimulationConfig();
107 enableThresholdPressure_ = simConfig.useThresholdPressure();
108 if (!enableThresholdPressure_) {
112 numEquilRegions_ = eclState_.getTableManager().getEqldims().getNumEquilRegions();
113 const decltype(numEquilRegions_) maxRegions =
114 std::numeric_limits<std::decay_t<
decltype(elemEquilRegion_[0])>>::max();
116 if (numEquilRegions_ > maxRegions) {
119 OPM_THROW(std::invalid_argument,
120 (fmt::format(
"The maximum number of supported "
121 "equilibration regions by OPM flow is {}, but "
123 maxRegions, numEquilRegions_)));
126 if (numEquilRegions_ > 2048) {
128 OpmLog::warning(fmt::format(
"Number of equilibration regions is {}, which is "
129 "rather large. Note, that this might "
130 "have a negative impact on performance "
131 "and memory consumption.", numEquilRegions_));
135 elemEquilRegion_ = lookUpData_.
136 template assignFieldPropsIntOnLeaf<short unsigned int>(eclState_.fieldProps(),
144 if (simConfig.getThresholdPressure().restart()) {
149 thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0);
150 thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0);
153template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
157 const SimulationConfig& simConfig = eclState_.getSimulationConfig();
158 const auto& thpres = simConfig.getThresholdPressure();
162 for (
const auto& elem : elements(gridView_, Dune::Partitions::interior)) {
163 for (
const auto& intersection : intersections(gridView_, elem)) {
164 if (intersection.boundary()) {
167 else if (!intersection.neighbor()) {
171 const auto& inside = intersection.inside();
172 const auto& outside = intersection.outside();
174 auto equilRegionInside = elemEquilRegion_[elementMapper_.index(inside)];
175 auto equilRegionOutside = elemEquilRegion_[elementMapper_.index(outside)];
177 if (thpres.hasRegionBarrier(equilRegionInside + 1, equilRegionOutside + 1)) {
179 if (thpres.hasThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1)) {
181 pth = thpres.getThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1);
185 unsigned offset = equilRegionInside*numEquilRegions_ + equilRegionOutside;
186 pth = thpresDefault_[offset];
189 unsigned offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside;
190 unsigned offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside;
192 thpres_[offset1] = pth;
193 thpres_[offset2] = pth;
199 if (thpres.ftSize() > 0) {
200 configureThpresft_();
204template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
205void GenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
209 const FaultCollection& faults = eclState_.getFaults();
210 const SimulationConfig& simConfig = eclState_.getSimulationConfig();
211 const auto& thpres = simConfig.getThresholdPressure();
214 int numFaults = faults.size();
215 int numCartesianElem = eclState_.getInputGrid().getCartesianSize();
216 thpresftValues_.resize(numFaults, -1.0);
217 cartElemFaultIdx_.resize(numCartesianElem, -1);
218 for (std::size_t faultIdx = 0; faultIdx < faults.size(); ++faultIdx) {
219 const auto& fault = faults.getFault(faultIdx);
220 thpresftValues_[faultIdx] = thpres.getThresholdPressureFault(faultIdx);
221 for (
const FaultFace& face : fault) {
224 for (std::size_t cartElemIdx : face) {
225 cartElemFaultIdx_[cartElemIdx] = faultIdx;
231template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
236 if (!enableThresholdPressure_) {
240 return this->thpres_;
254GenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
267 " ", units.from_si(UnitSystem::measure::pressure, val),
288 str += lineFormat(i, j, thpres.getThresholdPressure(j, i));