From 988e3ef9077c0972dcdcde3ac5839769ebbfaaf7 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 11 Aug 2026 15:39:36 +0200 Subject: [PATCH] Add a fused relperm + capillary pressure entry point krw(), krn(), krg(), pcnw() and pcgn() each re-read the fluid state and re-derive their own two-phase saturation argument, and krn() calls two helpers that derive sw and sg a third time. A caller that needs relperms and capillary pressures together -- the intensive quantities do -- pays for that five times over. relativePermeabilitiesAndCapillaryPressures() forms the three distinct arguments once and calls the two-phase laws directly. Same arithmetic, same table lookups, same results; the new test checks the fused and separate paths agree exactly over a saturation sweep, float and double. Measured on Norne (44 927 cells, 3-derivative AD, standalone harness, identical states and checksums): 89.9 -> 59.8 ns/cell for relperm + pc. Co-Authored-By: Claude Fable 5 --- .../EclDefaultMaterial.hpp | 68 +++++++++++++++++++ tests/material/test_eclmateriallawmanager.cpp | 45 ++++++++++++ 2 files changed, 113 insertions(+) diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp index f27084213f9..7cd694ff020 100644 --- a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp @@ -358,6 +358,74 @@ class EclDefaultMaterial : public TraitsT values[gasPhaseIdx] = krg(params, fluidState); } + /*! + * \brief Relative permeabilities and capillary pressures in one pass. + * + * Produces exactly the same values as relativePermeabilities() followed + * by capillaryPressures(), but derives each two-phase saturation + * argument once instead of letting krw(), krn(), krg(), pcnw() and + * pcgn() each re-read the fluid state and re-derive it. Callers that + * need both -- the intensive quantities do -- should prefer this. + */ + template + OPM_HOST_DEVICE static void relativePermeabilitiesAndCapillaryPressures(ContainerT& kr, + ContainerT& pc, + const Params& params, + const FluidState& fluidState) + { + OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps); + using Evaluation = typename std::remove_reference::type; + + const Scalar Swco = params.Swl(); + const Evaluation swRaw = decay(fluidState.saturation(waterPhaseIdx)); + const Evaluation sg = decay(fluidState.saturation(gasPhaseIdx)); + + // the three distinct two-phase arguments, each formed once + const Evaluation sw = max(Evaluation(Swco), swRaw); + const Evaluation Sw_ow = sw + sg; + const Evaluation So_go = 1.0 - Sw_ow; + const Evaluation Sg_arg = 1.0 - Swco - sg; + + const Evaluation krw_ = OilWaterMaterialLaw::template twoPhaseSatKrw( + params.oilWaterParams(), swRaw); + const Evaluation krg_ = GasOilMaterialLaw::template twoPhaseSatKrn( + params.gasOilParams(), Sg_arg); + const Evaluation kro_ow = OilWaterMaterialLaw::template twoPhaseSatKrn( + params.oilWaterParams(), Sw_ow); + const Evaluation kro_go = GasOilMaterialLaw::template twoPhaseSatKrw( + params.gasOilParams(), So_go); + const Evaluation pcnw_ = OilWaterMaterialLaw::template twoPhaseSatPcnw( + params.oilWaterParams(), swRaw); + const Evaluation pcgn_ = GasOilMaterialLaw::template twoPhaseSatPcnw( + params.gasOilParams(), Sg_arg); + + // same regularization as krn() near Sw_ow == Swco + Evaluation kro; + constexpr const Scalar epsilon = 1e-5; + if (scalarValue(Sw_ow) - Swco < epsilon) { + const Evaluation kro2 = (kro_ow + kro_go) / 2; + if (scalarValue(Sw_ow) - Swco > epsilon / 2) { + const Evaluation kro1 = (sg * kro_go + (sw - Swco) * kro_ow) / (Sw_ow - Swco); + const Evaluation alpha = (epsilon - (Sw_ow - Swco)) / (epsilon / 2); + kro = kro2 * alpha + kro1 * (1 - alpha); + } + else { + kro = kro2; + } + } + else { + kro = (sg * kro_go + (sw - Swco) * kro_ow) / (Sw_ow - Swco); + } + + kr[waterPhaseIdx] = krw_; + kr[oilPhaseIdx] = kro; + kr[gasPhaseIdx] = krg_; + + pc[gasPhaseIdx] = pcgn_; + pc[oilPhaseIdx] = 0; + pc[waterPhaseIdx] = -pcnw_; + } + /*! * \brief The relative permeability of the gas phase. */ diff --git a/tests/material/test_eclmateriallawmanager.cpp b/tests/material/test_eclmateriallawmanager.cpp index 5980e3d0b7e..59fab1db6de 100644 --- a/tests/material/test_eclmateriallawmanager.cpp +++ b/tests/material/test_eclmateriallawmanager.cpp @@ -965,3 +965,48 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Let, Scalar, Types) } } } + +// The fused entry point must reproduce the two separate calls exactly: it is +// the same arithmetic with each two-phase saturation argument formed once. +BOOST_AUTO_TEST_CASE_TEMPLATE(FusedRelpermAndCapillaryPressureMatch, Scalar, Types) +{ + using MaterialLawManager = typename Fixture::MaterialLawManager; + using MaterialLaw = typename Fixture::MaterialLaw; + using DefaultMaterial = typename MaterialLaw::DefaultMaterial; + using FluidState = typename Fixture::FluidState; + constexpr int numPhases = Fixture::numPhases; + + Opm::Parser parser; + const auto deck = parser.parseString(fam1DeckString); + const Opm::EclipseState eclState(deck); + const std::size_t n = eclState.getInputGrid().getNumActive(); + + MaterialLawManager manager; + manager.initFromState(eclState); + manager.initParamsForElements(eclState, n, doOldLookup, doNothing); + + const auto& mp = manager.materialLawParams(0); + const auto& dp = mp.template getRealParams(); + + for (int a = 0; a <= 20; ++a) { + for (int b = 0; a + b <= 20; ++b) { + const Scalar sw = Scalar(a) / 20; + const Scalar sg = Scalar(b) / 20; + + FluidState fs; + fs.setSaturation(Fixture::waterPhaseIdx, sw); + fs.setSaturation(Fixture::gasPhaseIdx, sg); + fs.setSaturation(Fixture::oilPhaseIdx, 1.0 - sw - sg); + + std::array krSep{}, pcSep{}, krFused{}, pcFused{}; + DefaultMaterial::relativePermeabilities(krSep, dp, fs); + DefaultMaterial::capillaryPressures(pcSep, dp, fs); + DefaultMaterial::relativePermeabilitiesAndCapillaryPressures(krFused, pcFused, dp, fs); + + for (int ph = 0; ph < numPhases; ++ph) { + BOOST_CHECK_EQUAL(krSep[ph], krFused[ph]); + BOOST_CHECK_EQUAL(pcSep[ph], pcFused[ph]); + } + } + } +}