From 504c8f62e52441f3009c0434ae0155b33e433ca3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Mon, 16 Feb 2026 11:04:08 +0100 Subject: [PATCH 1/2] Add tests for autochoke group target rate Add two test cases for getAutoChokeGroupProductionTargetRate(): - TestAutoChokeGroupProductionTargetRate: normal case with MANI (autochoke, FLD) under intermediate group PLAT that has a guide rate. Verifies the autochoke reduction loop produces the same result as applyReductionsAndFractions_() when both are applicable. - TestAutoChokeGroupTargetRateDoubleCountingBug: uses a NETWORK DATA file with NODEPROP (as_choke=YES) and lower well rates to naturally trigger underperformance detection (GCW=0). Documents the double-counting bug where missing local_reduction_level gating causes reductions to be subtracted at both FIELD and PLAT, producing a clamped- to-zero target instead of the correct value. Also adds GROUP_HIGHER_CONSTRAINTS_NETWORK.DATA with BRANPROP/NODEPROP keywords and parameterizes the test fixture with mani2_control, data_file, and well rates. --- CMakeLists_files.cmake | 1 + tests/GROUP_HIGHER_CONSTRAINTS_NETWORK.DATA | 278 +++++++++++++++++ tests/test_group_higher_constraints.cpp | 314 +++++++++++++++++++- 3 files changed, 589 insertions(+), 4 deletions(-) create mode 100644 tests/GROUP_HIGHER_CONSTRAINTS_NETWORK.DATA diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 57f85d8efc..60cd76a32b 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -679,6 +679,7 @@ list (APPEND TEST_DATA_FILES tests/GCONSUMP.DATA tests/GCONSUMP_COMPLEX.DATA tests/GROUP_HIGHER_CONSTRAINTS.DATA + tests/GROUP_HIGHER_CONSTRAINTS_NETWORK.DATA tests/GLIFT1.DATA tests/RC-01_MAST_PRED.DATA tests/include/flowl_b_vfp.ecl diff --git a/tests/GROUP_HIGHER_CONSTRAINTS_NETWORK.DATA b/tests/GROUP_HIGHER_CONSTRAINTS_NETWORK.DATA new file mode 100644 index 0000000000..4683840084 --- /dev/null +++ b/tests/GROUP_HIGHER_CONSTRAINTS_NETWORK.DATA @@ -0,0 +1,278 @@ +-- This reservoir simulation deck is made available under the Open Database +-- License: http://opendatacommons.org/licenses/odbl/1.0/. Any rights in +-- individual contents of the database are licensed under the Database Contents +-- License: http://opendatacommons.org/licenses/dbcl/1.0/ + +-- Extended network variant of GROUP_HIGHER_CONSTRAINTS.DATA +-- +-- Adds NETWORK, BRANPROP and NODEPROP keywords to enable the extended network +-- model. MANI is designated as a sub-sea manifold (autochoke group) via +-- NODEPROP, which sets as_choke()=true on the MANI group. +-- +-- Group hierarchy: +-- +-- FIELD (ORAT = 10000, terminal pressure = 200 barsa) +-- / \ +-- / \ +-- PLAT (GEFAC=0.9) PLAT-2 (GEFAC=0.85) +-- / \ | +-- / \ WELL-C +-- MANI (0.8,SSM) MANI-2 (0.75) +-- / \ | +-- WELL-A WELL-A2 WELL-B +-- +-- Network branches (all VFP=9999, no pressure drop): +-- MANI -> PLAT (autochoke branch) +-- MANI-2 -> PLAT +-- PLAT -> FIELD +-- PLAT-2 -> FIELD +-- +-- NODEPROP effects on MANI (as_choke=YES): +-- - group.as_choke() = true +-- - Wells switched to THP control (ProducerCMode::THP) +-- - Wells' guide rate availability disabled +-- +RUNSPEC + +TITLE +EFFICIENCY_FACTOR_NETWORK_TEST + +DIMENS + 5 5 3 / + +OIL +WATER +GAS +DISGAS + +METRIC + +START + 1 'JAN' 2020 / + +WELLDIMS +-- max wells, max cons/well, max groups, max wells/group + 6 5 6 4 / + +TABDIMS +/ + +EQLDIMS +/ + +UNIFIN +UNIFOUT + +NETWORK +-- max_nodes max_branches + 5 4 / + +-- ============================================================================= +GRID +-- ============================================================================= + +DXV +5*100.0 / + +DYV +5*100.0 / + +DZV +3*10.0 / + +TOPS +25*2000 / + +PORO +75*0.25 / + +PERMX +75*500 / + +COPY + PERMX PERMY / + PERMX PERMZ / +/ + +MULTIPLY + PERMZ 0.1 / +/ + +-- ============================================================================= +PROPS +-- ============================================================================= + +PVTW +-- ref.pres Bw Cw Vw Vw.visc + 250.0 1.02 4.0E-05 0.5 0.0 / + +ROCK +-- ref.pres compressibility + 250.0 4.0E-05 / + +DENSITY +-- oil water gas + 850 1014 1.0 / + +PVDO +-- P Bo Vo + 100.0 1.05 2.0 + 250.0 1.02 2.5 + 400.0 1.00 3.0 / + +PVDG +-- P Bg Vg + 100.0 0.010 0.015 + 250.0 0.005 0.020 + 400.0 0.003 0.025 / + +SWOF +-- Sw Krw Kro Pcow + 0.2 0.0 1.0 0.0 + 0.3 0.05 0.6 0.0 + 0.5 0.15 0.3 0.0 + 0.7 0.35 0.1 0.0 + 0.8 0.5 0.0 0.0 / + +SGOF +-- Sg Krg Kro Pcog + 0.0 0.0 1.0 0.0 + 0.1 0.02 0.6 0.0 + 0.3 0.10 0.3 0.0 + 0.5 0.25 0.1 0.0 + 0.7 0.5 0.0 0.0 / + +-- ============================================================================= +SOLUTION +-- ============================================================================= + +EQUIL +-- datum pressure OWC Pc@OWC GOC Pc@GOC RSVD RVVD SOLN + 2000 250.0 2100 0.0 1900 0.0 1 0 0 / + +RSVD + 1900 100.0 + 2100 100.0 / + +-- ============================================================================= +SUMMARY +-- ============================================================================= + +WOPR +/ +WGPR +/ +WWPR +/ +WBHP +/ +GOPR +/ +FOPR +FGPR +FWPR + +-- ============================================================================= +SCHEDULE +-- ============================================================================= + +GRUPTREE + 'PLAT' 'FIELD' / + 'PLAT-2' 'FIELD' / + 'MANI' 'PLAT' / + 'MANI-2' 'PLAT' / +/ + +-- Well specifications +WELSPECS + 'WELL-A' 'MANI' 3 3 2010 'OIL' / + 'WELL-A2' 'MANI' 4 4 2010 'OIL' / + 'WELL-B' 'MANI-2' 5 5 2010 'OIL' / + 'WELL-C' 'PLAT-2' 2 2 2010 'OIL' / +/ + +COMPDAT + 'WELL-A' 3 3 1 3 'OPEN' 1* 1* 0.2 / + 'WELL-A2' 4 4 1 3 'OPEN' 1* 1* 0.2 / + 'WELL-B' 5 5 1 3 'OPEN' 1* 1* 0.2 / + 'WELL-C' 2 2 1 3 'OPEN' 1* 1* 0.2 / +/ + +-- Group efficiency factors +GEFAC + 'PLAT' 0.9 / + 'PLAT-2' 0.85 / + 'MANI' 0.8 / + 'MANI-2' 0.75 / +/ + +-- Well efficiency factors +WEFAC + 'WELL-A' 1.0 / + 'WELL-A2' 1.0 / + 'WELL-B' 1.0 / + 'WELL-C' 1.0 / +/ + +-- Group controls +-- FIELD has overall oil rate target +GCONPROD +'FIELD' 'ORAT' 10000 1* 1* 1* 'RATE' / +/ + +-- PLAT is under FLD control with explicit guide rate +GCONPROD +'PLAT' 'FLD' 1* 1* 1* 1* 'RATE' 'YES' 9500 'OIL' / +/ + +-- PLAT-2 is under FLD control with explicit guide rate +GCONPROD +'PLAT-2' 'FLD' 1* 1* 1* 1* 'RATE' 'YES' 2550 'OIL' / +/ + +-- MANI is under FLD control with explicit guide rate (autochoke group) +-- Guide rate = E_MANI * (WELL-A + WELL-A2) = 0.8 * 10000 = 8000 +GCONPROD +'MANI' 'FLD' 1* 1* 1* 1* 'RATE' 'YES' 8000 'OIL' / +/ + +-- MANI-2 is under FLD control with explicit guide rate +GCONPROD +'MANI-2' 'FLD' 1* 1* 1* 1* 'RATE' 'YES' 1500 'OIL' / +/ + +-- Guide rate formula for wells (guide_rate = oil_potential) +GUIDERAT +-- min.int type A B C D E F allow.inc damp.factor use.free.gas + 0 'OIL' 1.0 1.0 0.0 0.0 0.0 0.0 'YES' 1.0 'NO' / + +-- Wells under group control (NODEPROP will switch MANI wells to THP) +WCONPROD + 'WELL-A' 'OPEN' 'GRUP' 1* 1* 1* 1* 1* 100.0 / + 'WELL-A2' 'OPEN' 'GRUP' 1* 1* 1* 1* 1* 100.0 / + 'WELL-B' 'OPEN' 'GRUP' 1* 1* 1* 1* 1* 100.0 / + 'WELL-C' 'OPEN' 'GRUP' 1* 1* 1* 1* 1* 100.0 / +/ + +-- Extended network branches (all VFP=9999: no pressure drop) +BRANPROP +-- Downtree Uptree VFP ALQ + 'MANI' 'PLAT' 9999 1* / + 'MANI-2' 'PLAT' 9999 1* / + 'PLAT' 'FIELD' 9999 1* / + 'PLAT-2' 'FIELD' 9999 1* / +/ + +-- Extended network node properties +-- FIELD: terminal node with fixed pressure (200 barsa) +-- MANI: sub-sea manifold (autochoke group, as_choke=YES) +NODEPROP +-- Node Pressure autoChoke addGasLift + 'FIELD' 200 'NO' 'NO' / + 'MANI' 1* 'YES' 'NO' / +/ + +TSTEP +1 / + +END diff --git a/tests/test_group_higher_constraints.cpp b/tests/test_group_higher_constraints.cpp index c7fa614b72..32561213b7 100644 --- a/tests/test_group_higher_constraints.cpp +++ b/tests/test_group_higher_constraints.cpp @@ -62,6 +62,7 @@ #include #include #include +#include #include @@ -92,12 +93,20 @@ struct GroupConstraintTestFixture : ToleranceAndUnitFixture bool plat_has_guide_rate, Opm::Group::ProductionCMode mani_control, Opm::Well::ProducerCMode well_a_control, - Opm::Well::ProducerCMode well_a2_control) + Opm::Well::ProducerCMode well_a2_control, + Opm::Group::ProductionCMode mani2_control = Opm::Group::ProductionCMode::FLD, + const std::string& data_file = "GROUP_HIGHER_CONSTRAINTS.DATA", + double well_a_rate_metric = 6000.0, + double well_a2_rate_metric = 4000.0) { plat_has_guide_rate_ = plat_has_guide_rate; mani_control_ = mani_control; + mani2_control_ = mani2_control; well_a_control_ = well_a_control; well_a2_control_ = well_a2_control; + data_file_ = data_file; + well_a_rate_metric_ = well_a_rate_metric; + well_a2_rate_metric_ = well_a2_rate_metric; setupSimulator_(); setupWellRates_(); @@ -208,6 +217,12 @@ private: //! Control mode for MANI group (ORAT or FLD) Opm::Group::ProductionCMode mani_control_{Opm::Group::ProductionCMode::ORAT}; + //! Control mode for MANI-2 group (FLD or ORAT) + Opm::Group::ProductionCMode mani2_control_{Opm::Group::ProductionCMode::FLD}; + + //! DATA file to load + std::string data_file_{"GROUP_HIGHER_CONSTRAINTS.DATA"}; + // Well rates (SM3/day, metric units for readability) double well_a_rate_metric_{6000.0}; double well_a2_rate_metric_{4000.0}; @@ -228,8 +243,7 @@ private: void setupSimulator_() { - const std::string filename = "GROUP_HIGHER_CONSTRAINTS.DATA"; - simulator_ = Opm::initSimulator(filename.data()); + simulator_ = Opm::initSimulator(data_file_.data()); simulator_->model().applyInitialSolution(); simulator_->setEpisodeIndex(-1); @@ -339,7 +353,7 @@ private: // Set group controls gs.production_control("MANI", mani_control_); - gs.production_control("MANI-2", Opm::Group::ProductionCMode::FLD); + gs.production_control("MANI-2", mani2_control_); gs.production_control("PLAT", Opm::Group::ProductionCMode::FLD); gs.production_control("PLAT-2", Opm::Group::ProductionCMode::FLD); gs.production_control("FIELD", Opm::Group::ProductionCMode::ORAT); @@ -971,4 +985,296 @@ BOOST_AUTO_TEST_CASE(TestWellGroupTargetProducerWithSiblingIndividualControl) checkAlgo(metric_rate(target), expected_target_metric); } +//! \brief Test getAutoChokeGroupProductionTargetRate() for autochoke group MANI +//! +//! This test exercises the autochoke target calculation which uses a simpler reduction loop +//! than applyReductionsAndFractions_(): it omits both local_reduction_level gating and addback. +//! With MANI-2 under individual ORAT control, a non-zero reduction exists at the PLAT level, +//! exposing the code-path difference. +//! +//! Code-path difference vs applyReductionsAndFractions_(): +//! - Standard algorithm: applies reduction only at ii >= local_reduction_level, plus addback +//! - Autochoke loop: applies reduction at ii==0 OR where guideRate->has(chain[ii]) is true +//! In this scenario both produce the same result because: +//! (a) The autochoke group (MANI) has no individual wells, so no addback is needed +//! (b) PLAT has a guide rate, so both algorithms apply the PLAT-level reduction +//! +//! \code +//! FIELD (ORAT = 10000) +//! ├── PLAT (E=0.9, FLD, GUIDERATE=9500, GCW=2) +//! │ ├── MANI (E=0.8, FLD, GUIDERATE=8000, GCW=2, autochoke group) +//! │ │ ├── WELL-A (RATE=6000, POT=6000, GRUP) +//! │ │ └── WELL-A2 (RATE=4000, POT=4000, GRUP) +//! │ └── MANI-2 (E=0.75, IDV, ORAT, GUIDERATE=1500) +//! │ └── WELL-B (RATE=2000, POT=2000, GRUP) +//! └── PLAT-2 (E=0.85, FLD, GUIDERATE=2550, GCW=1) +//! └── WELL-C (RATE=3000, POT=3000, GRUP) +//! \endcode +BOOST_AUTO_TEST_CASE(TestAutoChokeGroupProductionTargetRate) +{ + // Configuration: MANI under FLD control (autochoke group follows parent), + // both wells under GRUP control, MANI-2 under individual ORAT control. + // MANI-2 under ORAT creates a non-zero reduction at PLAT level (MANI-2's + // wells contribute) which is the key ingredient that differentiates the + // autochoke loop from applyReductionsAndFractions_(). + commonSetupForAllTests( + /*plat_has_guide_rate=*/true, + /*mani_control =*/Opm::Group::ProductionCMode::FLD, + /*well_a_control =*/Opm::Well::ProducerCMode::GRUP, + /*well_a2_control =*/Opm::Well::ProducerCMode::GRUP, + /*mani2_control =*/Opm::Group::ProductionCMode::ORAT + ); + + auto& gs = groupState(); + auto& gsh = groupStateHelper(); + + // Verify MANI-2 is under individual ORAT control + BOOST_CHECK(gs.production_control("MANI-2") == Opm::Group::ProductionCMode::ORAT); + BOOST_CHECK(gs.production_control("MANI") == Opm::Group::ProductionCMode::FLD); + + // Verify reduction rates: + // - PLAT reduction: MANI is FLD with guide rate → skip. + // MANI-2 is ORAT (individual) → add E_MANI2 * wellB = 0.75 * 2000 = 1500 + // - FIELD reduction: PLAT has guide rate and GCW > 0 → skip. PLAT-2 likewise → 0 + { + const auto& plat_red = gs.production_reduction_rates("PLAT"); + checkRate(metric_rate(plat_red[oilPhasePos()]), eMani2() * mani2TotalRateMetric()); // 1500 + const auto& field_red = gs.production_reduction_rates("FIELD"); + checkRate(metric_rate(field_red[oilPhasePos()]), 0.0); + } + + // Verify groupControlledWells: + // MANI-2 is under individual ORAT → its wells don't count as group-controlled + // for the parent, but MANI-2 itself is excluded from PLAT's FLD children + { + const auto always_included = std::string(""); + const auto is_prod = true; + const auto phase = Opm::Phase::OIL; + BOOST_CHECK_EQUAL(gsh.groupControlledWells("MANI", always_included, is_prod, phase), 2); + BOOST_CHECK_EQUAL(gsh.groupControlledWells("PLAT", always_included, is_prod, phase), 2); + BOOST_CHECK_EQUAL(gsh.groupControlledWells("FIELD", always_included, is_prod, phase), 3); + } + + // Call getAutoChokeGroupProductionTargetRate() + // This is how it's called in BlackoilWellModelNetwork: + // getAutoChokeGroupProductionTargetRate(bottom_group=MANI, group=MANI, gsh, resv_coeff, 1.0) + // The function recurses upward from MANI (FLD) → PLAT (FLD) → FIELD (ORAT=controlling). + // During recursion: efficiencyFactor accumulates E_MANI * E_PLAT = 0.8 * 0.9 = 0.72. + auto resv_coeff = resvCoeff(); + + using WGC = Opm::WellGroupControls; + auto [target_rate_si, control_mode] = WGC::getAutoChokeGroupProductionTargetRate( + maniGroup(), // bottom_group + maniGroup(), // group (starts recursion here, recurses up to controlling group) + gsh, + resv_coeff, + /*efficiencyFactor=*/1.0 + ); + + // Verify control mode is ORAT (from FIELD, the controlling group) + BOOST_CHECK(control_mode == Opm::Group::ProductionCMode::ORAT); + + // ======================================================================== + // EXPECTED COMPUTATION + // ======================================================================== + // + // The function recurses from MANI (FLD) → PLAT (FLD) → FIELD (ORAT). + // During recursion, efficiencyFactor accumulates: + // 1.0 * E_MANI * E_PLAT = 0.8 * 0.9 = 0.72 + // (Both MANI and PLAT multiply into efficiencyFactor before reaching FIELD.) + // + // chain = groupChainTopBot("MANI", "FIELD") = ["FIELD", "PLAT", "MANI"] + // num_ancestors = 2 + // orig_target = 10000 SM3/day (FIELD's ORAT target) + // + // Loop ii=0 (FIELD): + // (ii==0) → apply reduction + // localReduction("FIELD") = 0 (no reduction at FIELD level) + // target = 10000 - 0 = 10000 + // localFraction("PLAT"): PLAT=9500, PLAT-2=2550, both FLD → 9500/12050 + // target = 10000 * 9500/12050 ≈ 7884.65 + // + // Loop ii=1 (PLAT): + // guideRate->has("PLAT") = true → apply reduction + // localReduction("PLAT") = E_MANI2 * wellB_rate = 0.75 * 2000 = 1500 + // target = 7884.65 - 1500 = 6384.65 + // localFraction("MANI"): MANI is FLD with GCW=2, MANI-2 is ORAT (excluded) → 1.0 + // target = 6384.65 * 1.0 = 6384.65 + // + // Final: target / efficiencyFactor = 6384.65 / 0.72 ≈ 8867.57 SM3/day + // + // Code-path difference: in applyReductionsAndFractions_(), local_reduction_level + // gating and addback would apply at the intermediate PLAT level. Here, the autochoke + // loop applies reduction wherever guideRate->has() is true (PLAT), without addback. + // Both algorithms produce the same result here because MANI has no individual wells + // (reduction and addback would cancel if present). + // + const double eff = eMani() * ePlat(); // 0.72 + const double plat_fraction = platGuideRateMetric() + / (platGuideRateMetric() + plat2GuideRateMetric()); // 9500/12050 + const double plat_reduction_metric = eMani2() * mani2TotalRateMetric(); // 1500 + const double field_target_metric = 10000.0; + const double target_after_field = field_target_metric * plat_fraction; // ≈ 7884.65 + const double target_after_plat = target_after_field - plat_reduction_metric; // ≈ 6384.65 + const double mani_fraction = 1.0; // MANI is only FLD child at PLAT level + const double expected_target_metric = (target_after_plat * mani_fraction) + / eff; // ≈ 8867.57 + + const double target_rate_metric = metric_rate(target_rate_si); + checkAlgo(target_rate_metric, expected_target_metric); +} + +//! \brief Test getAutoChokeGroupProductionTargetRate() with GCW=0 (underperformance scenario) +//! +//! Uses the NETWORK DATA file where MANI is a sub-sea manifold (as_choke=true via NODEPROP). +//! NODEPROP switches MANI's wells to THP control and disables their guide rate availability. +//! +//! MANI's well rates are set low enough (total=5000) that isAutoChokeGroupUnderperforming_() +//! detects underperformance and updateGroupControlledWells() naturally sets MANI GCW=0. +//! +//! This test demonstrates the double-counting bug in the autochoke loop when an intermediate +//! group has a guide rate but GCW=0. The autochoke loop subtracts reductions unconditionally +//! at every level where guideRate->has() is true, without local_reduction_level gating. +//! When GCW=0 at MANI and PLAT, updateGroupTargetReductionRecursive_() propagates MANI's +//! total production into PLAT's reduction (case 1), which then propagates into FIELD's +//! reduction (case 1 again). The autochoke loop subtracts reduction(FIELD) at ii=0 and +//! reduction(PLAT) at ii=1, double-counting MANI's production. +//! +//! BUG: The buggy code gives target=0 (double-counted reductions make the target negative, +//! clamped to 0 by std::max at WellGroupControls.cpp:406). +//! FIX: With local_reduction_level gating, the correct target is ~7013.89 SM3/day. +//! +//! \code +//! FIELD (ORAT = 10000, GCW=1) +//! ├── PLAT (E=0.9, FLD, GUIDERATE=9500, GCW=0) +//! │ ├── MANI (E=0.8, FLD, GUIDERATE=8000, GCW=0, as_choke, underperforming) +//! │ │ ├── WELL-A (RATE=3000, POT=3000, THP) +//! │ │ └── WELL-A2 (RATE=2000, POT=2000, THP) +//! │ └── MANI-2 (E=0.75, IDV, ORAT, GUIDERATE=1500, GCW=1) +//! │ └── WELL-B (RATE=2000, POT=2000, GRUP) +//! └── PLAT-2 (E=0.85, FLD, GUIDERATE=2550, GCW=1) +//! └── WELL-C (RATE=3000, POT=3000, GRUP) +//! \endcode +BOOST_AUTO_TEST_CASE(TestAutoChokeGroupTargetRateDoubleCountingBug) +{ + // Uses NETWORK DATA file: MANI has as_choke()=true from NODEPROP, + // wells under MANI are switched to THP control by the NODEPROP handler. + // Lower MANI well rates (3000+2000=5000) trigger underperformance detection + // in isAutoChokeGroupUnderperforming_(), which naturally sets MANI GCW=0. + commonSetupForAllTests( + /*plat_has_guide_rate=*/true, + /*mani_control =*/Opm::Group::ProductionCMode::FLD, + /*well_a_control =*/Opm::Well::ProducerCMode::THP, + /*well_a2_control =*/Opm::Well::ProducerCMode::THP, + /*mani2_control =*/Opm::Group::ProductionCMode::ORAT, + /*data_file =*/"GROUP_HIGHER_CONSTRAINTS_NETWORK.DATA", + /*well_a_rate_metric =*/3000.0, + /*well_a2_rate_metric=*/2000.0 + ); + + // Verify NODEPROP correctly set as_choke() on MANI + BOOST_CHECK(maniGroup().as_choke()); + BOOST_CHECK(!platGroup().as_choke()); + + auto& gs = groupState(); + auto& gsh = groupStateHelper(); + + // finalizeSetup_() called updateGroupControlledWells() which runs + // isAutoChokeGroupUnderperforming_(): MANI (total=5000) is below its target + // share, so all wells are excluded → MANI GCW=0. + // With MANI GCW=0 and MANI-2 under ORAT (not FLD), PLAT GCW also becomes 0. + // + // Verify GCW values (naturally 0 from underperformance detection) + { + const auto always_included = std::string(""); + const auto is_prod = true; + const auto phase = Opm::Phase::OIL; + BOOST_CHECK_EQUAL(gsh.groupControlledWells("MANI", always_included, is_prod, phase), 0); + BOOST_CHECK_EQUAL(gsh.groupControlledWells("PLAT", always_included, is_prod, phase), 0); + BOOST_CHECK_EQUAL(gsh.groupControlledWells("FIELD", always_included, is_prod, phase), 1); + } + + // finalizeSetup_() called updateGroupTargetReduction() with + // GCW=0 at MANI and PLAT. updateGroupTargetReductionRecursive_() uses + // case 1 (total production propagated up) instead of case 3 (nothing): + // + // reduction(MANI) = 0 (wells are THP, and !group.as_choke() is false) + // + // reduction(PLAT): + // MANI (FLD, GCW=0): case 1 → E_MANI * sumWellSurfaceRates(MANI) = 0.8 * 5000 = 4000 + // MANI-2 (ORAT): case 1 → E_MANI2 * sumWellSurfaceRates(MANI-2) = 0.75 * 2000 = 1500 + // reduction(PLAT) = 5500 + // + // reduction(FIELD): + // PLAT (FLD, GCW=0): case 1 → E_PLAT * sumWellSurfaceRates(PLAT) = 0.9 * 5500 = 4950 + // PLAT-2 (FLD, GCW>0, guide rate): case 3 → nothing + // reduction(FIELD) = 4950 + // + // Verify reduction rates + { + const auto& plat_red = gs.production_reduction_rates("PLAT"); + const double plat_red_metric = eMani() * maniTotalRateMetric() + + eMani2() * mani2TotalRateMetric(); // 4000 + 1500 = 5500 + checkRate(metric_rate(plat_red[oilPhasePos()]), plat_red_metric); + + const auto& field_red = gs.production_reduction_rates("FIELD"); + const double field_red_metric = ePlat() * platTotalRateMetric(); // 0.9 * 5500 = 4950 + checkRate(metric_rate(field_red[oilPhasePos()]), field_red_metric); + } + + // Call getAutoChokeGroupProductionTargetRate + auto resv_coeff = resvCoeff(); + using WGC = Opm::WellGroupControls; + auto [target_rate_si, control_mode] = WGC::getAutoChokeGroupProductionTargetRate( + maniGroup(), // bottom_group + maniGroup(), // group (starts recursion here) + gsh, + resv_coeff, + /*efficiencyFactor=*/1.0 + ); + + BOOST_CHECK(control_mode == Opm::Group::ProductionCMode::ORAT); + + // ======================================================================== + // EXPECTED COMPUTATION (buggy code, no local_reduction_level gating) + // ======================================================================== + // + // chain = [FIELD, PLAT, MANI], efficiencyFactor = E_MANI * E_PLAT = 0.72 + // + // ii=0 (FIELD): + // reduction(FIELD) = 4950 + // target = 10000 - 4950 = 5050 + // localFraction("PLAT") = 1.0 (PLAT GCW=0 → excluded from sum, + // but PLAT-2 is the only active sibling → num_active=1 → return 1.0) + // target = 5050 * 1.0 = 5050 + // + // ii=1 (PLAT): + // guideRate->has("PLAT") = true → apply reduction [BUG: should be gated] + // reduction(PLAT) = 5500 + // target = 5050 - 5500 = -450 ← DOUBLE-COUNTED! + // localFraction("MANI") = 1.0 (MANI GCW=0, MANI-2 ORAT → num_active=0, + // total=0, potentials=0 → return 1.0) + // target = -450 * 1.0 = -450 + // + // Final: max(0, -450 / 0.72) = 0.0 + // + // The double-counting: reduction(FIELD) = 4950 already includes + // E_PLAT * sumWellSurfaceRates(PLAT) = 0.9 * 5500 = 4950 + // which accounts for ALL of PLAT's production. Then reduction(PLAT) = 5500 + // subtracts MANI's production AGAIN (via E_MANI * 5000 = 4000). + // + // With local_reduction_level gating (commit 2 fix): + // local_reduction_level = 0 (PLAT has GR but GCW=0 → doesn't qualify) + // ii=0: 0 >= 0 → true → target = 5050 + // ii=1: 0 >= 1 → false → skip reduction → target stays 5050 + // Final: max(0, 5050 / 0.72) ≈ 7013.89 + // + const double target_rate_metric = metric_rate(target_rate_si); + + // BUG: double-counting makes target negative, clamped to 0 by std::max + // After fix in commit 2, this should be updated to expect ~7013.89 + const double expected_buggy_target = 0.0; + BOOST_CHECK_EQUAL(target_rate_metric, expected_buggy_target); +} + BOOST_AUTO_TEST_SUITE_END() From 9b9864551d40c677a919243ce866bd267265f04d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Tue, 17 Feb 2026 10:40:11 +0100 Subject: [PATCH 2/2] Fix autochoke target rate and clean up guide rate Add local_reduction_level gating to the autochoke target loop in getAutoChokeGroupProductionTargetRate(), preventing double-counting of reductions at intermediate groups with group-controlled-wells count of zero. Remove getGuideRate() from GroupStateHelper: the group path is dead code (only caller was isAutoChokeGroupUnderperforming_ which now uses guide_rate_ directly), and the well path is inlined into FractionCalculator::guideRate(). Improve isAutoChokeGroupUnderperforming_() to compute the control group guide rate as the sum of its FLD children's guide rates (matching how targets are distributed), and accumulate efficiency factors from the autochoke group up to the control group. This makes the formula correct for the common case where the control group is the immediate parent with no reductions. Rename getLocalReductionLevel_ to getLocalReductionLevel and make it public, as it is now called from WellGroupControls.cpp. Update test to verify the corrected target rate and remove references to the now-fixed double-counting issue. --- opm/simulators/wells/FractionCalculator.cpp | 5 +- opm/simulators/wells/GroupStateHelper.cpp | 111 ++++++++------------ opm/simulators/wells/GroupStateHelper.hpp | 29 +++-- opm/simulators/wells/WellGroupControls.cpp | 14 ++- tests/test_group_higher_constraints.cpp | 51 +++------ 5 files changed, 90 insertions(+), 120 deletions(-) diff --git a/opm/simulators/wells/FractionCalculator.cpp b/opm/simulators/wells/FractionCalculator.cpp index eaa95fffb5..9b16f1cc04 100644 --- a/opm/simulators/wells/FractionCalculator.cpp +++ b/opm/simulators/wells/FractionCalculator.cpp @@ -168,7 +168,10 @@ guideRate(const std::string& name, const bool always_use_potentials) { if (schedule_.hasWell(name, report_step_)) { - return this->groupStateHelper().getGuideRate(name, target_); + if (guide_rate_->has(name) || guide_rate_->hasPotentials(name)) { + return guide_rate_->get(name, target_, this->groupStateHelper().getWellRateVector(name)); + } + return 0.0; } else { if (groupControlledWells(name, always_included_child) > 0) { if (is_producer_ && guide_rate_->has(name) && !always_use_potentials) { diff --git a/opm/simulators/wells/GroupStateHelper.cpp b/opm/simulators/wells/GroupStateHelper.cpp index 15a19423a4..cd428444e1 100644 --- a/opm/simulators/wells/GroupStateHelper.cpp +++ b/opm/simulators/wells/GroupStateHelper.cpp @@ -189,7 +189,7 @@ GroupStateHelper::checkGroupConstraintsInj(const std::strin const Scalar orig_target = tcalc.groupTarget(); const Scalar current_rate_available = tcalc.calcModeRateFromRates(rates); - const std::size_t local_reduction_level = this->getLocalReductionLevel_( + const std::size_t local_reduction_level = this->getLocalReductionLevel( chain, /*is_production_group=*/false, injection_phase); const Scalar target = this->applyReductionsAndFractions_(chain, @@ -330,7 +330,7 @@ GroupStateHelper::checkGroupConstraintsProd(const std::stri return std::make_pair(current_well_rate_available > group_target_rate_available, scale); } - const std::size_t local_reduction_level = this->getLocalReductionLevel_( + const std::size_t local_reduction_level = this->getLocalReductionLevel( chain, /*is_production_group=*/true, /*injection_phase=*/Phase::OIL); const Scalar target = this->applyReductionsAndFractions_(chain, @@ -353,62 +353,6 @@ GroupStateHelper::checkGroupConstraintsProd(const std::stri return std::make_pair(current_rate_available > target_rate_available, scale); } -template -Scalar -GroupStateHelper::getGuideRate(const std::string& name, - const GuideRateModel::Target target) const -{ - if (this->schedule_.hasWell(name, this->report_step_)) { - if (this->guide_rate_.has(name) || this->guide_rate_.hasPotentials(name)) { - return this->guide_rate_.get(name, target, this->getWellRateVector(name)); - } else { - return 0.0; - } - } - - if (this->guide_rate_.has(name)) { - return this->guide_rate_.get(name, target, this->getProductionGroupRateVector(name)); - } - - Scalar total_guide_rate = 0.0; - const Group& group = this->schedule_.getGroup(name, this->report_step_); - - for (const std::string& group_name : group.groups()) { - const Group::ProductionCMode& current_group_control - = this->groupState().production_control(group_name); - if (current_group_control == Group::ProductionCMode::FLD - || current_group_control == Group::ProductionCMode::NONE) { - // accumulate from sub wells/groups - total_guide_rate += this->getGuideRate(group_name, target); - } - } - - for (const std::string& well_name : group.wells()) { - const auto& well_tmp = this->schedule_.getWell(well_name, this->report_step_); - - if (well_tmp.isInjector()) - continue; - - const auto well_index = this->wellState().index(well_name); - if (!well_index.has_value()) - continue; - - const auto& ws = this->wellState().well(well_index.value()); - if (ws.status == Well::Status::SHUT) - continue; - - if (!this->wellState().isProductionGrup(well_name)) - continue; - - // Only count wells under group control or the ru - if (!this->wellState().isProductionGrup(well_name)) - continue; - - total_guide_rate += this->getGuideRate(well_name, target); - } - return total_guide_rate; -} - template Scalar GroupStateHelper:: @@ -671,7 +615,7 @@ GroupStateHelper::getWellGroupTargetInjector(const std::str const Scalar current_rate_available = tcalc.calcModeRateFromRates(rates); const auto chain = this->groupChainTopBot(name, group.name()); - const std::size_t local_reduction_level = this->getLocalReductionLevel_( + const std::size_t local_reduction_level = this->getLocalReductionLevel( chain, /*is_production_group=*/false, injection_phase); // Add-back is only performed when the well is under individual control (not GRUP) @@ -767,7 +711,7 @@ GroupStateHelper::getWellGroupTargetProducer(const std::str const Scalar current_rate_available = -tcalc.calcModeRateFromRates(rates); const auto chain = this->groupChainTopBot(name, group.name()); - const std::size_t local_reduction_level = this->getLocalReductionLevel_( + const std::size_t local_reduction_level = this->getLocalReductionLevel( chain, /*is_production_group=*/true, /*injection_phase=*/Phase::OIL); // Add-back is only performed when the well is under individual control (not GRUP) @@ -1596,7 +1540,7 @@ GroupStateHelper::getGuideRateVector_(const std::vector std::size_t -GroupStateHelper::getLocalReductionLevel_(const std::vector& chain, +GroupStateHelper::getLocalReductionLevel(const std::vector& chain, const bool is_production_group, const Phase injection_phase) const { @@ -1684,21 +1628,58 @@ GroupStateHelper::isAutoChokeGroupUnderperforming_(const Gr const auto& control_group_name = this->controlGroup_(group); const auto& control_group = this->schedule_.getGroup(control_group_name, this->report_step_); - // Calculate target rate for this auto choke group + // Calculate target rate for this auto choke group using a guide rate ratio + // formula. This is an approximation that ignores intermediate reductions. + // A more accurate calculation would use WellGroupControls::getAutoChokeGroupProductionTargetRate(), + // but that requires reduction rates not yet available during + // updateGroupControlledWells() (circular dependency). The efficiency factor is however accounted for + // below meaning that the current target rate will be correct for the special case where the control + // group has no reductions and the auto choke group is a direct child of the control group. + // TODO: The issue with missing handling of reductions should be fixed by refactoring the auto choke + // target calculation below to use the more accurate + // WellGroupControls::getAutoChokeGroupProductionTargetRate() function. const auto num_phases = this->numPhases(); std::vector resv_coeff(num_phases, 1.0); GroupStateHelpers::TargetCalculator tcalc{*this, resv_coeff, control_group}; const auto& control_group_target = tcalc.groupTarget(); - const auto& control_group_guide_rate - = this->getGuideRate(control_group_name, tcalc.guideTargetMode()); + + // Sum guide rates of the control group's FLD children. The control group's + // own guide rate (from GCONPROD) is NOT used here — it governs the control + // group's share of its parent's target, not how its own target is distributed. + Scalar control_group_guide_rate = 0.0; + for (const std::string& child : control_group.groups()) { + const auto child_ctrl = this->groupState().production_control(child); + if (child_ctrl == Group::ProductionCMode::FLD + || child_ctrl == Group::ProductionCMode::NONE) { + if (this->guide_rate_.has(child)) { + control_group_guide_rate += this->guide_rate_.get( + child, tcalc.guideTargetMode(), + this->getProductionGroupRateVector(child)); + } + } + } if (control_group_guide_rate <= 0) { return false; } - const Scalar target_rate = control_group_target * group_guide_rate / control_group_guide_rate; + // Accumulate group efficiency factors from the autochoke group up to (but not + // including) the control group. This matches what + // getAutoChokeGroupProductionTargetRate() does during its recursion. + Scalar accumulated_efficiency = 1.0; + { + std::string current = group.name(); + while (current != control_group_name) { + const auto& g = this->schedule_.getGroup(current, this->report_step_); + accumulated_efficiency *= g.getGroupEfficiencyFactor(); + current = g.parent(); + } + } + + const Scalar target_rate = control_group_target * group_guide_rate + / control_group_guide_rate / accumulated_efficiency; // Calculate current rate for this group std::vector rates(num_phases, 0.0); diff --git a/opm/simulators/wells/GroupStateHelper.hpp b/opm/simulators/wells/GroupStateHelper.hpp index e9cb45f533..773327fab5 100644 --- a/opm/simulators/wells/GroupStateHelper.hpp +++ b/opm/simulators/wells/GroupStateHelper.hpp @@ -234,8 +234,6 @@ public: std::vector getGroupRatesAvailableForHigherLevelControl(const Group& group, const bool is_injector) const; - Scalar getGuideRate(const std::string& name, const GuideRateModel::Target target) const; - Scalar getInjectionGroupTarget(const Group& group, const Phase& injection_phase, const std::vector& resv_coeff) const; @@ -245,6 +243,20 @@ public: /// @return The corresponding GuideRateModel::Target for the injection phase GuideRateModel::Target getInjectionGuideTargetMode(Phase injection_phase) const; + //! \brief Find the local reduction level in a group chain. + //! + //! The local reduction level is the deepest level in the chain (starting from level 1) + //! where a group has both a guide rate and group-controlled wells (GCW > 0). + //! This level determines where reductions are applied and add-back occurs. + //! + //! \param chain Group chain from control group (top) to bottom group + //! \param is_production_group True for production, false for injection + //! \param injection_phase Phase for injection groups (ignored for production) + //! \return The local reduction level (0 if no intermediate group qualifies) + std::size_t getLocalReductionLevel(const std::vector& chain, + bool is_production_group, + Phase injection_phase) const; + Scalar getProductionGroupTarget(const Group& group) const; /// @brief Get the guide rate target mode for a production group @@ -525,19 +537,6 @@ private: GuideRate::RateVector getGuideRateVector_(const std::vector& rates) const; - //! \brief Find the local reduction level in a group chain. - //! - //! The local reduction level is the deepest level in the chain (starting from level 1) - //! where a group has both a guide rate and group-controlled wells (GCW > 0). - //! This level determines where reductions are applied and add-back occurs. - //! - //! \param chain Group chain from control group (top) to bottom group - //! \param is_production_group True for production, false for injection - //! \param injection_phase Phase for injection groups (ignored for production) - //! \return The local reduction level (0 if no intermediate group qualifies) - std::size_t getLocalReductionLevel_(const std::vector& chain, - bool is_production_group, - Phase injection_phase) const; #ifdef RESERVOIR_COUPLING_ENABLED Scalar getReservoirCouplingMasterGroupRate_(const Group& group, diff --git a/opm/simulators/wells/WellGroupControls.cpp b/opm/simulators/wells/WellGroupControls.cpp index 437da22743..4636d9cec2 100644 --- a/opm/simulators/wells/WellGroupControls.cpp +++ b/opm/simulators/wells/WellGroupControls.cpp @@ -392,13 +392,19 @@ getAutoChokeGroupProductionTargetRate(const Group& bottom_group, const auto chain = groupStateHelper.groupChainTopBot(bottom_group.name(), group.name()); // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. const std::size_t num_ancestors = chain.size() - 1; + const std::size_t local_reduction_level = groupStateHelper.getLocalReductionLevel( + chain, /*is_production_group=*/true, Phase::OIL); Scalar target = orig_target; for (std::size_t ii = 0; ii < num_ancestors; ++ii) { if ((ii == 0) || guideRate->has(chain[ii])) { - // Apply local reductions only at the control level - // (top) and for levels where we have a specified - // group guide rate. - target -= localReduction(chain[ii]); + // Apply local reductions only at the control level (top) and for + // intermediate levels where we have both a guide rate and + // group-controlled wells (local_reduction_level >= ii). Without + // this gating, reductions from a group with GCW=0 would be + // double-counted. + if (local_reduction_level >= ii) { + target -= localReduction(chain[ii]); + } } target *= localFraction(chain[ii+1]); } diff --git a/tests/test_group_higher_constraints.cpp b/tests/test_group_higher_constraints.cpp index 32561213b7..fe535324c8 100644 --- a/tests/test_group_higher_constraints.cpp +++ b/tests/test_group_higher_constraints.cpp @@ -1131,18 +1131,12 @@ BOOST_AUTO_TEST_CASE(TestAutoChokeGroupProductionTargetRate) //! //! MANI's well rates are set low enough (total=5000) that isAutoChokeGroupUnderperforming_() //! detects underperformance and updateGroupControlledWells() naturally sets MANI GCW=0. -//! -//! This test demonstrates the double-counting bug in the autochoke loop when an intermediate -//! group has a guide rate but GCW=0. The autochoke loop subtracts reductions unconditionally -//! at every level where guideRate->has() is true, without local_reduction_level gating. -//! When GCW=0 at MANI and PLAT, updateGroupTargetReductionRecursive_() propagates MANI's -//! total production into PLAT's reduction (case 1), which then propagates into FIELD's -//! reduction (case 1 again). The autochoke loop subtracts reduction(FIELD) at ii=0 and -//! reduction(PLAT) at ii=1, double-counting MANI's production. -//! -//! BUG: The buggy code gives target=0 (double-counted reductions make the target negative, -//! clamped to 0 by std::max at WellGroupControls.cpp:406). -//! FIX: With local_reduction_level gating, the correct target is ~7013.89 SM3/day. +//! With GCW=0 at MANI and PLAT, updateGroupTargetReductionRecursive_() uses case 1 (total +//! production propagated up) instead of case 3 (guide rate distribution). The +//! local_reduction_level gating in getAutoChokeGroupProductionTargetRate() ensures that +//! the PLAT-level reduction is skipped (since PLAT has GCW=0 and therefore +//! local_reduction_level=0 < ii=1), avoiding double-counting of reductions that were +//! already propagated into FIELD's reduction via case 1. //! //! \code //! FIELD (ORAT = 10000, GCW=1) @@ -1155,7 +1149,7 @@ BOOST_AUTO_TEST_CASE(TestAutoChokeGroupProductionTargetRate) //! └── PLAT-2 (E=0.85, FLD, GUIDERATE=2550, GCW=1) //! └── WELL-C (RATE=3000, POT=3000, GRUP) //! \endcode -BOOST_AUTO_TEST_CASE(TestAutoChokeGroupTargetRateDoubleCountingBug) +BOOST_AUTO_TEST_CASE(TestAutoChokeGroupTargetRateWithUnderperformance) { // Uses NETWORK DATA file: MANI has as_choke()=true from NODEPROP, // wells under MANI are switched to THP control by the NODEPROP handler. @@ -1236,12 +1230,14 @@ BOOST_AUTO_TEST_CASE(TestAutoChokeGroupTargetRateDoubleCountingBug) BOOST_CHECK(control_mode == Opm::Group::ProductionCMode::ORAT); // ======================================================================== - // EXPECTED COMPUTATION (buggy code, no local_reduction_level gating) + // EXPECTED COMPUTATION // ======================================================================== // // chain = [FIELD, PLAT, MANI], efficiencyFactor = E_MANI * E_PLAT = 0.72 + // local_reduction_level = 0 (PLAT has guide rate but GCW=0 → doesn't qualify) // // ii=0 (FIELD): + // local_reduction_level(0) >= ii(0) → true → apply reduction // reduction(FIELD) = 4950 // target = 10000 - 4950 = 5050 // localFraction("PLAT") = 1.0 (PLAT GCW=0 → excluded from sum, @@ -1249,32 +1245,17 @@ BOOST_AUTO_TEST_CASE(TestAutoChokeGroupTargetRateDoubleCountingBug) // target = 5050 * 1.0 = 5050 // // ii=1 (PLAT): - // guideRate->has("PLAT") = true → apply reduction [BUG: should be gated] - // reduction(PLAT) = 5500 - // target = 5050 - 5500 = -450 ← DOUBLE-COUNTED! + // local_reduction_level(0) >= ii(1) → false → skip reduction + // (PLAT's reduction was already propagated into FIELD's via case 1) // localFraction("MANI") = 1.0 (MANI GCW=0, MANI-2 ORAT → num_active=0, // total=0, potentials=0 → return 1.0) - // target = -450 * 1.0 = -450 + // target = 5050 * 1.0 = 5050 // - // Final: max(0, -450 / 0.72) = 0.0 - // - // The double-counting: reduction(FIELD) = 4950 already includes - // E_PLAT * sumWellSurfaceRates(PLAT) = 0.9 * 5500 = 4950 - // which accounts for ALL of PLAT's production. Then reduction(PLAT) = 5500 - // subtracts MANI's production AGAIN (via E_MANI * 5000 = 4000). - // - // With local_reduction_level gating (commit 2 fix): - // local_reduction_level = 0 (PLAT has GR but GCW=0 → doesn't qualify) - // ii=0: 0 >= 0 → true → target = 5050 - // ii=1: 0 >= 1 → false → skip reduction → target stays 5050 - // Final: max(0, 5050 / 0.72) ≈ 7013.89 + // Final: max(0, 5050 / 0.72) ≈ 7013.89 // const double target_rate_metric = metric_rate(target_rate_si); - - // BUG: double-counting makes target negative, clamped to 0 by std::max - // After fix in commit 2, this should be updated to expect ~7013.89 - const double expected_buggy_target = 0.0; - BOOST_CHECK_EQUAL(target_rate_metric, expected_buggy_target); + const double expected_target = 5050.0 / 0.72; // ≈ 7013.89 + checkAlgo(target_rate_metric, expected_target); } BOOST_AUTO_TEST_SUITE_END()