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

Fix efficiency factor addback in group constraint checking

Browse Source 
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.
This commit is contained in:
Håkon Hægland committed 2025-12-15 08:28:39 +01:00
1 parent 172d62e93a
commit da720e1900
3 files changed
+74 -28
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

opm/simulators/wells/GroupStateHelper.cpp
+40 -4
View File
@@ -232,7 +232,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]);
@@ -404,7 +406,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]);
@@ -603,7 +607,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)) {
@@ -731,7 +737,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)) {
@@ -1381,6 +1389,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
@@ -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<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/test_group_higher_constraints.cpp
+21 -24
View File
@@ -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()