83 Scalar val1 = (fault1Idx >= 0) ? thpresftValues_[fault1Idx] : 0.0;
84 Scalar val2 = (fault2Idx >= 0) ? thpresftValues_[fault2Idx] : 0.0;
85 return std::max(val1, val2);
90 auto equilRegion1Idx = elemEquilRegion_[elem1Idx];
91 auto equilRegion2Idx = elemEquilRegion_[elem2Idx];
93 if (equilRegion1Idx == equilRegion2Idx)
96 return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx];
99template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
103 const auto& simConfig = eclState_.getSimulationConfig();
105 enableThresholdPressure_ = simConfig.useThresholdPressure();
106 if (!enableThresholdPressure_)
109 numEquilRegions_ = eclState_.getTableManager().getEqldims().getNumEquilRegions();
110 const decltype(numEquilRegions_) maxRegions =
111 std::numeric_limits<std::decay_t<
decltype(elemEquilRegion_[0])>>::max();
113 if (numEquilRegions_ > maxRegions) {
116 OPM_THROW(std::invalid_argument,
117 (fmt::format(
"The maximum number of supported "
118 "equilibration regions by OPM flow is {}, but "
120 maxRegions, numEquilRegions_)));
123 if (numEquilRegions_ > 2048) {
125 OpmLog::warning(fmt::format(
"Number of equilibration regions is {}, which is "
126 "rather large. Note, that this might "
127 "have a negative impact on performance "
128 "and memory consumption.", numEquilRegions_));
132 elemEquilRegion_ = lookUpData_.template assignFieldPropsIntOnLeaf<short unsigned int>(eclState_.fieldProps(),
140 if (simConfig.getThresholdPressure().restart())
144 thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0);
145 thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0);
148template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
152 const SimulationConfig& simConfig = eclState_.getSimulationConfig();
153 const auto& thpres = simConfig.getThresholdPressure();
157 for (
const auto& elem : elements(gridView_, Dune::Partitions::interior)) {
158 for (
const auto& intersection : intersections(gridView_, elem)) {
159 if (intersection.boundary())
161 else if (!intersection.neighbor())
164 const auto& inside = intersection.inside();
165 const auto& outside = intersection.outside();
167 auto equilRegionInside = elemEquilRegion_[elementMapper_.index(inside)];
168 auto equilRegionOutside = elemEquilRegion_[elementMapper_.index(outside)];
170 if (thpres.hasRegionBarrier(equilRegionInside + 1, equilRegionOutside + 1)) {
172 if (thpres.hasThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1)) {
174 pth = thpres.getThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1);
178 unsigned offset = equilRegionInside*numEquilRegions_ + equilRegionOutside;
179 pth = thpresDefault_[offset];
182 unsigned offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside;
183 unsigned offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside;
185 thpres_[offset1] = pth;
186 thpres_[offset2] = pth;
192 if (thpres.ftSize() > 0)
193 configureThpresft_();
196template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
197void GenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
201 const FaultCollection& faults = eclState_.getFaults();
202 const SimulationConfig& simConfig = eclState_.getSimulationConfig();
203 const auto& thpres = simConfig.getThresholdPressure();
206 int numFaults = faults.size();
207 int numCartesianElem = eclState_.getInputGrid().getCartesianSize();
208 thpresftValues_.resize(numFaults, -1.0);
209 cartElemFaultIdx_.resize(numCartesianElem, -1);
210 for (std::size_t faultIdx = 0; faultIdx < faults.size(); faultIdx++) {
211 auto& fault = faults.getFault(faultIdx);
212 thpresftValues_[faultIdx] = thpres.getThresholdPressureFault(faultIdx);
213 for (
const FaultFace& face : fault)
216 for (std::size_t cartElemIdx : face)
217 cartElemFaultIdx_[cartElemIdx] = faultIdx;
221template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
226 if (!enableThresholdPressure_)
229 std::vector<Scalar> result(numEquilRegions_ * numEquilRegions_, 0.0);
230 const auto& simConfig = eclState_.getSimulationConfig();
231 const auto& thpres = simConfig.getThresholdPressure();
234 for (
unsigned j = 1; j <= numEquilRegions_; ++j) {
235 for (
unsigned i = 1; i <= numEquilRegions_; ++i, ++idx) {
236 if (thpres.hasRegionBarrier(i, j)) {
237 if (thpres.hasThresholdPressure(i, j)) {
238 result[idx] = thpres.getThresholdPressure(i, j);
240 result[idx] = this->thpresDefault_[idx];
263 " ", units.from_si(UnitSystem::measure::pressure, val),
284 str += lineFormat(i, j, thpres.getThresholdPressure(j, i));