[WIP] Simultaneous network solve: injection and production networks, autochoke, gas-lift response

[hnil]

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 GNETINJE correct 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 make NETWORK_MODEL5_STDW_AUTOCHK run at all.

Needs OPM/opm-common#5314.

Trying it

Numbers are from this branch on an M4 Pro. --debug-verbosity-level=2 puts 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.DATA

flow $DECK --output-dir=/tmp/a
flow $DECK --output-dir=/tmp/b --network-solver=newton --network-analytic-jacobian=true

Compare GPRG:M5N and WGIR: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.DATA

flow $DECK --output-dir=/tmp/c
flow $DECK --output-dir=/tmp/d --network-solver=newton --network-autochoke=true --network-analytic-jacobian=true

	well solves	Newton	wall
legacy	475 143	790	15.0 s
simultaneous	3 774	123	2.1 s

B1 holds its 6000 sm³/d target both ways (GOPR:B1 at 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.DATA

flow $DECK --output-dir=/tmp/e
flow $DECK --output-dir=/tmp/f --network-solver=newton --gas-lift-network-response=true

	FOPT (sm³)	lift gas (sm³)	well solves	Newton	wall
legacy	1 029 703	1 120 162	72 919	2 733	53.7 s
simultaneous	1 049 109	85 430	7 456	302	6 s

Leave --network-analytic-jacobian off 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-01 and 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:

ctest -R networksolve

tests/test_networksolve.cpp covers the production formulation, the control rule and the group split against the rule-based allocation; --log_level=message on the binary prints the comparison tables.

[atgeirr]

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.

[atgeirr]

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.

[hnil]

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.

[atgeirr]

Why is a different convention used than in the existing code? Would prefer that conventions are the same unless for a good reason.

[hnil]

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.

[atgeirr]

This class reimplements functionality already available in Dune, prefer to use existing library functionality if equivalent.

[hnil]

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.

[atgeirr]

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.

[hnil]

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.

[atgeirr]

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.

[hnil]

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.

[atgeirr]

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.

[hnil]

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.

[hnil]

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.

[atgeirr]

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.

[hnil]

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

[atgeirr]

Please document this class, and rename it InjectionSystem (assuming that I am correct that this is the counterpart to ProductionSystem).

[hnil]

Done in d37d21e — InjectionSystem, with a class comment saying what the unknowns and equations are and that it is the counterpart to ProductionSystem.

[atgeirr]

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.

[hnil]

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.

[atgeirr]

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.

[hnil]

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.

[hnil]

Went through the rest of the network code the same way. Fourteen things worth doing, roughly: shareByGuide duplicates the recursion in #7220; the parallel gather is a hand-packed flat array with a manual stride; Node::efficiency and LineSearch use double in Scalar code; five OPM_NETWORK_* env vars switch behaviour next to a proper --network-dump-failures; the two Control enums disagree; Result carries production-only fields and guides_moving is dead; the rate/pressure scaling is duplicated in both systems; 1e6 is a dimensional sentinel for "no rate limit"; the solve memoises on float equality; a singular Jacobian is logged as ordinary non-convergence; and the injection system has no Shut. Happy to take them in any order you prefer — I would start with the last two and the sentinel.

Also swept for default arguments rather than only fixing the one you flagged: updatePressures(), thpPotential() and addNode() each had one this branch had added, all gone in fed7d15. update()'s relax_network_tolerance is left alone since it is upstream's.

[atgeirr]

This member, and its use in e.g. hasLeafNodeRates() does nothing. It should be removed, and also from functions such as computePressures() in BlackoilWellModelNetworkGeneric that only needed it to pass to this constructor.

[hnil]

Right, it does nothing — done in 5cda052. Once the leaf rates stopped being keyed by phase it was only handed to hasLeafNodeRate()/leafNodeRate(), which ignore it, since the VFPINJ table picks the phase its FLO type names. Member, constructor argument, the computePressures() parameter and five test call sites all gone; the injection overload now has the same signature as the production one. The local in the caller stays because newtonNodePressures() still needs it.

[atgeirr]

Can this function be replaced by, or rewritten to use the existing RegulaFalsi implementation in opm-common's RootFinders.hpp?

If necessary, suggest improvements for that code to make it usable here.

[atgeirr]

It is not well explained what a per-node update means here.

[hnil]

Rewritten in a9b77c3: one of these per node, held by the caller in a map and stepped once per node per network sub-iteration, each with its own bracket — which is also the limitation, since the nodes are coupled through the wells below them and a bracket end has to be dropped when a sibling moving contradicts it. Solving every node at once instead is what NetworkSolve::solve() does.

[atgeirr]

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

I think the arguments above are very much in favor of switching to inheritance and virtual function-based polymorphism for the NetworkInjectionSystem and NetworkProductionSystem classes and the solve() function that use them. Please also make that change.

[hnil]

Base class done in 6a4df6c — SystemBase<Scalar>, both systems derive, solve() takes SystemBase<Scalar>&. All five if constexpr (requires ...) are gone, and with them systemUsesAnalytic()/systemJacobian(), which existed only to ask a type whether it knew how. Twelve pure virtuals for what both implement, five virtual-with-default for the one-sided ones. It compiled first try with every override accepted, which is the check the old version could not make. Tests and all three decks unchanged.

Also in this push: the dead injection phase (5cda052) and the per-node comment (a9b77c3).

On RootFinders.hpp I would push back, and would like your view. RegulaFalsi::solve(f, a, b, ...) owns the iteration and calls f itself, and it wants a bracket up front. NodePressureUpdater cannot give it either: one evaluation of r(P) is a well-model and network pass over the whole tree on every rank, driven by the outer loop, so the class is a stepper — called once per node per sub-iteration, carrying its bracket between calls, building it as evidence arrives, dropping an end when the wells change state under it, and handling the no-VFP-solution and plateau cases. Wrapping opm-common's version would mean inverting that control flow. What would make it usable is splitting it: a state-holding stepper with a step(x, f(x)) -> next x, and solve() as a thin loop over that. I am happy to do that in opm-common if you think it is worth it — it is a bigger change than it looks, and it touches everything that uses RegulaFalsi today.

[atgeirr]

I think the changes to this file are ready to be merged to the master branch. Please extract the changes here, as well as follow-on changes in test_networkpressure.cpp and other files that adapt to the new return type NetworkBranchPressure and the invalid node concept to its own separate pull request.

[hnil]

Happy to — but the caller side does not stand alone on master, and I would rather ask than guess. In master computePressures() returns a pair and there is no network-domain concept at all; the invalid-node handling in this PR lives in the multi-domain BlackoilWellModelNetworkGeneric, which is your own base commits, not master. It also needs the two has_network_leaf_node_*_rates() accessors added to GroupState.

So either the new PR carries your nine base commits with it (which duplicates the first nine of this one), or I port the invalid-node handling onto master single-domain caller, which is a rewrite of that part rather than an extraction. Which would you prefer? The header itself and test_networkpressure.cpp are straightforward either way.

[atgeirr]

So either the new PR carries your nine base commits with it (which duplicates the first nine of this one), or I port the invalid-node handling onto master single-domain caller, which is a rewrite of that part rather than an extraction. Which would you prefer? The header itself and test_networkpressure.cpp are straightforward either way.

I see, I think then we should try to start from the nine base commits and include them. That means the first PR is a bit larger than ideal, but it should still be a consistent and understandable step forward. If there are parts of those commits that have been made obsolete or otherwise rewritten completely by later commits here, leave out those parts unless it makes the PR inconsistent or would lead to existing functionality (production network) failing where it succeeds in master.

[hnil]

Done — #7368. It came out cleaner than expected: your nine commits plus the VFP one are already a contiguous prefix of this branch, so it is a straight cut rather than a port, with NetworkBranchPressure, the invalid-node concept and the six test cases.

I checked your condition on production explicitly, by building master and diffing the runs rather than trusting the test verdicts: NETWORK-01 and NETWORK-01_STANDARD agree with master on all 369 summary vectors except TCPU. Worth knowing that the seven compareECLFiles_flow+NETWORK-01* tests fail for me on plain master too — my opm-tests checkout predates GWGR/FWGR being added to the summary — so they are not a usable gate here.

Nothing from the base nine looked obsolete enough to drop; the only thing later work undid was our own per-phase leaf-rate keying, which is already back to your node-keyed map.

[atgeirr]

Rename this function newtonInjectionNodePressures() to parallel the production version.

[hnil]

Renamed in 7a8b7fa.

[atgeirr]

The branch data will be stale if the newton solver was used: it is unchanged from the initial calls to computePressures(). Please update the data from the converged solution.

[hnil]

Good catch — fixed in 7a8b7fa. GPRB against GPR was out by 0.0037 bar on NETWORK_MODEL5_STDW_AUTOCHK at day 91; now consistent.

Two things I found doing it. The refresh has to go after the relaxed update rather than straight after the solve: that update moves the stored pressures again (it bounds the wells step towards what the solve placed), so refreshing any earlier just picks a different stale value — my first attempt did exactly that and was still 0.0037 bar out.

And it is not only the newton path. The default path has the same drift for the same reason — the damped update moves the pressure after computePressures() took the branch data — 0.00045 bar on that deck and step. I have left those nodes alone, since correcting them changes GPRB on every network deck; say if you would rather have it fixed everywhere and I will do it as its own commit. Default output here is byte-identical.

[hnil]

Picked up the other half of your first review — the split landed but not the network/ subdir. NetworkSolve.hpp, NetworkInjectionSystem.hpp, NetworkProductionSystem.hpp and NetworkNodePressureUpdater.hpp are now in opm/simulators/wells/network/, following wells/rescoup/: prefixes kept, listed after the flat wells headers. Paths only (a2f4989).

I took the node pressure updater along even though it belongs to the relaxed path rather than the simultaneous one — leaving the one remaining network header outside a network/ directory seemed worse than that small inconsistency. Easy to move back if you disagree.

Also: I had not been running the ctest regressions until now, only the unit tests — flow was not being built here. Both GNETINJE cases pass, as do the two network_cost tests and the parallel NETWORK-01 case. The seven NETWORK-01* comparisons fail in my checkout for an unrelated reason: my opm-tests is on a branch from 11 May, before GWGR/FWGR were added to the summary on 27 May, so the run emits keywords the references do not have. Nothing to do with this PR, but worth saying out loud since I had been quoting my own hand comparisons up to now.

[hnil]

Carried the three Copilot fixes from #7368 onto this branch too (71efe14), so the two do not drift apart: details::domainForWell() at all three call sites, the per-node pressure/inflow trace behind OPM_NETWORK_PRESSURE_TRACE, and the && indentation. Tests, both GNETINJE regressions and both cost tests pass.