From 8b1bcb36bce8d301a866e501fd7069fda0688602 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Fri, 1 Jul 2022 08:19:51 +0200 Subject: [PATCH 1/6] Adds support for KRNUM --- ebos/eclfluxmodule.hh | 12 ++++++++---- ebos/eclgenericproblem.cc | 9 +++++++++ ebos/eclgenericproblem.hh | 4 ++++ ebos/eclproblem.hh | 2 ++ .../utils/PartiallySupportedFlowKeywords.cpp | 1 - 5 files changed, 23 insertions(+), 5 deletions(-) diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index 19557e97a..943e22874 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -36,9 +36,12 @@ #include #include +#include +#include #include #include +#include namespace Opm { @@ -227,7 +230,7 @@ public: unsigned I = stencil.globalSpaceIndex(interiorDofIdx); unsigned J = stencil.globalSpaceIndex(exteriorDofIdx); - + auto facedir = scvf.dirId(); // direction (X, Y, or Z) of the face Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx); Scalar faceArea = scvf.area(); Scalar thpres = problem.thresholdPressure(I, J); @@ -265,7 +268,7 @@ public: intQuantsEx, phaseIdx,//input interiorDofIdx,//input - exteriorDofIdx,//intput + exteriorDofIdx,//input Vin, Vex, I, @@ -285,10 +288,11 @@ public: if (upwindIsInterior) volumeFlux[phaseIdx] = - pressureDifferences[phaseIdx]*up.mobility(phaseIdx)*transMult*(-trans/faceArea); + pressureDifferences[phaseIdx]*up.mobility(phaseIdx, facedir)*transMult*(-trans/faceArea); else volumeFlux[phaseIdx] = - pressureDifferences[phaseIdx]*(Toolbox::value(up.mobility(phaseIdx))*Toolbox::value(transMult)*(-trans/faceArea)); + pressureDifferences[phaseIdx]* + (Toolbox::value(up.mobility(phaseIdx, facedir))*Toolbox::value(transMult)*(-trans/faceArea)); } } diff --git a/ebos/eclgenericproblem.cc b/ebos/eclgenericproblem.cc index 4aa2709b6..c014208ab 100644 --- a/ebos/eclgenericproblem.cc +++ b/ebos/eclgenericproblem.cc @@ -322,6 +322,15 @@ updatePlmixnum_() updateNum("PLMIXNUM", plmixnum_); } +template +void EclGenericProblem:: +updateKrnum_() +{ + updateNum("KRNUMX", krnumx_); + updateNum("KRNUMY", krnumy_); + updateNum("KRNUMZ", krnumz_); +} + template bool EclGenericProblem:: vapparsActive(int episodeIdx) const diff --git a/ebos/eclgenericproblem.hh b/ebos/eclgenericproblem.hh index bfbdcfe9e..290c8ae89 100644 --- a/ebos/eclgenericproblem.hh +++ b/ebos/eclgenericproblem.hh @@ -306,6 +306,7 @@ protected: void updateSatnum_(); void updateMiscnum_(); void updatePlmixnum_(); + void updateKrnum_(); const EclipseState& eclState_; const Schedule& schedule_; @@ -318,6 +319,9 @@ protected: std::vector satnum_; std::vector miscnum_; std::vector plmixnum_; + std::vector krnumx_; + std::vector krnumy_; + std::vector krnumz_; std::vector rockParams_; std::vector rockTableIdx_; diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 1400f97d3..6e7b4439d 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -2223,6 +2223,8 @@ private: // the PLMIX region numbers (polymer model) this->updatePlmixnum_(); + // directional relative permeabilities + this->updateKrnum_(); //////////////////////////////// // porosity updateReferencePorosity_(); diff --git a/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp b/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp index f369266eb..bd54deaa0 100644 --- a/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp +++ b/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp @@ -196,7 +196,6 @@ partiallySupported() { "SATOPTS", { - {1,{false, allow_values {"HYSTER"}, "SATOPTS(DIRECT): directional relative permeability assignment option not supported - value ignored"}}, // options {2,{false, allow_values {"HYSTER"}, "SATOPTS(IRREVERS): reversible directional relative permeability option not supported – value ignored"}}, // IRREVERS {3,{false, allow_values {"HYSTER"}, "SATOPTS(HYSTER): hysteresis directional relative permeability option not supported - value ignored"}}, // HYSTER {4,{false, allow_values {"HYSTER"}, "SATOPTS(SURFTENS): capillary pressure surface tension pressure dependency option not supported – value ignored"}}, // SURFTENS From 7d0265dd44eb7fbb8d496ed1c3d9eb16038628c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Tue, 6 Sep 2022 19:27:20 +0200 Subject: [PATCH 2/6] Limit access to face direction. Only access the grid face direction information if directional relperms have been enabled. --- ebos/eclfluxmodule.hh | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index 943e22874..b0ee0f32f 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -230,7 +230,6 @@ public: unsigned I = stencil.globalSpaceIndex(interiorDofIdx); unsigned J = stencil.globalSpaceIndex(exteriorDofIdx); - auto facedir = scvf.dirId(); // direction (X, Y, or Z) of the face Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx); Scalar faceArea = scvf.area(); Scalar thpres = problem.thresholdPressure(I, J); @@ -286,13 +285,22 @@ public: // or averaged? all fluids should see the same compaction?! const Evaluation& transMult = up.rockCompTransMultiplier(); + const auto& materialLawManager = problem.materialLawManager(); + const Evaluation *mob = nullptr; + if (materialLawManager->hasDirectionalRelperms()) { + auto facedir = scvf.dirId(); // direction (X, Y, or Z) of the face + mob = &(up.mobility(phaseIdx, facedir)); + } + else { + mob = &(up.mobility(phaseIdx)); + } if (upwindIsInterior) volumeFlux[phaseIdx] = - pressureDifferences[phaseIdx]*up.mobility(phaseIdx, facedir)*transMult*(-trans/faceArea); + pressureDifferences[phaseIdx]*(*mob)*transMult*(-trans/faceArea); else volumeFlux[phaseIdx] = pressureDifferences[phaseIdx]* - (Toolbox::value(up.mobility(phaseIdx, facedir))*Toolbox::value(transMult)*(-trans/faceArea)); + (Toolbox::value(*mob)*Toolbox::value(transMult)*(-trans/faceArea)); } } From 9ed974fda39b59252ac87bfe0f7109878dd86487 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Wed, 14 Sep 2022 15:33:59 +0200 Subject: [PATCH 3/6] Minor simplification Simplify code using the unknown facedir value. --- ebos/eclfluxmodule.hh | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index b0ee0f32f..e21b39245 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -286,21 +286,17 @@ public: const Evaluation& transMult = up.rockCompTransMultiplier(); const auto& materialLawManager = problem.materialLawManager(); - const Evaluation *mob = nullptr; + FaceDir::DirEnum facedir = FaceDir::DirEnum::Unknown; if (materialLawManager->hasDirectionalRelperms()) { - auto facedir = scvf.dirId(); // direction (X, Y, or Z) of the face - mob = &(up.mobility(phaseIdx, facedir)); - } - else { - mob = &(up.mobility(phaseIdx)); + facedir = scvf.dirId(); // direction (X, Y, or Z) of the face } if (upwindIsInterior) volumeFlux[phaseIdx] = - pressureDifferences[phaseIdx]*(*mob)*transMult*(-trans/faceArea); + pressureDifferences[phaseIdx]*up.mobility(phaseIdx, facedir)*transMult*(-trans/faceArea); else volumeFlux[phaseIdx] = pressureDifferences[phaseIdx]* - (Toolbox::value(*mob)*Toolbox::value(transMult)*(-trans/faceArea)); + (Toolbox::value(up.mobility(phaseIdx, facedir))*Toolbox::value(transMult)*(-trans/faceArea)); } } From cbdec23e50dc3655818331c4fa1cfa65d2a2c608 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Fri, 16 Sep 2022 09:07:25 +0200 Subject: [PATCH 4/6] Also check for unknown face direction To avoid compiler warnings, we also need to check for the unknown face direction due to change in the FaceDir::DirEnum values in opm-common. --- ebos/eclproblem.hh | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 6e7b4439d..30b32456e 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -2851,6 +2851,8 @@ private: case FaceDir::ZPlus: data = &massratebcZ_; break; + case FaceDir::Unknown: + throw std::runtime_error("Unexpected unknown face direction"); } const Evaluation rate = bcface.rate; @@ -2885,6 +2887,8 @@ private: case FaceDir::ZPlus: data = &freebcZ_; break; + case FaceDir::Unknown: + throw std::runtime_error("Unexpected unknown face direction"); } for (int i = bcface.i1; i <= bcface.i2; ++i) { From c6129ad9a41194d8cc710a05270445fb9da637e2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Fri, 16 Sep 2022 09:10:42 +0200 Subject: [PATCH 5/6] Renamed function in ecfvstencil.hh Use the new name of the dirId() function in ecfvstencil.hh --- ebos/eclfluxmodule.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index e21b39245..138673586 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -288,7 +288,7 @@ public: const auto& materialLawManager = problem.materialLawManager(); FaceDir::DirEnum facedir = FaceDir::DirEnum::Unknown; if (materialLawManager->hasDirectionalRelperms()) { - facedir = scvf.dirId(); // direction (X, Y, or Z) of the face + facedir = scvf.faceDirFromDirId(); // direction (X, Y, or Z) of the face } if (upwindIsInterior) volumeFlux[phaseIdx] = From affc5cb53b27aaccb05192d9629b59ba1587344b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Thu, 22 Sep 2022 23:54:29 +0200 Subject: [PATCH 6/6] Return a pointer to the material manager This is used by blackoilintensivequantities.hh. Further, for other problems derived from MultiPhaseBaseProblem, that base class will return a null pointer such that problems that do not override the materialLawManagerPtr() method still can use the blackoil intensive quantitites (which in the case of directional relative permeabilities makes use of the problem reference to access the materialLawManager) --- ebos/eclproblem.hh | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 30b32456e..aa8645e86 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -1505,6 +1505,12 @@ public: std::shared_ptr materialLawManager() const { return materialLawManager_; } + // TODO: See discussion in multiphasebaseproblem.hh for the reason why we need this method + const EclMaterialLawManager* materialLawManagerPtr() const + { + return materialLawManager_.get(); + } + /*! * \copydoc materialLawManager() */