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()