23#ifndef OPM_GENERIC_THRESHOLD_PRESSURE_IMPL_HPP
24#define OPM_GENERIC_THRESHOLD_PRESSURE_IMPL_HPP
26#include <dune/grid/common/mcmgmapper.hh>
27#include <dune/grid/common/rangegenerators.hh>
29#include <opm/common/ErrorMacros.hpp>
31#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
32#include <opm/input/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
33#include <opm/input/eclipse/EclipseState/Tables/Eqldims.hpp>
34#include <opm/input/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
35#include <opm/input/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
39#include <fmt/format.h>
48template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
51 const GridView& gridView,
52 const ElementMapper& elementMapper,
53 const EclipseState& eclState)
54 : cartMapper_(cartMapper)
56 , elementMapper_(elementMapper)
57 , lookUpData_(gridView)
58 , lookUpCartesianData_(gridView, cartMapper_)
63template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
67 if (!enableThresholdPressure_)
71 if (!thpresftValues_.empty()) {
73 int fault1Idx = lookUpCartesianData_(elem1Idx, cartElemFaultIdx_);
74 int fault2Idx = lookUpCartesianData_(elem2Idx, cartElemFaultIdx_);
76 if (fault1Idx != -1 && fault1Idx == fault2Idx)
80 if (fault1Idx != fault2Idx) {
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>
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];
249template<
class Gr
id,
class Gr
idView,
class ElementMapper,
class Scalar>
254 if (!enableThresholdPressure_)
257 auto lineFormat = [
this](
unsigned i,
unsigned j,
double val)
259 const auto& units = eclState_.getUnits();
260 return fmt::format(
"{:4}{:>6}{:23}{:>6}{:24}{:>11.07}{:7}{}\n",
263 " ", units.from_si(UnitSystem::measure::pressure, val),
264 " ", units.name(UnitSystem::measure::pressure));
266 auto lineFormatS = [](
auto s1,
auto s2,
auto s3)
268 return fmt::format(
"{:4}{:^16}{:13}{:^9}{:21}{:^18}\n",
269 " ", s1,
" ", s2,
" ", s3);
272 std::string str =
"\nLIST OF ALL NON-ZERO THRESHOLD PRESSURES\n"
273 "----------------------------------------\n"
275 str += lineFormatS(
"FLOW FROM REGION",
"TO REGION",
"THRESHOLD PRESSURE");
276 str += lineFormatS(std::string(16,
'-'), std::string(9,
'-'), std::string(18,
'-'));
277 const auto& simConfig = eclState_.getSimulationConfig();
278 const auto& thpres = simConfig.getThresholdPressure();
280 for (
unsigned i = 1; i <= numEquilRegions_; ++i) {
281 for (
unsigned j = (thpres.irreversible() ? 1 : i); j <= numEquilRegions_; ++j) {
282 if (thpres.hasRegionBarrier(i, j)) {
283 if (thpres.hasThresholdPressure(i, j)) {
284 str += lineFormat(i, j, thpres.getThresholdPressure(j, i));
286 std::size_t idx = (j - 1) * numEquilRegions_ + (i - 1);
287 str += lineFormat(i, j, this->thpresDefault_[idx]);
292 str += lineFormatS(std::string(16,
'-'), std::string(9,
'-'), std::string(18,
'-'));
Definition: CollectDataOnIORank.hpp:49
void configureThpresft_()
Definition: GenericThresholdPressure_impl.hpp:198
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
GenericThresholdPressure(const CartesianIndexMapper &cartMapper, const GridView &gridView, const ElementMapper &elementMapper, const EclipseState &eclState)
Definition: GenericThresholdPressure_impl.hpp:50
Scalar thresholdPressure(int elem1Idx, int elem2Idx) const
Returns the theshold pressure [Pa] for the intersection between two elements.
Definition: GenericThresholdPressure_impl.hpp:65
std::vector< Scalar > getRestartVector() const
Returns a fully expanded vector for restart file writing.
Definition: GenericThresholdPressure_impl.hpp:224
void applyExplicitThresholdPressures_()
Definition: GenericThresholdPressure_impl.hpp:150
Definition: blackoilboundaryratevector.hh:37