From d90107841c9aac729421467c8990a220528b6b87 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 1 Oct 2026 00:10:25 +0200 Subject: [PATCH 01/10] Accept a COMPVD phase change away from the type-1 gas-oil contact With EQUIL item 10 = 1 the COMPVD rows give the composition versus depth, so they decide where the gas zone ends. A contact outside the gap between the last vapour row and the first liquid row now moves to the nearest edge of that gap with a warning instead of stopping the run. Type 3 keeps the error, as it uses the contact as reference depth. --- .../flow/equil/InitStateEquilComp.hpp | 59 +++++++++++---- tests/test_compequil.cpp | 71 +++++++++++++++++++ 2 files changed, 116 insertions(+), 14 deletions(-) diff --git a/opm/simulators/flow/equil/InitStateEquilComp.hpp b/opm/simulators/flow/equil/InitStateEquilComp.hpp index d5a250e3262..1dff1b80133 100644 --- a/opm/simulators/flow/equil/InitStateEquilComp.hpp +++ b/opm/simulators/flow/equil/InitStateEquilComp.hpp @@ -201,7 +201,9 @@ class WaterDensityODE * * The COMPVD phase column selects the EOS root. If its rows name both phases, * the vapour rows define the gas zone above the gas-oil contact and the liquid - * rows define the liquid zone below it. + * rows define the liquid zone below it. For type 1 a contact outside the gap + * between the last vapour row and the first liquid row moves to the nearest + * edge of that gap. * * Only cell-centre initialization is supported (EQUIL item 9 = 0). * Gas-oil contact capillary pressure must be zero: the downstream flash uses @@ -325,6 +327,7 @@ class InitialStateComputer /// The equilibrated vertical distributions within one region. struct Region { int initType{1}; // EQUIL item 10 + /// Gas-oil contact. A two-zone COMPVD table may move it, see zoneBoundary(). Scalar zgoc{}; /// Name of the selected composition keyword, used in diagnostics. std::string_view compositionKeyword; @@ -455,13 +458,17 @@ class InitialStateComputer return rows; } - /// Verify that vapour rows are at or above the gas-oil contact and liquid - /// rows are at or below it. - static void checkZonesStraddleContact(const CompvdTable& compvd, - const std::vector& vaporRows, - const std::vector& liquidRows, - const Scalar zgoc, - const std::size_t regionIdx) + /// Depth at which the gas zone of a COMPVD table naming both phases meets + /// its liquid zone: the gas-oil contact when it lies between the last + /// vapour row and the first liquid row. With EQUIL item 10 = 1 the rows + /// give the composition versus depth, so a contact outside that gap moves + /// to its nearest edge and every row keeps the phase it names. Type 3 uses + /// the contact as its reference depth, so there the rows must agree with it. + static Scalar zoneBoundary(const CompvdTable& compvd, + const std::vector& vaporRows, + const std::vector& liquidRows, + const Region& reg, + const std::size_t regionIdx) { if (vaporRows.empty() || liquidRows.empty()) { OPM_THROW(std::runtime_error, @@ -470,14 +477,39 @@ class InitialStateComputer } const auto& depth = compvd.getDepthColumn(); - if ((depth[vaporRows.back()] > zgoc) || (depth[liquidRows.front()] < zgoc)) { + const Scalar lastVapor = depth[vaporRows.back()]; + const Scalar firstLiquid = depth[liquidRows.front()]; + if (lastVapor > firstLiquid) { + OPM_THROW(std::runtime_error, + fmt::format("The COMPVD table of region {} has a vapour row at {} m " + "below its liquid row at {} m. All vapour rows must lie " + "above the liquid rows.", + regionIdx + 1, lastVapor, firstLiquid)); + } + + if ((reg.zgoc >= lastVapor) && (reg.zgoc <= firstLiquid)) { + return reg.zgoc; + } + + if (reg.initType != 1) { OPM_THROW(std::runtime_error, fmt::format("The COMPVD table of region {} puts its vapour rows down to " "{} m and its liquid rows from {} m, which do not meet at " - "the gas-oil contact at {} m.", - regionIdx + 1, depth[vaporRows.back()], - depth[liquidRows.front()], zgoc)); + "the gas-oil contact at {} m. EQUIL item 10 = {} takes the " + "contact as the reference depth, so the rows must agree " + "with it.", + regionIdx + 1, lastVapor, firstLiquid, reg.zgoc, + reg.initType)); } + + const Scalar boundary = std::clamp(reg.zgoc, lastVapor, firstLiquid); + OpmLog::warning(fmt::format("Equilibration region {}: the gas-oil contact at {} m lies " + "outside the COMPVD phase change between the last vapour " + "row at {} m and the first liquid row at {} m. With EQUIL " + "item 10 = 1 the rows give the composition versus depth, " + "so the gas zone ends at {} m instead.", + regionIdx + 1, reg.zgoc, lastVapor, firstLiquid, boundary)); + return boundary; } /// The COMPVD rows carrying \p phase. @@ -645,8 +677,7 @@ class InitialStateComputer // the liquid rows the one below the contact. const auto vaporRows = rowsOfPhase(compvd, CompvdTable::Phase::Vapor); const auto liquidRows = rowsOfPhase(compvd, CompvdTable::Phase::Liquid); - checkZonesStraddleContact(compvd, vaporRows, liquidRows, - record.gasOilContactDepth(), regionIdx); + reg.zgoc = zoneBoundary(compvd, vaporRows, liquidRows, reg, regionIdx); reg.twoZone = true; setupComposition(reg.vaporVdTable, compvd, vaporRows); setupComposition(reg.compositionVdTable, compvd, liquidRows); diff --git a/tests/test_compequil.cpp b/tests/test_compequil.cpp index b8d3ac6e788..8fdb02d7388 100644 --- a/tests/test_compequil.cpp +++ b/tests/test_compequil.cpp @@ -45,6 +45,8 @@ #include #include #include +#include +#include #include namespace { @@ -848,6 +850,75 @@ BOOST_AUTO_TEST_CASE(CompvdTwoZoneClampsCompositionOutsideRows) } } +BOOST_AUTO_TEST_CASE(CompvdTwoZoneContactOutsideRowGap) +{ + // With EQUIL item 10 = 1 the COMPVD rows give the composition versus depth, + // so a gas-oil contact outside the gap between the last vapour row (2049 m) + // and the first liquid row (2051 m) moves to the nearest edge of the gap. + const auto states = [](const std::string& goc) { + const EquilFixture fix(deckString( + "EQUIL\n 2072.5 200 2300 0 " + goc + " 0 /\n", "EQLDIMS\n/\n", "", + "COMPVD\n" + " 2000 0 0.95 0.05 0 150.0\n" + " 2049 0 0.95 0.05 0 150.0\n" + " 2051 0 0.60 0.40 1 150.0\n" + " 2100 0 0.40 0.60 1 150.0 /\n")); + return fix.compute(std::vector(20, 0)).fluidStates(); + }; + + // The cell at 2057.5 m lies above a contact at 2060 m but below the first + // liquid row, and the one at 2042.5 m below a contact at 2030 m but above + // the last vapour row. + for (const auto& [goc, edge, cell, phaseIdx] : + {std::tuple{"2060", "2051", std::size_t{11}, FluidSystem::oilPhaseIdx}, + std::tuple{"2030", "2049", std::size_t{8}, FluidSystem::gasPhaseIdx}}) { + BOOST_TEST_CONTEXT("Contact at " << goc << " m") { + const auto moved = states(goc); + BOOST_CHECK_CLOSE(moved[cell].saturation(phaseIdx), 1.0, 1e-10); + + const auto expected = states(edge); + BOOST_REQUIRE_EQUAL(moved.size(), expected.size()); + for (std::size_t c = 0; c < moved.size(); ++c) { + BOOST_TEST_CONTEXT("Cell " << c) { + for (const auto p : {FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}) { + BOOST_CHECK_EQUAL(moved[c].pressure(p), expected[c].pressure(p)); + BOOST_CHECK_EQUAL(moved[c].saturation(p), expected[c].saturation(p)); + } + for (int comp = 0; comp < 3; ++comp) { + BOOST_CHECK_EQUAL(moved[c].moleFraction(comp), + expected[c].moleFraction(comp)); + } + } + } + } + } +} + +BOOST_AUTO_TEST_CASE(CompvdTwoZoneRejectsRowsAgainstContact) +{ + // Type 3 takes the gas-oil contact as its reference depth, so the rows may + // not put the phase change elsewhere. No item 10 accepts a vapour row below + // a liquid one. + const std::string gapAbove = "COMPVD\n" + " 2000 0 0.95 0.05 0 150.0\n" + " 2049 0 0.95 0.05 0 150.0\n" + " 2051 0 0.60 0.40 1 150.0\n" + " 2100 0 0.40 0.60 1 150.0 /\n"; + const std::string interleaved = "COMPVD\n" + " 2000 0 0.95 0.05 0 150.0\n" + " 2040 0 0.60 0.40 1 150.0\n" + " 2060 0 0.95 0.05 0 150.0\n" + " 2100 0 0.40 0.60 1 150.0 /\n"; + for (const auto& [equil, compvd] : + {std::pair{"EQUIL\n 2060 200 2300 0 2060 0 3* 3 /\n", gapAbove}, + std::pair{"EQUIL\n 2072.5 200 2300 0 2050 0 /\n", interleaved}}) { + BOOST_TEST_CONTEXT(equil << compvd) { + const EquilFixture fix(deckString(equil, "EQLDIMS\n/\n", "", compvd)); + BOOST_CHECK_THROW(fix.compute(std::vector(20, 0)), std::runtime_error); + } + } +} + BOOST_AUTO_TEST_CASE(WaterZoneBelowContact) { // The water-oil contact at 2050 m lies inside the column. Above it the From 841fd81fac842574144a7e4a5d5bde0f19092d29 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Sat, 3 Oct 2026 00:47:36 +0200 Subject: [PATCH 02/10] Register the COMPVD contact-mismatch regression Its deck comes from OPM/opm-tests#1638. --- regressionTests.cmake | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/regressionTests.cmake b/regressionTests.cmake index 29f184d7f4d..1f9ba0d241c 100644 --- a/regressionTests.cmake +++ b/regressionTests.cmake @@ -218,6 +218,23 @@ add_test_compareECLFiles( compositional/equilibration ) +add_test_compareECLFiles( + CASENAME + equil_1d_compvd_water_gascap_contact_mismatch + FILENAME + EQUIL_1D_COMPVD_WATER_GASCAP_CONTACT_MISMATCH + SIMULATOR + flow_comp + REFERENCE_SIMULATOR + flow_comp + ABS_TOL + ${abs_tol} + REL_TOL + ${rel_tol} + DIR + compositional/equilibration +) + add_test_compareECLFiles( CASENAME equil_1d_compvd_oil From c0b31b426d7b973453c5d7a13e6198b100c302fa Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Sat, 3 Oct 2026 11:39:23 +0200 Subject: [PATCH 03/10] Keep the last vapour row of a two-zone COMPVD table in the gas zone A cell on the contact belongs to the liquid zone, so a contact moved onto the last vapour row gave that row's depth the liquid composition. With EQUIL item 10 = 1 the rows decide the phase. --- .../flow/equil/InitStateEquilComp.hpp | 19 ++++++++++++- tests/test_compequil.cpp | 27 +++++++++++++++++++ 2 files changed, 45 insertions(+), 1 deletion(-) diff --git a/opm/simulators/flow/equil/InitStateEquilComp.hpp b/opm/simulators/flow/equil/InitStateEquilComp.hpp index 1dff1b80133..5424a66e831 100644 --- a/opm/simulators/flow/equil/InitStateEquilComp.hpp +++ b/opm/simulators/flow/equil/InitStateEquilComp.hpp @@ -338,6 +338,10 @@ class InitialStateComputer /// COMPVD naming both phases describes a gas zone over a liquid one, /// each with its own composition and its own hydrostatic column. bool twoZone{false}; + /// Depth of the last vapour row of such a table. With EQUIL item 10 = 1 + /// the rows decide the phase, so a cell at this depth stays in the gas + /// zone even when the contact lies on it. + Scalar lastVaporDepth{}; /// Gas-zone composition for a COMPVD table that names both phases. std::vector vaporVdTable; /// Equilibrium vapour at the contact for type 3 without gas-zone rows. @@ -679,6 +683,7 @@ class InitialStateComputer const auto liquidRows = rowsOfPhase(compvd, CompvdTable::Phase::Liquid); reg.zgoc = zoneBoundary(compvd, vaporRows, liquidRows, reg, regionIdx); reg.twoZone = true; + reg.lastVaporDepth = compvd.getDepthColumn()[vaporRows.back()]; setupComposition(reg.vaporVdTable, compvd, vaporRows); setupComposition(reg.compositionVdTable, compvd, liquidRows); } @@ -1040,10 +1045,22 @@ class InitialStateComputer numSamplePoints, waterContactSpan(reg, span)); } + /// Whether \p depth lies in the gas zone of a region that has one. A depth + /// on the gas-oil contact belongs to the liquid zone, except the last vapour + /// row of a two-zone COMPVD table with EQUIL item 10 = 1, whose rows decide + /// the phase. + static bool isInGasZone(const Region& reg, const Scalar depth) + { + if (reg.twoZone && (reg.initType == 1) && (depth <= reg.lastVaporDepth)) { + return true; + } + return ((reg.initType == 3) || reg.twoZone) && (depth < reg.zgoc); + } + Scalar assignCell(FluidState& fs, const Region& reg, const Scalar depth, const std::size_t cell) const { - const bool inGasZone = ((reg.initType == 3) || reg.twoZone) && (depth < reg.zgoc); + const bool inGasZone = isInGasZone(reg, depth); const CompVec z = [®, depth, inGasZone]() { if (!inGasZone) { diff --git a/tests/test_compequil.cpp b/tests/test_compequil.cpp index 8fdb02d7388..5e0d1ef7bd4 100644 --- a/tests/test_compequil.cpp +++ b/tests/test_compequil.cpp @@ -894,6 +894,33 @@ BOOST_AUTO_TEST_CASE(CompvdTwoZoneContactOutsideRowGap) } } +BOOST_AUTO_TEST_CASE(CompvdTwoZoneKeepsTheLastVapourRowGas) +{ + // A cell on the gas-oil contact belongs to the liquid zone, but with EQUIL + // item 10 = 1 the rows decide the phase. The cell at 2032.5 m sits on the + // last vapour row, both when the contact lies on that row and when it moves + // there from above. + for (const std::string goc : {"2032.5", "2020"}) { + BOOST_TEST_CONTEXT("Contact at " << goc << " m") { + const EquilFixture fix(deckString( + "EQUIL\n 2072.5 200 2300 0 " + goc + " 0 /\n", "EQLDIMS\n/\n", "", + "COMPVD\n" + " 2000 0 0.95 0.05 0 150.0\n" + " 2032.5 0 0.95 0.05 0 150.0\n" + " 2051 0 0.60 0.40 1 150.0\n" + " 2100 0 0.40 0.60 1 150.0 /\n")); + const auto states = fix.compute(std::vector(20, 0)).fluidStates(); + BOOST_REQUIRE_EQUAL(fix.depths[6], 2032.5); + BOOST_CHECK_CLOSE(states[6].saturation(FluidSystem::gasPhaseIdx), 1.0, 1e-10); + BOOST_CHECK_SMALL(states[6].moleFraction(1) - 0.95, 1e-10); + + // The cells below that row belong to the liquid zone. + BOOST_CHECK_CLOSE(states[7].saturation(FluidSystem::oilPhaseIdx), 1.0, 1e-10); + BOOST_CHECK_SMALL(states[7].moleFraction(1) - 0.60, 1e-10); + } + } +} + BOOST_AUTO_TEST_CASE(CompvdTwoZoneRejectsRowsAgainstContact) { // Type 3 takes the gas-oil contact as its reference depth, so the rows may From 05492eced43f47a66c700b2c5a066cee219fe277 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Sat, 3 Oct 2026 11:39:23 +0200 Subject: [PATCH 04/10] Anchor the water to the hydrocarbon zone at its contact The water took the liquid column's pressure even when the gas zone reaches the water-oil contact, as it does once the gas-oil contact lies below it. --- .../flow/equil/InitStateEquilComp.hpp | 6 ++++- tests/test_compequil.cpp | 27 +++++++++++++++++++ 2 files changed, 32 insertions(+), 1 deletion(-) diff --git a/opm/simulators/flow/equil/InitStateEquilComp.hpp b/opm/simulators/flow/equil/InitStateEquilComp.hpp index 5424a66e831..071be33c729 100644 --- a/opm/simulators/flow/equil/InitStateEquilComp.hpp +++ b/opm/simulators/flow/equil/InitStateEquilComp.hpp @@ -810,7 +810,11 @@ class InitialStateComputer if (!reg.oilPressure.has_value() && !reg.gasPressure.has_value()) { return; } - const auto& hcPressure = reg.oilPressure.has_value() ? reg.oilPressure : reg.gasPressure; + // The water meets the hydrocarbon zone just above its contact: the gas + // zone when the gas-oil contact lies at or below the water-oil contact. + const bool gasAtContact = reg.gasPressure.has_value() + && (!reg.oilPressure.has_value() || (reg.zwoc <= reg.zgoc)); + const auto& hcPressure = gasAtContact ? reg.gasPressure : reg.oilPressure; const Scalar pcow = record.waterOilContactCapillaryPressure(); const Scalar pContact = hcPressure->value(reg.zwoc) - pcow; diff --git a/tests/test_compequil.cpp b/tests/test_compequil.cpp index 5e0d1ef7bd4..c0a24739515 100644 --- a/tests/test_compequil.cpp +++ b/tests/test_compequil.cpp @@ -982,6 +982,33 @@ BOOST_AUTO_TEST_CASE(WaterZoneBelowContact) BOOST_CHECK_GT(rhoWater, rhoHc); } +BOOST_AUTO_TEST_CASE(GasZoneMeetingTheWaterAnchorsIt) +{ + // The vapour rows reach 2060 m, so the gas-oil contact at 2000 m moves down + // to 2060 m, below the water-oil contact at 2052.5 m. The gas zone then + // meets the water, which takes its pressure from the gas at the contact + // rather than from a liquid column extended above its own zone. + const WaterEquilFixture fix(waterDeckString( + "EQUIL\n 2010 150 2052.5 0 2000 0 /\n", + "COMPVD\n" + " 2000 0 0.95 0.05 0 150.0\n" + " 2060 0 0.95 0.05 0 150.0\n" + " 2070 0 0.60 0.40 1 150.0\n" + " 2100 0 0.40 0.60 1 150.0 /\n")); + const auto states = fix.compute(std::vector(20, 0), + std::vector(20, connateSw), + std::vector(20, 1.0)).fluidStates(); + + // Cell 10 is centred on the water-oil contact, where the zero capillary + // pressure leaves the water and the gas at one pressure. + BOOST_REQUIRE_EQUAL(fix.depths[10], 2052.5); + BOOST_CHECK_CLOSE(states[10].pressure(WaterFluidSystem::waterPhaseIdx), + states[10].pressure(WaterFluidSystem::gasPhaseIdx), 1e-8); + BOOST_CHECK_CLOSE(states[10].saturation(WaterFluidSystem::gasPhaseIdx), + 1.0 - connateSw, 1e-10); + BOOST_CHECK_CLOSE(states[11].saturation(WaterFluidSystem::waterPhaseIdx), 1.0, 1e-10); +} + BOOST_AUTO_TEST_CASE(CoincidentContactsKeepTheGasRoot) { // Contacts that coincide leave no liquid: every hydrocarbon cell sits above From 6437cc523f30fc53b0e47b3895f2a1e64b0b25db Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Sat, 3 Oct 2026 20:21:17 +0200 Subject: [PATCH 05/10] Tighten the COMPVD zone-boundary checks and messages Clamp once and compare, so the row gap is defined in one place. Once the boundary moves it is no longer the EQUIL contact, so the later messages stop calling it that. The rejection test checks which error it gets. --- .../flow/equil/InitStateEquilComp.hpp | 21 +++++++++++-------- tests/test_compequil.cpp | 18 ++++++++++++---- 2 files changed, 26 insertions(+), 13 deletions(-) diff --git a/opm/simulators/flow/equil/InitStateEquilComp.hpp b/opm/simulators/flow/equil/InitStateEquilComp.hpp index 071be33c729..e911853ee9d 100644 --- a/opm/simulators/flow/equil/InitStateEquilComp.hpp +++ b/opm/simulators/flow/equil/InitStateEquilComp.hpp @@ -491,8 +491,9 @@ class InitialStateComputer regionIdx + 1, lastVapor, firstLiquid)); } - if ((reg.zgoc >= lastVapor) && (reg.zgoc <= firstLiquid)) { - return reg.zgoc; + const Scalar boundary = std::clamp(reg.zgoc, lastVapor, firstLiquid); + if (boundary == reg.zgoc) { + return boundary; } if (reg.initType != 1) { @@ -506,7 +507,6 @@ class InitialStateComputer reg.initType)); } - const Scalar boundary = std::clamp(reg.zgoc, lastVapor, firstLiquid); OpmLog::warning(fmt::format("Equilibration region {}: the gas-oil contact at {} m lies " "outside the COMPVD phase change between the last vapour " "row at {} m and the first liquid row at {} m. With EQUIL " @@ -934,10 +934,11 @@ class InitialStateComputer const std::array datumSpan{std::min(pressureSpan[0], reg.zgoc), std::max(pressureSpan[1], reg.zgoc)}; if ((reg.zgoc < span[0]) || (reg.zgoc > span[1])) { - OpmLog::warning(fmt::format("Equilibration region {}: the gas-oil contact at {} m " - "lies outside the cells of the region, so the COMPVD " + OpmLog::warning(fmt::format("Equilibration region {}: the gas zone ends at {} m, " + "outside the cells of the region, so the COMPVD " "rows of one phase describe no cell.", - regionIdx + 1, reg.zgoc)); + regionIdx + 1, + reg.zgoc)); } if (datum < reg.zgoc) { @@ -959,9 +960,11 @@ class InitialStateComputer numSamplePoints, pressureSpan); } - OpmLog::info(fmt::format("Equilibration region {}: COMPVD gives a gas zone above the " - "contact at {} m and a liquid one below it " - "(EQUIL item 10 = 1).", regionIdx + 1, reg.zgoc)); + OpmLog::info(fmt::format("Equilibration region {}: COMPVD gives a gas zone above " + "{} m and a liquid one below it " + "(EQUIL item 10 = 1).", + regionIdx + 1, + reg.zgoc)); } /// EQUIL item 10 type 3: the selected table supplies the liquid composition, diff --git a/tests/test_compequil.cpp b/tests/test_compequil.cpp index c0a24739515..5c282f39cb6 100644 --- a/tests/test_compequil.cpp +++ b/tests/test_compequil.cpp @@ -45,6 +45,7 @@ #include #include #include +#include #include #include #include @@ -936,12 +937,21 @@ BOOST_AUTO_TEST_CASE(CompvdTwoZoneRejectsRowsAgainstContact) " 2040 0 0.60 0.40 1 150.0\n" " 2060 0 0.95 0.05 0 150.0\n" " 2100 0 0.40 0.60 1 150.0 /\n"; - for (const auto& [equil, compvd] : - {std::pair{"EQUIL\n 2060 200 2300 0 2060 0 3* 3 /\n", gapAbove}, - std::pair{"EQUIL\n 2072.5 200 2300 0 2050 0 /\n", interleaved}}) { + for (const auto& [equil, compvd, reason] : + {std::tuple {"EQUIL\n 2060 200 2300 0 2060 0 3* 3 /\n", + gapAbove, + "so the rows must agree with it"}, + std::tuple {"EQUIL\n 2072.5 200 2300 0 2050 0 /\n", + interleaved, + "must lie above the liquid rows"}}) { BOOST_TEST_CONTEXT(equil << compvd) { const EquilFixture fix(deckString(equil, "EQLDIMS\n/\n", "", compvd)); - BOOST_CHECK_THROW(fix.compute(std::vector(20, 0)), std::runtime_error); + BOOST_CHECK_EXCEPTION(fix.compute(std::vector(20, 0)), + std::runtime_error, + [reason](const std::runtime_error& e) { + return std::string_view {e.what()}.find(reason) + != std::string_view::npos; + }); } } } From 139f34f95224d4cebb4877b8ec150057548fa303 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 28 Sep 2026 14:34:12 +0200 Subject: [PATCH 06/10] Apply rock compressibility in the compositional model The flash intensive quantities took the reference porosity, so ROCK had no effect on compositional runs. Scale the porosity with the pressure as the black-oil model does; the rock table follows PVTNUM or ROCKNUM. ROCKCOMP is rejected, as its compaction tables do not enter the porosity. --- examples/problems/co2ptflashproblem.hh | 12 +++++++++ .../ptflash/flashintensivequantities.hh | 9 ++++++- tests/test_flashphasepresence.cpp | 26 +++++++++++++++++++ 3 files changed, 46 insertions(+), 1 deletion(-) diff --git a/examples/problems/co2ptflashproblem.hh b/examples/problems/co2ptflashproblem.hh index 8a9a25083e4..cb41d20c765 100644 --- a/examples/problems/co2ptflashproblem.hh +++ b/examples/problems/co2ptflashproblem.hh @@ -285,6 +285,18 @@ public: return Opm::CompositionalConfig::EOSType::PR; } + template + Scalar rockCompressibility(const Context&, unsigned, unsigned) const + { + return 0.0; + } + + template + Scalar rockReferencePressure(const Context&, unsigned, unsigned) const + { + return 1.0e5; + } + /*! * \copydoc FvBaseProblem::finishInit */ diff --git a/opm/models/ptflash/flashintensivequantities.hh b/opm/models/ptflash/flashintensivequantities.hh index 9b97252426d..12f1f73ef34 100644 --- a/opm/models/ptflash/flashintensivequantities.hh +++ b/opm/models/ptflash/flashintensivequantities.hh @@ -323,8 +323,15 @@ public: // Compute the remaining quantities ///////////// - // porosity + // porosity, scaled with the pressure by the rock compressibility of the + // cell's region as in the black-oil model porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx); + const Scalar rockCompressibility = problem.rockCompressibility(elemCtx, dofIdx, timeIdx); + if (rockCompressibility > 0.0) { + const Scalar rockRefPressure = problem.rockReferencePressure(elemCtx, dofIdx, timeIdx); + const Evaluation x = rockCompressibility * (p - rockRefPressure); + porosity_ *= 1.0 + x + 0.5 * x * x; + } Valgrind::CheckDefined(porosity_); // intrinsic permeability diff --git a/tests/test_flashphasepresence.cpp b/tests/test_flashphasepresence.cpp index ec616677f46..3cf0573050c 100644 --- a/tests/test_flashphasepresence.cpp +++ b/tests/test_flashphasepresence.cpp @@ -176,6 +176,12 @@ struct TestProblem template double porosity(const Context&, unsigned, unsigned) const { return 0.2; } template + double rockCompressibility(const Context&, unsigned, unsigned) const { return compressibility; } + template + double rockReferencePressure(const Context&, unsigned, unsigned) const { return 5.0e6; } + + double compressibility = 0.0; // 1/Pa + template Dune::FieldMatrix intrinsicPermeability(const Context&, unsigned, unsigned) const { return Dune::FieldMatrix{1.0}; } }; @@ -283,3 +289,23 @@ BOOST_AUTO_TEST_CASE(PresenceIsRefreshedWhenHydrocarbonEntersAndLeaves) BOOST_CHECK(!iq.phaseIsPresent(1)); BOOST_CHECK_EQUAL(iq.saturationForOutput(2), 1.0); } + +BOOST_AUTO_TEST_CASE(PorosityFollowsTheRockCompressibility) +{ + // The pore volume grows with the pressure over the reference pressure, + // with the second-order expansion of exp(x) the black-oil model uses. + TestContext context; + context.problem_.compressibility = 1.0e-9; + context.priVars.values[2] = 0.5; + IntensiveQuantities iq; + iq.update(context, 0, 0); + + const double x = 1.0e-9 * (1.0e7 - 5.0e6); + BOOST_CHECK_CLOSE(iq.porosity().value(), 0.2 * (1.0 + x + 0.5 * x * x), 1e-12); + BOOST_CHECK_CLOSE(iq.porosity().derivative(0), 0.2 * 1.0e-9 * (1.0 + x), 1e-10); + + context.problem_.compressibility = 0.0; + iq.update(context, 0, 0); + BOOST_CHECK_EQUAL(iq.porosity().value(), 0.2); + BOOST_CHECK_EQUAL(iq.porosity().derivative(0), 0.0); +} From 4ab87a651e15cd7920c71b1409959a96b3fdba85 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 28 Sep 2026 14:36:49 +0200 Subject: [PATCH 07/10] Take the water properties of each cell's PVT region in compositional Flow The cells took the water PVT and the water surface density of the first PVT region. Select the cell's PVT region in the parameter cache. --- examples/problems/co2ptflashproblem.hh | 6 ++++ .../ptflash/flashintensivequantities.hh | 2 ++ tests/test_flashphasepresence.cpp | 28 +++++++++++++++++-- 3 files changed, 34 insertions(+), 2 deletions(-) diff --git a/examples/problems/co2ptflashproblem.hh b/examples/problems/co2ptflashproblem.hh index cb41d20c765..0d2c77d5fa5 100644 --- a/examples/problems/co2ptflashproblem.hh +++ b/examples/problems/co2ptflashproblem.hh @@ -285,6 +285,12 @@ public: return Opm::CompositionalConfig::EOSType::PR; } + template + unsigned pvtRegionIndex(const Context&, unsigned, unsigned) const + { + return 0; + } + template Scalar rockCompressibility(const Context&, unsigned, unsigned) const { diff --git a/opm/models/ptflash/flashintensivequantities.hh b/opm/models/ptflash/flashintensivequantities.hh index 12f1f73ef34..1f2ed0d409c 100644 --- a/opm/models/ptflash/flashintensivequantities.hh +++ b/opm/models/ptflash/flashintensivequantities.hh @@ -228,6 +228,8 @@ public: // Update phases typename FluidSystem::template ParameterCache paramCache(eos_type); + // the water properties of the cell's PVT region + paramCache.setRegionIndex(problem.pvtRegionIndex(elemCtx, dofIdx, timeIdx)); paramCache.updatePhase(fluidState_, FluidSystem::oilPhaseIdx); paramCache.updatePhase(fluidState_, FluidSystem::gasPhaseIdx); diff --git a/tests/test_flashphasepresence.cpp b/tests/test_flashphasepresence.cpp index 3cf0573050c..85fe8208bd2 100644 --- a/tests/test_flashphasepresence.cpp +++ b/tests/test_flashphasepresence.cpp @@ -59,11 +59,18 @@ struct TestFluidSystem void updatePhase(const FluidState&, unsigned) {} Eval molarVolume(unsigned) const { return Eval{1.0}; } Eval correctedMolarVolume(unsigned) const { return Eval{1.0}; } + void setRegionIndex(unsigned region) { regionIdx = region; } + unsigned regionIndex() const { return regionIdx; } + unsigned regionIdx = 0; }; + // Water is 100 kg/m3 denser in each further PVT region. template - static auto density(const FluidState&, const Cache&, unsigned) - { return typename FluidState::ValueType{1000.0}; } + static auto density(const FluidState&, const Cache& cache, unsigned phase) + { + return typename FluidState::ValueType{ + phase == waterPhaseIdx ? 1000.0 + 100.0 * cache.regionIndex() : 1000.0}; + } template static auto viscosity(const FluidState&, const Cache&, unsigned) @@ -179,8 +186,11 @@ struct TestProblem double rockCompressibility(const Context&, unsigned, unsigned) const { return compressibility; } template double rockReferencePressure(const Context&, unsigned, unsigned) const { return 5.0e6; } + template + unsigned pvtRegionIndex(const Context&, unsigned, unsigned) const { return pvtRegion; } double compressibility = 0.0; // 1/Pa + unsigned pvtRegion = 0; template Dune::FieldMatrix intrinsicPermeability(const Context&, unsigned, unsigned) const { return Dune::FieldMatrix{1.0}; } @@ -309,3 +319,17 @@ BOOST_AUTO_TEST_CASE(PorosityFollowsTheRockCompressibility) BOOST_CHECK_EQUAL(iq.porosity().value(), 0.2); BOOST_CHECK_EQUAL(iq.porosity().derivative(0), 0.0); } + +BOOST_AUTO_TEST_CASE(WaterPropertiesTakeTheCellsPvtRegion) +{ + TestContext context; + context.priVars.values[2] = 0.5; + IntensiveQuantities iq; + iq.update(context, 0, 0); + BOOST_CHECK_EQUAL(iq.fluidState().density(2).value(), 1000.0); + + context.problem_.pvtRegion = 1; + iq.update(context, 0, 0); + BOOST_CHECK_EQUAL(iq.fluidState().density(2).value(), 1100.0); + BOOST_CHECK_EQUAL(iq.fluidState().density(0).value(), 1000.0); +} From 8528fb2a3c687805a3888b438417d5dcd0a3833c Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 28 Sep 2026 14:41:05 +0200 Subject: [PATCH 08/10] Equilibrate each water column with the water of its PVT region The water column of every equilibration region took the water PVT of the first PVT region. Integrate a water column for each PVT region among the cells of an equilibration region and give each cell the column of its own region. A datum below the water-oil contact gives the pressure of all of these waters, so there they must reach the contact at the same pressure. --- opm/simulators/flow/FlowProblemComp.hpp | 1 + .../flow/equil/InitStateEquilComp.hpp | 140 ++++++++++-- tests/test_compequil.cpp | 214 +++++++++++++++++- 3 files changed, 336 insertions(+), 19 deletions(-) diff --git a/opm/simulators/flow/FlowProblemComp.hpp b/opm/simulators/flow/FlowProblemComp.hpp index f8eb41b4daa..40f49afeb8b 100644 --- a/opm/simulators/flow/FlowProblemComp.hpp +++ b/opm/simulators/flow/FlowProblemComp.hpp @@ -437,6 +437,7 @@ class FlowProblemComp : public FlowProblem getEosType(), vanguard.cellCenterDepths(), eqlnum, + this->pvtnum_, vanguard.gridView().comm(), this->gravity()[dimWorld - 1], this->numPressurePointsEquil(), diff --git a/opm/simulators/flow/equil/InitStateEquilComp.hpp b/opm/simulators/flow/equil/InitStateEquilComp.hpp index e911853ee9d..27b184cb159 100644 --- a/opm/simulators/flow/equil/InitStateEquilComp.hpp +++ b/opm/simulators/flow/equil/InitStateEquilComp.hpp @@ -62,6 +62,7 @@ #include #include #include +#include #include #include #include @@ -148,9 +149,11 @@ class WaterDensityODE WaterDensityODE(const TabulatedFunction& tempVdTable, const CompositionalConfig::EOSType eosType, + const unsigned pvtRegion, const Scalar normGrav) : tempVdTable_(tempVdTable) , eosType_(eosType) + , pvtRegion_(pvtRegion) , g_(normGrav) {} @@ -162,6 +165,7 @@ class WaterDensityODE fs.setPressure(FluidSystem::waterPhaseIdx, press); typename FluidSystem::template ParameterCache paramCache(eosType_); + paramCache.setRegionIndex(pvtRegion_); return FluidSystem::density(fs, paramCache, FluidSystem::waterPhaseIdx) * g_; } @@ -169,6 +173,7 @@ class WaterDensityODE private: const TabulatedFunction& tempVdTable_; CompositionalConfig::EOSType eosType_; + unsigned pvtRegion_; Scalar g_; }; @@ -208,6 +213,11 @@ class WaterDensityODE * Only cell-centre initialization is supported (EQUIL item 9 = 0). * Gas-oil contact capillary pressure must be zero: the downstream flash uses * a single pressure for all phases. + * + * Each PVT region among the cells of an equilibration region has a water column + * of its own, integrated with its water, and each cell takes the column of its + * PVT region. A datum below the water-oil contact gives the pressure of all of + * these waters, so there they must reach the contact at the same pressure. */ template class InitialStateComputer @@ -221,6 +231,8 @@ class InitialStateComputer /// \param[in] eosType Equation of state used by the fluid system. /// \param[in] cellCenterDepth Depth of each cell centre. /// \param[in] eqlnum Zero-based equilibration region of each cell. + /// \param[in] pvtnum Zero-based PVT region of each cell, or empty + /// for a single PVT region. /// \param[in] comm Communicator for parallel runs. /// \param[in] gravity Norm of the gravity vector. /// \param[in] numSamplePoints Sample points in each pressure integration. @@ -232,6 +244,7 @@ class InitialStateComputer const CompositionalConfig::EOSType eosType, const std::vector& cellCenterDepth, const std::vector& eqlnum, + const std::vector& pvtnum, const Parallel::Communication& comm, const Scalar gravity, const int numSamplePoints, @@ -250,6 +263,7 @@ class InitialStateComputer "versus depth from the ZMFVD or the COMPVD keyword."); } + const auto numPvtRegions = tables.getTabdims().getNumPVTTables(); OPM_BEGIN_PARALLEL_TRY_CATCH(); if (eqlnum.size() != cellCenterDepth.size()) { OPM_THROW(std::runtime_error, @@ -265,6 +279,20 @@ class InitialStateComputer cell, region + 1, records.size())); } } + if (!pvtnum.empty() && (pvtnum.size() != cellCenterDepth.size())) { + OPM_THROW(std::runtime_error, + fmt::format("PVTNUM contains {} entries for {} cell depths.", + pvtnum.size(), cellCenterDepth.size())); + } + for (std::size_t cell = 0; cell < pvtnum.size(); ++cell) { + const auto region = pvtnum[cell]; + if (region < 0 || std::cmp_greater_equal(region, numPvtRegions)) { + OPM_THROW(std::runtime_error, + fmt::format("Cell {} has PVTNUM {} outside the {} " + "PVT regions.", + cell, region + 1, numPvtRegions)); + } + } // The endpoint vectors are optional, but a non-empty one is indexed for // every cell. for (const auto& [name, limits] : {std::pair{"connate water", std::cref(connateWater)}, @@ -278,18 +306,22 @@ class InitialStateComputer } OPM_END_PARALLEL_TRY_CATCH("Invalid equilibration input: ", comm); + const auto pvtRegions = + pvtRegionsOfEachRegion(eqlnum, pvtnum, records.size(), numPvtRegions, comm); std::vector regions; regions.reserve(records.size()); for (std::size_t r = 0; r < records.size(); ++r) { - regions.push_back(setupRegion(records.getRecord(r), tables, cellCenterDepth, - eqlnum, comm, gravity, numSamplePoints, r)); + regions.push_back(setupRegion(records.getRecord(r), tables, pvtRegions[r], + cellCenterDepth, eqlnum, comm, gravity, + numSamplePoints, r)); } fluidStates_.resize(cellCenterDepth.size()); referencePressures_.resize(cellCenterDepth.size()); for (std::size_t cell = 0; cell < cellCenterDepth.size(); ++cell) { - referencePressures_[cell] = - assignCell(fluidStates_[cell], regions[eqlnum[cell]], cellCenterDepth[cell], cell); + const auto pvtRegion = pvtnum.empty() ? 0u : static_cast(pvtnum[cell]); + referencePressures_[cell] = assignCell(fluidStates_[cell], regions[eqlnum[cell]], + pvtRegion, cellCenterDepth[cell], cell); } } @@ -327,6 +359,8 @@ class InitialStateComputer /// The equilibrated vertical distributions within one region. struct Region { int initType{1}; // EQUIL item 10 + /// Zero-based PVT regions of the region's cells on all processes. + std::vector pvtRegions; /// Gas-oil contact. A two-zone COMPVD table may move it, see zoneBoundary(). Scalar zgoc{}; /// Name of the selected composition keyword, used in diagnostics. @@ -358,7 +392,8 @@ class InitialStateComputer /// anchored at the water-oil contact instead. The fluid there is the one /// just above the contact, which the EOS root has to follow. bool anchoredAtWaterContact{false}; - std::optional waterPressure; + /// The water column of each PVT region in pvtRegions. + std::map waterPressure; }; /// Each pressure column must reach the water-oil contact before another @@ -568,8 +603,39 @@ class InitialStateComputer }); } + /// The zero-based PVT regions of the cells of each equilibration region, on + /// all processes. A region without cells has none. + static std::vector> + pvtRegionsOfEachRegion(const std::vector& eqlnum, + const std::vector& pvtnum, + const std::size_t numRegions, + const std::size_t numPvtRegions, + const Parallel::Communication& comm) + { + // Whether equilibration region r holds a cell of PVT region p, at + // r * numPvtRegions + p. + std::vector present(numRegions * numPvtRegions, 0); + for (std::size_t cell = 0; cell < eqlnum.size(); ++cell) { + const auto pvtRegion = pvtnum.empty() ? 0 : pvtnum[cell]; + present[static_cast(eqlnum[cell]) * numPvtRegions + + static_cast(pvtRegion)] = 1; + } + comm.max(present.data(), static_cast(present.size())); + + std::vector> pvtRegions(numRegions); + for (std::size_t r = 0; r < numRegions; ++r) { + for (std::size_t p = 0; p < numPvtRegions; ++p) { + if (present[r * numPvtRegions + p] != 0) { + pvtRegions[r].push_back(static_cast(p)); + } + } + } + return pvtRegions; + } + Region setupRegion(const EquilRecord& record, const TableManager& tables, + const std::vector& pvtRegions, const std::vector& cellCenterDepth, const std::vector& eqlnum, const Parallel::Communication& comm, @@ -579,6 +645,7 @@ class InitialStateComputer { Region reg; + reg.pvtRegions = pvtRegions; reg.initType = record.compositionalInitType(); if (reg.initType != 1 && reg.initType != 3) { OPM_THROW(std::runtime_error, @@ -751,7 +818,7 @@ class InitialStateComputer integrateWaterPressure(reg, span, gravity, numSamplePoints, record.datumDepth(), record.datumDepthPressure()); hcDatum = reg.zwoc; - hcPressure = reg.waterPressure->value(reg.zwoc) + hcPressure = waterContactPressure(reg, record.datumDepth(), regionIdx) + record.waterOilContactCapillaryPressure(); OpmLog::info(fmt::format("Equilibration region {}: the datum at {} m lies below the " "water-oil contact at {} m, so it gives the water pressure; " @@ -824,8 +891,8 @@ class InitialStateComputer "is at {} m.", regionIdx + 1, reg.zwoc)); } - /// Integrates the water pressure over \p span from \p depth, where it is - /// \p pressure. + /// Integrates the water pressure of each PVT region of \p reg over \p span + /// from \p depth, where it is \p pressure. void integrateWaterPressure(Region& reg, const std::array& span, const Scalar gravity, @@ -833,10 +900,43 @@ class InitialStateComputer const Scalar depth, const Scalar pressure) const { - const WaterODE ode(reg.tempVdTable, eosType_, gravity); - reg.waterPressure.emplace(ode, - typename WaterPressFunc::InitCond{depth, pressure}, - numSamplePoints, waterContactSpan(reg, span)); + for (const auto pvtRegion : reg.pvtRegions) { + const WaterODE ode(reg.tempVdTable, eosType_, pvtRegion, gravity); + reg.waterPressure.try_emplace(pvtRegion, ode, + typename WaterPressFunc::InitCond{depth, pressure}, + numSamplePoints, waterContactSpan(reg, span)); + } + } + + /// The water pressure at the water-oil contact when the datum at \p datum, + /// below the contact, gives it. The datum states the pressure of the water + /// of every PVT region in the region, so their columns must reach the + /// contact at one pressure for the hydrocarbon to be anchored there. + static Scalar waterContactPressure(const Region& reg, + const Scalar datum, + const std::size_t regionIdx) + { + // Waters that differ physically part by more than this round-off + // tolerance between the datum and the contact. + constexpr Scalar samePressure{1.0e-10}; + + const auto& [firstPvtRegion, firstColumn] = *reg.waterPressure.begin(); + const Scalar pressure = firstColumn.value(reg.zwoc); + for (const auto& [pvtRegion, column] : reg.waterPressure) { + const Scalar other = column.value(reg.zwoc); + if (std::abs(other - pressure) > samePressure * pressure) { + OPM_THROW(std::runtime_error, + fmt::format("Equilibration region {} places the datum at {} m, below " + "the water-oil contact at {} m, but the water of its PVT " + "regions {} and {} reaches the contact at {:.6g} and " + "{:.6g} bar. Put the datum in the hydrocarbon column.", + regionIdx + 1, datum, reg.zwoc, + firstPvtRegion + 1, pvtRegion + 1, + unit::convert::to(pressure, unit::barsa), + unit::convert::to(other, unit::barsa))); + } + } + return pressure; } /// A per-cell saturation endpoint, or \p fallback when the caller supplied @@ -1064,8 +1164,8 @@ class InitialStateComputer return ((reg.initType == 3) || reg.twoZone) && (depth < reg.zgoc); } - Scalar assignCell(FluidState& fs, const Region& reg, const Scalar depth, - const std::size_t cell) const + Scalar assignCell(FluidState& fs, const Region& reg, const unsigned pvtRegion, + const Scalar depth, const std::size_t cell) const { const bool inGasZone = isInGasZone(reg, depth); @@ -1088,8 +1188,12 @@ class InitialStateComputer // saturation leaves some residual hydrocarbon. Keep this separate from // the phase pressures so equilibration retains the capillary offset. const Scalar hydrocarbonPressure = pressFunc->value(depth); - const bool inWaterZone = (depth > reg.zwoc) && reg.waterPressure.has_value(); - const Scalar press = inWaterZone ? reg.waterPressure->value(depth) + // The water column of the cell's own PVT region. + const auto water = reg.waterPressure.find(pvtRegion); + const WaterPressFunc* waterPressure = + (water != reg.waterPressure.end()) ? &water->second : nullptr; + const bool inWaterZone = (depth > reg.zwoc) && (waterPressure != nullptr); + const Scalar press = inWaterZone ? waterPressure->value(depth) : hydrocarbonPressure; fs.setTemperature(Details::evalDepthTable(reg.tempVdTable, depth)); @@ -1110,8 +1214,8 @@ class InitialStateComputer sWat = (depth > reg.zwoc) ? waterLimit(maxWater_, cell, Scalar{1}) : waterLimit(connateWater_, cell, Scalar{0}); fs.setSaturation(FluidSystem::waterPhaseIdx, sWat); - if (reg.waterPressure.has_value()) { - fs.setPressure(FluidSystem::waterPhaseIdx, reg.waterPressure->value(depth)); + if (waterPressure != nullptr) { + fs.setPressure(FluidSystem::waterPhaseIdx, waterPressure->value(depth)); } } diff --git a/tests/test_compequil.cpp b/tests/test_compequil.cpp index 5c282f39cb6..d029e1a3193 100644 --- a/tests/test_compequil.cpp +++ b/tests/test_compequil.cpp @@ -134,14 +134,17 @@ struct BasicEquilFixture /// \param connateWater, maxWater Per-cell water saturation endpoints, as /// the simulator reads them from the scaled saturation functions. + /// \param pvtnum Zero-based PVT region of each cell, empty for one region. Computer compute(const std::vector& eqlnum, const std::vector& connateWater = {}, - const std::vector& maxWater = {}) const + const std::vector& maxWater = {}, + const std::vector& pvtnum = {}) const { return Computer(eclState, eclState.compositionalConfig().eosType(0), {depths.begin(), depths.end()}, eqlnum, + pvtnum, Opm::Parallel::Communication{}, gravity, /*numSamplePoints=*/100, @@ -191,6 +194,66 @@ std::string waterDeckString(const std::string& equil, "SOF3\n 0.00 0.0 0.0\n 0.99 1.0 1.0 /\n"); } +// A deck from deckString() declaring two PVT regions, both with its rock. +std::string withTwoPvtRegions(std::string deck) +{ + for (const auto& [from, to] : {std::pair{std::string{"TABDIMS\n/\n"}, + std::string{"TABDIMS\n1* 2 /\n"}}, + std::pair{std::string{"ROCK\n68.9476 0 /\n"}, + std::string{"ROCK\n68.9476 0 /\n68.9476 0 /\n"}}}) { + const auto pos = deck.find(from); + BOOST_REQUIRE(pos != std::string::npos); + deck.replace(pos, from.size(), to); + } + return deck; +} + +// A water zone at the base of each block of ten cells. +constexpr auto twoBlockEquil = "EQUIL\n" + " 2010 150 2035 0 2000 0 /\n" + " 2060 155 2085 0 2050 0 /\n"; + +// The fresh water of PVT region 1 and the brine of PVT region 2. +constexpr auto freshWaterAndBrine = "PVTW\n" + " 75.0 1.03 4.0E-5 0.3 0.0 /\n" + " 75.0 1.01 2.5E-5 0.9 0.0 /\n" + "DENSITY\n" + " 800.0 1000.0 1.0 /\n" + " 800.0 1150.0 1.0 /\n"; + +// Two equilibration regions and two PVT regions over the column. The arguments +// give the EQLNUM and PVTNUM records, the two EQUIL records, and the PVTW and +// DENSITY of the two PVT regions. +std::string twoWaterRegionsDeck(const std::string& eqlnum, const std::string& pvtnum, + const std::string& equil = twoBlockEquil, + const std::string& water = freshWaterAndBrine) +{ + return withTwoPvtRegions(deckString(equil, + "EQLDIMS\n2 /\n", + "REGIONS\n" + "EQLNUM\n" + eqlnum + " /\n" + "PVTNUM\n" + pvtnum + " /\n", + "ZMFVD\n" + " 2000 0 0.7 0.3\n" + " 2100 0 0.3 0.7 /\n" + " 2000 0 0.7 0.3\n" + " 2100 0 0.3 0.7 /\n", + "RTEMP\n100\n/\n", + "OIL\nGAS\nWATER\n", + "SWFN\n 0.01 0.0 0.0\n 1.00 1.0 0.0 /\n" + "SGFN\n 0.00 0.0 0.0\n 0.99 1.0 0.0 /\n" + "SOF3\n 0.00 0.0 0.0\n 0.99 1.0 1.0 /\n" + + water)); +} + +// The water density of a PVT region above at a pressure in bar. +Scalar waterDensity(const Scalar surfaceDensity, const Scalar bw, const Scalar cw, + const Scalar pressure) +{ + const Scalar x = cw * (pressure - 75.0); + return surfaceDensity * (1.0 + x + 0.5 * x * x) / bw; +} + } // Anonymous namespace BOOST_AUTO_TEST_CASE(Type1LiquidRootPressureIntegration) @@ -1019,6 +1082,155 @@ BOOST_AUTO_TEST_CASE(GasZoneMeetingTheWaterAnchorsIt) BOOST_CHECK_CLOSE(states[11].saturation(WaterFluidSystem::waterPhaseIdx), 1.0, 1e-10); } +BOOST_AUTO_TEST_CASE(EachWaterColumnTakesThePvtRegionOfItsCells) +{ + // The lower block holds brine: its water column is steeper than the fresh + // water column of the upper block, and than its own column would be in the + // first PVT region. + std::vector regions(20, 0); + std::fill(regions.begin() + 10, regions.end(), 1); + const std::vector connate(20, connateSw); + const std::vector maxWater(20, 1.0); + const WaterEquilFixture fix(twoWaterRegionsDeck("10*1 10*2", "10*1 10*2")); + const auto states = fix.compute(regions, connate, maxWater, regions).fluidStates(); + const auto first = fix.compute(regions, connate, maxWater).fluidStates(); + + const auto pressure = [](const auto& fs) { + return Opm::getValue(fs.pressure(WaterFluidSystem::waterPhaseIdx)); + }; + // Cells 8 and 9 lie below the upper contact at 2035 m, 18 and 19 below the + // lower one at 2085 m. + const Scalar fresh = impliedDensity(pressure(states[8]), pressure(states[9])); + const Scalar brine = impliedDensity(pressure(states[18]), pressure(states[19])); + const Scalar pFresh = 0.5 * (pressure(states[8]) + pressure(states[9])) / barsa; + const Scalar pBrine = 0.5 * (pressure(states[18]) + pressure(states[19])) / barsa; + BOOST_CHECK_CLOSE(fresh, waterDensity(1000.0, 1.03, 4.0e-5, pFresh), 0.01); + BOOST_CHECK_CLOSE(brine, waterDensity(1150.0, 1.01, 2.5e-5, pBrine), 0.01); + BOOST_CHECK_CLOSE(impliedDensity(pressure(first[18]), pressure(first[19])), + waterDensity(1000.0, 1.03, 4.0e-5, pBrine), 0.01); +} + +BOOST_AUTO_TEST_CASE(OneEquilibrationRegionHoldsTheWaterOfEachPvtRegion) +{ + // One equilibration region over both PVT regions, with its contact at 2035 m + // among the fresh-water cells. The hydrocarbon and the fresh water are as + // without the second PVT region, and the brine column starts from the same + // contact. + const std::vector eqlnum(20, 0); + std::vector pvtnum(20, 0); + std::fill(pvtnum.begin() + 10, pvtnum.end(), 1); + const std::vector connate(20, connateSw); + const std::vector maxWater(20, 1.0); + const WaterEquilFixture fix(twoWaterRegionsDeck("20*1", "10*1 10*2")); + const auto states = fix.compute(eqlnum, connate, maxWater, pvtnum).fluidStates(); + const auto first = fix.compute(eqlnum, connate, maxWater).fluidStates(); + + const auto pressure = [](const auto& fs, const unsigned phaseIdx) { + return Opm::getValue(fs.pressure(phaseIdx)); + }; + constexpr auto water = WaterFluidSystem::waterPhaseIdx; + for (std::size_t c = 0; c < states.size(); ++c) { + BOOST_TEST_CONTEXT("Cell " << c) { + BOOST_CHECK_EQUAL(pressure(states[c], WaterFluidSystem::oilPhaseIdx), + pressure(first[c], WaterFluidSystem::oilPhaseIdx)); + if (c < 10) { + BOOST_CHECK_EQUAL(pressure(states[c], water), pressure(first[c], water)); + } + } + } + + // Cells 18 and 19 hold brine. Cell 10 lies 17.5 m below the contact, where + // the brine exceeds the fresh water by the weight of the denser column. + const Scalar pBrine = 0.5 * (pressure(states[18], water) + pressure(states[19], water)); + BOOST_CHECK_CLOSE(impliedDensity(pressure(states[18], water), pressure(states[19], water)), + waterDensity(1150.0, 1.01, 2.5e-5, pBrine / barsa), 0.01); + const Scalar pCell10 = 0.5 * (pressure(states[10], water) + pressure(first[10], water)); + BOOST_CHECK_CLOSE(pressure(states[10], water) - pressure(first[10], water), + (waterDensity(1150.0, 1.01, 2.5e-5, pCell10 / barsa) + - waterDensity(1000.0, 1.03, 4.0e-5, pCell10 / barsa)) + * gravity * (fix.depths[10] - 2035.0), + 0.1); +} + +BOOST_AUTO_TEST_CASE(DatumInDifferentWatersIsRejected) +{ + // A datum below the contact gives the pressure of both the fresh water and + // the brine, which cannot both reach the contact at one pressure. + const WaterEquilFixture fix(twoWaterRegionsDeck("20*1", "10*1 10*2", + "EQUIL\n" + " 2045 150 2035 0 2000 0 /\n" + " 2060 155 2085 0 2050 0 /\n")); + std::vector pvtnum(20, 0); + std::fill(pvtnum.begin() + 10, pvtnum.end(), 1); + BOOST_CHECK_EXCEPTION(fix.compute(std::vector(20, 0), + std::vector(20, connateSw), + std::vector(20, 1.0), pvtnum), + std::runtime_error, + [](const std::runtime_error& error) { + return std::string_view{error.what()}.find( + "the water of its PVT regions 1 and 2 reaches the contact") + != std::string_view::npos; + }); +} + +BOOST_AUTO_TEST_CASE(DatumInTheSameWaterIsAccepted) +{ + // PVT regions that share their water, as when they differ in the rock only, + // equilibrate as one. + const WaterEquilFixture fix(twoWaterRegionsDeck("20*1", "10*1 10*2", + "EQUIL\n" + " 2045 150 2035 0 2000 0 /\n" + " 2060 155 2085 0 2050 0 /\n", + "PVTW\n" + " 75.0 1.03 4.0E-5 0.3 0.0 /\n" + " 75.0 1.03 4.0E-5 0.3 0.0 /\n" + "DENSITY\n" + " 800.0 1000.0 1.0 /\n" + " 800.0 1000.0 1.0 /\n")); + const std::vector eqlnum(20, 0); + std::vector pvtnum(20, 0); + std::fill(pvtnum.begin() + 10, pvtnum.end(), 1); + const std::vector connate(20, connateSw); + const std::vector maxWater(20, 1.0); + const auto states = fix.compute(eqlnum, connate, maxWater, pvtnum).fluidStates(); + const auto single = fix.compute(eqlnum, connate, maxWater).fluidStates(); + for (std::size_t c = 0; c < states.size(); ++c) { + for (const auto phaseIdx : {WaterFluidSystem::oilPhaseIdx, + WaterFluidSystem::waterPhaseIdx}) { + BOOST_CHECK_EQUAL(Opm::getValue(states[c].pressure(phaseIdx)), + Opm::getValue(single[c].pressure(phaseIdx))); + } + } +} + +BOOST_AUTO_TEST_CASE(PvtRegionsMayVaryWithoutWater) +{ + // Without water the PVT regions take no part in the equilibration. + const EquilFixture fix(withTwoPvtRegions(deckString( + "EQUIL\n 2010 150 2300 0 2000 0 /\n", "EQLDIMS\n/\n", + "REGIONS\nPVTNUM\n10*1 10*2 /\n"))); + std::vector pvtnum(20, 0); + std::fill(pvtnum.begin() + 10, pvtnum.end(), 1); + const auto states = fix.compute(std::vector(20, 0), {}, {}, pvtnum).fluidStates(); + const auto single = fix.compute(std::vector(20, 0)).fluidStates(); + for (std::size_t c = 0; c < states.size(); ++c) { + BOOST_CHECK_EQUAL(Opm::getValue(states[c].pressure(FluidSystem::oilPhaseIdx)), + Opm::getValue(single[c].pressure(FluidSystem::oilPhaseIdx))); + } +} + +BOOST_AUTO_TEST_CASE(InvalidPvtnumIsRejected) +{ + const WaterEquilFixture fix(twoWaterRegionsDeck("10*1 10*2", "10*1 10*2")); + std::vector regions(20, 0); + std::fill(regions.begin() + 10, regions.end(), 1); + auto pvtnum = regions; + pvtnum.back() = 2; + BOOST_CHECK_THROW(fix.compute(regions, std::vector(20, connateSw), + std::vector(20, 1.0), pvtnum), + std::runtime_error); +} + BOOST_AUTO_TEST_CASE(CoincidentContactsKeepTheGasRoot) { // Contacts that coincide leave no liquid: every hydrocarbon cell sits above From 023e083623f43456240e6b321e31038fe093eaf3 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 29 Sep 2026 13:14:02 +0200 Subject: [PATCH 09/10] Register compositional PVT region regressions --- regressionTests.cmake | 51 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/regressionTests.cmake b/regressionTests.cmake index 1f9ba0d241c..db011a27ffe 100644 --- a/regressionTests.cmake +++ b/regressionTests.cmake @@ -252,6 +252,57 @@ add_test_compareECLFiles( compositional/equilibration ) +add_test_compareECLFiles( + CASENAME + pvtnum_water + FILENAME + PVTNUM_WATER + SIMULATOR + flow_comp + REFERENCE_SIMULATOR + flow_comp + ABS_TOL + ${abs_tol} + REL_TOL + ${rel_tol} + DIR + compositional/pvt_regions +) + +add_test_compareECLFiles( + CASENAME + pvtnum_equil + FILENAME + PVTNUM_EQUIL + SIMULATOR + flow_comp + REFERENCE_SIMULATOR + flow_comp + ABS_TOL + ${abs_tol} + REL_TOL + ${rel_tol} + DIR + compositional/pvt_regions +) + +add_test_compareECLFiles( + CASENAME + pvtnum_rock + FILENAME + PVTNUM_ROCK + SIMULATOR + flow_comp + REFERENCE_SIMULATOR + flow_comp + ABS_TOL + ${abs_tol} + REL_TOL + ${rel_tol} + DIR + compositional/pvt_regions +) + add_test_compareECLFiles( CASENAME spe12 From 474a2e9d146e4e7e2740a7cb5a6b87490909abee Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Sun, 4 Oct 2026 23:48:56 +0200 Subject: [PATCH 10/10] Register shared equilibration PVT region regression --- regressionTests.cmake | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/regressionTests.cmake b/regressionTests.cmake index db011a27ffe..32fbb0c686b 100644 --- a/regressionTests.cmake +++ b/regressionTests.cmake @@ -286,6 +286,23 @@ add_test_compareECLFiles( compositional/pvt_regions ) +add_test_compareECLFiles( + CASENAME + pvtnum_equil_shared + FILENAME + PVTNUM_EQUIL_SHARED + SIMULATOR + flow_comp + REFERENCE_SIMULATOR + flow_comp + ABS_TOL + ${abs_tol} + REL_TOL + ${rel_tol} + DIR + compositional/pvt_regions +) + add_test_compareECLFiles( CASENAME pvtnum_rock