OilfieldToolsNet/opm-simulators
4
0 Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-10-26 15:40:09 -05:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #6681 from hakonhagland/fix_eff_factor

Browse Source 
Bugfix: Use correct efficiency factor when adding back local reduction contribution in group constraint checking
This commit is contained in:
Markus Blatt authored and GitHub committed 2025-12-19 17:31:16 +01:00
2 parents 17e35dc14e da720e1900
commit 3eee2885a2
8 files changed
+830 -86
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

CMakeLists.txt
+1
View File
@@ -509,6 +509,7 @@ endif()
target_sources(test_outputdir PRIVATE $<TARGET_OBJECTS:moduleVersion>)
target_sources(test_equil PRIVATE $<TARGET_OBJECTS:moduleVersion>)
target_sources(test_group_higher_constraints PRIVATE $<TARGET_OBJECTS:moduleVersion>)
target_sources(test_RestartSerialization PRIVATE $<TARGET_OBJECTS:moduleVersion>)
target_sources(test_glift1 PRIVATE $<TARGET_OBJECTS:moduleVersion>)
if(MPI_FOUND)
CMakeLists_files.cmake
+2
View File
@@ -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
opm/simulators/wells/GroupStateHelper.cpp
+40 -4
View File
@@ -216,7 +216,9 @@ GroupStateHelper<Scalar, IndexTraits>::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]);
@@ -374,7 +376,9 @@ GroupStateHelper<Scalar, IndexTraits>::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]);
@@ -671,7 +675,9 @@ GroupStateHelper<Scalar, IndexTraits>::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)) {
@@ -785,7 +791,9 @@ GroupStateHelper<Scalar, IndexTraits>::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)) {
@@ -1464,6 +1472,34 @@ GroupStateHelper<Scalar, IndexTraits>::worstOffendingWell(const Group& group,
// Private methods
// ---------------------------------------------------------------------
template <typename Scalar, typename IndexTraits>
Scalar
GroupStateHelper<Scalar, IndexTraits>::computeAddbackEfficiency_(
const std::vector<std::string>& 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 <typename Scalar, typename IndexTraits>
std::string
GroupStateHelper<Scalar, IndexTraits>::controlGroup_(const Group& group) const
opm/simulators/wells/GroupStateHelper.hpp
+13
View File
@@ -321,6 +321,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<std::string>& chain,
std::size_t local_reduction_level) const;
std::string controlGroup_(const Group& group) const;
GuideRate::RateVector getGuideRateVector_(const std::vector<Scalar>& rates) const;
tests/GROUP_HIGHER_CONSTRAINTS.DATA
+230
View File
@@ -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
tests/SimulatorFixture.hpp
+125
View File
@@ -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 <http://www.gnu.org/licenses/>.
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 <opm/models/utils/propertysystem.hh>
#include <opm/models/utils/parametersystem.hpp>
#include <opm/models/utils/start.hh>
#include <opm/simulators/flow/FlowGenericVanguard.hpp>
#include <opm/simulators/flow/BlackoilModelParameters.hpp>
#include <opm/simulators/timestepping/EclTimeSteppingParams.hpp>
#if HAVE_DUNE_FEM
#include <dune/fem/misc/mpimanager.hh>
#else
#include <dune/common/parallel/mpihelper.hh>
#endif
#include <boost/test/unit_test.hpp>
#include <boost/version.hpp>
#if BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 < 71
#include <boost/test/floating_point_comparison.hpp>
#else
#include <boost/test/tools/floating_point_comparison.hpp>
#endif
#include <memory>
#include <string>
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<Parallel::Communication>());
}
};
/*!
* \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 <class TypeTag>
std::unique_ptr<GetPropType<TypeTag, Properties::Simulator>>
initSimulator(const char* filename,
const char* test_name = "test_simulator",
int threads_per_process = 1)
{
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
std::string filename_arg = "--ecl-deck-file-name=";
filename_arg += filename;
const char* argv[] = {
test_name,
filename_arg.c_str()
};
Parameters::reset();
registerAllParameters_<TypeTag>(false);
registerEclTimeSteppingParameters<double>();
BlackoilModelParameters<double>::registerParameters();
Parameters::Register<Parameters::EnableTerminalOutput>("Do *NOT* use!");
Parameters::SetDefault<Parameters::ThreadsPerProcess>(threads_per_process);
Parameters::endRegistration();
setupParameters_<TypeTag>(/*argc=*/sizeof(argv) / sizeof(argv[0]),
argv,
/*registerParams=*/false,
/*allowUnused=*/false,
/*handleHelp=*/true,
/*myRank=*/0);
FlowGenericVanguard::readDeck(filename);
return std::make_unique<Simulator>();
}
} // namespace Opm
#endif // OPM_SIMULATOR_FIXTURE_HPP
tests/test_glift1.cpp
+8 -82
View File
@@ -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 <opm/models/utils/propertysystem.hh>
#include <opm/models/utils/parametersystem.hpp>
#include <opm/models/utils/start.hh>
#include "SimulatorFixture.hpp"
#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/Schedule/Well/Well.hpp>
#include <opm/simulators/utils/DeferredLogger.hpp>
#include <opm/simulators/flow/BlackoilModel.hpp>
#include <opm/simulators/flow/FlowProblemBlackoil.hpp>
#include <opm/simulators/flow/equil/EquilibrationHelpers.hpp>
#include <opm/simulators/wells/BlackoilWellModel.hpp>
@@ -45,87 +40,18 @@
#include <opm/simulators/wells/BlackoilWellModelGasLift.hpp>
#include <opm/simulators/wells/WellState.hpp>
#if HAVE_DUNE_FEM
#include <dune/fem/misc/mpimanager.hh>
#else
#include <dune/common/parallel/mpihelper.hh>
#endif
#include <memory>
#include <optional>
#include <string>
#include <vector>
#include <boost/test/unit_test.hpp>
#include <boost/version.hpp>
#if BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 < 71
#include <boost/test/floating_point_comparison.hpp>
#else
#include <boost/test/tools/floating_point_comparison.hpp>
#endif
namespace Opm::Properties {
namespace TTag {
struct TestGliftTypeTag {
using InheritsFrom = std::tuple<TestTypeTag>;
};
}
}
template <class TypeTag>
std::unique_ptr<Opm::GetPropType<TypeTag, Opm::Properties::Simulator>>
initSimulator(const char *filename)
{
using namespace Opm;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
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<TestTypeTag>;
};
Parameters::reset();
registerAllParameters_<TypeTag>(false);
registerEclTimeSteppingParameters<double>();
BlackoilModelParameters<double>::registerParameters();
Parameters::Register<Parameters::EnableTerminalOutput>("Do *NOT* use!");
Opm::Parameters::SetDefault<Opm::Parameters::ThreadsPerProcess>(2);
Parameters::endRegistration();
setupParameters_<TypeTag>(/*argc=*/sizeof(argv) / sizeof(argv[0]),
argv,
/*registerParams=*/false,
/*allowUnused*/false,
/*handleHelp*/true,
/*myRank*/0);
FlowGenericVanguard::readDeck(filename);
return std::make_unique<Simulator>();
}
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<Opm::Parallel::Communication>());
}
};
}
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<TypeTag>(filename.data());
auto simulator = Opm::initSimulator<TypeTag>(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<TypeTag>(filename.data());
auto simulator = Opm::initSimulator<TypeTag>(filename.data(), "test_glift1", /*threads=*/2);
simulator->model().applyInitialSolution();
simulator->setEpisodeIndex(-1);
tests/test_group_higher_constraints.cpp
+411
View File
@@ -0,0 +1,411 @@
// -*- 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 <http://www.gnu.org/licenses/>.
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 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
* - 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 <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/Schedule/Group/Group.hpp>
#include <opm/input/eclipse/Schedule/Group/GuideRate.hpp>
#include <opm/input/eclipse/Schedule/Group/GuideRateModel.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/simulators/flow/FlowProblemBlackoil.hpp>
#include <opm/simulators/wells/BlackoilWellModel.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <opm/simulators/wells/GroupState.hpp>
#include <boost/test/unit_test.hpp>
namespace Opm::Properties::TTag {
struct TestEfficiencyFactorTypeTag {
using InheritsFrom = std::tuple<TestTypeTag>;
};
}
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<TypeTag>;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
using IndexTraits = typename FluidSystem::IndexTraitsType;
const std::string filename = "GROUP_HIGHER_CONSTRAINTS.DATA";
auto simulator = Opm::initSimulator<TypeTag>(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<Opm::Well::Status> 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<double> mani_prod_rates(num_phases, 0.0);
mani_prod_rates[oil_phase_pos] = -well_a_oil_rate; // 10000 SM3/day (positive)
std::vector<double> mani2_prod_rates(num_phases, 0.0);
mani2_prod_rates[oil_phase_pos] = -well_b_oil_rate; // 4000 SM3/day (positive)
std::vector<double> 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<double> 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<double> 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<double> 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 fixed 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
//
// Addback uses partial efficiency E_MANI (not accumulated 0.72):
// addback = current_rate * E_MANI = 10000 * 0.8 = 8000
// target = 2000 + 8000 = 10000
//
// local_fraction("MANI") = 8000/10800 = 20/27
// target = 10000 * 20/27 = 200000/27
//
// Final calculation:
// target_rate_available = target / efficiency_factor
// = (200000/27) / 0.72 = 200000/27 * 25/18 = 5000000/486
// = 10288.0658... SM3/day
//
// Scale factor:
// scale = target_rate_available / current_rate_available
// = (5000000/486) / 10000 = 500/486 = 250/243 = 1.02880658...
//
// Because scale > 1, no constraint violation - MANI can produce its full rate.
//
// 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_fix = scale * current_rate;
checkAlgo(target_rate_available_with_fix, 5000000.0 / 486.0); // 10288.0658...
}
BOOST_AUTO_TEST_SUITE_END()