Skip to content

Implement injection networks - #7345

Closed
atgeirr wants to merge 34 commits into
OPM:masterfrom
atgeirr:injection-network
Closed

Implement injection networks#7345
atgeirr wants to merge 34 commits into
OPM:masterfrom
atgeirr:injection-network

Conversation

@atgeirr

@atgeirr atgeirr commented Aug 21, 2026

Copy link
Copy Markdown
Member

This is not ready for review or user testing yet, but is made a draft PR to ease developer testing and improvement.

Atgeirr Flø Rasmussen and others added 30 commits August 20, 2026 16:20
A network branch lookup clamped the rate and upstream pressure to nothing:
VFPHelpers::findInterpData extrapolates linearly past the axis ends, and the
flow-line tables are zero-filled where the line cannot deliver the rate, so an
injection network with rates beyond the table axis produced node pressures of
-260 bar, and a THP below the axis start gave 0. Both were then handed to the
wells as THP limits.

The branch calculator now clamps the lookup point to the table axes (rates
scaled uniformly so prod WFR/GFR are kept) and treats a result <= 1 atm as
'no solution'. Such nodes, and their descendants, are reported by
NetworkPressureComputation::invalidNodes(); updatePressures() keeps their
previous pressure, counts the network as unbalanced (max update bound), warns
once per report step, and does not push a dynamic THP to wells under them.

Also make updateActiveStateImpl accumulate instead of clearing active_ on an
inactive domain (unreachable today, but wrong).

Tests: gas_injection_rate_beyond_flow_axis, gas_injection_zero_cell_region and
gas_injection_thp_below_axis in test_networkpressure.cpp; on 0996672 they
give -213 bar, 0 bar and 28.9 bar respectively.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The explicit node-pressure iteration P <- P + w (P_computed - P) cannot
converge an injection network: with BHP pinned by the reservoir the injectors'
THP-mode rate is very steep in THP (~6e4 sm3/d/bar on GNETINJE_GAS-01) and the
flow-line pressure falls ~0.5 bar per 1000 sm3/d, a loop gain of ~30, so with
w=0.1 the iteration falls into a period-2 limit cycle (leaf inflow 0 <-> 1.5e6,
pressure 500 <-> 0). Wells stopped at the wrong moment were then shut for good.

NodePressureUpdater (new header) replaces the update per node: keep a
sign-change bracket on r(P) = P_computed(P) - P and take Illinois regula-falsi
steps inside it; before a bracket exists take a secant step clipped to
[P, P_computed], capped at max(25% P, 10 bar), and, when the leaf's wells are on
group/rate control (flat response down to their own THP), go straight to just
below that kink. Stale bracket ends (from wells that were momentarily stopped,
or from the other leaves moving) are dropped when contradicted, when the bracket
collapses, or when their Illinois weight is halved away. Convergence is judged
on the remaining pressure uncertainty (bracket width), not on the residual,
which is amplified by the well response. Enabled by default; the old damped
update is the fallback and is kept behind --network-pressure-update-secant=false.

Well side: the injector's dynamic THP is clamped to its VFPINJ THP axis instead
of being silently left stale; injectors get the previous step's node pressure at
report-step start like producers (starting at the WCONINJE THP made them inject
at their rate limit and starve the other leaves); the 'Option B' potentials
refresh is removed (it wrote well_potentials = surface_rates, so a well that had
been stopped restarted every THP solve at zero rate and was stopped again).

Result on opm-tests/network/GNETINJE_GAS-01 and _WAT-01: 0 unconverged network
messages, no wells shut, and rates/THP/control modes match the E100 reference
(opm-tests/eclref) to within 1% / 0.2 bar at every report step.

Tests: NodePressureUpdaterTests in test_networkpressure.cpp (stiff response
converges where the damped update does not, stale end dropped, plateau floor,
step cap, invalid evaluation moves down).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The physical-shut path is meant for wells that cannot operate under their own
THP/BHP limits. An injector whose THP is a network node pressure may be unable
to inject at this iteration's pressure and yet inject fine once the network is
balanced (the pressure depends on the other wells' rates); shutting it makes the
network converge to a wrong state with the remaining wells carrying the whole
group target. Keep such wells stopped for the step; the local solve reopens
them when the pressure allows. Injectors only for now.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
GAS and WATER injection networks were written into the same leaf-rate map, so
a group that is a leaf of both would see the second network's rates, and the
injection phase threaded into NetworkPressureComputation was ignored. Key the
map by (Phase, group) and use the phase in the calculator; the mock in
test_networkpressure.cpp now returns only the network's own phase and its
phase-less lookup returns nothing, so the gas/water cases fail without this.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
assignNodeAndBranchValues() only exported the production network. Fill the
gas/water injection maps of data::GroupAndNetworkValues (opm-common
c4b9da0a4) from the per-domain node pressures and branch data as well.

On GNETINJE_GAS-01, GPRG now matches the E100 reference: 340.0 / 209.4 /
209.4 / 204.3 / 204.3 bar (PLAT-A / M5S / G1 / M5N / F1) at day 31.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
… default

--network-pressure-update-secant becomes a string: injection (default), all,
none. With the bracketing update also on production networks NETWORK-01 gives
a slightly different (better converged) solution than the reference (BPR
0.04 bar at one step); keep production networks on the damped update until the
references are regenerated. Production-network results are now identical to
the pre-change scheme.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…upported

Take GNETINJE off the unsupported-keyword list (it made the decks fail at the
default parsing strictness) and register the two opm-tests injection network
decks as regression tests (--solver-max-time-step-in-days=1). The references
are generated with this branch; both cases agree with the E100 runs in
opm-tests/eclref within 2% on node pressures and 3% on well rates at every
report step (tools/plot_injection_network.py --check).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
An injection group target (GCONINJE VREP/REIN) is a function of the produced
voidage, so within a network sub-iteration the injectors must see this
iteration's production, not the previous one's. The well loop ran in deck
order, so that was a matter of luck.

Producers are now re-solved first, then the injectors. The group data is
refreshed in between only when a deck has both a production and an injection
network -- the only case where the producers re-solved here can feed an
injection target in the same sub-iteration -- so nothing else pays for the
extra update. details::activeNetworks() already returns the production
network first, so the pressure computation was in the right order already.

No deck in opm-tests has both network kinds, so this is a correctness guard
rather than a fix for an observed failure: NETWORK-01, NETWORK-01_STANDARD,
NETWORK-01-WTEST, GNETINJE_GAS-01 and GNETINJE_WAT-01 are all bit-identical
(compareECL -t SMRY at 1e-9) before and after.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…t help yet)

An experiment against the per-node bracketing update: Anderson uses the last few
(P, G(P)) pairs of the whole pressure vector of a network, so it sees the coupling
between nodes that a per-node update cannot, without needing any derivative.

Kept deliberately isolated: NetworkAndersonAcceleration.hpp is self-contained (no
OPM dependencies, own tiny least-squares solve), there is one call site in
updatePressures(), and it is off unless
--network-pressure-update-acceleration=anderson is given
(--network-anderson-depth sets the history length, default 4). With the default
the results are bit-identical to the previous commit on all five network decks.

Measured on the opm-tests decks, against the E100 reference:

  GNETINJE_WAT-01   bracketing: 3 unconverged steps   anderson(4): 0
  GNETINJE_GAS-01   bracketing: 0 unconverged steps   anderson(4): 18
  GNETINJE_GAS-01 WGIR:G-3H vs E100, over all reported steps
                    bracketing: median 0.24 %, p90 3.4 %, 86 % within 2 %
                    anderson(4): median 0.02 %, p90 100 %, 63 % within 2 %

So it is better where the response is smooth and much worse where it is not: the
gas case fails the report-step comparison with Eclipse that the bracketing update
passes. That is the expected failure of unsafeguarded extrapolation across the
kinks where wells switch control -- Anderson replaces the bracket, the step caps
and the plateau rule that the per-node updater relies on. It should not be
enabled as it stands; a variant that keeps those safeguards (limit the Anderson
step to the same bound, reject it when the residual grows) is the obvious next
thing to try.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
… yet)

Balance the injection networks against each injector's linearised rate response
instead of re-solving the well equations, so a residual evaluation costs a VFP
lookup. --network-well-proxy=ipr, default none; nothing changes when it is off.

Two things had to be fixed before this could even be measured.

The explicit IPR (ipr_a_/ipr_b_, WellInterface::updateIPR) is identically zero
for injectors: its crossflow guard drops every connection with
pressure_diff = p_r - h_perf > 0, which for an injector is the ordinary case.
Nothing noticed because the only consumers are the producer operability checks.

The implicit IPR is the usable one -- it differentiates the converged well
equation through the well Jacobian -- but updateIPRImplicit hard-coded the
producer control swap, while the rhs picks out the control equation and so
requires that equation to be the bhp one. Made symmetric here; the producer path
is unchanged.

With that, the response is exact at the point of linearisation (it reproduces
each injector's solved rate to four figures), and the proxy is still not good
enough: on GNETINJE_GAS-01 and _WAT-01 it takes unconverged network steps from
0/3 to 0/1 but the deviation from the E100 reference goes from 44/1 violations
to 620/646, because from about two thirds through both runs it drives every
injector to zero rate. dq/dbhp is large enough here that the intersection bhp
moves 0.04 bar for a 292 bar change in applied thp, so any node-pressure error
lands on the steep part, the linear rate goes negative and the well shuts. A
Newton on this Jacobian will need globalisation and the well's own rate/group
limits, not a full step.

Kept because the linearisation itself is sound and is the building block for
that Newton. See injection_network_findings.md.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
tests/test_networksolve.cpp extracts the GNETINJE_GAS-01 network with the wells
replaced by their inflow performance plus the control logic, so solution methods
can be compared without a reservoir or a well solve. The three VFPINJ tables are
the deck's, verbatim; the wells are calibrated to the Eclipse 100 operating point
at day 31, and the bench then reproduces it to 0.1 bar.

Table 9999 means "no table", so G1 carries M5S's pressure and F1 carries M5N's
and the whole problem has two unknowns. That is small enough to look at a
Jacobian by hand.

From the wells' WCONINJE THP, to 0.01 bar:

  damped (omega 0.1)     FAILED  (limit cycle -- the original branch's method)
  bracketing (shipped)   42
  anderson (depth 4)     21
  newton (full step)     FAILED  (overshoots off the plateau to 450/506 bar)
  newton + line search   14

and over dq/dbhp = 1e4 .. 1e6 sm3/d/bar, bracketing 29/42/30/38 against
newton+ls 10/14/21/22.

So a Newton on this problem is worth doing and must be globalised: the full step
diverges on its own, and the same line search that fixes it also keeps the
advantage as the wells stiffen.

Note for anyone extending this: VFPInjProperties::addTable stores a
reference_wrapper, so the tables have to outlive the properties object and must
not be moved.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The Newton direction was already there; what was missing was any measurement of
what to do with it. Four strategies now share one Newton driver -- full step,
step cap, backtracking line search (plain or Armijo), trust region -- so the only
thing that varies between them is how the step is accepted.

Iteration count from one good start says little about a globalisation, so the
bench sweeps a 23x23 grid of starting pressures across the tables' THP axis and
counts how many reach the right answer:

  bracketing (shipped)   529/529   mean 31 iterations
  newton, full step        3/529   mean  5
  newton, capped step     54/529   mean  8
  newton, line search    529/529   mean 12
  newton, trust region   529/529   mean 13

Two things worth having measured. Capping the step -- which is what
--network-max-pressure-update-in-bars already does -- is not a globalisation:
it recovers from a tenth of the space. And a real one costs nothing in
robustness, so the choice between the line search and the trust region is a
choice about code, and the line search is the smaller piece of code.

The trust region needed its collapse case handled to get there: when the
Jacobian is taken across a control switch the radius halves away to nothing, and
returning the unchanged point stalls the solver. Taking the smallest step and
reopening the region fixes it (83% -> 100%).

Also covers the group-target branch of the well control logic, which the earlier
tests left dead.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
… one

The bench solved the network in the node pressures alone, with every rate
recovered from them by an inner solve. This adds the same network without the
eliminations -- node pressures, branch rates and each well's (rate, bhp) as
unknowns, with the eliminations restated as node balances, branch drops, inflow
performance and one control equation per well -- and lets the same Newton solve
either. 16 unknowns against 2.

Over the grid of 529 starting pressures, at the measured stiffness:

  eliminated, plain newton     3/529   mean  5 iterations
  eliminated, line search    529/529   mean 12
  full,       plain newton   529/529   mean  9
  full,       line search    529/529   mean  9

The eliminated residual needs globalising because the control clamps put kinks
in it. The full system holds its controls fixed while the step is taken, so
there are no kinks and an unglobalised Newton is already fully robust -- and
faster. That is the case for carrying the rates.

Three things had to be right for that to hold, and each is a finding in its own
right:

Controls must be chosen by most-restrictive-wins, matching the clamp the
eliminated form applies. A fixed priority chain makes the active set chatter and
the Newton never terminates.

Globalisation must not veto an active-set change. The residual jumps when a
control switches and that is not a failure to make progress; letting the line
search reject it stalls the switch instead of resolving it (363/529 -> 529/529).

The table limits want to be bounds on the unknowns, not a clamp on the lookup.
Outside the tables' box the zero-filled cells extrapolate away and the full
system has a spurious root there -- at the softest wells a plain Newton reaches
it from a third of the grid, ending with every well at its rate limit and node
pressures of -683 and -10975 bar. Clamping the lookup to the axes, which is what
the simulator's pressure computation does, removes that root and replaces it
with a flat residual the Newton cannot descend (383/529 -> 23/529). Neither is
right; the limits belong in the active set beside the well controls.
table_bounds_want_to_be_constraints pins both halves down.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
… build

Three things, of which the first two were the open questions.

The table limits belong on the unknowns. Clamping the lookups to the axes, which
is what the simulator's pressure computation does, flattens the residual and
leaves a Newton nothing to descend (22/529 starts). Leaving the tables to
extrapolate keeps the derivatives but admits a root outside them (412/529).
Holding the branch flows inside the box does best (421/529), and only by
projecting the offending components -- a scalar fraction-to-boundary lets one
binding rate throttle the pressure updates as well, which is worse than doing
nothing (299/529).

The group target is now an equation with a multiplier, so wells the group does
not bind stay on their own controls. It converges in 4 iterations and lands on
the target exactly, with the network at higher pressure than it runs free, which
is what curtailment should do. One trap: the activation test has to include
equality, because at the solution a well's rate is its share exactly and a strict
test flips the control every iteration.

The case is now data rather than code. A NetworkCase is nodes with a parent and a
table, wells hanging off nodes, a terminal pressure and an optional group target;
gnetinjeGas() is one instance of it, assembled the way a deck reader would --
topology from GRUPTREE/GNETINJE, tables from the VFPINJ includes, limits from
WCONINJE, and the calibration point from the reference summary. Both problems and
all the solvers work off that, with no node names or unknown indices written into
them, so a second case is a builder function rather than a rewrite.

The refactor reproduces every previously measured number.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A case needs a reference operating point per well -- the rate it takes and the
pressure at its node -- because that is what its IPR is a linearisation about.
Both ways of supplying one are now there.

Reference::set() writes it down. That is what the built-in GNETINJE_GAS-01 case
uses, so the bench still needs no files and runs everywhere, and it is also how
you would pose a case no reference run covers.

Reference::fromSummary() reads it out of a reference run at a given time, paired
with fromDeck(), which takes the topology from GRUPTREE/GNETINJE, the tables from
the VFPINJ keywords and the limits from WCONINJE. It reads the keywords off the
parsed Deck rather than building a Schedule, so it stays usable from a test.

built_from_deck_matches_the_builtin_case checks the two agree, node for node and
well for well, and that the deck-built case solves to the same place. It skips
when opm-tests is absent.

Two things that cost time and are worth knowing:

WCONINJE's RATE is a UDA item whose unit depends on the injected phase, so the
parser leaves it dimensionless and getSI() hands back the raw deck number. Taken
at face value the rate limit came out 86400x too large, which silently removes
the rate control -- the case still looked right and simply stopped converging.

The operating point has to be read at a report step. Off the ministep vectors the
same nominal time gives a slightly different state, which is a different case.

Correcting the built-in reference against the summary at day 31 moved it 0.6% on
the G wells, and the bench now lands on (209.405, 204.240) against E100's
(209.409, 204.244) -- previously (209.30, 204.19).

Rates are carried as an (aqua, liquid, vapour) triple rather than one number,
which is what a production network will need: there the branch rate splits and
VFPPROD wants WFR and GFR fractions and an ALQ, and those fractions are extra
unknowns per branch with their own mixing equations at the nodes.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Until now the bench solved one hand-picked operating point. This runs the full
formulation at every report step of both GNETINJE decks, calibrated at each step
to what the Eclipse reference says the wells were doing, with the field target
applied where the reference has them on group control.

  GNETINJE_GAS-01   4/4 report steps, worst 0.05 bar, 4-16 iterations
  GNETINJE_WAT-01   6/6 report steps, worst 0.39 bar, 4-8 iterations

Both group-controlled steps land on the reference exactly. Those are the steps
where the method we ship is weakest -- day 91 of the gas case is where it spends
two reported steps at the validity floor.

Getting there needed two corrections, both of which matter beyond the bench.

A well must be calibrated against its own bhp, not the pressure at its node. The
two agree only while it is on THP control. At day 91 the wells sit 162 bar below
their node (WTHP 202 against GPRG 365), and calibrating there against the node
pressure builds a well that does not exist -- which is why that step would not
converge at all before.

The stiffness knob only transfers between cases as a fraction of the well's own
rate. 6e4 sm3/d/bar is reasonable for a gas injector taking 5e5 sm3/d and
meaningless for a water injector taking 700; carrying it across is what made
three of the water steps diverge. setRelativeStiffness(0.12) reproduces the gas
figure and works for both.

fromDeck also had to learn that GNETINJE spells the phase WAT where WCONINJE
spells it WATER.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
opm/simulators/wells/NetworkSystem.hpp now holds the injection network solved
simultaneously in its pressures and its rates: the branch drops, node balances,
inflow performance and per-well control equation, an optional group multiplier,
the active-set control selection, the branch-flow bounds and the Newton.

The bench is now a caller of it rather than a separate implementation. All it
still owns is where a well's inflow performance comes from -- a linearisation
about a reference operating point, where the simulator will use the well
Jacobian -- plus the eliminated formulation and the fixed-point methods it is
compared against. Every measured number is unchanged, which is the check that
the extraction is faithful.

setClampToAxes is on the system only so the comparison that rules it out can
still be made, and says so.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…newton

The simulator now fills in the same NetworkSolve::System the bench does. The
only thing that differs is where a well's rate response comes from: the bench
linearises about a reference operating point, the simulator reads the implicit
IPR out of the well state. Everything downstream -- the equations, the
active-set control selection, the branch-flow bounds, the Newton -- is shared.

Default is unchanged (fixedpoint), and both GNETINJE references still compare
clean. With --network-solver=newton, both decks at 1-day steps:

                  unconverged network steps    deviations from E100
  gas   fixedpoint          0                          44
  gas   newton              0                          10
  water fixedpoint          3                           4
  water newton              0                           4

so the gas case lands four times closer to Eclipse and the water case stops
taking unconverged steps, in slightly fewer Newton iterations either way.

The node pressures come back already at the fixed point, so the relaxation that
follows sees no imbalance and stops; the branch data from the ordinary
evaluation is kept for the output. Anything the solve cannot handle -- more than
one root, a well with no usable inflow performance, a network that does not
converge -- returns nullopt and leaves the fixed point in charge.

The group multiplier the system can carry is not used here. The simulator's own
group machinery has already decided each well's share by the time the network is
solved, so a well on GRUP enters as one held at that rate.

updateIPRImplicit is refreshed for injectors before the solve; the well solve
only maintains it for producers.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
--network-analytic-jacobian assembles the Jacobian from the derivatives the VFP
interpolation already computes and throws away, instead of differencing the
residual n+1 times. Everything but the two table lookups is constant. Checked
against the difference quotient in four states -- the ordinary THP one, a
converged one, one with the group holding the wells, and one driven onto the
bhp and rate limits -- and agreeing to about 1e-7 relative. On both decks it
gives the same answer as differencing, which is the point: it is a cost, not an
accuracy, change.

--network-group-control hands the network the group's injection total and lets
it place the split, so a well that runs into its own limit is taken up by the
others rather than the total quietly dropping. On these two decks it changes
nothing, because the group machinery has already made the split and no well is
near a limit, so the multiplier has nothing to redistribute. It is off by
default and unproven; a case where a group-controlled well hits its bhp limit is
what would show whether it earns its place.

Both decks, 1-day steps, against the Eclipse reference:

                            unconverged   deviations
  gas   fixedpoint               0            44
  gas   newton                   0            10
  gas   newton + analytic        0            10
  gas   newton + analytic + grp  0            10
  water fixedpoint               3             4
  water newton                   0             4
  water newton + analytic        0             4
  water newton + analytic + grp  0             4

Also fixes two things in the previous commit's well setup. A group-controlled
well is held at the rate the group gave it, which had been lost; without it the
gas case went from 10 deviations to 252. And a well with neither a rate nor a
target now takes the whole network back to the fixed point instead of entering
with an invented limit -- std::max(rate, 1) is 86400 sm3/d in SI, which is no
limit at all for a water injector taking 700, and cost the water case 4
deviations against 240.

The solve is now started from the wells' current rates where the caller knows
them, so the first control selection is made on the state the well is actually
in.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A well's share of a group total should be proportional to what it can inject at
the pressure the network gives it. The previous version used its current rate,
which is the split being decided -- so the allocation reproduced whatever it
already was, which is why --network-group-control changed nothing.

System::thpPotential() meets a well's inflow performance with its tubing curve
at a given node pressure. The guides are refreshed from it each iteration and
have to settle before the solve is called converged, the same way the controls
do.

It fixes the node pressures at the gas case's day 61 -- GPRG:M5N/F1 were 2.6 %
high and are now inside tolerance, taking the gas deviations from 10 to 8 -- and
leaves the well rates bit-identical.

That last part is the useful half of the result. Under group control the network
writes nothing to a well but a dynamic THP limit, and for a group-held well that
limit does not bind, so the network cannot move the split however good its guide
rates are. The rates come from the group allocation, which runs before the
network and does not see it. Closing that needs the allocation to become part of
the same solve, not a better guide handed to it.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The simultaneous solve read the rank-local well container and communicated
nothing, so under MPI each rank would assemble a network out of its own wells
and get its own answer. It was latent because the solver is opt-in and had only
been run serially.

The relaxed computation stays right by working from group and node quantities,
which are the same everywhere by the time it runs. This does the same thing one
level down. The well list and everything static about a well -- its node, its
vfp table, its bhp and rate limits -- now come from the schedule, which is
replicated. Only what a well is currently doing is rank-local: its inflow
performance, its rate and its control mode are contributed once by the rank that
owns it and summed, after which every rank holds the same numbers and reaches
the same decisions, including the decision to give up and leave it to the fixed
point.

Contributed by the owner rather than by every holder because a distributed well
appears on several ranks carrying the same values; summing those would count it
twice.

On GNETINJE_GAS-01 the node pressures now agree across 1, 2 and 4 ranks to about
1e-5 relative, and are identical to six figures on the group-controlled steps.
The deviation count against the Eclipse reference is 8 on all three, and 4 on
the water case serial and on two ranks. What is left between serial and parallel
is a connection total on a shut well that the default solver differs on in
exactly the same way, so it is not this.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
hnil and others added 4 commits August 20, 2026 16:20
Every way out of newtonNodePressures() returned nullopt without a word, so there
was no way to tell whether the simultaneous solve had run or the relaxed update
had quietly done the work. Each exit now names its cause -- more than one root,
no terminal pressure, the phase not active, no injectors on the network, a well
with no usable inflow performance, a well with neither a rate nor a target, or
simply not converging -- and a successful solve reports its iteration count.

Debug level, because this runs once per network sub-iteration; it lands in the
.DBG file next to the network trace.

It paid for itself immediately. On the two decks the solve is used for most of
the run but not all of it: gas 282 solved against 21 handed back (6.9 %), water
344 against 37 (9.7 %), nearly all of them non-convergence in 50 iterations.
Neither run's answer changes -- the relaxed update covers those steps, which is
what the fallback is for -- but a tenth of the water case was taking a path
nobody could see.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A non-converged solve now carries its residual, whether a control was still
switching, whether the group shares were still moving, and one letter per well
per iteration for the last eight iterations. The caller prints all of it.

That was enough to identify what the fallbacks actually are. Every one of them
is an active set that will not settle, and the trace shows a clean period-2
cycle:

  water   TTTT GGGG TTTT GGGG TTTT GGGG TTTT GGGG
  gas     TTRR TTGG TTRR TTGG   and   RRTR GGTG RRTR GGTG

so the wells are flipping between group and thp control, not converging slowly.
The share a well is tested against moves with the guides and with the multiplier
while its rate is chasing that share, and the two never meet.

Two candidate fixes were measured and neither is right yet, so both are written
down rather than shipped.

Making the bhp and rate activation tests inclusive, as the group one already is,
looks correct -- a well at a limit sits exactly on it -- and is not: start()
clamps the opening rate to the limit, so an inclusive test latches every well
onto rate control at the first iteration and holds it there, which sends the
bench to the out-of-table root at -683 bar. The tests stay strict, with a comment
saying why.

Deciding group membership from the group rather than the well -- it binds when
the wells could between them take more than the target -- does cure the cycling,
taking the water case from 37 fallbacks to 2 and the gas case from 21 to 7. It
also moves the gas answer away from the reference, 8 deviations to 26. Not a
trade worth making blind, so it is described in the code where the next person
will find it.

Nothing here changes an answer: gas 8 deviations and water 4, as before.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The simulator refreshes each well's share of a group total from what it can
inject at the network pressure. The bench never turned that on, so the one code
path that fails in the simulator was the one path with no coverage.

It has coverage now, in the simulator's own group configuration -- the target is
the total the wells are already injecting, the guides are their current rates, so
the multiplier starts at one and every well starts exactly on its share, right on
the activation boundary. It converges, in three iterations, on the same answer as
holding the guides fixed.

So the GRUP/THP cycling is not the guide logic by itself. The half the bench
cannot supply is the inflow performance: the simulator's comes from the well
Jacobian at a start-up or post-control-change state, and every one of its
failures falls at exactly those. Reproducing one needs the failing system dumped
and replayed here, which is the next thing to build rather than the next thing to
guess.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
--network-dump-failures=<prefix> writes every network system that fails to
converge: nodes, wells with the inflow performance the well Jacobian gave them,
limits, group target and starting pressures. Point OPM_NETWORK_DUMP at the
directory and tests/test_networksolve.cpp replays each one against the same
tables, with no simulator. The failures reproduce exactly, cycle for cycle, in
about a second.

It found its first bug immediately. The group residual summed every well's rate
while the target was accumulated only over the wells the group was actually
holding, so the group was being asked to account for rates it does not control --
a constraint nothing can satisfy, which the active set then thrashes against
instead of settling. 26 of the 146 captured systems converge once the residual
counts only what the group holds.

Guide rates are now taken once per solve rather than refreshed each iteration.
That is what they are in the simulator -- explicit, set per timestep -- and
refreshing them inside the Newton makes each well's share a moving target while
its rate is chasing it. It takes the water case with group control from 37 solves
handed back to 2.

The solve also opens just inside a well's limits rather than exactly on them, so
the first control selection is not made on a rate that clamping put there.

No answer changes: gas 8 deviations from the reference with group control and 10
without, water 4 either way, both references still comparing clean by default.

One thing measured and not taken: making the bhp and rate activation tests
inclusive, as the group one is, takes the gas case from 70 solves handed back to
4 and the water case from 39 to 2 -- and sends the bench to the out-of-table root
from two thirds of its starting points. The simulator and the bench disagree
about it, which is worth understanding before either is believed.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@atgeirr

atgeirr commented Aug 24, 2026

Copy link
Copy Markdown
Member Author

Closing in favour of #7348 which is a revised version.

@atgeirr atgeirr closed this Aug 24, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

manual:new-feature This is a new feature and should be described in the manual

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants