27 #ifndef OPM_ECL_EPS_SCALING_POINTS_HPP 28 #define OPM_ECL_EPS_SCALING_POINTS_HPP 33 #include <opm/parser/eclipse/Deck/Deck.hpp> 34 #include <opm/parser/eclipse/Deck/DeckRecord.hpp> 35 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp> 36 #include <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp> 37 #include <opm/parser/eclipse/EclipseState/Tables/SgfnTable.hpp> 38 #include <opm/parser/eclipse/EclipseState/Tables/SgofTable.hpp> 39 #include <opm/parser/eclipse/EclipseState/Tables/SlgofTable.hpp> 40 #include <opm/parser/eclipse/EclipseState/Tables/Sof3Table.hpp> 41 #include <opm/parser/eclipse/EclipseState/Tables/SwfnTable.hpp> 42 #include <opm/parser/eclipse/EclipseState/Tables/SwofTable.hpp> 43 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp> 63 typedef std::vector<int> IntData;
64 typedef std::vector<double> DoubleData;
68 void initFromDeck(
const Opm::Deck& ,
69 const Opm::EclipseState& eclState,
72 std::string kwPrefix = useImbibition?
"I":
"";
75 satnum = &eclState.get3DProperties().getIntGridProperty(
"IMBNUM").getData();
77 satnum = &eclState.get3DProperties().getIntGridProperty(
"SATNUM").getData();
79 retrieveGridPropertyData_(&swl, eclState, kwPrefix+
"SWL");
80 retrieveGridPropertyData_(&sgl, eclState, kwPrefix+
"SGL");
81 retrieveGridPropertyData_(&swcr, eclState, kwPrefix+
"SWCR");
82 retrieveGridPropertyData_(&sgcr, eclState, kwPrefix+
"SGCR");
83 retrieveGridPropertyData_(&sowcr, eclState, kwPrefix+
"SOWCR");
84 retrieveGridPropertyData_(&sogcr, eclState, kwPrefix+
"SOGCR");
85 retrieveGridPropertyData_(&swu, eclState, kwPrefix+
"SWU");
86 retrieveGridPropertyData_(&sgu, eclState, kwPrefix+
"SGU");
87 retrieveGridPropertyData_(&pcw, eclState, kwPrefix+
"PCW");
88 retrieveGridPropertyData_(&pcg, eclState, kwPrefix+
"PCG");
89 retrieveGridPropertyData_(&krw, eclState, kwPrefix+
"KRW");
90 retrieveGridPropertyData_(&kro, eclState, kwPrefix+
"KRO");
91 retrieveGridPropertyData_(&krg, eclState, kwPrefix+
"KRG");
94 const auto& ecl3dProps = eclState.get3DProperties();
95 poro = &ecl3dProps.getDoubleGridProperty(
"PORO").getData();
97 if (ecl3dProps.hasDeckDoubleGridProperty(
"PERMX")) {
98 permx = &ecl3dProps.getDoubleGridProperty(
"PERMX").getData();
103 if (ecl3dProps.hasDeckDoubleGridProperty(
"PERMY"))
104 permy = &ecl3dProps.getDoubleGridProperty(
"PERMY").getData();
106 if (ecl3dProps.hasDeckDoubleGridProperty(
"PERMZ"))
107 permz = &ecl3dProps.getDoubleGridProperty(
"PERMZ").getData();
111 const IntData* satnum;
113 const DoubleData* swl;
114 const DoubleData* sgl;
115 const DoubleData* swcr;
116 const DoubleData* sgcr;
117 const DoubleData* sowcr;
118 const DoubleData* sogcr;
119 const DoubleData* swu;
120 const DoubleData* sgu;
121 const DoubleData* pcw;
122 const DoubleData* pcg;
123 const DoubleData* krw;
124 const DoubleData* kro;
125 const DoubleData* krg;
126 const DoubleData* poro;
127 const DoubleData* permx;
128 const DoubleData* permy;
129 const DoubleData* permz;
135 void retrieveGridPropertyData_(
const DoubleData **data,
136 const Opm::EclipseState& eclState,
137 const std::string& properyName)
140 if (eclState.get3DProperties().hasDeckDoubleGridProperty(properyName))
141 (*data) = &eclState.get3DProperties().getDoubleGridProperty(properyName).getData();
153 template <
class Scalar>
180 Scalar pcowLeverettFactor;
181 Scalar pcgoLeverettFactor;
191 std::cout <<
" Swl: " << Swl <<
"\n" 192 <<
" Sgl: " << Sgl <<
"\n" 193 <<
" Sowl: " << Sowl <<
"\n" 194 <<
" Sogl: " << Sogl <<
"\n" 195 <<
" Swcr: " << Swcr <<
"\n" 196 <<
" Sgcr: " << Sgcr <<
"\n" 197 <<
" Sowcr: " << Sowcr <<
"\n" 198 <<
" Sogcr: " << Sogcr <<
"\n" 199 <<
" Swu: " << Swu <<
"\n" 200 <<
" Sgu: " << Sgu <<
"\n" 201 <<
" Sowu: " << Sowu <<
"\n" 202 <<
" Sogu: " << Sogu <<
"\n" 203 <<
" maxPcow: " << maxPcow <<
"\n" 204 <<
" maxPcgo: " << maxPcgo <<
"\n" 205 <<
" pcowLeverettFactor: " << pcowLeverettFactor <<
"\n" 206 <<
" pcgoLeverettFactor: " << pcgoLeverettFactor <<
"\n" 207 <<
" maxKrw: " << maxKrw <<
"\n" 208 <<
" maxKrg: " << maxKrg <<
"\n" 209 <<
" maxKrow: " << maxKrow <<
"\n" 210 <<
" maxKrog: " << maxKrog <<
"\n";
220 void extractUnscaled(
const Opm::Deck& deck,
221 const Opm::EclipseState& eclState,
222 unsigned satRegionIdx)
225 const auto& tables = eclState.getTableManager();
226 const TableContainer& swofTables = tables.getSwofTables();
227 const TableContainer& sgofTables = tables.getSgofTables();
228 const TableContainer& slgofTables = tables.getSlgofTables();
229 const TableContainer& swfnTables = tables.getSwfnTables();
230 const TableContainer& sgfnTables = tables.getSgfnTables();
231 const TableContainer& sof3Tables = tables.getSof3Tables();
233 bool hasWater = deck.hasKeyword(
"WATER");
234 bool hasGas = deck.hasKeyword(
"GAS");
235 bool hasOil = deck.hasKeyword(
"OIL");
241 if (!sgofTables.empty())
242 extractUnscaledSgof_(sgofTables.getTable<SgofTable>(satRegionIdx));
244 assert(!slgofTables.empty());
245 extractUnscaledSlgof_(slgofTables.getTable<SlgofTable>(satRegionIdx));
250 assert(!swofTables.empty());
254 extractUnscaledSwof_(swofTables.getTable<SwofTable>(satRegionIdx));
260 if (!hasWater || !hasGas || !hasOil)
261 throw std::domain_error(
"The specified phase configuration is not suppored");
263 bool family1 = (!sgofTables.empty() || !slgofTables.empty()) && !swofTables.empty();
264 bool family2 = !swfnTables.empty() && !sgfnTables.empty() && !sof3Tables.empty();
267 extractUnscaledSwof_(swofTables.getTable<SwofTable>(satRegionIdx));
269 if (!sgofTables.empty()) {
271 extractUnscaledSgof_(sgofTables.getTable<SgofTable>(satRegionIdx));
275 assert(!slgofTables.empty());
277 extractUnscaledSlgof_(slgofTables.getTable<SlgofTable>(satRegionIdx));
281 extractUnscaledSwfn_(swfnTables.getTable<SwfnTable>(satRegionIdx));
282 extractUnscaledSgfn_(sgfnTables.getTable<SgfnTable>(satRegionIdx));
283 extractUnscaledSof3_(sof3Tables.getTable<Sof3Table>(satRegionIdx));
286 throw std::domain_error(
"No valid saturation keyword family specified");
290 pcowLeverettFactor = 1.0;
291 pcgoLeverettFactor = 1.0;
299 void extractScaled(
const Opm::EclipseState& eclState,
301 unsigned cartesianCellIdx)
305 extractGridPropertyValue_(Swl, epsProperties.swl, cartesianCellIdx);
306 extractGridPropertyValue_(Sgl, epsProperties.sgl, cartesianCellIdx);
307 extractGridPropertyValue_(Swcr, epsProperties.swcr, cartesianCellIdx);
308 extractGridPropertyValue_(Sgcr, epsProperties.sgcr, cartesianCellIdx);
309 extractGridPropertyValue_(Sowcr, epsProperties.sowcr, cartesianCellIdx);
310 extractGridPropertyValue_(Sogcr, epsProperties.sogcr, cartesianCellIdx);
311 extractGridPropertyValue_(Swu, epsProperties.swu, cartesianCellIdx);
312 extractGridPropertyValue_(Sgu, epsProperties.sgu, cartesianCellIdx);
314 extractGridPropertyValue_(maxPcow, epsProperties.pcw, cartesianCellIdx);
315 extractGridPropertyValue_(maxPcgo, epsProperties.pcg, cartesianCellIdx);
317 extractGridPropertyValue_(maxKrw, epsProperties.krw, cartesianCellIdx);
318 extractGridPropertyValue_(maxKrg, epsProperties.krg, cartesianCellIdx);
321 extractGridPropertyValue_(maxKrow, epsProperties.kro, cartesianCellIdx);
322 extractGridPropertyValue_(maxKrog, epsProperties.kro, cartesianCellIdx);
327 pcowLeverettFactor = 1.0;
328 pcgoLeverettFactor = 1.0;
329 if (eclState.getTableManager().useJFunc()) {
330 const auto& jfunc = eclState.getTableManager().getJFunc();
331 const auto& jfuncDir = jfunc.direction();
334 if (jfuncDir == Opm::JFunc::Direction::X)
336 (*epsProperties.permx)[cartesianCellIdx];
337 else if (jfuncDir == Opm::JFunc::Direction::Y)
339 (*epsProperties.permy)[cartesianCellIdx];
340 else if (jfuncDir == Opm::JFunc::Direction::Z)
342 (*epsProperties.permz)[cartesianCellIdx];
343 else if (jfuncDir == Opm::JFunc::Direction::XY)
350 (*epsProperties.permy)[cartesianCellIdx]);
352 OPM_THROW(std::runtime_error,
"Illegal direction indicator for the JFUNC " 353 "keyword ("<<static_cast<int>(jfuncDir)<<
")");
358 Scalar poro = (*epsProperties.poro)[cartesianCellIdx];
359 Scalar alpha = jfunc.alphaFactor();
360 Scalar beta = jfunc.betaFactor();
364 Scalar commonFactor = std::pow(poro, alpha)/std::pow(perm, beta);
368 const Scalar Uconst = 0.318316 * 1e5;
371 const auto& jfuncFlag = jfunc.flag();
372 if (jfuncFlag == Opm::JFunc::Flag::WATER || jfuncFlag == Opm::JFunc::Flag::BOTH) {
375 jfunc.owSurfaceTension();
376 pcowLeverettFactor = commonFactor*gamma*Uconst;
380 if (jfuncFlag == Opm::JFunc::Flag::GAS || jfuncFlag == Opm::JFunc::Flag::BOTH) {
383 jfunc.goSurfaceTension();
384 pcgoLeverettFactor = commonFactor*gamma*Uconst;
392 void extractUnscaledSgof_(
const Opm::SgofTable& sgofTable)
395 Sgl = sgofTable.getSgColumn().front();
396 Sogl = 1.0 - sgofTable.getSgColumn().back();
399 Sgu = sgofTable.getSgColumn().back();
400 Sogu = 1.0 - sgofTable.getSgColumn().front();
404 for (
size_t rowIdx = 0; rowIdx < sgofTable.numRows(); ++ rowIdx) {
405 if (sgofTable.getKrgColumn()[rowIdx] > 0.0)
408 Sgcr = sgofTable.getSgColumn()[rowIdx];
413 for (
int rowIdx = static_cast<int>(sgofTable.numRows() - 1);
417 if (sgofTable.getKrogColumn()[
static_cast<size_t>(rowIdx)] > 0.0)
420 Sogcr = 1.0 - sgofTable.getSgColumn()[
static_cast<size_t>(rowIdx)];
424 maxPcgo = sgofTable.getPcogColumn().back();
427 maxKrg = sgofTable.getKrgColumn().back();
428 maxKrog = sgofTable.getKrogColumn().front();
431 void extractUnscaledSlgof_(
const Opm::SlgofTable& slgofTable)
434 Sgl = 1.0 - slgofTable.getSlColumn().back();
435 Sogl = slgofTable.getSlColumn().front();
438 Sgu = 1.0 - slgofTable.getSlColumn().front();
439 Sogu = slgofTable.getSlColumn().back();
443 for (
int rowIdx = static_cast<int>(slgofTable.numRows()) - 1;
447 if (slgofTable.getKrgColumn()[
static_cast<size_t>(rowIdx)] > 0.0)
450 Sgcr = 1 - slgofTable.getSlColumn()[
static_cast<size_t>(rowIdx)];
455 for (
size_t rowIdx = 0; rowIdx < slgofTable.numRows(); ++ rowIdx) {
456 if (slgofTable.getKrogColumn()[rowIdx] > 0.0)
459 Sogcr = slgofTable.getSlColumn()[rowIdx];
463 maxPcgo = slgofTable.getPcogColumn().front();
466 maxKrg = slgofTable.getKrgColumn().front();
467 maxKrog = slgofTable.getKrogColumn().back();
470 void extractUnscaledSwof_(
const Opm::SwofTable& swofTable)
473 Swl = swofTable.getSwColumn().front();
474 Sowl = 1.0 - swofTable.getSwColumn().back();
477 Swu = swofTable.getSwColumn().back();
478 Sowu = 1.0 - swofTable.getSwColumn().front();
482 for (
size_t rowIdx = 0; rowIdx < swofTable.numRows(); ++ rowIdx) {
483 if (swofTable.getKrwColumn()[rowIdx] > 0.0)
486 Swcr = swofTable.getSwColumn()[rowIdx];
491 for (
int rowIdx = static_cast<int>(swofTable.numRows()) - 1;
495 if (swofTable.getKrowColumn()[
static_cast<size_t>(rowIdx)] > 0.0)
498 Sowcr = 1.0 - swofTable.getSwColumn()[
static_cast<size_t>(rowIdx)];
502 maxPcow = swofTable.getPcowColumn().front();
505 maxKrw = swofTable.getKrwColumn().back();
506 maxKrow = swofTable.getKrowColumn().front();
509 void extractUnscaledSwfn_(
const Opm::SwfnTable& swfnTable)
512 Swl = swfnTable.getSwColumn().front();
515 Swu = swfnTable.getSwColumn().back();
519 for (
size_t rowIdx = 0; rowIdx < swfnTable.numRows(); ++ rowIdx) {
520 if (swfnTable.getKrwColumn()[rowIdx] > 0)
523 Swcr = swfnTable.getSwColumn()[rowIdx];
527 maxPcow = swfnTable.getPcowColumn().front();
530 maxKrw = swfnTable.getKrwColumn().back();
533 void extractUnscaledSgfn_(
const Opm::SgfnTable& sgfnTable)
536 Sgl = sgfnTable.getSgColumn().front();
539 Sgu = sgfnTable.getSgColumn().back();
540 Sogu = 1 - sgfnTable.getSgColumn().front();
544 for (
size_t rowIdx = 0; rowIdx < sgfnTable.numRows(); ++ rowIdx) {
545 if (sgfnTable.getKrgColumn()[rowIdx] > 0)
548 Sgcr = sgfnTable.getSgColumn()[rowIdx];
552 maxPcgo = sgfnTable.getPcogColumn().back();
555 maxKrg = sgfnTable.getKrgColumn().back();
558 void extractUnscaledSof3_(
const Opm::Sof3Table& sof3Table)
561 Sowl = sof3Table.getSoColumn().front() + Sgl;
562 Sogl = sof3Table.getSoColumn().front() + Swl;
565 Sowu = sof3Table.getSoColumn().back();
569 for (
size_t rowIdx = 0 ; rowIdx < sof3Table.numRows(); ++ rowIdx) {
570 if (sof3Table.getKrowColumn()[rowIdx] > 0) {
574 Sowcr = sof3Table.getSoColumn()[rowIdx];
579 for (
size_t rowIdx = 0 ; rowIdx < sof3Table.numRows(); ++ rowIdx) {
580 if (sof3Table.getKrogColumn()[rowIdx] > 0)
583 Sogcr = sof3Table.getSoColumn()[rowIdx];
587 maxKrow = sof3Table.getKrowColumn().back();
588 maxKrog = sof3Table.getKrogColumn().back();
590 #endif // HAVE_OPM_PARSER 592 void extractGridPropertyValue_(Scalar& targetValue,
593 const std::vector<double>* propData,
594 unsigned cartesianCellIdx)
599 targetValue = (*propData)[cartesianCellIdx];
609 template <
class Scalar>
620 if (epsSystemType == EclOilWaterSystem) {
622 saturationPcPoints_[0] = epsInfo.Swl;
623 saturationPcPoints_[1] = epsInfo.Swu;
627 saturationKrwPoints_[0] = epsInfo.Swcr;
628 saturationKrwPoints_[1] = 1.0 - epsInfo.Sowcr - epsInfo.Sgl;
629 saturationKrwPoints_[2] = epsInfo.Swu;
632 saturationKrwPoints_[0] = epsInfo.Swcr;
633 saturationKrwPoints_[1] = epsInfo.Swu;
641 saturationKrnPoints_[2] = 1.0 - epsInfo.Sowcr;
642 saturationKrnPoints_[1] = epsInfo.Swcr + epsInfo.Sgl;
643 saturationKrnPoints_[0] = epsInfo.Swl + epsInfo.Sgl;
646 saturationKrnPoints_[1] = 1 - epsInfo.Sowcr;
647 saturationKrnPoints_[0] = epsInfo.Swl + epsInfo.Sgl;
651 maxPcnwOrLeverettFactor_ = epsInfo.pcowLeverettFactor;
653 maxPcnwOrLeverettFactor_ = epsInfo.maxPcow;
654 maxKrw_ = epsInfo.maxKrw;
655 maxKrn_ = epsInfo.maxKrow;
658 assert(epsSystemType == EclGasOilSystem);
661 saturationPcPoints_[0] = 1.0 - epsInfo.Sgu;
662 saturationPcPoints_[1] = 1.0 - epsInfo.Sgl;
666 saturationKrwPoints_[0] = epsInfo.Sogcr;
667 saturationKrwPoints_[1] = 1 - epsInfo.Sgcr - epsInfo.Swl;
668 saturationKrwPoints_[2] = 1 - epsInfo.Swl - epsInfo.Sgl;
671 saturationKrwPoints_[0] = epsInfo.Sogcr;
672 saturationKrwPoints_[1] = 1 - epsInfo.Swl - epsInfo.Sgl;
680 saturationKrnPoints_[2] = 1.0 - epsInfo.Sgcr;
681 saturationKrnPoints_[1] = epsInfo.Sogcr + epsInfo.Swl;
682 saturationKrnPoints_[0] = 1.0 - epsInfo.Sgu;
685 saturationKrnPoints_[1] = 1.0 - epsInfo.Sgcr;
686 saturationKrnPoints_[0] = 1.0 - epsInfo.Sgu;
690 maxPcnwOrLeverettFactor_ = epsInfo.pcgoLeverettFactor;
692 maxPcnwOrLeverettFactor_ = epsInfo.maxPcgo;
694 maxKrw_ = epsInfo.maxKrog;
695 maxKrn_ = epsInfo.maxKrg;
703 { saturationPcPoints_[pointIdx] = value; }
709 {
return saturationPcPoints_; }
715 { saturationKrwPoints_[pointIdx] = value; }
721 {
return saturationKrwPoints_; }
727 { saturationKrnPoints_[pointIdx] = value; }
733 {
return saturationKrnPoints_; }
739 { maxPcnwOrLeverettFactor_ = value; }
745 {
return maxPcnwOrLeverettFactor_; }
751 { maxPcnwOrLeverettFactor_ = value; }
757 {
return maxPcnwOrLeverettFactor_; }
785 std::cout <<
" saturationKrnPoints_[0]: " << saturationKrnPoints_[0] <<
"\n" 786 <<
" saturationKrnPoints_[1]: " << saturationKrnPoints_[1] <<
"\n" 787 <<
" saturationKrnPoints_[2]: " << saturationKrnPoints_[2] <<
"\n";
792 Scalar maxPcnwOrLeverettFactor_;
801 std::array<Scalar, 2> saturationPcPoints_;
804 std::array<Scalar, 3> saturationKrwPoints_;
807 std::array<Scalar, 3> saturationKrnPoints_;
Implements some common averages.
void setMaxKrw(Scalar value)
Sets the maximum wetting phase relative permeability.
Definition: EclEpsScalingPoints.hpp:762
void setSaturationPcPoint(unsigned pointIdx, Scalar value)
Sets an saturation value for capillary pressure saturation scaling.
Definition: EclEpsScalingPoints.hpp:702
void setMaxPcnw(Scalar value)
Sets the maximum capillary pressure.
Definition: EclEpsScalingPoints.hpp:738
Scalar maxKrw() const
Returns the maximum wetting phase relative permeability.
Definition: EclEpsScalingPoints.hpp:768
void setSaturationKrnPoint(unsigned pointIdx, Scalar value)
Sets an saturation value for non-wetting phase relperm saturation scaling.
Definition: EclEpsScalingPoints.hpp:726
Collects all grid properties which are relevant for end point scaling.
Definition: EclEpsScalingPoints.hpp:61
void setSaturationKrwPoint(unsigned pointIdx, Scalar value)
Sets an saturation value for wetting-phase relperm saturation scaling.
Definition: EclEpsScalingPoints.hpp:714
Definition: Air_Mesitylene.hpp:33
Specifies the configuration used by the endpoint scaling code.
Definition: EclEpsConfig.hpp:63
Scalar leverettFactor() const
Returns the Leverett scaling factor for capillary pressure.
Definition: EclEpsScalingPoints.hpp:756
This structure represents all values which can be possibly used as scaling points in the endpoint sca...
Definition: EclEpsScalingPoints.hpp:154
Scalar maxKrn() const
Returns the maximum wetting phase relative permeability.
Definition: EclEpsScalingPoints.hpp:780
Represents the points on the X and Y axis to be scaled if endpoint scaling is used.
Definition: EclEpsScalingPoints.hpp:610
void setLeverettFactor(Scalar value)
Sets the Leverett scaling factor for capillary pressure.
Definition: EclEpsScalingPoints.hpp:750
Scalar arithmeticMean(Scalar x, Scalar y)
Computes the arithmetic average of two values.
Definition: Means.hpp:45
bool enableLeverettScaling() const
Returns whether the Leverett capillary pressure scaling is enabled.
Definition: EclEpsConfig.hpp:156
const std::array< Scalar, 3 > & saturationKrwPoints() const
Returns the points used for wetting phase relperm saturation scaling.
Definition: EclEpsScalingPoints.hpp:720
Specifies the configuration used by the endpoint scaling code.
void setMaxKrn(Scalar value)
Sets the maximum wetting phase relative permeability.
Definition: EclEpsScalingPoints.hpp:774
Scalar maxPcnw() const
Returns the maximum capillary pressure.
Definition: EclEpsScalingPoints.hpp:744
bool enableThreePointKrSatScaling() const
Returns whether three point saturation scaling is enabled for the relative permeabilities.
Definition: EclEpsConfig.hpp:100
const std::array< Scalar, 3 > & saturationKrnPoints() const
Returns the points used for non-wetting phase relperm saturation scaling.
Definition: EclEpsScalingPoints.hpp:732
EclTwoPhaseSystemType
Specified which fluids are involved in a given twophase material law for endpoint scaling...
Definition: EclEpsConfig.hpp:51
const std::array< Scalar, 2 > & saturationPcPoints() const
Returns the points used for capillary pressure saturation scaling.
Definition: EclEpsScalingPoints.hpp:708
void init(const EclEpsScalingPointsInfo< Scalar > &epsInfo, const EclEpsConfig &config, EclTwoPhaseSystemType epsSystemType)
Assigns the scaling points which actually ought to be used.
Definition: EclEpsScalingPoints.hpp:616