From 172d62e93a0665f2e1e9251a621566c8e80d76b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Wed, 10 Dec 2025 18:46:57 +0100 Subject: [PATCH 1/2] Add test for checkGroupConstraintsProd efficiency factor bug This test exposes a bug in GroupStateHelper::checkGroupConstraintsProd() where the addback calculation uses the accumulated efficiency factor instead of the partial efficiency factor at local_reduction_level. Test setup: - Group hierarchy: FIELD -> PLAT (E=0.9) -> MANI (E=0.8) / MANI2 (E=0.7) - MANI has individual ORAT control, PLAT has guide rate - WELL-A produces 10000 SM3/day under MANI The bug causes: - Addback uses 0.72 (E_MANI * E_PLAT) instead of 0.8 (E_MANI) - This results in scale = 0.9465 (false constraint violation) - Correct calculation would give scale = 1.0288 (no violation) The test currently passes with the buggy expected values. A subsequent commit will fix the bug and update the test expectations. --- CMakeLists.txt | 1 + CMakeLists_files.cmake | 2 + tests/GROUP_HIGHER_CONSTRAINTS.DATA | 230 +++++++++++++ tests/SimulatorFixture.hpp | 125 +++++++ tests/test_glift1.cpp | 90 +----- tests/test_group_higher_constraints.cpp | 414 ++++++++++++++++++++++++ 6 files changed, 780 insertions(+), 82 deletions(-) create mode 100644 tests/GROUP_HIGHER_CONSTRAINTS.DATA create mode 100644 tests/SimulatorFixture.hpp create mode 100644 tests/test_group_higher_constraints.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index c1eb0b1737..3f78ddf9f5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -509,6 +509,7 @@ endif() target_sources(test_outputdir PRIVATE $) target_sources(test_equil PRIVATE $) +target_sources(test_group_higher_constraints PRIVATE $) target_sources(test_RestartSerialization PRIVATE $) target_sources(test_glift1 PRIVATE $) if(MPI_FOUND) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 979502e40f..f87e313d1d 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -476,6 +476,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_convergencereport.cpp tests/test_deferredlogger.cpp tests/test_dilu.cpp + tests/test_group_higher_constraints.cpp tests/test_equil.cpp tests/test_extractMatrix.cpp tests/test_flexiblesolver.cpp @@ -669,6 +670,7 @@ list (APPEND TEST_DATA_FILES tests/rhs3rep.txt tests/options_flexiblesolver.json tests/options_flexiblesolver_simple.json + tests/GROUP_HIGHER_CONSTRAINTS.DATA tests/GLIFT1.DATA tests/RC-01_MAST_PRED.DATA tests/include/flowl_b_vfp.ecl diff --git a/tests/GROUP_HIGHER_CONSTRAINTS.DATA b/tests/GROUP_HIGHER_CONSTRAINTS.DATA new file mode 100644 index 0000000000..aac49e666f --- /dev/null +++ b/tests/GROUP_HIGHER_CONSTRAINTS.DATA @@ -0,0 +1,230 @@ +-- 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/ + +-- Test case for efficiency factor handling in checkGroupConstraintsProd() + +-- Group hierarchy: +-- +-- FIELD (ORAT = 10000) +-- | +-- PLAT (GEFAC=0.9, FLD, has GUIDERAT) +-- / \ +-- / \ +-- MANI (GEFAC=0.8) MANI2 (GEFAC=0.7) +-- ORAT control FLD control +-- | | +-- WELL-A WELL-B +-- +-- This setup triggers checkGroupConstraintsProd() for MANI with: +-- - local_reduction_level at PLAT (because PLAT has guide rate) +-- - efficiency_factor = 0.9 * 0.8 = 0.72 + +RUNSPEC + +TITLE +EFFICIENCY_FACTOR_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 + 4 5 5 4 / + +TABDIMS +/ + +EQLDIMS +/ + +UNIFIN +UNIFOUT + +-- ============================================================================= +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' / + 'MANI' 'PLAT' / + 'MANI2' 'PLAT' / +/ + +-- Well specifications +WELSPECS + 'WELL-A' 'MANI' 3 3 2010 'OIL' / + 'WELL-B' 'MANI2' 2 2 2010 'OIL' / +/ + +COMPDAT + 'WELL-A' 3 3 1 3 'OPEN' 1* 1* 0.2 / + 'WELL-B' 2 2 1 3 'OPEN' 1* 1* 0.2 / +/ + +-- Group efficiency factors +GEFAC + 'PLAT' 0.9 / + 'MANI' 0.8 / + 'MANI2' 0.7 / +/ + +-- Well efficiency factors +WEFAC + 'WELL-A' 1.0 / + 'WELL-B' 1.0 / +/ + +-- Group controls +-- FIELD has overall oil rate target +GCONPROD +'FIELD' 'ORAT' 10000 1* 1* 1* 'RATE' / +/ + +-- PLAT is under FLD control with guide rate enabled ('YES' and 'FORM') +-- Having a guide rate on PLAT triggers local_reduction_level at PLAT in checkGroupConstraintsProd() +GCONPROD +'PLAT' 'FLD' 1* 1* 1* 1* 'RATE' 'YES' 1* 'FORM' / +/ + +-- Guide rate formula (applies to all groups/wells using 'FORM' target) +GUIDERAT +-- min.int type A B C D E F allow.inc damp.factor use.free.gas + 0 'OIL' 1.0 0.0 1.0 1.0 0.0 0.0 'YES' 1.0 'NO' / + +-- MANI has individual ORAT control (this triggers checkGroupConstraintsProd) +GCONPROD +'MANI' 'ORAT' 5000 1* 1* 1* 'RATE' / +/ + +-- MANI2 has FLD control - this creates a sibling for MANI to share PLAT's allocation via guide rates +GCONPROD +'MANI2' 'FLD' 1* 1* 1* 1* 'RATE' / +/ + +-- Wells under group control +WCONPROD + 'WELL-A' 'OPEN' 'GRUP' 1* 1* 1* 1* 1* 100.0 / + 'WELL-B' 'OPEN' 'GRUP' 1* 1* 1* 1* 1* 100.0 / +/ + +TSTEP +1 / + +END diff --git a/tests/SimulatorFixture.hpp b/tests/SimulatorFixture.hpp new file mode 100644 index 0000000000..54d60bc5a5 --- /dev/null +++ b/tests/SimulatorFixture.hpp @@ -0,0 +1,125 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \brief Common fixtures and utilities for simulator-based tests + */ +#ifndef OPM_SIMULATOR_FIXTURE_HPP +#define OPM_SIMULATOR_FIXTURE_HPP + +#include "config.h" +#include "TestTypeTag.hpp" + +#include +#include +#include + +#include +#include +#include + +#if HAVE_DUNE_FEM +#include +#else +#include +#endif + +#include +#include +#if BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 < 71 +#include +#else +#include +#endif + +#include +#include + +namespace Opm { + +/*! + * \brief Global fixture for simulator-based tests + * + * Initializes MPI and sets up communication. Use with BOOST_GLOBAL_FIXTURE. + */ +struct SimulatorFixture +{ + SimulatorFixture() + { + int argc = boost::unit_test::framework::master_test_suite().argc; + char** argv = boost::unit_test::framework::master_test_suite().argv; +#if HAVE_DUNE_FEM + Dune::Fem::MPIManager::initialize(argc, argv); +#else + Dune::MPIHelper::instance(argc, argv); +#endif + FlowGenericVanguard::setCommunication(std::make_unique()); + } +}; + +/*! + * \brief Initialize a simulator from an Eclipse deck file + * + * \tparam TypeTag The type tag for the simulator + * \param filename Path to the Eclipse deck file + * \param test_name Name of the test (used in argv[0]) + * \param threads_per_process Number of threads per process (default: 1) + * \return Unique pointer to the initialized simulator + */ +template +std::unique_ptr> +initSimulator(const char* filename, + const char* test_name = "test_simulator", + int threads_per_process = 1) +{ + using Simulator = GetPropType; + + std::string filename_arg = "--ecl-deck-file-name="; + filename_arg += filename; + + const char* argv[] = { + test_name, + filename_arg.c_str() + }; + + Parameters::reset(); + registerAllParameters_(false); + registerEclTimeSteppingParameters(); + BlackoilModelParameters::registerParameters(); + Parameters::Register("Do *NOT* use!"); + Parameters::SetDefault(threads_per_process); + Parameters::endRegistration(); + setupParameters_(/*argc=*/sizeof(argv) / sizeof(argv[0]), + argv, + /*registerParams=*/false, + /*allowUnused=*/false, + /*handleHelp=*/true, + /*myRank=*/0); + + FlowGenericVanguard::readDeck(filename); + return std::make_unique(); +} + +} // namespace Opm + +#endif // OPM_SIMULATOR_FIXTURE_HPP diff --git a/tests/test_glift1.cpp b/tests/test_glift1.cpp index 2b92703278..3d645f23fa 100644 --- a/tests/test_glift1.cpp +++ b/tests/test_glift1.cpp @@ -20,20 +20,15 @@ module for the precise wording of the license and the list of copyright holders. */ -#include "config.h" -#include "TestTypeTag.hpp" #define BOOST_TEST_MODULE Glift1 -#include -#include -#include +#include "SimulatorFixture.hpp" #include #include #include #include -#include #include #include #include @@ -45,87 +40,18 @@ #include #include -#if HAVE_DUNE_FEM -#include -#else -#include -#endif - -#include #include #include #include -#include -#include -#if BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 < 71 -#include -#else -#include -#endif - -namespace Opm::Properties { - namespace TTag { - struct TestGliftTypeTag { - using InheritsFrom = std::tuple; - }; - } -} - -template -std::unique_ptr> -initSimulator(const char *filename) -{ - using namespace Opm; - using Simulator = GetPropType; - - std::string filename_arg = "--ecl-deck-file-name="; - filename_arg += filename; - - const char* argv[] = { - "test_equil", - filename_arg.c_str() +namespace Opm::Properties::TTag { + struct TestGliftTypeTag { + using InheritsFrom = std::tuple; }; - - Parameters::reset(); - registerAllParameters_(false); - registerEclTimeSteppingParameters(); - BlackoilModelParameters::registerParameters(); - Parameters::Register("Do *NOT* use!"); - Opm::Parameters::SetDefault(2); - Parameters::endRegistration(); - setupParameters_(/*argc=*/sizeof(argv) / sizeof(argv[0]), - argv, - /*registerParams=*/false, - /*allowUnused*/false, - /*handleHelp*/true, - /*myRank*/0); - - FlowGenericVanguard::readDeck(filename); - return std::make_unique(); } - -namespace { - -struct GliftFixture -{ - GliftFixture() - { - int argc = boost::unit_test::framework::master_test_suite().argc; - char** argv = boost::unit_test::framework::master_test_suite().argv; - #if HAVE_DUNE_FEM - Dune::Fem::MPIManager::initialize(argc, argv); - #else - Dune::MPIHelper::instance(argc, argv); -#endif - Opm::FlowGenericVanguard::setCommunication(std::make_unique()); - } -}; - -} - -BOOST_GLOBAL_FIXTURE(GliftFixture); +using SimulatorFixture = Opm::SimulatorFixture; +BOOST_GLOBAL_FIXTURE(SimulatorFixture); BOOST_AUTO_TEST_CASE(G1) { @@ -145,7 +71,7 @@ BOOST_AUTO_TEST_CASE(G1) const std::string filename = "GLIFT1.DATA"; using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups; - auto simulator = initSimulator(filename.data()); + auto simulator = Opm::initSimulator(filename.data(), "test_glift1", /*threads=*/2); simulator->model().applyInitialSolution(); simulator->setEpisodeIndex(-1); @@ -296,7 +222,7 @@ BOOST_AUTO_TEST_CASE(G2) const std::string filename = "GLIFT1.DATA"; using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups; - auto simulator = initSimulator(filename.data()); + auto simulator = Opm::initSimulator(filename.data(), "test_glift1", /*threads=*/2); simulator->model().applyInitialSolution(); simulator->setEpisodeIndex(-1); diff --git a/tests/test_group_higher_constraints.cpp b/tests/test_group_higher_constraints.cpp new file mode 100644 index 0000000000..3f92d2fa40 --- /dev/null +++ b/tests/test_group_higher_constraints.cpp @@ -0,0 +1,414 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \brief Test for efficiency factor handling in checkGroupConstraintsProd() + * + * This test verifies the behavior of GroupStateHelper::checkGroupConstraintsProd() + * when checking higher-level group constraints with efficiency factors (GEFAC). + * + * Group hierarchy used in the test: + * \code + * FIELD (ORAT = 10000) + * | + * PLAT (GEFAC=0.9, FLD, has guide rate) + * / \ + * / \ + * MANI (GEFAC=0.8) MANI2 (GEFAC=0.7) + * ORAT control FLD control + * | | + * WELL-A WELL-B + * \endcode + * + * The test exposes a bug where the "addback" calculation at local_reduction_level + * uses the accumulated efficiency factor (E_MANI * E_PLAT = 0.72) instead of + * the partial efficiency factor (E_MANI = 0.8). This causes incorrect constraint + * violation detection when MANI's production should actually be within limits. + * + * Key concepts tested: + * - local_reduction_level: Set at PLAT because it has a guide rate + * - Reduction rates: PLAT's reduction = E_MANI * MANI_rate (groups with individual control) + * - Addback: Should restore MANI's contribution using E_MANI, not accumulated efficiency + * - Guide rate fractions: How MANI and MANI2 share PLAT's allocation + */ + +#define BOOST_TEST_MODULE GroupHigherConstraints +#include "SimulatorFixture.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace Opm::Properties::TTag { + struct TestEfficiencyFactorTypeTag { + using InheritsFrom = std::tuple; + }; +} + +using SimulatorFixture = Opm::SimulatorFixture; +BOOST_GLOBAL_FIXTURE(SimulatorFixture); + +namespace { + +//! \brief Test fixture with unit conversion helpers and tolerance checks +struct ToleranceAndUnitFixture +{ + //! Relative tolerance for exact mathematical relationships (e.g., efficiency factors) + static constexpr double tight_tol = 1e-12; + + //! Relative tolerance for rate values with minor numerical differences + static constexpr double rate_tol = 1e-8; + + //! Relative tolerance for algorithm results (e.g., scale factors) + static constexpr double algo_tol = 1e-3; + + //! Convert SI rate (m³/s) to metric rate (SM3/day) + static double metric_rate(double si_rate) + { + using namespace Opm::unit; + return convert::to(si_rate, cubic(meter) / day); + } + + //! Convert metric rate (SM3/day) to SI rate (m³/s) + static double si_rate(double metric_rate) + { + using namespace Opm::unit; + return convert::from(metric_rate, cubic(meter) / day); + } + + //! Check that two values are close within tight tolerance (for exact relationships) + static void checkClose(double actual, double expected) + { + BOOST_CHECK_CLOSE(actual, expected, tight_tol); + } + + //! Check that two rate values are close + static void checkRate(double actual, double expected) + { + BOOST_CHECK_CLOSE(actual, expected, rate_tol); + } + + //! Check that two algorithm results are close (looser tolerance) + static void checkAlgo(double actual, double expected) + { + BOOST_CHECK_CLOSE(actual, expected, algo_tol); + } +}; + +} // anonymous namespace + +BOOST_FIXTURE_TEST_SUITE(GroupHigherConstraintsTests, ToleranceAndUnitFixture) + +BOOST_AUTO_TEST_CASE(TestGroupHigherConstraints) +{ + using TypeTag = Opm::Properties::TTag::TestEfficiencyFactorTypeTag; + using WellModel = Opm::BlackoilWellModel; + using FluidSystem = Opm::GetPropType; + using IndexTraits = typename FluidSystem::IndexTraitsType; + + const std::string filename = "GROUP_HIGHER_CONSTRAINTS.DATA"; + auto simulator = Opm::initSimulator(filename.data()); + + simulator->model().applyInitialSolution(); + simulator->setEpisodeIndex(-1); + simulator->setEpisodeLength(0.0); + simulator->startNextEpisode(/*episodeStartTime=*/0.0, /*episodeLength=*/1e30); + simulator->setTimeStepSize(Opm::unit::day); + simulator->model().newtonMethod().setIterationIndex(0); + + WellModel& well_model = simulator->problem().wellModel(); + const auto& schedule = simulator->vanguard().schedule(); + const int report_step_idx = 0; + const auto& pu = well_model.phaseUsage(); + + const int oil_phase_pos = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); + + well_model.beginReportStep(report_step_idx); + + // We want predictable rates that exceed FIELD's target to trigger constraint violation: + // WELL-A: 10000 SM3/day oil (negative for production) + // WELL-B: 4000 SM3/day oil (negative for production) + // Total at FIELD level = E_PLAT * (E_MANI * 10000 + E_MANI2 * 4000) + // = 0.9 * (0.8 * 10000 + 0.7 * 4000) = 0.9 * 10800 = 9720 SM3/day + // This is close to FIELD's ORAT of 10000, so MANI's share will be constrained. + const double well_a_oil_rate = si_rate(-10000.0); // m³/s (negative = production) + const double well_b_oil_rate = si_rate(-4000.0); // m³/s + + auto& ws = well_model.wellState(); + + // Set WELL-A rates + { + auto well_idx = ws.index("WELL-A"); + BOOST_REQUIRE(well_idx.has_value()); + auto& well_a = ws.well(well_idx.value()); + std::fill(well_a.surface_rates.begin(), well_a.surface_rates.end(), 0.0); + well_a.surface_rates[oil_phase_pos] = well_a_oil_rate; + well_a.production_cmode = Opm::Well::ProducerCMode::GRUP; // Under group control + } + + // Set WELL-B rates + { + auto well_idx = ws.index("WELL-B"); + BOOST_REQUIRE(well_idx.has_value()); + auto& well_b = ws.well(well_idx.value()); + std::fill(well_b.surface_rates.begin(), well_b.surface_rates.end(), 0.0); + well_b.surface_rates[oil_phase_pos] = well_b_oil_rate; + well_b.production_cmode = Opm::Well::ProducerCMode::GRUP; // Under group control + } + + // After setting well control modes, we need to update the global well info + // so that isProductionGrup() returns true for wells under group control. + // This is needed for FractionCalculator::guideRateSum() to include wells. + std::vector well_status(ws.size(), Opm::Well::Status::OPEN); + const auto& comm = simulator->vanguard().grid().comm(); + ws.updateGlobalIsGrup(comm, well_status); + + // Set up guide rates for wells to control how rates are distributed. + // Without this, localFraction() falls back to potentials or returns 1.0. + // We use guide rates proportional to well rates for simplicity. + // + // Guide rate values are in output units (SM3/day for metric), not SI units. + // The init_grvalue() method stores values directly without conversion. + auto& guide_rate = well_model.guideRate(); + const double sim_time = 0.0; + + // Set guide rates for wells (values in SM3/day, same as their production rates) + guide_rate.init_grvalue(report_step_idx, "WELL-A", + Opm::GuideRate::GuideRateValue{sim_time, 10000.0, Opm::GuideRateModel::Target::OIL}); + guide_rate.init_grvalue(report_step_idx, "WELL-B", + Opm::GuideRate::GuideRateValue{sim_time, 4000.0, Opm::GuideRateModel::Target::OIL}); + + // Set guide rate for PLAT group - this is needed to trigger local_reduction_level + // at PLAT level in checkGroupConstraintsProd(). Without this, the algorithm + // doesn't set local_reduction_level and the addback logic at line checkGroupConstraintsProd() is skipped. + guide_rate.init_grvalue(report_step_idx, "PLAT", + Opm::GuideRate::GuideRateValue{sim_time, 9000.0, Opm::GuideRateModel::Target::OIL}); + + // The checkGroupConstraintsProd function uses production_rates() to compute + // current_rate_fraction when checking guide rate violations. We need to + // populate these for the groups in the chain. + auto& gs = well_model.groupState(); + const int num_phases = pu.numActivePhases(); + + // Get group objects and their efficiency factors + const auto& mani_group = schedule.getGroup("MANI", report_step_idx); + const auto& mani2_group = schedule.getGroup("MANI2", report_step_idx); + const auto& plat_group = schedule.getGroup("PLAT", report_step_idx); + const auto& field_group = schedule.getGroup("FIELD", report_step_idx); + + const double E_MANI = mani_group.getGroupEfficiencyFactor(); + const double E_MANI2 = mani2_group.getGroupEfficiencyFactor(); + const double E_PLAT = plat_group.getGroupEfficiencyFactor(); + + // Set production rates for groups (these are the "current" group rates) + // MANI produces 10000 SM3/day, but with efficiency factor this shows up as: + // MANI: 10000 SM3/day (direct well rate) + // PLAT: E_MANI * 10000 + E_MANI2 * 4000 = 8000 + 2800 = 10800 SM3/day + // FIELD: E_PLAT * 10800 = 0.9 * 10800 = 9720 SM3/day + std::vector mani_prod_rates(num_phases, 0.0); + mani_prod_rates[oil_phase_pos] = -well_a_oil_rate; // 10000 SM3/day (positive) + + std::vector mani2_prod_rates(num_phases, 0.0); + mani2_prod_rates[oil_phase_pos] = -well_b_oil_rate; // 4000 SM3/day (positive) + + std::vector plat_prod_rates(num_phases, 0.0); + plat_prod_rates[oil_phase_pos] = E_MANI * (-well_a_oil_rate) + E_MANI2 * (-well_b_oil_rate); // 10800 SM3/day + + std::vector field_prod_rates(num_phases, 0.0); + field_prod_rates[oil_phase_pos] = E_PLAT * plat_prod_rates[oil_phase_pos]; // 9720 SM3/day + + gs.update_production_rates("MANI", mani_prod_rates); + gs.update_production_rates("MANI2", mani2_prod_rates); + gs.update_production_rates("PLAT", plat_prod_rates); + gs.update_production_rates("FIELD", field_prod_rates); + + // MANI has individual ORAT control (this triggers checkGroupConstraintsProd) + gs.production_control("MANI", Opm::Group::ProductionCMode::ORAT); + + // MANI2 has FLD control (follows parent) + gs.production_control("MANI2", Opm::Group::ProductionCMode::FLD); + + // PLAT has FLD control (follows parent FIELD) + gs.production_control("PLAT", Opm::Group::ProductionCMode::FLD); + + // FIELD has ORAT control + gs.production_control("FIELD", Opm::Group::ProductionCMode::ORAT); + + Opm::DeferredLogger deferred_logger; + auto& gsh = well_model.groupStateHelper(); + + // Update the groupControlledWells count based on well control modes + gsh.updateGroupControlledWells(/*is_production_group=*/true, Opm::Phase::OIL, deferred_logger); + + gsh.updateGroupTargetReduction(field_group, /*is_injector=*/false); + + // ======================================================================== + // VERIFY: Check the test setup + // ======================================================================== + + // Verify efficiency factors from schedule + checkClose(E_MANI, 0.8); + checkClose(E_MANI2, 0.7); + checkClose(E_PLAT, 0.9); + const double accumulated_eff = E_MANI * E_PLAT; + checkClose(accumulated_eff, 0.72); + + // Verify well rates in WellState + checkRate(metric_rate(ws.well("WELL-A").surface_rates[oil_phase_pos]), -10000.0); + checkRate(metric_rate(ws.well("WELL-B").surface_rates[oil_phase_pos]), -4000.0); + + // Verify wells are marked as under group control + BOOST_CHECK(ws.isProductionGrup("WELL-A")); + BOOST_CHECK(ws.isProductionGrup("WELL-B")); + + // Verify guide rates are set + BOOST_CHECK(guide_rate.has("WELL-A")); + BOOST_CHECK(guide_rate.has("WELL-B")); + BOOST_CHECK(guide_rate.has("PLAT")); + + // Verify groupControlledWells count (1 well per leaf group) + BOOST_CHECK_EQUAL(gsh.groupControlledWells("MANI", "", true, Opm::Phase::OIL), 1); + BOOST_CHECK_EQUAL(gsh.groupControlledWells("MANI2", "", true, Opm::Phase::OIL), 1); + BOOST_CHECK_EQUAL(gsh.groupControlledWells("PLAT", "", true, Opm::Phase::OIL), 1); + + // Verify sumWellSurfaceRates for each group + // MANI: direct well rate = 10000 SM3/day + // MANI2: direct well rate = 4000 SM3/day + // PLAT: E_MANI * 10000 + E_MANI2 * 4000 = 0.8 * 10000 + 0.7 * 4000 = 10800 SM3/day + const double sum_mani = gsh.sumWellSurfaceRates(mani_group, oil_phase_pos, /*injector=*/false); + const double sum_mani2 = gsh.sumWellSurfaceRates(mani2_group, oil_phase_pos, /*injector=*/false); + const double sum_plat = gsh.sumWellSurfaceRates(plat_group, oil_phase_pos, /*injector=*/false); + checkRate(metric_rate(sum_mani), 10000.0); + checkRate(metric_rate(sum_mani2), 4000.0); + checkRate(metric_rate(sum_plat), E_MANI * 10000.0 + E_MANI2 * 4000.0); // 10800 + + // Verify FIELD's ORAT target + const auto& field_group_ctrl = field_group.productionControls(gsh.summaryState()); + checkRate(metric_rate(field_group_ctrl.oil_target), 10000.0); + + // Verify group reduction rates + // FIELD: 0 (PLAT has guide rate, so doesn't contribute to FIELD's reduction) + // PLAT: E_MANI * 10000 = 8000 (MANI has individual ORAT control) + // MANI, MANI2: 0 (leaf groups with no sub-groups having individual control) + const auto& field_red = gs.production_reduction_rates("FIELD"); + const auto& plat_red = gs.production_reduction_rates("PLAT"); + BOOST_CHECK_SMALL(metric_rate(field_red[oil_phase_pos]), tight_tol); + checkRate(metric_rate(plat_red[oil_phase_pos]), E_MANI * 10000.0); // 8000 + + // Verify guide rate fraction for MANI + // MANI's guide rate = E_MANI * well_guide_rate = 0.8 * 10000 = 8000 + // MANI2's guide rate = E_MANI2 * well_guide_rate = 0.7 * 4000 = 2800 + // Total = 10800, MANI fraction = 8000/10800 = 20/27 + const double expected_mani_fraction = (E_MANI * 10000.0) / (E_MANI * 10000.0 + E_MANI2 * 4000.0); + checkRate(expected_mani_fraction, 20.0 / 27.0); // 0.740740740... + + // Prepare rates array for MANI (negative for production) + // This represents the "current rate" that will be compared against the target + std::vector mani_rates(num_phases, 0.0); + mani_rates[oil_phase_pos] = well_a_oil_rate; // -10000 SM3/day in SI units + + // Prepare resv_coeff (reservoir volume coefficients) - use 1.0 for simplicity + std::vector resv_coeff(num_phases, 1.0); + + // Initial efficiency factor is E_MANI (will accumulate E_PLAT during recursion) + const double initial_efficiency = E_MANI; + + // Call checkGroupConstraintsProd + // Arguments: name, parent, group, rates, efficiency_factor, resv_coeff, check_guide_rate, deferred_logger + auto [constraint_violated, scale] = gsh.checkGroupConstraintsProd( + "MANI", // name (the group/well being checked) + "PLAT", // parent (immediate parent group name) + plat_group, // group (PLAT - will recurse to FIELD since PLAT has FLD) + mani_rates.data(), // rates (MANI's production rates) + initial_efficiency, // efficiency_factor (E_MANI, will accumulate E_PLAT) + resv_coeff, // resv_coeff + /*check_guide_rate=*/true, // check_guide_rate + deferred_logger // deferred_logger + ); + + // ======================================================================== + // EXPECTED RESULT (with current buggy implementation) + // ======================================================================== + // + // The checkGroupConstraintsProd algorithm proceeds as follows: + // + // Setup: + // - Chain: [FIELD, PLAT, MANI] + // - efficiency_factor = E_MANI * E_PLAT = 0.8 * 0.9 = 0.72 + // - current_rate_available = 10000 SM3/day (MANI's production) + // - orig_target = 10000 SM3/day (FIELD's ORAT) + // - local_reduction_level = 1 (PLAT has guide rate) + // + // Loop iteration ii=0 (FIELD): + // target = 10000 + // local_reduction("FIELD") = 0 (PLAT has guide rate) + // target = 10000 - 0 = 10000 + // local_fraction("PLAT") = 1.0 (only child of FIELD) + // target = 10000 * 1.0 = 10000 + // + // Loop iteration ii=1 (PLAT): + // target = 10000 + // local_reduction("PLAT") = E_MANI * 10000 = 8000 (MANI has ORAT control) + // target = 10000 - 8000 = 2000 + // + // BUG: Addback uses accumulated efficiency (0.72) instead of E_MANI (0.8): + // addback = current_rate * efficiency_factor = 10000 * 0.72 = 7200 + // target = 2000 + 7200 = 9200 + // + // (Correct addback would be: 10000 * 0.8 = 8000, giving target = 10000) + // + // local_fraction("MANI") = 8000/10800 = 0.740741 + // target = 9200 * 0.740741 = 6814.81 + // + // Final calculation (line 413): + // target_rate_available = target / efficiency_factor + // = (184000/27) / 0.72 = 9465.0205761... SM3/day + // + // Scale factor: + // scale = target_rate_available / current_rate_available + // = 9465.0205761... / 10000 = 0.94650205761... + // + // Because scale < 1, the buggy implementation incorrectly constrains MANI, + // reducing its production even though the correct calculation would give + // scale = 1.0288 (no constraint violation). + // + // Exact value: 4600000 / 486 / 10000 = 2300/243 / 5 = 460/243 + const double expected_scale_with_bug = 460.0 / 486.0; // 0.9465020576... + BOOST_CHECK(constraint_violated); // Bug causes false constraint violation + checkAlgo(scale, expected_scale_with_bug); + + // Derived values for verification + const double current_rate = 10000.0; // SM3/day + const double target_rate_available_with_bug = scale * current_rate; + checkAlgo(target_rate_available_with_bug, 4600000.0 / 486.0); // 9465.0205761... +} + +BOOST_AUTO_TEST_SUITE_END() From da720e1900e169d4cf64b7e115db653fe4178397 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Sat, 13 Dec 2025 10:02:02 +0100 Subject: [PATCH 2/2] Fix efficiency factor addback in group constraint checking The addback calculation at local_reduction_level was using the accumulated efficiency factor (from entity to control group) instead of the partial efficiency factor (from entity to local_reduction_level). This caused incorrect constraint violations when checking higher-level group constraints. The fix computes the partial efficiency by multiplying efficiency factors from local_reduction_level down to the entity. --- opm/simulators/wells/GroupStateHelper.cpp | 44 ++++++++++++++++++++-- opm/simulators/wells/GroupStateHelper.hpp | 13 +++++++ tests/test_group_higher_constraints.cpp | 45 +++++++++++------------ 3 files changed, 74 insertions(+), 28 deletions(-) diff --git a/opm/simulators/wells/GroupStateHelper.cpp b/opm/simulators/wells/GroupStateHelper.cpp index cddb2cc15f..2893574490 100644 --- a/opm/simulators/wells/GroupStateHelper.cpp +++ b/opm/simulators/wells/GroupStateHelper.cpp @@ -232,7 +232,9 @@ GroupStateHelper::checkGroupConstraintsInj(const std::strin // Add my reduction back at the level where it is included in the local reduction if (local_reduction_level == ii) { - target += current_rate_available * efficiency_factor; + const Scalar addback_efficiency + = this->computeAddbackEfficiency_(chain, local_reduction_level); + target += current_rate_available * addback_efficiency; } } target *= local_fraction_lambda(chain[ii + 1]); @@ -404,7 +406,9 @@ GroupStateHelper::checkGroupConstraintsProd(const std::stri } // Add my reduction back at the level where it is included in the local reduction if (local_reduction_level == ii) { - target += current_rate_available * efficiency_factor; + const Scalar addback_efficiency + = this->computeAddbackEfficiency_(chain, local_reduction_level); + target += current_rate_available * addback_efficiency; } } target *= local_fraction_lambda(chain[ii + 1]); @@ -603,7 +607,9 @@ GroupStateHelper::getWellGroupTargetInjector(const std::str // If we are under individual control we need to add the wells rate back at the level where it is // included in the local reduction if (local_reduction_level == ii && !this->wellState().isInjectionGrup(name)) { - target += current_rate_available * efficiency_factor; + const Scalar addback_efficiency + = this->computeAddbackEfficiency_(chain, local_reduction_level); + target += current_rate_available * addback_efficiency; } } if (!this->wellState().isInjectionGrup(name)) { @@ -731,7 +737,9 @@ GroupStateHelper::getWellGroupTargetProducer(const std::str // If we are under individual control we need to add the wells rate back at the level where it is // included in the local reduction if (local_reduction_level == ii && !this->wellState().isProductionGrup(name)) { - target += current_rate_available * efficiency_factor; + const Scalar addback_efficiency + = this->computeAddbackEfficiency_(chain, local_reduction_level); + target += current_rate_available * addback_efficiency; } } if (this->wellState().isProductionGrup(name)) { @@ -1381,6 +1389,34 @@ GroupStateHelper::worstOffendingWell(const Group& group, // Private methods // --------------------------------------------------------------------- +template +Scalar +GroupStateHelper::computeAddbackEfficiency_( + const std::vector& chain, + const std::size_t local_reduction_level) const +{ + // Compute partial efficiency factor from local_reduction_level down to the entity. + // Chain is ordered [control_group, ..., local_reduction_level, ..., entity]. + // We multiply efficiency factors from index (local_reduction_level + 1) to the entity. + // Note: The entity at chain[num_ancestors] may be a well or a group. + // Both wells (WEFAC) and groups (GEFAC) have efficiency factors that must be included, + // since the local_reduction rates include these efficiency factors. + const std::size_t num_ancestors = chain.size() - 1; + Scalar efficiency = 1.0; + for (std::size_t jj = local_reduction_level + 1; jj <= num_ancestors; ++jj) { + const std::string& name = chain[jj]; + if (this->schedule_.hasGroup(name, this->report_step_)) { + const auto& grp = this->schedule_.getGroup(name, this->report_step_); + efficiency *= grp.getGroupEfficiencyFactor(); + } else if (this->schedule_.hasWell(name, this->report_step_)) { + const auto& well = this->schedule_.getWell(name, this->report_step_); + efficiency *= well.getEfficiencyFactor() + * this->wellState()[name].efficiency_scaling_factor; + } + } + return efficiency; +} + template std::string GroupStateHelper::controlGroup_(const Group& group) const diff --git a/opm/simulators/wells/GroupStateHelper.hpp b/opm/simulators/wells/GroupStateHelper.hpp index d9b1664588..8d912b299a 100644 --- a/opm/simulators/wells/GroupStateHelper.hpp +++ b/opm/simulators/wells/GroupStateHelper.hpp @@ -320,6 +320,19 @@ public: DeferredLogger& deferred_logger) const; private: + //! \brief Compute partial efficiency factor for addback calculation. + //! + //! The addback in constraint checking must use the partial efficiency factor + //! (from local_reduction_level down to the entity), not the accumulated efficiency + //! (from the entity to the control group). This helper multiplies efficiency + //! factors from local_reduction_level+1 to the entity. + //! + //! \param chain The group chain from control group (top) to entity (bottom) + //! \param local_reduction_level The level at which addback is applied + //! \return The partial efficiency factor for addback + Scalar computeAddbackEfficiency_(const std::vector& chain, + std::size_t local_reduction_level) const; + std::string controlGroup_(const Group& group) const; GuideRate::RateVector getGuideRateVector_(const std::vector& rates) const; diff --git a/tests/test_group_higher_constraints.cpp b/tests/test_group_higher_constraints.cpp index 3f92d2fa40..ca569bb8b3 100644 --- a/tests/test_group_higher_constraints.cpp +++ b/tests/test_group_higher_constraints.cpp @@ -40,10 +40,10 @@ * WELL-A WELL-B * \endcode * - * The test exposes a bug where the "addback" calculation at local_reduction_level - * uses the accumulated efficiency factor (E_MANI * E_PLAT = 0.72) instead of - * the partial efficiency factor (E_MANI = 0.8). This causes incorrect constraint - * violation detection when MANI's production should actually be within limits. + * The test verifies that the "addback" calculation at local_reduction_level + * correctly uses the partial efficiency factor (E_MANI = 0.8) instead of + * the accumulated efficiency factor (E_MANI * E_PLAT = 0.72). This ensures + * correct constraint violation detection. * * Key concepts tested: * - local_reduction_level: Set at PLAT because it has a guide rate @@ -355,7 +355,7 @@ BOOST_AUTO_TEST_CASE(TestGroupHigherConstraints) ); // ======================================================================== - // EXPECTED RESULT (with current buggy implementation) + // EXPECTED RESULT (with fixed implementation) // ======================================================================== // // The checkGroupConstraintsProd algorithm proceeds as follows: @@ -379,36 +379,33 @@ BOOST_AUTO_TEST_CASE(TestGroupHigherConstraints) // local_reduction("PLAT") = E_MANI * 10000 = 8000 (MANI has ORAT control) // target = 10000 - 8000 = 2000 // - // BUG: Addback uses accumulated efficiency (0.72) instead of E_MANI (0.8): - // addback = current_rate * efficiency_factor = 10000 * 0.72 = 7200 - // target = 2000 + 7200 = 9200 + // Addback uses partial efficiency E_MANI (not accumulated 0.72): + // addback = current_rate * E_MANI = 10000 * 0.8 = 8000 + // target = 2000 + 8000 = 10000 // - // (Correct addback would be: 10000 * 0.8 = 8000, giving target = 10000) + // local_fraction("MANI") = 8000/10800 = 20/27 + // target = 10000 * 20/27 = 200000/27 // - // local_fraction("MANI") = 8000/10800 = 0.740741 - // target = 9200 * 0.740741 = 6814.81 - // - // Final calculation (line 413): + // Final calculation: // target_rate_available = target / efficiency_factor - // = (184000/27) / 0.72 = 9465.0205761... SM3/day + // = (200000/27) / 0.72 = 200000/27 * 25/18 = 5000000/486 + // = 10288.0658... SM3/day // // Scale factor: // scale = target_rate_available / current_rate_available - // = 9465.0205761... / 10000 = 0.94650205761... + // = (5000000/486) / 10000 = 500/486 = 250/243 = 1.02880658... // - // Because scale < 1, the buggy implementation incorrectly constrains MANI, - // reducing its production even though the correct calculation would give - // scale = 1.0288 (no constraint violation). + // Because scale > 1, no constraint violation - MANI can produce its full rate. // - // Exact value: 4600000 / 486 / 10000 = 2300/243 / 5 = 460/243 - const double expected_scale_with_bug = 460.0 / 486.0; // 0.9465020576... - BOOST_CHECK(constraint_violated); // Bug causes false constraint violation - checkAlgo(scale, expected_scale_with_bug); + // Exact value: 250/243 + const double expected_scale_with_fix = 250.0 / 243.0; // 1.02880658... + BOOST_CHECK(!constraint_violated); // No constraint violation with fix + checkAlgo(scale, expected_scale_with_fix); // Derived values for verification const double current_rate = 10000.0; // SM3/day - const double target_rate_available_with_bug = scale * current_rate; - checkAlgo(target_rate_available_with_bug, 4600000.0 / 486.0); // 9465.0205761... + const double target_rate_available_with_fix = scale * current_rate; + checkAlgo(target_rate_available_with_fix, 5000000.0 / 486.0); // 10288.0658... } BOOST_AUTO_TEST_SUITE_END()