[WIP] Simultaneous network solve: injection and production networks, autochoke, gas-lift response - #7348
Draft
hnil wants to merge 68 commits into
Draft
[WIP] Simultaneous network solve: injection and production networks, autochoke, gas-lift response#7348hnil wants to merge 68 commits into
hnil wants to merge 68 commits into
Conversation
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>
newtonProductionNodePressures() is the production counterpart of the injection
adapter: build a ProductionSystem from the replicated schedule, take each well's
per-phase implicit IPR off the well Jacobian summed once by the owning rank, and
hand back node pressures. The well state holds q = b*bhp - a in opm's signed
rates where production is negative; the system wants production positive and
falling with bhp, which is the same line negated.
Three things had to be fixed to make it agree with the relaxed update it now
runs alongside.
The implicit IPR was refreshed only for injectors, on the grounds that the well
solve maintains it for producers. It maintains it only where its own control
logic happens to need it: 103 of 376 production network solves arrived with a
well they could not linearise and handed the whole network back. Refreshing it
for every prediction well takes that to zero.
thpPotential() bracketed between the bhp limit and the bhp at which *oil* stops.
Water and gas have their own zero crossings and are still flowing there, so the
tubing still has something to lift, the bracket did not contain the crossing, and
the well read as unable to produce at all. It now brackets to the last phase to
stop.
A zero from thpPotential() means the well cannot lift against this node pressure
-- its table does not reach that high. That makes thp *unavailable*, not a
control allowing nothing: as an allowance of zero it won "most restrictive" every
time, and the thp row it then imposed says nothing about a rate, so the well
produced whatever the tubing crossing happened to be. On NETWORK-01, whose wells
are tabulated for a thp of 20 bar while the node sits near 90, that was every
well on every step: 3730 sm3/d against a 1000 limit, and a node 5 bar high.
Wells not under a network-held group target are pinned at the rate they already
have rather than at the deck's WCONPROD limit. Re-deriving the operating point
would overwrite one the group and the well solve have already agreed on, and a
limit that does not bind in the simulator would become one here.
NETWORK-01 365 solved, 0 fell back, no deviation from the relaxed
NETWORK-01_STANDARD 365 / 0, no deviation
NETWORK-01-REROUTE 212 / 346, no deviation; every fallback is BRANPROP
re-routing giving the network more than one
root, which this formulation does not model
injection gas 151 / 0, water 364 / 0, unchanged
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ot model A network with several roots is a forest of independent trees -- every node has one parent, so the trees share nothing. The adapters now take a root and solve each tree; whichever converges replaces the relaxed answer for its own nodes. That was the only cause of fallback on seven decks: both REROUTE variants (the BRANPROP re-route leaves the old root behind), both MULTIROOT variants, WTEST, 6_UDA_MODEL5 and both GSATPROD decks all reported "the network has more than one root" and nothing else. Being able to converge is not the same as being right, so every deck with a network was then compared against the relaxed update, which found four node features the system silently mishandles. Each is now refused by name instead: - a fixed-pressure node below the root (MULTIROOT's GRPB at 82 bar): a boundary the extended network holds fixed, which this would compute from the branch table -- worth 7 % on field rates; - gas lift added at a node (NODEPROP item 4) and wells under WLIFTOPT: the alq is the wells' own, not the branch constant -- GASLIFT-13 got a simulation 128 % off with 4x fewer Newton iterations, which is not a faster solve but a different problem; - satellite production (GSATPROD): rates that arrive without wells -- 9-22 % off; - a well efficiency factor: scales the branch contribution but not the well's own rate, and the system has one rate per well. With the guards, every deck in the regression suite that the solve accepts agrees with the relaxed update to 0.5 % at every report time -- the seven production decks exactly, water within its known day-92 transient, gas closer to the E100 reference than the relaxed update is (8 deviating pairs against 44, which include its node-pressure collapse to 1 bar). The guarded decks run unchanged, every fallback naming its reason. Of the ten regression-covered network decks: seven solve outright with zero fallbacks, two are refused entirely (gas lift at a node, satellites), one (GASLIFT-13-style) partially. Nothing is silently wrong. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A review before taking on gas lift and autochoke, which turned up that most of
what the adapter was refusing is a constant inside one network solve and only
needed a slot:
- WEFAC and NEFAC. A well's own rate is q; the branch above it sees
efficiency * q, and a node passes on efficiency * what it collects. Both
systems carry the factors in their balance rows, the injection Jacobian in
its, and the dump format. The injection adapter had been ignoring them
silently; the production one refused. The well factor is the network one
(WEFAC item 3) times the per-rank scaling, summed once by the owner.
- Lift gas. Added to the gas stream at the leaf when NODEPROP item 4 says so,
scaled by the well's efficiency, and not part of q -- it goes up the tubing,
it is not produced. Each well carries its own; a satellite's is the node's.
- Satellite production. A constant triple at the node; the group sum returns it
instead of the group's wells, so those wells are left out to match.
- A well's alq is what its own tubing table sees, which is alq_state, not the
WCONPROD item -- the same number unless an optimiser has moved it.
And one defect in the production control rule: thpPotential() reported exactly
what the bhp limit allows when thp did not bind, which tied, and the tie went to
thp -- whose row then settled the bhp below the limit the deck set. It now
reports more than the limit allows, and the well lands on bhp control at its
limit.
before after vs relaxed
GSATPROD6 (regression) 0 / 1156 refused 1194 / 0 no deviation
GSATPROD5 0 / 1146 1184 / 0 no deviation
6_UDA_MODEL5 (regr.) 0 / 440 528 / 0 FOPR identical,
FWPR 0.015 sm3/d
NETWORK-01-WEFAC 303 / 0 (factors 1) 303 / 0 no deviation
Wells under gas lift optimisation stay refused. With the alq taken as a
constant the solves converge, and the runs then diverge from the relaxed update
by 100 % and more over the schedule (GASLIFT-13/14): the optimiser reacts to the
node pressures it sees between passes, and that coupling is the real gas lift
work, not the constant. The guard names it.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
NETWORK_MODEL5_STDW_AUTOCHK aborted at day 3.2, and NETWORK_MODEL5_STDW_AUTOCHK with GROUPGUIDERATES with it. The autochoke block decides whether the target is the group's own or an ancestor's from the *deck* control (ORAT here), then asks the group *state* for the target -- and the state says NONE for a group under its limit, which getProductionGroupTargetForMode_ throws on. One layer down the TargetCalculator was built from the same state cmode and asserted on NONE in calcModeRateFromRates. A group under its limit is simply a choke that ends up open; read the target off the control already in hand and build the calculator on it. Also guard the upstream-node lookup on the first pass of a run, when the node pressure map is still empty: dereferencing end() there handed the group a garbage pressure. The deck now runs to the end (790 Newton iterations); the GROUPGUIDERATES variant gets much further and then loses its timestep, which is a separate problem. Both fail identically on upstream master without this. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
--network-autochoke=true, with --network-solver=newton, hands autochoke nodes to
the simultaneous network solve instead of the legacy search. Defaults stay as
they were; without the option a tree with a choke node is refused by name and
the legacy block keeps the node.
In the system a choke is one row and one flag. The node's pass-through row
becomes "oil collected here equals the target", and its pressure is the group's
common thp -- the multiplier, with no lambda. Closed when the wells behind it
could deliver more than the target with the valve open, at the upstream
pressure; open otherwise. Because the wells' controls have to be chosen at the
pressure the choke settles at and not at the iterate's, the selection first
finds that pressure from the wells' allowances alone -- the cheap model the
control rule already runs on -- and the Newton only polishes it. No well solves
anywhere: the legacy search calls iterateWellEqWithSwitching on every well in
the group per sample of a brute-force bracket, up to 1300 samples a pass.
Getting there exposed four things the production system lacked, each found by
comparing with the legacy answer at a converged point:
- the hydrostatic correction between a tubing table's datum and the well's
reference depth, which every well's own thp evaluation applies -- 13 bar on
these wells, and the reason the tubing looked twice as permeable as it is.
Both systems carry it now, computed on the typed side with the well's density
and shared through the owner-summed entries;
- the liquid-loading hump: at low rates a tubing table needs *more* pressure
than at moderate rates, so h(bhp) is negative at both ends of the bracket and
positive between, and a bisection sees "cannot lift" for a well that lifts
fine. thpPotential scans for the crossing where h turns positive with rising
bhp, the one the well settles on; the other is the loading point;
- a start consistent with the control rule: opening every well at the bhp limit
put a well whose limit is the 1 atm default at a rate off every table;
- a margin in the control choice, so the marginal well a choke leaves exactly
where its tubing passes its own limit does not flip every iteration.
Two rules about what the network may decide. A well it is not deciding for --
not under a choke, not under a network-held group target -- is pinned at its
current rate and not offered thp; with every well pinned a tree has nothing to
decide and the solve is skipped, because it would only reproduce the relaxed
evaluation and walk the wells down a different path to it, which a well test
downstream can turn into a different decision. And a well the well model has
at zero rate is dead at its thp: the linearised inflow with the table still
finds a flowing crossing for it, and more back-pressure cannot revive it.
Nodes a tree solve has placed go through the bracketing update rather than the
damped one: the damped update creeps ten per cent per sub-iteration toward an
answer that does not change and runs out the cap -- a hundred solves a Newton
iteration -- while handing the wells the whole jump at once is how one gets
shut as inoperable in a transient.
NETWORK_MODEL5_STDW_AUTOCHK, new path vs legacy
B1 oil day 1 5929 / 6019 day 31 5997 / 6001
day 60 6001 / 6016 day 91 4904 / 4846
node pressure within 0.4 bar throughout; 4118 solves, 34 fell back, 43 s
against 790 Newton iterations and 895 brute-force brackets in 17 s.
Open: B-1H is shut by the well model's operability test at day 91 (day 152
on the legacy path), after which the two diverge. Slower than legacy for now:
the scan costs 96 table lookups per well per control selection.
Every other deck is unchanged to the digit -- the seven NETWORK-01 variants,
6_UDA_MODEL5, the MSW deck -- and the injection decks keep their E100 match.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Inside the network sub-loop the wells are frozen, so the system the
simultaneous solve is handed comes back unchanged sub-iteration after
sub-iteration -- and was being solved again each time, from a different guess:
at best a repeat (1490 of 4118 solves converged in two iterations, the guess
already the answer), at worst a different root, which then kept the bracketing
update busy for a dozen sub-iterations. Keep the inputs and the answer per tree
root and hand the answer back while the inputs match.
NETWORK_MODEL5_STDW_AUTOCHK, same build, before and after, against legacy:
legacy before after
well solves 467054 11636 7612
network solves - 4118 870
reservoir Newton its 790 313 288
fell back - 34 0
wall clock 15 s 43 s 12 s
Within 1 % of the legacy answer to day 91; B-1H is shut by the well model at
day 106 against day 152, after which they diverge -- still the open item.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
--gas-lift-network-response=true, with --network-solver=newton, answers the gas
lift optimiser's trial evaluations from the simultaneous network solve instead
of a well solve at a fixed thp. The optimiser is untouched: same increments,
same economic gradient, same weights and limits, same stage-2 redistribution.
Only the oracle changes -- which is the whole point, since every gradient it
takes is currently at a node pressure that does not respond to the lift gas
being tried.
A trial copies the last solved tree, sets the well's alq, re-solves, and reads
the well's potential off the node pressure that state gives. The well keeps its
limits while the network is solved, so the state is one the field could be in;
the potential is then taken at that pressure, because the optimiser asks for a
potential and applies the limits itself. Freeing the well instead -- the first
thing I tried -- lets it flow at a rate it is never allowed, moves the node
pressure somewhere it never sits, and the optimiser then applies its limit on
top of a potential taken at the wrong pressure.
What it costs, GASLIFT-13, against the legacy oracle:
legacy network-answered
well solves 72 920 10 773
reservoir Newton its 2 733 446
wall clock 61 s 6 s
What it is not: validated. The two oracles disagree, and I have not established
which is right. Within one optimiser pass (same well, same ipr) C-1H reads
2492 / 2982 sm3/d oil at alq 12500 / 25000 from the network, against 2203 / 897
from a well solve -- the legacy answer collapsing by 60 % as the lift gas
doubles. It is not the alq axis running out: the table's axis is 0 to 219000 and
both trials sit inside its first interval. It could be the well genuinely going
unstable, or the legacy bracket landing on another root; judging it needs an
independent measurement -- fixed-alq runs, no trial machinery on either side --
which is the next step and wants no more guessing before it.
So: default off, and the flag's help says experimental. The simulation it
produces differs from the legacy one by more than a rounding, and until the
oracle is judged that difference is not evidence of anything.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
--network-autochoke-bracket-samples, default 300, which is the value the search
has always used and so changes nothing. Each sample solves every well in the
group, and the search runs on every outer network iteration, so this count is
where the legacy autochoke's cost lives: on NETWORK_MODEL5_STDW_AUTOCHK,
samples well solves reservoir Newton its wall
300 475 143 790 15 s
12 151 033 492 11 s
6 30 528 271 3 s
The answers differ -- by 29 % and 40 % in late-time field oil against the 300
run. That is the finding: the legacy autochoke's result depends on how finely
its bracket is sampled, so it is not a converged quantity and 300 is not more
right than 12, only the value everyone has been running. The new simultaneous
path tracked the 300 run to day 91 and then parted from it at B-1H's shut; the
legacy variants part from each other in the same region. The late-time
behaviour of this deck is fragile under every method, including legacy against
itself.
Also drops a temporary gas lift diagnostic.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The choke's root-find asks for each well's thp allowance at dozens of pressures per control selection, and each answer is a 96-sample scan of the tubing table. Keep the answers on a quarter-bar grid and interpolate where the potential is smooth; where it is not -- one side cannot lift, or is unbounded, or the two grid points differ by more than five per cent, which is the stable crossing jumping along the loading hump -- evaluate exactly. Without that last guard the interpolation invented allowances across the jump, the choke pressure ran to 41 bar and every B well was shut by day 97. With the guard the autochoke deck is unchanged to the digit (288 iterations, 7613 well solves, the same shut day) and no faster: with inflow slopes near 1000 sm3/d per bar the potential moves ten per cent per cell and the guard is almost always taken. GASLIFT-13, whose wells are gentler, goes from 6 s to 4 s at the same well-solve count. The remaining cost is the scan itself. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The three rate derivatives of every table lookup come from differentiating the lookup itself -- VFPProdProperties::bhp on DenseAd::Evaluation<Scalar, 3>, the instantiation that exists -- so the chain rule through FLO, WFR and GFR types is the library's, not a hand-written one per type. The thp derivative is one more lookup a hundredth of a bar up: the table is piecewise linear in thp, so inside an interval that is the derivative. Balance, inflow, control, group and choke rows are structural. Checked against differences on a plain network, under a group target and with a closed choke: largest mismatch 0.0056 on entries of 14400. Opt-in through --network-analytic-jacobian, which already existed for the injection side; the default still differences, so nothing changes by itself. Measured on NETWORK_MODEL5_STDW_AUTOCHK: 265 reservoir Newton iterations and 6624 well solves against 288 and 7613 with differences, the same answer -- and 50 network solves that fail where none did, every one of them a well flipping between thp and its own rate limit, TTTO / TOTO. The exact Jacobian lands the well precisely on the tie; the differenced one was inaccurate enough to nudge it off. Locking a tied well on its rate limit removed every fallback and gave the best counts on both decks, and shut every B well by day 121; locking it on the tubing gave the closest late-time answer with this Jacobian and doubled the gas lift deck's iterations. Neither lock is in. A tie is a complementarity condition -- both rows hold there -- and it wants a formulation, not a margin; that is the next thing to do on the control rule. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Three things, in the order the bench forced them. A Fischer-Burmeister row for a well on the tie between its rate limit and its tubing, Control::Tied: a + b - sqrt(a^2 + b^2) on the scaled rate slack and tubing slack, with the generalised derivative at the origin. Only with the assembled Jacobian -- the row is not smooth there, and a difference straddling the kink collapsed the differenced path from 0 fallbacks to 678. A tied well whose lift gas is changed for a trial is released: a different alq is a different tubing curve, and the tie went with it. Production dump and replay. write()/readProduction() carry everything a solve depends on -- sources, choke targets, efficiencies, the hydrostatic term, lift gas, pinning, dead wells -- and the bench replays a directory of them against the MODEL5 tables with both Jacobians and compares the answers. The adapter writes tree failures, trial failures (which used to fall back to the well solve without a trace), and with OPM_NETWORK_DUMP_ALL=N the first N solved systems as well, because a converged answer can still be the wrong root. What the replay said about GASLIFT-13, where the analytic Jacobian turns a 302-iteration run into a 592-iteration one with 51 % different late-time oil: of 241 systems the analytic run solved, 240 replay identical with both Jacobians (gaps under 0.02 bar) and the differenced one fails one; the differenced run's four tree failures are the tie, which the tied row converges; its three trial failures fail both ways because two wells tie at once, which one row cannot resolve. So the two Jacobians do not reach different network answers. They differ in seven solves that one converges and the other does not, and seven trial answers are enough to send a greedy, discrete optimiser down another path -- deterministically (592 twice). That is the optimiser's sensitivity, not the Jacobian's error, and it is now measurable standalone. Also a bench case proving both Jacobians reach the same solution on plain, grouped, choked and tied systems from fourteen starts each (gaps 1e-10 bar). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
--network-complementarity=true (with --network-analytic-jacobian=true) closes
every production well that has limits of its own with one row instead of an
active set: the box complementarity q in [0, limit] against the tubing slack
bhp - tubing(p, q) and the bhp slack bhp - bhp_limit, as nested
Fischer-Burmeister functions. q = 0 when the tubing needs more than the
reservoir gives, q = limit when both slacks are positive, a slack at zero in
between; nothing switches, and two wells tying at once need no treatment.
Eligibility is decided once per solve from the starting point -- deciding it
from the iterate re-introduced the switching the row exists to remove -- and
complementarity wells start from the rate the well model has, because a tubing
table with a loading hump gives a well two branches and the start chooses.
With it, bounds by projection for the production step (no node below an
atmosphere, no bhp below its limit, no pressure step over 50 bar), which the
active-set rows had not needed and the complementarity row did: its first full
step from a poor start ran a choke node to minus three thousand bar.
Standalone, in the bench:
complementarity_agrees_with_the_active_set 5 shapes x 14 starts: converges
from every start, agrees with the active set to 3e-2 % or better
replay_production_failures, 7 GASLIFT-13 dumps: differenced 0/7, analytic
active set 4/7, complementarity 7/7 -- the three two-well ties included
40 autochoke failure dumps (earlier run): 33/40 in 9 iterations, to the
physically right answer where the active set had three wells producing
2500 sm3/d through tubing that cannot lift them.
In the simulator it is not yet better:
NETWORK_MODEL5_STDW_AUTOCHK 134 Newton, 6208 well solves, 156 fell back
(active set + analytic: 223 / 5097 / 0)
GASLIFT-13 350 Newton, 8402 well solves, oil per lift gas
4.4 (active set, differenced: 302 / 7456 / 12.3)
Every remaining failure is one class, and the bench reproduces it: the choke
is decided closed by the allowance test while the complementarity wells sit on
their dead branch and cannot reach the target. The choke has to become a
complementarity too -- p_node - p_up >= 0 against target - sum q >= 0 -- rather
than a pre-solved state. That is the next item. Default off.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Under --network-complementarity the choke node is one row, fb(p_n - p_open, target - q_n), with p_open the branch's own table value when the choke node has one; no open/closed pre-solve. The outer box on the well row is dropped (degenerate when the bhp limit binds); whether a well flows or is dead is the scan's answer at the current node pressure, inside the row. Control::Shut (q = 0) for a bhp limit at or above shut-in, in both formulations: the active set used to report injection there, and production_thp_that_does_not_bind_leaves_the_well_on_bhp had asserted it. limitStep caps a complementarity well's oil step and re-seeds a rate at or below zero to its allowance. OPM_NETWORK_CM_TRACE prints the iteration. Bench: complementarity_agrees_with_the_active_set, 7 shapes x 14 starts all agree (both chokes, ties, bhp above reservoir). AUTOCHK replay still 26/1396: the real systems hold a well at its hump tangency, which the FB row cannot linearise. Decks not better; numbers in the findings note.
Four faults found on the AUTOCHK failure dumps with the new per-iteration trace (OPM_NETWORK_CM_TRACE): the re-seed set the rate but not the bhp, so the well never sat on its IPR; a well reviving as the pressure it left behind fell made a system with no fixed point, so a well that loses its crossing as the pressure rises stays shut for the solve; wells the well model has dead were handed to the active set, which produced 2500 m3/d from each of them under a 6000 choke target, so under the flag they are Shut; and one alpha for the whole step let a choke row's pressure demand cut the shut rows to 0.1 % of their rate per iteration, so pressures are clamped componentwise and rates move freely. The scan starts where the table's flow axis ends, and a choke over wells that cannot respond has no row. Replay of 1396 AUTOCHK dumps: 1396/1396 in 2-5 iterations (was 26). AUTOCHK deck: 0 fallbacks, 201 Newton, 8 % more oil than the active set, but 2.5x the well solves on the one dying well -- the shut has to be carried across the well model's iterations next.
A well the well model has at zero rate has no usable IPR, and the production adapter refused the whole system for it -- every one of GASLIFT-13's 280 fallbacks. Under --network-complementarity such a well is Shut in the system (Well::shut; q = 0 with an IPR, bhp = limit without). GASLIFT-13: 298 Newton / 6583 well solves / 0 fallbacks against 424 / 8885 / 108 for the active set; the 4.6x lift gas meets the deck's LIFTOPT gradient (0.036 vs 5e-3), the active set's does not. Carrying a network shut across the well model's iterations was tried per report step and per time step; both chop the step (AUTOCHK 27 -> 74-88 time steps). Kept behind OPM_NETWORK_CARRY_SHUT, off. Tests: complementarity_shuts_dead_wells (a dying well, a choke over dead wells) and the_dumps_behind_the_complementarity_fixes_converge (four AUTOCHK dumps under tests/network_dumps, OPM_VFP_INCLUDE); both fail on 2e318a256.
The default path gives a well the well model has at zero rate bhp control and produces from it. Asserted as it is, so a fix flips the test instead of passing by.
Anderson acceleration (--network-pressure-update-acceleration, --network-anderson-depth) and the well-index proxy balance (--network-well-proxy, --network-well-proxy-max-iterations) were both added measured and off by default. Neither is worth carrying further. Anderson is better than the per-node bracketing update where the response is smooth and much worse where it is not: on GNETINJE_WAT-01 it takes 3 unconverged network steps to 0, on GNETINJE_GAS-01 it takes 0 to 18 and fails the report-step comparison with E100 that bracketing passes. It replaces the bracket, the step cap and the plateau rule with unsafeguarded extrapolation, which is exactly what the control switches punish. A safeguarded variant would be new code, not this code. The proxy balance is superseded: it linearised each injector against its implicit IPR and relaxed the node pressures against that, and the same linearisation is now solved simultaneously with the pressures in NetworkSystem.hpp, which does converge. On its own it drove every injector to zero rate two thirds of the way through both runs (E100 violations 44/1 -> 620/646). What that commit fixed in updateIPRImplicit stays: it was producer-only by an accident of the hard-coded control swap, and the simultaneous solve needs the injector side of it. Numbers and the full write-up are in injection_network_findings.md; the code is on the wip/network-2026-08-24 tag. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
atgeirr
reviewed
Aug 24, 2026
atgeirr
left a comment
atgeirr
left a comment
Member
There was a problem hiding this comment.
Starting with a few quick comments. Nowhere near finished of course!
I will also add here that the NetworkSystem.hpp file is huge, maybe better to split it into more files if a logical organization can be found, consider putting it in a network/ subdir.
| std::map<std::string, std::vector<Scalar>> m_production_rates; | ||
| std::map<std::string, std::vector<Scalar>> m_network_leaf_node_injection_rates; | ||
| // Injection networks are per phase (GNETINJE GAS / WAT); a group can be a leaf of both. | ||
| std::map<std::pair<Phase, std::string>, std::vector<Scalar>> m_network_leaf_node_injection_rates; |
Member
There was a problem hiding this comment.
It seems unnecessary to have separate storage per phase, since the value objects in the map are small vectors they can store all phase injection rates in a single such object.
Member
Author
There was a problem hiding this comment.
Agreed, done in b5a50cb. The value already covers every phase and the writer fills it the same way whichever network asks, so the pair key stored the same vector twice; the per-phase masking only existed in the test mock. Back to your node-keyed map, both GNETINJE decks unchanged to the digit.
| int node = 0; | ||
| int vfp_table = 0; | ||
| /// Inflow performance, q = ipr_a + ipr_b * bhp. The simulator's implicit | ||
| /// IPR stores it as q = b*bhp - a, so ipr_a is the negated one. |
Member
There was a problem hiding this comment.
Why is a different convention used than in the existing code? Would prefer that conventions are the same unless for a good reason.
Member
Author
There was a problem hiding this comment.
Agreed, done in 1b32904: q = ipr_b*bhp - ipr_a as in WellState, and the negation in the adapter is gone. Sign-only, results identical.
| /// Dense square system. The networks this solves have tens of unknowns, so | ||
| /// Gaussian elimination with partial pivoting is the whole story. | ||
| template<class Scalar> | ||
| class DenseMatrix |
Member
There was a problem hiding this comment.
This class reimplements functionality already available in Dune, prefer to use existing library functionality if equivalent.
Member
Author
There was a problem hiding this comment.
Done in ace7891 — it now wraps Dune::DynamicMatrix. I kept a thin wrapper for (i,j) indexing and to return false on a singular matrix rather than throw, since that case just hands the network back to the relaxed update. Counts identical everywhere; AUTOCHK FOPT moves 1.7e-5 relative over 182 days from the different elimination order.
The injection leaf rates were keyed on (phase, node), but the value is already a vector over all phases and the writer fills it the same way whichever network asks -- so a node that is a leaf of both the gas and the water network stored the same vector twice. Nothing needed the split. The VFP injection table picks the phase its FLO type names, so handing it every phase's rate is what the lookup already expects; the per-phase masking only ever existed in the test's mock, which is now a plain pass-through. Restores atgeirr's original node-keyed map, keeping the has_() accessor added since. Both GNETINJE decks are unchanged to the digit. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The network system stored q = ipr_a + ipr_b * bhp while WellState stores q = implicit_ipr_b * bhp - implicit_ipr_a, so the adapter negated one term on the way in and the struct carried a comment explaining it. No reason for the difference; use OPM's convention throughout and drop the negation. Sign-only, exact in floating point. The dumped-system fixture in the bench carries the same column, so its ipr_a values are negated with it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
DenseMatrix reimplemented Gaussian elimination with partial pivoting, which Dune::DynamicMatrix already provides. Keep the class as a thin wrapper for two reasons only: (i,j) indexing at the assembly sites, and a bool for a singular matrix instead of an exception -- a singular Jacobian means this network hands back to the relaxed update, and throwing out of the well model would be worse. GNETINJE_GAS-01 and GASLIFT-13 are unchanged to the digit, iteration and well-solve counts identical everywhere. AUTOCHK's FOPT moves 1.7e-5 relative over 182 days at identical counts -- Dune's elimination orders the arithmetic differently, and this is the only change in the series that can move a bit. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
atgeirr
reviewed
Aug 24, 2026
| /// Solve any of the systems in this file. They differ in what a rate is -- one | ||
| /// number for an injection network, three for a production one -- but not in how | ||
| /// the Newton, the active set or the bounds work. | ||
| template<class Sys, class Globalisation = FullStep> |
Member
There was a problem hiding this comment.
Please remove default template and function arguments. Avoid them, except as a means to reduce the need to change existing code, and even then avoid unless it leads to very large changes.
Member
Author
There was a problem hiding this comment.
Done in 4bad252. Tolerance and iteration cap are now a Parameters struct and the globalisation is always named; two overloads instead of defaults, the two-argument one spelling out FullStep{}. The bench is the only caller that varies any of it.
| } | ||
| } | ||
|
|
||
| /// Solve any of the systems in this file. They differ in what a rate is -- one |
Member
There was a problem hiding this comment.
The comment should explain that this solves using a Newton-Raphson method. Do not refer to "this file". The rate comment is unneeded as nowhere in this function does the rate size matter.
Member
Author
There was a problem hiding this comment.
Done in d37d21e: it now says Newton-Raphson with an active set over the well controls, and the file reference and the rate remark are gone.
| return out; | ||
| }; | ||
|
|
||
| // Guide rates are explicit: the simulator sets them once per timestep, and |
Member
There was a problem hiding this comment.
Check this comment and the call to refreshGuides(). Is it necessary, and if it is, when and why? It is my belief that as of now guide rates are explicit and there should not be a need to refresh.
Member
Author
There was a problem hiding this comment.
You are right that OPM's guide rates are explicit and this does not touch them. This is a different quantity that unfortunately borrows the name: it only does anything under --network-group-control, where the network places the group split itself and a well's share is its own potential at the node pressure — not knowable until the starting pressures are, hence once here rather than at setup. Inside the Newton it cycles, which is why it is outside. Comment rewritten in d37d21e; happy to rename it if you have a better word.
Member
Author
There was a problem hiding this comment.
Added a note in ca143b4 that explicit is not meant as the end state: solving for the share — either as an unknown of the system, or by iterating an outer loop to a fixed point — is the way past the cycling, and it is also what would let the network correct a group split rather than only report pressures. Not in this PR.
| // is moving, and a capped allowance ties with it. Now that the | ||
| // iterate is not transient, drop the cap and carry on from here; | ||
| // whoever is over the line goes on rate control and the rest take | ||
| // it up. Once round only. |
Member
There was a problem hiding this comment.
I do not understand "Once round only", and I do not see any code that should make the enforcement happen only once. Please explain better.
Member
Author
There was a problem hiding this comment.
It is the enforcing flag two lines down: set on the one pass that drops the cap, and tested in the same condition, so a well still over its limit afterwards converges as it is. Comment now names it (d37d21e).
|
|
||
|
|
||
| template<class Scalar> | ||
| class System |
Member
There was a problem hiding this comment.
Please document this class, and rename it InjectionSystem (assuming that I am correct that this is the counterpart to ProductionSystem).
Member
Author
There was a problem hiding this comment.
Done in d37d21e — InjectionSystem, with a class comment saying what the unknowns and equations are and that it is the counterpart to ProductionSystem.
`System` is the counterpart to `ProductionSystem`, so call it `InjectionSystem` and document what it holds and what closes it. solve()'s comment now says it is a Newton-Raphson with an active set over the well controls, rather than describing the file. Two comments answered rather than changed, because the code was already right and the comment was not: - refreshGuides() is a no-op unless the network is placing a group's split itself (--network-group-control). It is not OPM's guide rates, which stay explicit and set once per timestep; it is each well's own potential at the node pressure, which cannot be known before the starting pressures are. Once before the Newton, not inside it. - "Once round only" is enforced by the `enforcing` flag two lines below; the comment now names it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Tolerance and iteration cap move into a `Parameters` struct whose members carry
the values, and the globalisation is always named. Two overloads instead of
defaults: the two-argument one is what the simulator uses and spells out
`FullStep{}` itself; the bench, which is the only caller that varies any of
this, passes them explicitly.
Results unchanged -- both GNETINJE decks, AUTOCHK and GASLIFT-13 identical to
the digit, GASLIFT-13's summary byte-identical.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Comment only. The share a group-controlled well is tested against is explicit, for the same reason OPM's guide rates are, and the comment now says that solving for it -- as an unknown of the system, or by iterating an outer loop to a fixed point -- is the way past the cycling rather than something ruled out. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
atgeirr
reviewed
Aug 25, 2026
| // targets are not handled, and re-routing (BRANPROP changing the tree) is not | ||
| // modelled -- the tree is taken as given. | ||
| template<class Scalar> | ||
| class ProductionSystem |
Member
There was a problem hiding this comment.
Please move ProductionSystem and InjectionSystem to separate files. Also consider if it makes sense to use inheritance-based polymorphism with virtual functions instead of taking these classes as template arguments to solve(), I am not convinced it would be better, but would like to explore the option.
Member
Author
There was a problem hiding this comment.
Split in 04a77a3: NetworkSolve.hpp (Node, Result, DenseMatrix, shareByGuide, globalisations, Parameters, solve — 380 lines), NetworkInjectionSystem.hpp and NetworkProductionSystem.hpp. No code changed; your line comments above will have moved.
On virtual functions: I think it would be an improvement and I am happy to do it. The two systems already share ten methods with identical signatures (size, start, residual, jacobian, updateControls, limitStep, columnScale, pressures, wellRates, controlLetter) and State is std::vector in both, so a base is straightforward. What it would buy is that the five one-sided ones stop being duck-typed: solve() currently discovers refreshGuides / setEnforceRateLimits / rateLimitsViolated (injection) and wellPhaseRates / wellBhps (production) with if constexpr (requires ...), so renaming one silently disables it. Dispatch cost is irrelevant at tens of unknowns against VFP lookups. Say the word and I will do it as its own commit.
| /// outside the bench wants. | ||
| template<class Sys> | ||
| Result<typename Sys::ScalarType> | ||
| solve(Sys& system, const std::vector<typename Sys::ScalarType>& node_pressure_guess) |
Member
There was a problem hiding this comment.
I am not happy about this overload. This is functionally very similar to the default parameter approach. I want a single overload with no default arguments, not hiding the default parameters inside a Parameters structure.
Member
Author
There was a problem hiding this comment.
Fair — done properly in 76dcc46. One solve(), no overload, and Parameters has no in-class initialisers, so nothing supplies a value a caller did not ask for. The simulator names its three call sites' settings once as kNetworkSolveParams beside them, and the bench has its own kParams for the cases not varying them; both are ordinary constants in the file that uses them, not defaults.
Drops the two-argument overload and the in-class initialisers on Parameters, so there is one solve() and nothing supplies a value a caller did not ask for. The simulator's three call sites share one named constant, kNetworkSolveParams, defined next to them with the values in plain sight; the bench has its own kParams for the cases that are not varying the settings, and the handful that are pass their own. Neither is a default: both are spelled out at namespace scope in the file that uses them. Tests and all three decks unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
NetworkSystem.hpp was 2527 lines holding both systems and the Newton around
them. Split three ways, no code changed:
NetworkSolve.hpp Node, Result, DenseMatrix, shareByGuide, the
globalisations, Parameters and solve() -- what
neither system owns. 380 lines.
NetworkInjectionSystem.hpp Well, Control, InjectionSystem, dump/replay.
NetworkProductionSystem.hpp ProductionSystem, dump/replay.
solve() is a template over the system, so the solver header knows nothing about
either; both system headers include it. A stale comment above systemUsesAnalytic
that had been left behind by an earlier edit is dropped.
Tests and all three decks unchanged.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Draft, not for merge as one piece — pushed so the work can be seen and tested. Builds on #7345, whose commits are included here rebased on master.
Everything new is opt-in; with default flags the answers are unchanged, except for the injection-network robustness fixes at the bottom of the branch, which are what make
GNETINJEcorrect at all.It splits into roughly: injection-network robustness (VFP axis clamping, bracketing node-pressure update, per-phase leaf rates,
GPRG/GPRW) · the simultaneous solve--network-solver=newton, injection then production ·--network-autochoke·--gas-lift-network-response· two legacy autochoke bugfixes that makeNETWORK_MODEL5_STDW_AUTOCHKrun at all.Needs OPM/opm-common#5314.
Trying it
Numbers are from this branch on an M4 Pro.
--debug-verbosity-level=2puts the counts in the.DBG:grep -c 'inner iterations'= well solves,grep -c 'Network: solved the'= simultaneous solves,grep -c 'is not possible'= fallbacks to the relaxed update.Injection network —
opm-tests/network/GNETINJE_GAS-01.DATACompare
GPRG:M5NandWGIR:G-3H. The relaxed update collapses M5N to 1.1 bar at day 121 and takes both G injectors to zero; the simultaneous solve holds 202.7 bar with all four injectors on.Autochoke —
opm-tests/network/NETWORK_MODEL5_STDW_AUTOCHK.DATAB1 holds its 6000 sm³/d target both ways (
GOPR:B1at day 1/31/60/91: 6019/6001/6016/4846 legacy, 6000/6003/5900/4866 new). They part late, where B-1H is shut at day 106 instead of 152 — open, and legacy's own answer there moves with its bracket sample count. On master this deck aborts at day 3.2 (Invalid Group::ProductionCMode ... B1 cmode: NONE); the two legacy bugfixes in this branch are what make it run at all.Gas lift —
opm-tests/gaslift/GASLIFT-13.DATALeave
--network-analytic-jacobianoff here: on this deck it changes the optimiser's path, not just its cost (302 -> 1151 Newton and a different lift-gas allocation), which is an open question rather than a recommendation.Only the optimiser's trial evaluation changes; its algorithm, increments, weights and limits are untouched. The trial is answered by the network solve instead of at a frozen node pressure, which is why it stops paying for lift gas the network takes straight back.
Production network
No deck in opm-tests exercises this on its own — in
NETWORK-01and its variants every well is a pinned source, so the solve declines the tree and the relaxed evaluation runs unchanged. It is exercised through the two recipes above (236 and 451 simultaneous solves respectively). Directly:tests/test_networksolve.cppcovers the production formulation, the control rule and the group split against the rule-based allocation;--log_level=messageon the binary prints the comparison tables.