From 4f21916b004513ca98b24dc397654d4a7c07e6ac Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 28 Feb 2026 21:44:34 +0000 Subject: [PATCH 1/8] Initial plan From 581745485930bcea001d7741b13a8f5d4536e7c7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 28 Feb 2026 21:54:35 +0000 Subject: [PATCH 2/8] Add standalone GroupTreeRates module for group-tree rate distribution Implements the group-tree rate distribution algorithm as standalone functions: - distributeGroupTreeRates: Main iterative rate balancing - setSubRates: Top-down guide-rate based distribution - findWorstOffendingChild: Find worst limit violation in subtree - updateParentStatus: Bottom-up status propagation Includes comprehensive unit tests covering unconstrained trees, single-well limits, group limit violations, cascading violations, and deep tree hierarchies. Co-authored-by: steink <6746854+steink@users.noreply.github.com> --- CMakeLists_files.cmake | 3 + opm/simulators/wells/GroupTreeRates.cpp | 196 +++++++++++ opm/simulators/wells/GroupTreeRates.hpp | 107 ++++++ tests/test_GroupTreeRates.cpp | 427 ++++++++++++++++++++++++ 4 files changed, 733 insertions(+) create mode 100644 opm/simulators/wells/GroupTreeRates.cpp create mode 100644 opm/simulators/wells/GroupTreeRates.hpp create mode 100644 tests/test_GroupTreeRates.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 32dfb6ab408..ec54e0342d2 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -230,6 +230,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/wells/GroupEconomicLimitsChecker.cpp opm/simulators/wells/GroupState.cpp opm/simulators/wells/GroupStateHelper.cpp + opm/simulators/wells/GroupTreeRates.cpp opm/simulators/wells/MSWellHelpers.cpp opm/simulators/wells/MultisegmentWellAssemble.cpp opm/simulators/wells/MultisegmentWellEquations.cpp @@ -489,6 +490,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_glift1.cpp tests/test_graphcoloring.cpp tests/test_GroupState.cpp + tests/test_GroupTreeRates.cpp tests/test_interregflows.cpp tests/test_invert.cpp tests/test_keyword_validator.cpp @@ -1179,6 +1181,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/wells/GroupEconomicLimitsChecker.hpp opm/simulators/wells/GroupState.hpp opm/simulators/wells/GroupStateHelper.hpp + opm/simulators/wells/GroupTreeRates.hpp opm/simulators/wells/GuideRateHandler.hpp opm/simulators/wells/MSWellHelpers.hpp opm/simulators/wells/MultisegmentWell.hpp diff --git a/opm/simulators/wells/GroupTreeRates.cpp b/opm/simulators/wells/GroupTreeRates.cpp new file mode 100644 index 00000000000..dc704bbf918 --- /dev/null +++ b/opm/simulators/wells/GroupTreeRates.cpp @@ -0,0 +1,196 @@ +/* + Copyright 2025 Equinor ASA + + 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 . +*/ + +#include +#include + +#include +#include + +namespace Opm { + +template +int distributeGroupTreeRates(std::vector>& tree, + int rootIndex, + int maxIter) +{ + if (rootIndex == 0) { + // Initialise: all nodes group-controlled except root + for (auto& node : tree) { + node.status = -1; + } + tree[rootIndex].status = 1; + tree[rootIndex].rate = tree[rootIndex].limit; + } + + int iter = 0; + for (; iter < maxIter; ++iter) { + setSubRates(tree, rootIndex); + + Scalar worstExcess{0}; + const int ix = findWorstOffendingChild(tree, rootIndex, worstExcess); + if (ix < 0) { + break; // converged + } + + // Fix the worst-offending node at its limit + tree[ix].rate = tree[ix].limit; + tree[ix].status = 1; + + // If it has children, recursively solve its subtree + if (!tree[ix].children.empty()) { + distributeGroupTreeRates(tree, ix, maxIter); + } + + updateParentStatus(tree, ix); + } + return iter; +} + +template +void setSubRates(std::vector>& tree, + int nodeIndex) +{ + const auto& children = tree[nodeIndex].children; + if (children.empty()) { + return; // leaf (well) node + } + + // Separate children into fixed (status != -1) and group-controlled + Scalar fixedRate{0}; + Scalar guideSum{0}; + bool allFixed = true; + + for (const int ci : children) { + if (tree[ci].status != -1) { + fixedRate += tree[ci].rate; + } else { + guideSum += tree[ci].guideRate; + allFixed = false; + } + } + + if (allFixed) { + tree[nodeIndex].status = 0; + tree[nodeIndex].rate = fixedRate; + return; + } + + const Scalar availRate = tree[nodeIndex].rate - fixedRate; + assert(availRate > Scalar{0} || guideSum == Scalar{0}); + + for (const int ci : children) { + if (tree[ci].status == -1) { + tree[ci].rate = (guideSum > Scalar{0}) + ? availRate * tree[ci].guideRate / guideSum + : Scalar{0}; + setSubRates(tree, ci); + } + } +} + +template +int findWorstOffendingChild(const std::vector>& tree, + int nodeIndex, + Scalar& worstExcess) +{ + const auto& children = tree[nodeIndex].children; + int worstIndex = -1; + + for (const int ci : children) { + const Scalar excess = tree[ci].rate - tree[ci].limit; + if (excess > worstExcess) { + worstExcess = excess; + worstIndex = ci; + } + } + + // Recurse into children that have subtrees + for (const int ci : children) { + if (!tree[ci].children.empty()) { + Scalar childExcess = worstExcess; + const int childWorst = findWorstOffendingChild(tree, ci, childExcess); + if (childExcess > worstExcess) { + worstExcess = childExcess; + worstIndex = childWorst; + } + } + } + + if (worstExcess > Scalar{0} && worstIndex >= 0) { + assert(tree[worstIndex].status == -1); + } else { + worstIndex = -1; + } + return worstIndex; +} + +template +void updateParentStatus(std::vector>& tree, + int nodeIndex) +{ + const int pix = tree[nodeIndex].parent; + if (pix < 0) { + return; + } + + const auto& siblings = tree[pix].children; + bool allFixed = true; + Scalar totalRate{0}; + + for (const int ci : siblings) { + if (tree[ci].status == -1) { + allFixed = false; + break; + } + totalRate += tree[ci].rate; + } + + if (allFixed) { + tree[pix].status = 0; + tree[pix].rate = totalRate; + updateParentStatus(tree, pix); + } +} + +// Explicit template instantiations +template struct GroupTreeNode; +template struct GroupTreeNode; + +template int distributeGroupTreeRates( + std::vector>&, int, int); +template int distributeGroupTreeRates( + std::vector>&, int, int); + +template void setSubRates( + std::vector>&, int); +template void setSubRates( + std::vector>&, int); + +template int findWorstOffendingChild( + const std::vector>&, int, double&); +template int findWorstOffendingChild( + const std::vector>&, int, float&); + +template void updateParentStatus( + std::vector>&, int); +template void updateParentStatus( + std::vector>&, int); + +} // namespace Opm diff --git a/opm/simulators/wells/GroupTreeRates.hpp b/opm/simulators/wells/GroupTreeRates.hpp new file mode 100644 index 00000000000..780a5f8adf5 --- /dev/null +++ b/opm/simulators/wells/GroupTreeRates.hpp @@ -0,0 +1,107 @@ +/* + Copyright 2025 Equinor ASA + + 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 . +*/ + +#ifndef OPM_GROUP_TREE_RATES_HPP +#define OPM_GROUP_TREE_RATES_HPP + +#include +#include +#include + +namespace Opm { + +/// Simplified node in the group/well hierarchy used for rate distribution. +/// +/// Each node represents either a group or a well. Wells are leaf nodes +/// (children is empty). The root node typically represents the FIELD group. +/// +/// The status field encodes the control state: +/// -1 : group-controlled (rate determined by parent allocation) +/// 0 : fully determined by children (all children are fixed) +/// 1 : individually rate-limited +template +struct GroupTreeNode { + int index{-1}; ///< Unique node index in the flat vector + int parent{-1}; ///< Parent node index (-1 for root) + Scalar rate{0}; ///< Current allocated rate + Scalar limit{0}; ///< Rate limit for this node + int status{-1}; ///< Control status (-1, 0, or 1) + Scalar guideRate{0}; ///< Guide rate used for allocation + std::string name; ///< Group or well name + std::vector children; ///< Indices of child nodes +}; + +/// Distribute rates through a group tree respecting individual node limits. +/// +/// Starting from the root node (at @p rootIndex), rates are distributed +/// top-down to children according to their guide-rate fractions. Whenever +/// a child's allocated rate exceeds its limit, that child is fixed at its +/// limit and the distribution is recomputed for the remaining +/// group-controlled children. The process repeats until no limit +/// violations remain. +/// +/// @param tree Flat vector of tree nodes (modified in-place). +/// @param rootIndex Index of the subtree root in @p tree. +/// @param maxIter Maximum number of outer iterations (safety limit). +/// @return Number of outer iterations used. +template +int distributeGroupTreeRates(std::vector>& tree, + int rootIndex = 0, + int maxIter = 1000); + +/// Distribute the rate of node @p nodeIndex to its children. +/// +/// Children that are already individually limited (status != -1) keep +/// their current rate. The remaining (available) rate is split among +/// group-controlled children proportionally to their guide rates. +/// If all children are fixed, the parent's status becomes 0 and its +/// rate is set to the sum of the children's rates. +/// +/// @param tree Flat vector of tree nodes (modified in-place). +/// @param nodeIndex Index of the node whose children are updated. +template +void setSubRates(std::vector>& tree, + int nodeIndex); + +/// Find the child (in the subtree rooted at @p nodeIndex) whose rate +/// most exceeds its limit. +/// +/// @param tree Flat vector of tree nodes (read-only access). +/// @param nodeIndex Root of the subtree to search. +/// @param[out] worstExcess The largest (rate - limit) value found. +/// @return Index of the worst offending node, or -1 if no +/// violation exists. +template +int findWorstOffendingChild(const std::vector>& tree, + int nodeIndex, + Scalar& worstExcess); + +/// After fixing a child node, propagate upward: if all siblings are +/// also fixed, the parent becomes fully determined (status 0) with +/// rate equal to the sum of its children. +/// +/// @param tree Flat vector of tree nodes (modified in-place). +/// @param nodeIndex Index of the node whose parent is updated. +template +void updateParentStatus(std::vector>& tree, + int nodeIndex); + +} // namespace Opm + +#endif // OPM_GROUP_TREE_RATES_HPP diff --git a/tests/test_GroupTreeRates.cpp b/tests/test_GroupTreeRates.cpp new file mode 100644 index 00000000000..3a76b2fe9ee --- /dev/null +++ b/tests/test_GroupTreeRates.cpp @@ -0,0 +1,427 @@ +/* + Copyright 2025 Equinor ASA + + 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 . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include + +#define BOOST_TEST_MODULE GroupTreeRatesTest +#include + +#include +#include + +using namespace Opm; + +namespace { + +// Helper: build a simple tree +// +// [0] FIELD (limit=10000) +// / \ +// [1] G1 [2] G2 +// (lim=6000) (lim=7000) +// / \ / \ +// [3]W1 [4]W2 [5]W3 [6]W4 +// (4000) (3000) (4500) (3500) +// +// Guide rates: G1=0.5, G2=0.5, W1=0.6, W2=0.4, W3=0.55, W4=0.45 +std::vector> makeSimpleTree() +{ + std::vector> tree(7); + + // FIELD (root) + tree[0] = {0, -1, 0.0, 10000.0, -1, 1.0, "FIELD", {1, 2}}; + // G1 + tree[1] = {1, 0, 0.0, 6000.0, -1, 0.5, "G1", {3, 4}}; + // G2 + tree[2] = {2, 0, 0.0, 7000.0, -1, 0.5, "G2", {5, 6}}; + // W1 + tree[3] = {3, 1, 0.0, 4000.0, -1, 0.6, "W1", {}}; + // W2 + tree[4] = {4, 1, 0.0, 3000.0, -1, 0.4, "W2", {}}; + // W3 + tree[5] = {5, 2, 0.0, 4500.0, -1, 0.55, "W3", {}}; + // W4 + tree[6] = {6, 2, 0.0, 3500.0, -1, 0.45, "W4", {}}; + + return tree; +} + +// Helper: build a tree where no limits are binding (all limits > field rate) +std::vector> makeUnconstrainedTree() +{ + std::vector> tree(5); + + // FIELD + tree[0] = {0, -1, 0.0, 1000.0, -1, 1.0, "FIELD", {1, 2}}; + // G1 + tree[1] = {1, 0, 0.0, 9000.0, -1, 0.6, "G1", {3, 4}}; + // G2 (leaf well acting as a group-less well) + tree[2] = {2, 0, 0.0, 9000.0, -1, 0.4, "W_SOLO", {}}; + // W1 + tree[3] = {3, 1, 0.0, 9000.0, -1, 0.5, "W1", {}}; + // W2 + tree[4] = {4, 1, 0.0, 9000.0, -1, 0.5, "W2", {}}; + + return tree; +} + +// Helper: build a single-well tree (FIELD -> W1) +std::vector> makeSingleWellTree() +{ + std::vector> tree(2); + tree[0] = {0, -1, 0.0, 500.0, -1, 1.0, "FIELD", {1}}; + tree[1] = {1, 0, 0.0, 300.0, -1, 1.0, "W1", {}}; + return tree; +} + +double totalLeafRate(const std::vector>& tree) +{ + double total = 0.0; + for (const auto& node : tree) { + if (node.children.empty()) { + total += node.rate; + } + } + return total; +} + +} // anonymous namespace + + +BOOST_AUTO_TEST_CASE(SingleWell_WellLimitBinding) +{ + auto tree = makeSingleWellTree(); + // FIELD limit = 500, W1 limit = 300 => W1 should be limited to 300 + const int iter = distributeGroupTreeRates(tree); + + BOOST_CHECK(iter > 0); + // W1 should be at its limit + BOOST_CHECK_CLOSE(tree[1].rate, 300.0, 1e-10); + // FIELD should reflect the well rate + BOOST_CHECK_CLOSE(tree[0].rate, 300.0, 1e-10); + // W1 individually limited + BOOST_CHECK_EQUAL(tree[1].status, 1); + // FIELD determined by children + BOOST_CHECK_EQUAL(tree[0].status, 0); +} + +BOOST_AUTO_TEST_CASE(SingleWell_FieldLimitBinding) +{ + auto tree = makeSingleWellTree(); + tree[1].limit = 800.0; // W1 limit > FIELD limit + const int iter = distributeGroupTreeRates(tree); + + BOOST_CHECK(iter >= 0); + // W1 rate should equal FIELD limit since it is the only well + BOOST_CHECK_CLOSE(tree[1].rate, 500.0, 1e-10); +} + +BOOST_AUTO_TEST_CASE(Unconstrained_RatesFollowGuideRates) +{ + auto tree = makeUnconstrainedTree(); + distributeGroupTreeRates(tree); + + // FIELD rate = 1000 (its limit) + const double fieldRate = 1000.0; + // G1 gets 60%, W_SOLO gets 40% + BOOST_CHECK_CLOSE(tree[1].rate, fieldRate * 0.6, 1e-10); + BOOST_CHECK_CLOSE(tree[2].rate, fieldRate * 0.4, 1e-10); + // W1, W2 each get 50% of G1 + BOOST_CHECK_CLOSE(tree[3].rate, fieldRate * 0.6 * 0.5, 1e-10); + BOOST_CHECK_CLOSE(tree[4].rate, fieldRate * 0.6 * 0.5, 1e-10); + + // Total leaf rate should equal field rate + BOOST_CHECK_CLOSE(totalLeafRate(tree), fieldRate, 1e-10); +} + +BOOST_AUTO_TEST_CASE(SimpleTree_GroupAndWellLimits) +{ + auto tree = makeSimpleTree(); + distributeGroupTreeRates(tree); + + // All well rates must not exceed their limits + for (const auto& node : tree) { + BOOST_CHECK_LE(node.rate, node.limit + 1e-10); + } + + // Total leaf rate must not exceed FIELD limit + BOOST_CHECK_LE(totalLeafRate(tree), tree[0].limit + 1e-10); + + // Group rates must not exceed their limits + BOOST_CHECK_LE(tree[1].rate, tree[1].limit + 1e-10); // G1 + BOOST_CHECK_LE(tree[2].rate, tree[2].limit + 1e-10); // G2 +} + +BOOST_AUTO_TEST_CASE(SimpleTree_RateConsistency) +{ + auto tree = makeSimpleTree(); + distributeGroupTreeRates(tree); + + // G1 rate == W1 rate + W2 rate + BOOST_CHECK_CLOSE(tree[1].rate, tree[3].rate + tree[4].rate, 1e-10); + // G2 rate == W3 rate + W4 rate + BOOST_CHECK_CLOSE(tree[2].rate, tree[5].rate + tree[6].rate, 1e-10); + // FIELD rate == G1 rate + G2 rate + BOOST_CHECK_CLOSE(tree[0].rate, tree[1].rate + tree[2].rate, 1e-10); +} + +BOOST_AUTO_TEST_CASE(SimpleTree_StatusConsistency) +{ + auto tree = makeSimpleTree(); + distributeGroupTreeRates(tree); + + // Root should be individually limited (status 1) + BOOST_CHECK_EQUAL(tree[0].status, 1); + + // If a node has status 0, its rate must equal the sum of children's rates + // If a node has status 1, its rate must equal its limit + for (const auto& node : tree) { + if (node.status == 1) { + BOOST_CHECK_CLOSE(node.rate, node.limit, 1e-10); + } + if (node.status == 0) { + double childSum = 0.0; + for (int ci : node.children) { + childSum += tree[ci].rate; + } + BOOST_CHECK_CLOSE(node.rate, childSum, 1e-10); + } + } +} + +BOOST_AUTO_TEST_CASE(TightWellLimit) +{ + // One well has a very tight limit + auto tree = makeUnconstrainedTree(); + tree[3].limit = 10.0; // W1 has tight limit + + distributeGroupTreeRates(tree); + + // W1 must be at its limit + BOOST_CHECK_CLOSE(tree[3].rate, 10.0, 1e-10); + BOOST_CHECK_EQUAL(tree[3].status, 1); + + // Total still respects FIELD + BOOST_CHECK_LE(totalLeafRate(tree), tree[0].limit + 1e-10); + + // Rate consistency + BOOST_CHECK_CLOSE(tree[1].rate, tree[3].rate + tree[4].rate, 1e-10); +} + +BOOST_AUTO_TEST_CASE(AllWellsTight) +{ + // All wells have very tight limits (sum < FIELD limit) + auto tree = makeSimpleTree(); + tree[3].limit = 100.0; // W1 + tree[4].limit = 100.0; // W2 + tree[5].limit = 100.0; // W3 + tree[6].limit = 100.0; // W4 + + distributeGroupTreeRates(tree); + + // Each well should be at its limit + BOOST_CHECK_CLOSE(tree[3].rate, 100.0, 1e-10); + BOOST_CHECK_CLOSE(tree[4].rate, 100.0, 1e-10); + BOOST_CHECK_CLOSE(tree[5].rate, 100.0, 1e-10); + BOOST_CHECK_CLOSE(tree[6].rate, 100.0, 1e-10); + + // Total = 400 < FIELD limit of 10000 + BOOST_CHECK_CLOSE(totalLeafRate(tree), 400.0, 1e-10); +} + +BOOST_AUTO_TEST_CASE(GroupLimitViolation) +{ + // Tree where initial guide-rate distribution violates a group limit: + // FIELD (limit=10000), G1 (lim=3000), G2 (lim=8000) + // Guide rates 0.5/0.5 => initial: G1=5000 > 3000 (violation!) + // After fixing G1=3000, G2 gets 7000 + std::vector> tree(5); + tree[0] = {0, -1, 0.0, 10000.0, -1, 1.0, "FIELD", {1, 2}}; + tree[1] = {1, 0, 0.0, 3000.0, -1, 0.5, "G1", {3}}; + tree[2] = {2, 0, 0.0, 8000.0, -1, 0.5, "G2", {4}}; + tree[3] = {3, 1, 0.0, 5000.0, -1, 1.0, "W1", {}}; + tree[4] = {4, 2, 0.0, 9000.0, -1, 1.0, "W2", {}}; + + const int iter = distributeGroupTreeRates(tree); + + // Should take at least 1 iteration to fix the G1 violation + BOOST_CHECK(iter >= 1); + // G1 should be at its limit + BOOST_CHECK_CLOSE(tree[1].rate, 3000.0, 1e-10); + BOOST_CHECK_EQUAL(tree[1].status, 1); + // W1 = G1's rate = 3000 (within its 5000 limit) + BOOST_CHECK_CLOSE(tree[3].rate, 3000.0, 1e-10); + // G2 gets the remaining 7000 + BOOST_CHECK_CLOSE(tree[2].rate, 7000.0, 1e-10); + // W2 = G2's rate = 7000 (within its 9000 limit) + BOOST_CHECK_CLOSE(tree[4].rate, 7000.0, 1e-10); + // Total + BOOST_CHECK_CLOSE(totalLeafRate(tree), 10000.0, 1e-10); +} + +BOOST_AUTO_TEST_CASE(CascadingViolations) +{ + // Multiple cascading violations: + // FIELD (10000) -> G1 (lim=4000, guide=0.5), G2 (lim=9000, guide=0.5) + // G1 -> W1 (lim=1000, guide=0.6), W2 (lim=5000, guide=0.4) + // G2 -> W3 (lim=3000, guide=0.5), W4 (lim=9000, guide=0.5) + // + // Initial: G1=5000 (>4000 viol), G2=5000 + // G1 children: W1=3000 (>1000 viol), W2=2000 + // G2 children: W3=2500, W4=2500 + // After fixing: multiple iterations needed + std::vector> tree(7); + tree[0] = {0, -1, 0.0, 10000.0, -1, 1.0, "FIELD", {1, 2}}; + tree[1] = {1, 0, 0.0, 4000.0, -1, 0.5, "G1", {3, 4}}; + tree[2] = {2, 0, 0.0, 9000.0, -1, 0.5, "G2", {5, 6}}; + tree[3] = {3, 1, 0.0, 1000.0, -1, 0.6, "W1", {}}; + tree[4] = {4, 1, 0.0, 5000.0, -1, 0.4, "W2", {}}; + tree[5] = {5, 2, 0.0, 3000.0, -1, 0.5, "W3", {}}; + tree[6] = {6, 2, 0.0, 9000.0, -1, 0.5, "W4", {}}; + + const int iter = distributeGroupTreeRates(tree); + BOOST_CHECK(iter >= 1); + + // All limits respected + for (const auto& node : tree) { + BOOST_CHECK_LE(node.rate, node.limit + 1e-10); + } + + // W1 must be at its limit (1000) + BOOST_CHECK_CLOSE(tree[3].rate, 1000.0, 1e-10); + BOOST_CHECK_EQUAL(tree[3].status, 1); + + // Rate consistency + BOOST_CHECK_CLOSE(tree[1].rate, tree[3].rate + tree[4].rate, 1e-10); + BOOST_CHECK_CLOSE(tree[2].rate, tree[5].rate + tree[6].rate, 1e-10); + BOOST_CHECK_CLOSE(tree[0].rate, tree[1].rate + tree[2].rate, 1e-10); +} + +BOOST_AUTO_TEST_CASE(DeepTree) +{ + // FIELD -> G1 -> G2 -> W1 + std::vector> tree(4); + tree[0] = {0, -1, 0.0, 5000.0, -1, 1.0, "FIELD", {1}}; + tree[1] = {1, 0, 0.0, 3000.0, -1, 1.0, "G1", {2}}; + tree[2] = {2, 1, 0.0, 2000.0, -1, 1.0, "G2", {3}}; + tree[3] = {3, 2, 0.0, 1500.0, -1, 1.0, "W1", {}}; + + distributeGroupTreeRates(tree); + + // W1 is the tightest at 1500 + BOOST_CHECK_CLOSE(tree[3].rate, 1500.0, 1e-10); + BOOST_CHECK_CLOSE(tree[2].rate, 1500.0, 1e-10); + BOOST_CHECK_CLOSE(tree[1].rate, 1500.0, 1e-10); + BOOST_CHECK_CLOSE(tree[0].rate, 1500.0, 1e-10); +} + +BOOST_AUTO_TEST_CASE(ZeroIterations_NoViolation) +{ + // If root limit is 0, no distribution needed + std::vector> tree(2); + tree[0] = {0, -1, 0.0, 0.0, -1, 1.0, "FIELD", {1}}; + tree[1] = {1, 0, 0.0, 100.0, -1, 1.0, "W1", {}}; + + const int iter = distributeGroupTreeRates(tree); + + // Should converge immediately (no violations since rates are 0) + BOOST_CHECK_EQUAL(iter, 0); + BOOST_CHECK_CLOSE(tree[1].rate, 0.0, 1e-10); +} + +BOOST_AUTO_TEST_CASE(EqualGuideRates) +{ + // Three wells with equal guide rates + std::vector> tree(4); + tree[0] = {0, -1, 0.0, 900.0, -1, 1.0, "FIELD", {1, 2, 3}}; + tree[1] = {1, 0, 0.0, 500.0, -1, 1.0, "W1", {}}; + tree[2] = {2, 0, 0.0, 500.0, -1, 1.0, "W2", {}}; + tree[3] = {3, 0, 0.0, 500.0, -1, 1.0, "W3", {}}; + + distributeGroupTreeRates(tree); + + // Each well gets 300 + BOOST_CHECK_CLOSE(tree[1].rate, 300.0, 1e-10); + BOOST_CHECK_CLOSE(tree[2].rate, 300.0, 1e-10); + BOOST_CHECK_CLOSE(tree[3].rate, 300.0, 1e-10); +} + +BOOST_AUTO_TEST_CASE(SetSubRates_Basic) +{ + auto tree = makeUnconstrainedTree(); + // Manually set root as limited with known rate + tree[0].status = 1; + tree[0].rate = 1000.0; + + setSubRates(tree, 0); + + BOOST_CHECK_CLOSE(tree[1].rate, 600.0, 1e-10); + BOOST_CHECK_CLOSE(tree[2].rate, 400.0, 1e-10); + BOOST_CHECK_CLOSE(tree[3].rate, 300.0, 1e-10); + BOOST_CHECK_CLOSE(tree[4].rate, 300.0, 1e-10); +} + +BOOST_AUTO_TEST_CASE(FindWorstOffending_Basic) +{ + auto tree = makeSimpleTree(); + // Set FIELD rate and distribute + tree[0].status = 1; + tree[0].rate = 10000.0; + setSubRates(tree, 0); + + double worstExcess = 0.0; + const int ix = findWorstOffendingChild(tree, 0, worstExcess); + + // At least one child should be offending since rates > limits for some + BOOST_CHECK(ix >= 0 || worstExcess <= 0.0); +} + +BOOST_AUTO_TEST_CASE(UpdateParentStatus_AllFixed) +{ + auto tree = makeUnconstrainedTree(); + // Fix all children of FIELD + tree[1].status = 1; + tree[1].rate = 500.0; + tree[2].status = 1; + tree[2].rate = 300.0; + tree[0].status = -1; + + updateParentStatus(tree, 1); + + BOOST_CHECK_EQUAL(tree[0].status, 0); + BOOST_CHECK_CLOSE(tree[0].rate, 800.0, 1e-10); +} + +BOOST_AUTO_TEST_CASE(UpdateParentStatus_NotAllFixed) +{ + auto tree = makeUnconstrainedTree(); + tree[1].status = 1; + tree[1].rate = 500.0; + tree[2].status = -1; // still group-controlled + tree[0].status = -1; + + updateParentStatus(tree, 1); + + // Parent should NOT be updated since not all children are fixed + BOOST_CHECK_EQUAL(tree[0].status, -1); +} From e625ead47cd1753b9e06b340556143fdda9fb0fd Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 28 Feb 2026 21:55:12 +0000 Subject: [PATCH 3/8] Address code review: improve assertion clarity and rename misleading variable Co-authored-by: steink <6746854+steink@users.noreply.github.com> --- opm/simulators/wells/GroupTreeRates.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/opm/simulators/wells/GroupTreeRates.cpp b/opm/simulators/wells/GroupTreeRates.cpp index dc704bbf918..43c856d513e 100644 --- a/opm/simulators/wells/GroupTreeRates.cpp +++ b/opm/simulators/wells/GroupTreeRates.cpp @@ -93,6 +93,7 @@ void setSubRates(std::vector>& tree, } const Scalar availRate = tree[nodeIndex].rate - fixedRate; + assert(availRate >= Scalar{0}); assert(availRate > Scalar{0} || guideSum == Scalar{0}); for (const int ci : children) { @@ -124,10 +125,10 @@ int findWorstOffendingChild(const std::vector>& tree, // Recurse into children that have subtrees for (const int ci : children) { if (!tree[ci].children.empty()) { - Scalar childExcess = worstExcess; - const int childWorst = findWorstOffendingChild(tree, ci, childExcess); - if (childExcess > worstExcess) { - worstExcess = childExcess; + Scalar subtreeExcess = worstExcess; + const int childWorst = findWorstOffendingChild(tree, ci, subtreeExcess); + if (subtreeExcess > worstExcess) { + worstExcess = subtreeExcess; worstIndex = childWorst; } } From c5363dc40072a8e309bd207b64e5fa2d783e1af3 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 28 Feb 2026 21:55:41 +0000 Subject: [PATCH 4/8] Fix BOOST_CHECK_CLOSE with zero value to use BOOST_CHECK_SMALL Co-authored-by: steink <6746854+steink@users.noreply.github.com> --- tests/test_GroupTreeRates.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_GroupTreeRates.cpp b/tests/test_GroupTreeRates.cpp index 3a76b2fe9ee..424c948bde8 100644 --- a/tests/test_GroupTreeRates.cpp +++ b/tests/test_GroupTreeRates.cpp @@ -346,7 +346,7 @@ BOOST_AUTO_TEST_CASE(ZeroIterations_NoViolation) // Should converge immediately (no violations since rates are 0) BOOST_CHECK_EQUAL(iter, 0); - BOOST_CHECK_CLOSE(tree[1].rate, 0.0, 1e-10); + BOOST_CHECK_SMALL(tree[1].rate, 1e-10); } BOOST_AUTO_TEST_CASE(EqualGuideRates) From 255e0232c5fc4268466f5993c4a889535aa53703 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 2 Mar 2026 13:11:09 +0000 Subject: [PATCH 5/8] Refactor GroupTreeRates to OPM class-based style and integrate with updateWellControls - Wrap free functions into GroupTreeRates class with static methods - Rename to snake_case: distributeGroupTreeRates -> distribute, setSubRates -> set_sub_rates, findWorstOffendingChild -> find_worst_offending_child, updateParentStatus -> update_parent_status - Rename struct member: guideRate -> guide_rate - Use OPM-style include guard and Doxygen documentation - Add include and integration comment in BlackoilWellModel_impl.hpp showing how distribute() fits into updateWellControls - Update all tests to use new class-based API Co-authored-by: steink <6746854+steink@users.noreply.github.com> --- .../wells/BlackoilWellModel_impl.hpp | 11 ++ opm/simulators/wells/GroupTreeRates.cpp | 142 ++++++++---------- opm/simulators/wells/GroupTreeRates.hpp | 122 ++++++++------- tests/test_GroupTreeRates.cpp | 40 ++--- 4 files changed, 162 insertions(+), 153 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 4cf608a9563..860387c5008 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -42,6 +42,7 @@ #include #include +#include #include #include #include @@ -1723,6 +1724,16 @@ namespace Opm { // update wsolvent fraction for REIN wells this->updateWsolvent(fieldGroup, episodeIdx, this->nupcolWellState()); + // TODO: Consider using GroupTreeRates::distribute() here to solve + // the group-tree rates to a more converged state under the simplifying + // assumption of constant well-rate fractions. This would involve: + // 1. Building a flat GroupTreeNode vector from the current group hierarchy + // (using fieldGroup, schedule, guide rates, and current well/group rates) + // 2. Calling GroupTreeRates::distribute(tree) + // 3. Applying the resulting rates back to the well/group state + // This approach could improve convergence of group-tree rate balancing + // without requiring network pressure recalculation. + return changed_well_group; } diff --git a/opm/simulators/wells/GroupTreeRates.cpp b/opm/simulators/wells/GroupTreeRates.cpp index 43c856d513e..2186f9c8746 100644 --- a/opm/simulators/wells/GroupTreeRates.cpp +++ b/opm/simulators/wells/GroupTreeRates.cpp @@ -21,30 +21,30 @@ #include #include -#include namespace Opm { template -int distributeGroupTreeRates(std::vector>& tree, - int rootIndex, - int maxIter) +int GroupTreeRates:: +distribute(std::vector>& tree, + int root_index, + int max_iter) { - if (rootIndex == 0) { + if (root_index == 0) { // Initialise: all nodes group-controlled except root for (auto& node : tree) { node.status = -1; } - tree[rootIndex].status = 1; - tree[rootIndex].rate = tree[rootIndex].limit; + tree[root_index].status = 1; + tree[root_index].rate = tree[root_index].limit; } int iter = 0; - for (; iter < maxIter; ++iter) { - setSubRates(tree, rootIndex); + for (; iter < max_iter; ++iter) { + set_sub_rates(tree, root_index); - Scalar worstExcess{0}; - const int ix = findWorstOffendingChild(tree, rootIndex, worstExcess); + Scalar worst_excess{0}; + const int ix = find_worst_offending_child(tree, root_index, worst_excess); if (ix < 0) { break; // converged } @@ -55,143 +55,129 @@ int distributeGroupTreeRates(std::vector>& tree, // If it has children, recursively solve its subtree if (!tree[ix].children.empty()) { - distributeGroupTreeRates(tree, ix, maxIter); + distribute(tree, ix, max_iter); } - updateParentStatus(tree, ix); + update_parent_status(tree, ix); } return iter; } template -void setSubRates(std::vector>& tree, - int nodeIndex) +void GroupTreeRates:: +set_sub_rates(std::vector>& tree, + int node_index) { - const auto& children = tree[nodeIndex].children; + const auto& children = tree[node_index].children; if (children.empty()) { return; // leaf (well) node } // Separate children into fixed (status != -1) and group-controlled - Scalar fixedRate{0}; - Scalar guideSum{0}; - bool allFixed = true; + Scalar fixed_rate{0}; + Scalar guide_sum{0}; + bool all_fixed = true; for (const int ci : children) { if (tree[ci].status != -1) { - fixedRate += tree[ci].rate; + fixed_rate += tree[ci].rate; } else { - guideSum += tree[ci].guideRate; - allFixed = false; + guide_sum += tree[ci].guide_rate; + all_fixed = false; } } - if (allFixed) { - tree[nodeIndex].status = 0; - tree[nodeIndex].rate = fixedRate; + if (all_fixed) { + tree[node_index].status = 0; + tree[node_index].rate = fixed_rate; return; } - const Scalar availRate = tree[nodeIndex].rate - fixedRate; - assert(availRate >= Scalar{0}); - assert(availRate > Scalar{0} || guideSum == Scalar{0}); + const Scalar avail_rate = tree[node_index].rate - fixed_rate; + assert(avail_rate >= Scalar{0}); + assert(avail_rate > Scalar{0} || guide_sum == Scalar{0}); for (const int ci : children) { if (tree[ci].status == -1) { - tree[ci].rate = (guideSum > Scalar{0}) - ? availRate * tree[ci].guideRate / guideSum + tree[ci].rate = (guide_sum > Scalar{0}) + ? avail_rate * tree[ci].guide_rate / guide_sum : Scalar{0}; - setSubRates(tree, ci); + set_sub_rates(tree, ci); } } } template -int findWorstOffendingChild(const std::vector>& tree, - int nodeIndex, - Scalar& worstExcess) +int GroupTreeRates:: +find_worst_offending_child(const std::vector>& tree, + int node_index, + Scalar& worst_excess) { - const auto& children = tree[nodeIndex].children; - int worstIndex = -1; + const auto& children = tree[node_index].children; + int worst_index = -1; for (const int ci : children) { const Scalar excess = tree[ci].rate - tree[ci].limit; - if (excess > worstExcess) { - worstExcess = excess; - worstIndex = ci; + if (excess > worst_excess) { + worst_excess = excess; + worst_index = ci; } } // Recurse into children that have subtrees for (const int ci : children) { if (!tree[ci].children.empty()) { - Scalar subtreeExcess = worstExcess; - const int childWorst = findWorstOffendingChild(tree, ci, subtreeExcess); - if (subtreeExcess > worstExcess) { - worstExcess = subtreeExcess; - worstIndex = childWorst; + Scalar subtree_excess = worst_excess; + const int child_worst = find_worst_offending_child(tree, ci, subtree_excess); + if (subtree_excess > worst_excess) { + worst_excess = subtree_excess; + worst_index = child_worst; } } } - if (worstExcess > Scalar{0} && worstIndex >= 0) { - assert(tree[worstIndex].status == -1); + if (worst_excess > Scalar{0} && worst_index >= 0) { + assert(tree[worst_index].status == -1); } else { - worstIndex = -1; + worst_index = -1; } - return worstIndex; + return worst_index; } template -void updateParentStatus(std::vector>& tree, - int nodeIndex) +void GroupTreeRates:: +update_parent_status(std::vector>& tree, + int node_index) { - const int pix = tree[nodeIndex].parent; + const int pix = tree[node_index].parent; if (pix < 0) { return; } const auto& siblings = tree[pix].children; - bool allFixed = true; - Scalar totalRate{0}; + bool all_fixed = true; + Scalar total_rate{0}; for (const int ci : siblings) { if (tree[ci].status == -1) { - allFixed = false; + all_fixed = false; break; } - totalRate += tree[ci].rate; + total_rate += tree[ci].rate; } - if (allFixed) { + if (all_fixed) { tree[pix].status = 0; - tree[pix].rate = totalRate; - updateParentStatus(tree, pix); + tree[pix].rate = total_rate; + update_parent_status(tree, pix); } } // Explicit template instantiations +template class GroupTreeRates; +template class GroupTreeRates; + template struct GroupTreeNode; template struct GroupTreeNode; -template int distributeGroupTreeRates( - std::vector>&, int, int); -template int distributeGroupTreeRates( - std::vector>&, int, int); - -template void setSubRates( - std::vector>&, int); -template void setSubRates( - std::vector>&, int); - -template int findWorstOffendingChild( - const std::vector>&, int, double&); -template int findWorstOffendingChild( - const std::vector>&, int, float&); - -template void updateParentStatus( - std::vector>&, int); -template void updateParentStatus( - std::vector>&, int); - } // namespace Opm diff --git a/opm/simulators/wells/GroupTreeRates.hpp b/opm/simulators/wells/GroupTreeRates.hpp index 780a5f8adf5..d328f57400f 100644 --- a/opm/simulators/wells/GroupTreeRates.hpp +++ b/opm/simulators/wells/GroupTreeRates.hpp @@ -17,16 +17,15 @@ along with OPM. If not, see . */ -#ifndef OPM_GROUP_TREE_RATES_HPP -#define OPM_GROUP_TREE_RATES_HPP +#ifndef OPM_GROUP_TREE_RATES_HEADER_INCLUDED +#define OPM_GROUP_TREE_RATES_HEADER_INCLUDED -#include #include #include namespace Opm { -/// Simplified node in the group/well hierarchy used for rate distribution. +/// \brief Simplified node in the group/well hierarchy used for rate distribution. /// /// Each node represents either a group or a well. Wells are leaf nodes /// (children is empty). The root node typically represents the FIELD group. @@ -42,66 +41,79 @@ struct GroupTreeNode { Scalar rate{0}; ///< Current allocated rate Scalar limit{0}; ///< Rate limit for this node int status{-1}; ///< Control status (-1, 0, or 1) - Scalar guideRate{0}; ///< Guide rate used for allocation + Scalar guide_rate{0}; ///< Guide rate used for allocation std::string name; ///< Group or well name std::vector children; ///< Indices of child nodes }; -/// Distribute rates through a group tree respecting individual node limits. +/// \brief Distribute rates through a group tree respecting individual node limits. /// -/// Starting from the root node (at @p rootIndex), rates are distributed -/// top-down to children according to their guide-rate fractions. Whenever -/// a child's allocated rate exceeds its limit, that child is fixed at its -/// limit and the distribution is recomputed for the remaining -/// group-controlled children. The process repeats until no limit -/// violations remain. -/// -/// @param tree Flat vector of tree nodes (modified in-place). -/// @param rootIndex Index of the subtree root in @p tree. -/// @param maxIter Maximum number of outer iterations (safety limit). -/// @return Number of outer iterations used. +/// Under the simplifying assumption that well-rate fractions are constant, +/// this class iteratively solves the group tree to a converged state. +/// Starting from the root node, rates are distributed top-down to children +/// according to their guide-rate fractions. Whenever a child's allocated +/// rate exceeds its limit, that child is fixed at its limit and the +/// distribution is recomputed for the remaining group-controlled children. +/// The process repeats until no limit violations remain. template -int distributeGroupTreeRates(std::vector>& tree, - int rootIndex = 0, - int maxIter = 1000); +class GroupTreeRates +{ +public: + /// \brief Distribute rates through a group tree respecting individual node limits. + /// + /// Starting from the root node (at \p root_index), rates are distributed + /// top-down to children according to their guide-rate fractions. Whenever + /// a child's allocated rate exceeds its limit, that child is fixed at its + /// limit and the distribution is recomputed for the remaining + /// group-controlled children. The process repeats until no limit + /// violations remain. + /// + /// \param[in,out] tree Flat vector of tree nodes (modified in-place). + /// \param[in] root_index Index of the subtree root in \p tree. + /// \param[in] max_iter Maximum number of outer iterations (safety limit). + /// \return Number of outer iterations used. + static int distribute(std::vector>& tree, + int root_index = 0, + int max_iter = 1000); -/// Distribute the rate of node @p nodeIndex to its children. -/// -/// Children that are already individually limited (status != -1) keep -/// their current rate. The remaining (available) rate is split among -/// group-controlled children proportionally to their guide rates. -/// If all children are fixed, the parent's status becomes 0 and its -/// rate is set to the sum of the children's rates. -/// -/// @param tree Flat vector of tree nodes (modified in-place). -/// @param nodeIndex Index of the node whose children are updated. -template -void setSubRates(std::vector>& tree, - int nodeIndex); + /// \brief Distribute the rate of a node to its children. + /// + /// Children that are already individually limited (status != -1) keep + /// their current rate. The remaining (available) rate is split among + /// group-controlled children proportionally to their guide rates. + /// If all children are fixed, the parent's status becomes 0 and its + /// rate is set to the sum of the children's rates. + /// + /// \param[in,out] tree Flat vector of tree nodes (modified in-place). + /// \param[in] node_index Index of the node whose children are updated. + static void set_sub_rates(std::vector>& tree, + int node_index); -/// Find the child (in the subtree rooted at @p nodeIndex) whose rate -/// most exceeds its limit. -/// -/// @param tree Flat vector of tree nodes (read-only access). -/// @param nodeIndex Root of the subtree to search. -/// @param[out] worstExcess The largest (rate - limit) value found. -/// @return Index of the worst offending node, or -1 if no -/// violation exists. -template -int findWorstOffendingChild(const std::vector>& tree, - int nodeIndex, - Scalar& worstExcess); + /// \brief Find the child whose rate most exceeds its limit. + /// + /// Searches the subtree rooted at \p node_index for the node whose + /// rate most exceeds its limit. + /// + /// \param[in] tree Flat vector of tree nodes (read-only access). + /// \param[in] node_index Root of the subtree to search. + /// \param[in,out] worst_excess The largest (rate - limit) value found. + /// \return Index of the worst offending node, or -1 + /// if no violation exists. + static int find_worst_offending_child(const std::vector>& tree, + int node_index, + Scalar& worst_excess); -/// After fixing a child node, propagate upward: if all siblings are -/// also fixed, the parent becomes fully determined (status 0) with -/// rate equal to the sum of its children. -/// -/// @param tree Flat vector of tree nodes (modified in-place). -/// @param nodeIndex Index of the node whose parent is updated. -template -void updateParentStatus(std::vector>& tree, - int nodeIndex); + /// \brief Update parent status after fixing a child node. + /// + /// If all siblings are also fixed, the parent becomes fully determined + /// (status 0) with rate equal to the sum of its children. + /// + /// \param[in,out] tree Flat vector of tree nodes (modified in-place). + /// \param[in] node_index Index of the node whose parent is updated. + static void update_parent_status(std::vector>& tree, + int node_index); +}; } // namespace Opm -#endif // OPM_GROUP_TREE_RATES_HPP +#endif // OPM_GROUP_TREE_RATES_HEADER_INCLUDED diff --git a/tests/test_GroupTreeRates.cpp b/tests/test_GroupTreeRates.cpp index 424c948bde8..a7c384477b3 100644 --- a/tests/test_GroupTreeRates.cpp +++ b/tests/test_GroupTreeRates.cpp @@ -112,7 +112,7 @@ BOOST_AUTO_TEST_CASE(SingleWell_WellLimitBinding) { auto tree = makeSingleWellTree(); // FIELD limit = 500, W1 limit = 300 => W1 should be limited to 300 - const int iter = distributeGroupTreeRates(tree); + const int iter = GroupTreeRates::distribute(tree); BOOST_CHECK(iter > 0); // W1 should be at its limit @@ -129,7 +129,7 @@ BOOST_AUTO_TEST_CASE(SingleWell_FieldLimitBinding) { auto tree = makeSingleWellTree(); tree[1].limit = 800.0; // W1 limit > FIELD limit - const int iter = distributeGroupTreeRates(tree); + const int iter = GroupTreeRates::distribute(tree); BOOST_CHECK(iter >= 0); // W1 rate should equal FIELD limit since it is the only well @@ -139,7 +139,7 @@ BOOST_AUTO_TEST_CASE(SingleWell_FieldLimitBinding) BOOST_AUTO_TEST_CASE(Unconstrained_RatesFollowGuideRates) { auto tree = makeUnconstrainedTree(); - distributeGroupTreeRates(tree); + GroupTreeRates::distribute(tree); // FIELD rate = 1000 (its limit) const double fieldRate = 1000.0; @@ -157,7 +157,7 @@ BOOST_AUTO_TEST_CASE(Unconstrained_RatesFollowGuideRates) BOOST_AUTO_TEST_CASE(SimpleTree_GroupAndWellLimits) { auto tree = makeSimpleTree(); - distributeGroupTreeRates(tree); + GroupTreeRates::distribute(tree); // All well rates must not exceed their limits for (const auto& node : tree) { @@ -175,7 +175,7 @@ BOOST_AUTO_TEST_CASE(SimpleTree_GroupAndWellLimits) BOOST_AUTO_TEST_CASE(SimpleTree_RateConsistency) { auto tree = makeSimpleTree(); - distributeGroupTreeRates(tree); + GroupTreeRates::distribute(tree); // G1 rate == W1 rate + W2 rate BOOST_CHECK_CLOSE(tree[1].rate, tree[3].rate + tree[4].rate, 1e-10); @@ -188,7 +188,7 @@ BOOST_AUTO_TEST_CASE(SimpleTree_RateConsistency) BOOST_AUTO_TEST_CASE(SimpleTree_StatusConsistency) { auto tree = makeSimpleTree(); - distributeGroupTreeRates(tree); + GroupTreeRates::distribute(tree); // Root should be individually limited (status 1) BOOST_CHECK_EQUAL(tree[0].status, 1); @@ -215,7 +215,7 @@ BOOST_AUTO_TEST_CASE(TightWellLimit) auto tree = makeUnconstrainedTree(); tree[3].limit = 10.0; // W1 has tight limit - distributeGroupTreeRates(tree); + GroupTreeRates::distribute(tree); // W1 must be at its limit BOOST_CHECK_CLOSE(tree[3].rate, 10.0, 1e-10); @@ -237,7 +237,7 @@ BOOST_AUTO_TEST_CASE(AllWellsTight) tree[5].limit = 100.0; // W3 tree[6].limit = 100.0; // W4 - distributeGroupTreeRates(tree); + GroupTreeRates::distribute(tree); // Each well should be at its limit BOOST_CHECK_CLOSE(tree[3].rate, 100.0, 1e-10); @@ -262,7 +262,7 @@ BOOST_AUTO_TEST_CASE(GroupLimitViolation) tree[3] = {3, 1, 0.0, 5000.0, -1, 1.0, "W1", {}}; tree[4] = {4, 2, 0.0, 9000.0, -1, 1.0, "W2", {}}; - const int iter = distributeGroupTreeRates(tree); + const int iter = GroupTreeRates::distribute(tree); // Should take at least 1 iteration to fix the G1 violation BOOST_CHECK(iter >= 1); @@ -299,7 +299,7 @@ BOOST_AUTO_TEST_CASE(CascadingViolations) tree[5] = {5, 2, 0.0, 3000.0, -1, 0.5, "W3", {}}; tree[6] = {6, 2, 0.0, 9000.0, -1, 0.5, "W4", {}}; - const int iter = distributeGroupTreeRates(tree); + const int iter = GroupTreeRates::distribute(tree); BOOST_CHECK(iter >= 1); // All limits respected @@ -326,7 +326,7 @@ BOOST_AUTO_TEST_CASE(DeepTree) tree[2] = {2, 1, 0.0, 2000.0, -1, 1.0, "G2", {3}}; tree[3] = {3, 2, 0.0, 1500.0, -1, 1.0, "W1", {}}; - distributeGroupTreeRates(tree); + GroupTreeRates::distribute(tree); // W1 is the tightest at 1500 BOOST_CHECK_CLOSE(tree[3].rate, 1500.0, 1e-10); @@ -342,7 +342,7 @@ BOOST_AUTO_TEST_CASE(ZeroIterations_NoViolation) tree[0] = {0, -1, 0.0, 0.0, -1, 1.0, "FIELD", {1}}; tree[1] = {1, 0, 0.0, 100.0, -1, 1.0, "W1", {}}; - const int iter = distributeGroupTreeRates(tree); + const int iter = GroupTreeRates::distribute(tree); // Should converge immediately (no violations since rates are 0) BOOST_CHECK_EQUAL(iter, 0); @@ -358,7 +358,7 @@ BOOST_AUTO_TEST_CASE(EqualGuideRates) tree[2] = {2, 0, 0.0, 500.0, -1, 1.0, "W2", {}}; tree[3] = {3, 0, 0.0, 500.0, -1, 1.0, "W3", {}}; - distributeGroupTreeRates(tree); + GroupTreeRates::distribute(tree); // Each well gets 300 BOOST_CHECK_CLOSE(tree[1].rate, 300.0, 1e-10); @@ -373,7 +373,7 @@ BOOST_AUTO_TEST_CASE(SetSubRates_Basic) tree[0].status = 1; tree[0].rate = 1000.0; - setSubRates(tree, 0); + GroupTreeRates::set_sub_rates(tree, 0); BOOST_CHECK_CLOSE(tree[1].rate, 600.0, 1e-10); BOOST_CHECK_CLOSE(tree[2].rate, 400.0, 1e-10); @@ -387,13 +387,13 @@ BOOST_AUTO_TEST_CASE(FindWorstOffending_Basic) // Set FIELD rate and distribute tree[0].status = 1; tree[0].rate = 10000.0; - setSubRates(tree, 0); + GroupTreeRates::set_sub_rates(tree, 0); - double worstExcess = 0.0; - const int ix = findWorstOffendingChild(tree, 0, worstExcess); + double worst_excess = 0.0; + const int ix = GroupTreeRates::find_worst_offending_child(tree, 0, worst_excess); // At least one child should be offending since rates > limits for some - BOOST_CHECK(ix >= 0 || worstExcess <= 0.0); + BOOST_CHECK(ix >= 0 || worst_excess <= 0.0); } BOOST_AUTO_TEST_CASE(UpdateParentStatus_AllFixed) @@ -406,7 +406,7 @@ BOOST_AUTO_TEST_CASE(UpdateParentStatus_AllFixed) tree[2].rate = 300.0; tree[0].status = -1; - updateParentStatus(tree, 1); + GroupTreeRates::update_parent_status(tree, 1); BOOST_CHECK_EQUAL(tree[0].status, 0); BOOST_CHECK_CLOSE(tree[0].rate, 800.0, 1e-10); @@ -420,7 +420,7 @@ BOOST_AUTO_TEST_CASE(UpdateParentStatus_NotAllFixed) tree[2].status = -1; // still group-controlled tree[0].status = -1; - updateParentStatus(tree, 1); + GroupTreeRates::update_parent_status(tree, 1); // Parent should NOT be updated since not all children are fixed BOOST_CHECK_EQUAL(tree[0].status, -1); From 04b1686e97f94fc78ba96bb49d9257f68e5aad2d Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 2 Mar 2026 13:11:56 +0000 Subject: [PATCH 6/8] Reorder template instantiations: struct before class for dependency order Co-authored-by: steink <6746854+steink@users.noreply.github.com> --- opm/simulators/wells/GroupTreeRates.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/simulators/wells/GroupTreeRates.cpp b/opm/simulators/wells/GroupTreeRates.cpp index 2186f9c8746..cfd636a78ed 100644 --- a/opm/simulators/wells/GroupTreeRates.cpp +++ b/opm/simulators/wells/GroupTreeRates.cpp @@ -174,10 +174,10 @@ update_parent_status(std::vector>& tree, } // Explicit template instantiations -template class GroupTreeRates; -template class GroupTreeRates; - template struct GroupTreeNode; template struct GroupTreeNode; +template class GroupTreeRates; +template class GroupTreeRates; + } // namespace Opm From 556f24386fd380de70f38e6e3b1516fee0999537 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 2 Mar 2026 14:15:39 +0000 Subject: [PATCH 7/8] Replace TODO with actual GroupTreeRates integration code; keep original as updateWellControlsOriginal - Rename existing updateWellControls to updateWellControlsOriginal - New updateWellControls calls updateWellControlsOriginal, then: 1. Builds flat GroupTreeNode vector from the group hierarchy 2. Calls GroupTreeRates::distribute(tree) 3. Applies resulting rates back to well/group state - Add updateWellControlsOriginal declaration to BlackoilWellModel.hpp Co-authored-by: steink <6746854+steink@users.noreply.github.com> --- opm/simulators/wells/BlackoilWellModel.hpp | 3 + .../wells/BlackoilWellModel_impl.hpp | 189 +++++++++++++++++- 2 files changed, 182 insertions(+), 10 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 77938a5c028..1591391b1bd 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -274,6 +274,9 @@ template class WellContributions; bool updateWellControls(DeferredLogger& deferred_logger); + bool + updateWellControlsOriginal(DeferredLogger& deferred_logger); + void updateAndCommunicate(const int reportStepIdx); bool updateGroupControls(const Group& group, diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 860387c5008..7c35a928b02 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -65,7 +65,9 @@ #include #include +#include #include +#include #include #include @@ -1652,6 +1654,183 @@ namespace Opm { bool BlackoilWellModel:: updateWellControls(DeferredLogger& deferred_logger) + { + OPM_TIMEFUNCTION(); + const bool changed = updateWellControlsOriginal(deferred_logger); + + if (!this->wellsActive()) { + return changed; + } + + const int episodeIdx = simulator_.episodeIndex(); + const Group& fieldGroup = this->schedule().getGroup("FIELD", episodeIdx); + const auto& pu = this->phaseUsage(); + + // Determine the target phase based on the FIELD group production control + const auto cmode = this->groupState().production_control(fieldGroup.name()); + int target_phase_idx = -1; + switch (cmode) { + case Group::ProductionCMode::ORAT: + if (pu.phaseIsActive(IndexTraits::oilPhaseIdx)) + target_phase_idx = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); + break; + case Group::ProductionCMode::WRAT: + if (pu.phaseIsActive(IndexTraits::waterPhaseIdx)) + target_phase_idx = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx); + break; + case Group::ProductionCMode::GRAT: + if (pu.phaseIsActive(IndexTraits::gasPhaseIdx)) + target_phase_idx = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); + break; + default: + // For LRAT, RESV, FLD or other modes the simplified tree + // distribution is not applicable. + break; + } + + if (target_phase_idx < 0) { + return changed; + } + + // Determine the guide-rate target matching the production control mode + const auto guide_target = this->groupStateHelper().getProductionGuideTargetMode(fieldGroup); + if (guide_target == GuideRateModel::Target::NONE) { + return changed; + } + + // Step 1: Build a flat GroupTreeNode vector from the current group hierarchy + std::vector> tree; + + // Recursive lambda to traverse the group hierarchy and populate the flat vector + std::function buildTree = + [&](const Group& group, int parent_idx) -> int + { + const int node_idx = static_cast(tree.size()); + tree.push_back(GroupTreeNode{}); + auto& node = tree[node_idx]; + node.index = node_idx; + node.parent = parent_idx; + node.name = group.name(); + + // Get group rate limit from production controls + const auto controls = group.productionControls(this->summaryState()); + switch (cmode) { + case Group::ProductionCMode::ORAT: node.limit = controls.oil_target; break; + case Group::ProductionCMode::WRAT: node.limit = controls.water_target; break; + case Group::ProductionCMode::GRAT: node.limit = controls.gas_target; break; + default: node.limit = std::numeric_limits::max(); break; + } + + // Get current rate from group state + if (this->groupState().has_production_rates(group.name())) { + const auto& rates = this->groupState().production_rates(group.name()); + if (target_phase_idx < static_cast(rates.size())) { + node.rate = rates[target_phase_idx]; + } + } + + // Get guide rate for this group + const auto grv = this->groupStateHelper().getProductionGroupRateVector(group.name()); + if (this->guideRate().has(group.name())) { + node.guide_rate = this->guideRate().get(group.name(), guide_target, grv); + } else { + node.guide_rate = Scalar{1}; + } + + // Process child groups + for (const std::string& child_group_name : group.groups()) { + const auto& child_group = this->schedule().getGroup(child_group_name, episodeIdx); + const int child_idx = buildTree(child_group, node_idx); + tree[node_idx].children.push_back(child_idx); + } + + // Process child wells (leaf nodes) + for (const std::string& well_name : group.wells()) { + const auto& well_ecl = this->schedule().getWell(well_name, episodeIdx); + if (!well_ecl.isProducer()) { + continue; + } + + const int well_node_idx = static_cast(tree.size()); + tree.push_back(GroupTreeNode{}); + auto& well_node = tree[well_node_idx]; + well_node.index = well_node_idx; + well_node.parent = node_idx; + well_node.name = well_name; + well_node.children = {}; + + // Get well rate limit from production controls + const auto well_controls = well_ecl.productionControls(this->summaryState()); + switch (cmode) { + case Group::ProductionCMode::ORAT: well_node.limit = well_controls.oil_rate; break; + case Group::ProductionCMode::WRAT: well_node.limit = well_controls.water_rate; break; + case Group::ProductionCMode::GRAT: well_node.limit = well_controls.gas_rate; break; + default: well_node.limit = std::numeric_limits::max(); break; + } + + // Get current well rate from well state + if (this->wellState().has(well_name)) { + const auto& ws = this->wellState().well(well_name); + if (target_phase_idx < static_cast(ws.surface_rates.size())) { + // Production rates are negative in OPM convention + well_node.rate = -ws.surface_rates[target_phase_idx]; + } + } + + // Get guide rate for this well + const auto wrv = this->groupStateHelper().getWellRateVector(well_name); + if (this->guideRate().has(well_name) || this->guideRate().hasPotentials(well_name)) { + well_node.guide_rate = this->guideRate().get(well_name, guide_target, wrv); + } else { + well_node.guide_rate = Scalar{1}; + } + + tree[node_idx].children.push_back(well_node_idx); + } + + return node_idx; + }; + + buildTree(fieldGroup, -1); + + if (tree.empty()) { + return changed; + } + + // Step 2: Distribute rates through the group tree + GroupTreeRates::distribute(tree); + + // Step 3: Apply the resulting rates back to the well/group state + for (const auto& node : tree) { + if (node.children.empty()) { + // Leaf node = well: apply rate back to well state + if (this->wellState().has(node.name)) { + auto& ws = this->wellState().well(node.name); + if (target_phase_idx < static_cast(ws.surface_rates.size())) { + // Store as negative (production convention) + ws.surface_rates[target_phase_idx] = -node.rate; + } + } + } else { + // Group node: update group state production rates + if (this->groupState().has_production_rates(node.name)) { + auto rates = this->groupState().production_rates(node.name); + if (target_phase_idx < static_cast(rates.size())) { + rates[target_phase_idx] = node.rate; + this->groupState().update_production_rates(node.name, rates); + } + } + } + } + + return changed; + } + + + template + bool + BlackoilWellModel:: + updateWellControlsOriginal(DeferredLogger& deferred_logger) { OPM_TIMEFUNCTION(); if (!this->wellsActive()) { @@ -1724,16 +1903,6 @@ namespace Opm { // update wsolvent fraction for REIN wells this->updateWsolvent(fieldGroup, episodeIdx, this->nupcolWellState()); - // TODO: Consider using GroupTreeRates::distribute() here to solve - // the group-tree rates to a more converged state under the simplifying - // assumption of constant well-rate fractions. This would involve: - // 1. Building a flat GroupTreeNode vector from the current group hierarchy - // (using fieldGroup, schedule, guide rates, and current well/group rates) - // 2. Calling GroupTreeRates::distribute(tree) - // 3. Applying the resulting rates back to the well/group state - // This approach could improve convergence of group-tree rate balancing - // without requiring network pressure recalculation. - return changed_well_group; } From 30f133d4ba1c5a07ea1b64132975a08530cda6ea Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 2 Mar 2026 14:16:46 +0000 Subject: [PATCH 8/8] Add comments explaining default guide rate fallback value of 1.0 Co-authored-by: steink <6746854+steink@users.noreply.github.com> --- opm/simulators/wells/BlackoilWellModel_impl.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 7c35a928b02..909c6a5cacb 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -1734,6 +1734,8 @@ namespace Opm { if (this->guideRate().has(group.name())) { node.guide_rate = this->guideRate().get(group.name(), guide_target, grv); } else { + // Default to 1.0 so that children without explicit guide rates + // share the parent's rate equally. node.guide_rate = Scalar{1}; } @@ -1782,6 +1784,8 @@ namespace Opm { if (this->guideRate().has(well_name) || this->guideRate().hasPotentials(well_name)) { well_node.guide_rate = this->guideRate().get(well_name, guide_target, wrv); } else { + // Default to 1.0 so that wells without explicit guide rates + // share the group's rate equally. well_node.guide_rate = Scalar{1}; }