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

CMakeLists_files.cmake

@@ -491,6 +491,7 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_milu.cpp
   tests/test_multmatrixtransposed.cpp
   tests/test_networkpressure.cpp
+  tests/test_networksolve.cpp
   tests/test_nonnc.cpp
   tests/test_norne_pvt.cpp
   tests/test_OilSatfuncConsistencyChecks.cpp
@@ -1315,6 +1316,10 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/WellTest.hpp
   opm/simulators/wells/WellTracerRate.hpp
   opm/simulators/wells/WGState.hpp
+  opm/simulators/wells/network/NetworkInjectionSystem.hpp
+  opm/simulators/wells/network/NetworkNodePressureUpdater.hpp
+  opm/simulators/wells/network/NetworkProductionSystem.hpp
+  opm/simulators/wells/network/NetworkSolve.hpp
   opm/simulators/wells/rescoup/RescoupProxy.hpp
 )
 if (USE_GPU_BRIDGE)

compareECLFiles.cmake

@@ -755,6 +755,7 @@ endif ()
 include (${CMAKE_CURRENT_SOURCE_DIR}/regressionTests.cmake)
 include (${CMAKE_CURRENT_SOURCE_DIR}/comparisonTests.cmake)
 include (${CMAKE_CURRENT_SOURCE_DIR}/restartTests.cmake)
+include (${CMAKE_CURRENT_SOURCE_DIR}/networkCostTests.cmake)
 
 # PORV test
 opm_set_test_driver(${PROJECT_SOURCE_DIR}/tests/run-porv-acceptanceTest.sh "")

networkCostTests.cmake

@@ -0,0 +1,131 @@
+# What the network solve costs, as opposed to what it computes.
+#
+# The regression tests compare the answer and would pass a run that reaches it
+# by falling back to the relaxed update on every step, or by re-solving one
+# well thousands of times because the network and the well model disagree
+# about whether it flows. Both have happened on these two decks. The bounds
+# are roughly 25 % above what the simultaneous solve does today; they are
+# there to catch a regression, not to pin an exact count.
+
+function(add_network_cost_test)
+  set(oneValueArgs CASENAME FILENAME DIR MAX_FALLBACKS MAX_NEWTON MAX_WELL_SOLVES MAX_WELL_SHARE)
+  set(multiValueArgs TEST_ARGS)
+  cmake_parse_arguments(PARAM "" "${oneValueArgs}" "${multiValueArgs}" ${ARGN})
+  add_test(NAME network_cost_${PARAM_CASENAME}
+           COMMAND ${PROJECT_SOURCE_DIR}/tests/run-network-cost-test.sh
+                   -i ${OPM_TESTS_ROOT}/${PARAM_DIR}
+                   -f ${PARAM_FILENAME}
+                   -r ${PROJECT_BINARY_DIR}/tests/results/network_cost_${PARAM_CASENAME}
+                   -e ${PROJECT_BINARY_DIR}/bin/flow_blackoil
+                   -F ${PARAM_MAX_FALLBACKS}
+                   -N ${PARAM_MAX_NEWTON}
+                   -W ${PARAM_MAX_WELL_SOLVES}
+                   -S ${PARAM_MAX_WELL_SHARE}
+                   -- ${PARAM_TEST_ARGS})
+  set_tests_properties(network_cost_${PARAM_CASENAME} PROPERTIES LABELS "network_cost")
+endfunction()
+
+# An autochoke under the simultaneous solve. Measured 2026-08-24: 0 fallbacks,
+# 199 Newton, 4603 well solves, B-1H 45 % of them. Legacy on the same deck
+# needs 746 Newton and 70040 well solves for the same answer.
+add_network_cost_test(
+  CASENAME autochoke_complementarity
+  FILENAME NETWORK_MODEL5_STDW_AUTOCHK
+  DIR network
+  MAX_FALLBACKS 0
+  MAX_NEWTON 260
+  MAX_WELL_SOLVES 6000
+  MAX_WELL_SHARE 60
+  TEST_ARGS --network-solver=newton --network-analytic-jacobian=true
+            --network-group-control=true --network-autochoke=true
+            --network-complementarity=true
+)
+
+# Gas lift answered by the network. Measured 2026-08-24: 0 fallbacks, 521
+# Newton, 9051 well solves, B-3H 32 % of them. Legacy: 2607 and 49320.
+add_network_cost_test(
+  CASENAME gaslift_complementarity
+  FILENAME GASLIFT-13
+  DIR gaslift
+  MAX_FALLBACKS 0
+  MAX_NEWTON 680
+  MAX_WELL_SOLVES 11500
+  MAX_WELL_SHARE 45
+  TEST_ARGS --network-solver=newton --network-analytic-jacobian=true
+            --network-group-control=true --gas-lift-network-response=true
+            --network-complementarity=true
+)
+
+# The simultaneous solve reaches the same balance as the relaxed update by a
+# different route, so on a deck where the relaxed update converges the two must
+# agree. Checkable without any external reference, which matters because these
+# decks have none.
+#
+# The tolerance is looser than the regression suite's (which compares a run
+# against its own stored reference): two solution paths through the same
+# schedule differ at the sub-percent level. It still discriminates -- where the
+# formulations genuinely disagree the gap is 4-6 % (the active set against
+# legacy on autochoke, and the gas-lift decks, which is why those are covered
+# by cost tests instead).
+set(network_agreement_abs_tol 2e-2)
+set(network_agreement_rel_tol 2e-2)
+
+function(add_network_agreement_test)
+  set(oneValueArgs CASENAME FILENAME DIR)
+  set(multiValueArgs BOTH_ARGS NEWTON_ARGS)
+  cmake_parse_arguments(PARAM "" "${oneValueArgs}" "${multiValueArgs}" ${ARGN})
+  string(REPLACE ";" " " _both "${PARAM_BOTH_ARGS}")
+  add_test(NAME network_agreement_${PARAM_CASENAME}
+           COMMAND ${PROJECT_SOURCE_DIR}/tests/run-network-agreement-test.sh
+                   -i ${OPM_TESTS_ROOT}/${PARAM_DIR}
+                   -f ${PARAM_FILENAME}
+                   -r ${PROJECT_BINARY_DIR}/tests/results/network_agreement_${PARAM_CASENAME}
+                   -e ${PROJECT_BINARY_DIR}/bin/flow_blackoil
+                   -c $<TARGET_FILE:compareECL>
+                   -a ${network_agreement_abs_tol}
+                   -t ${network_agreement_rel_tol}
+                   -b "${_both}"
+                   -- ${PARAM_NEWTON_ARGS})
+  set_tests_properties(network_agreement_${PARAM_CASENAME}
+                       PROPERTIES LABELS "network_agreement")
+endfunction()
+
+# Every deck where the production solve engages and the relaxed update also
+# converges. Measured 2026-08-24 (legacy FOPT vs simultaneous, well solves):
+#   NETWORK-01              538142 / 449   vs 538842 / 379
+#   NETWORK-01_STANDARD     538142 / 449   vs 538842 / 379
+#   NETWORK-01-REROUTE      543737 / 756   vs 544171 / 493
+#   NETWORK-01-WTEST        674921 / 1016  vs 675875 / 549
+#   NETWORK-01-WEFAC-...    128788 / 430   vs 128772 / 365
+#   NETWORK-01-MULTIROOT    778152 / 407   vs 778188 / 380
+#   6_UDA_MODEL5_STDW       317225 / 1781  vs 317225 / 1788  (exact)
+set(_agree_newton --network-solver=newton --network-analytic-jacobian=true)
+
+foreach(case NETWORK-01 NETWORK-01_STANDARD NETWORK-01-WTEST
+             NETWORK-01-WEFAC-GEFAC-ITEM3NO NETWORK-01-MULTIROOT)
+  string(TOLOWER ${case} _lc)
+  string(REPLACE "-" "_" _lc ${_lc})
+  add_network_agreement_test(
+    CASENAME ${_lc}
+    FILENAME ${case}
+    DIR network
+    BOTH_ARGS --enable-tuning=true
+    NEWTON_ARGS ${_agree_newton}
+  )
+endforeach()
+
+add_network_agreement_test(
+  CASENAME network_01_reroute
+  FILENAME NETWORK-01-REROUTE
+  DIR network
+  BOTH_ARGS --enable-tuning=true --local-well-solve-control-switching=true
+  NEWTON_ARGS ${_agree_newton}
+)
+
+add_network_agreement_test(
+  CASENAME uda_model5_stdw
+  FILENAME 6_UDA_MODEL5_STDW
+  DIR model5
+  BOTH_ARGS --enable-tuning=true
+  NEWTON_ARGS ${_agree_newton} --network-group-control=true
+)

opm/simulators/flow/BlackoilModelParameters.cpp

@@ -116,6 +116,15 @@ BlackoilModelParameters<Scalar>::BlackoilModelParameters()
     rc_network_loose_coupling_ = Parameters::Get<Parameters::RcNetworkLooseCoupling>();
     network_pressure_update_damping_factor_ = Parameters::Get<Parameters::NetworkPressureUpdateDampingFactor<Scalar>>();
     network_max_pressure_update_in_bars_ = Parameters::Get<Parameters::NetworkMaxPressureUpdateInBars<Scalar>>();
+    network_pressure_update_secant_ = Parameters::Get<Parameters::NetworkPressureUpdateSecant>();
+    network_solver_ = Parameters::Get<Parameters::NetworkSolver>();
+    network_analytic_jacobian_ = Parameters::Get<Parameters::NetworkAnalyticJacobian>();
+    network_group_control_ = Parameters::Get<Parameters::NetworkGroupControl>();
+    network_autochoke_ = Parameters::Get<Parameters::NetworkAutochoke>();
+    network_autochoke_bracket_samples_ = Parameters::Get<Parameters::NetworkAutochokeBracketSamples>();
+    network_complementarity_ = Parameters::Get<Parameters::NetworkComplementarity>();
+    gaslift_network_response_ = Parameters::Get<Parameters::GasLiftNetworkResponse>();
+    network_dump_failures_ = Parameters::Get<Parameters::NetworkDumpFailures>();
     local_domains_ordering_ = domainOrderingMeasureFromString(Parameters::Get<Parameters::LocalDomainsOrderingMeasure>());
     write_partitions_ = Parameters::Get<Parameters::DebugEmitCellPartition>();
 
@@ -276,6 +285,39 @@ void BlackoilModelParameters<Scalar>::registerParameters()
         ("Damping factor in the inner network pressure update iterations");
     Parameters::Register<Parameters::NetworkMaxPressureUpdateInBars<Scalar>>
         ("Maximum pressure update in the inner network pressure update iterations");
+    Parameters::Register<Parameters::NetworkSolver>
+        ("How the injection networks are solved: fixedpoint relaxes the node pressures against "
+         "the wells, newton solves pressures and rates simultaneously and falls back to the "
+         "fixed point when it does not converge");
+    Parameters::Register<Parameters::NetworkAnalyticJacobian>
+        ("Assemble the network Jacobian from the VFP table derivatives instead of differencing "
+         "the residual (--network-solver=newton only)");
+    Parameters::Register<Parameters::NetworkGroupControl>
+        ("Let the network hold a group's injection total and place the split itself, so a well "
+         "that hits its own limit is taken up by the others (--network-solver=newton only)");
+    Parameters::Register<Parameters::NetworkAutochoke>
+        ("Solve autochoke nodes inside the simultaneous network solve: the node pressure "
+         "becomes the group's common thp and is raised until the oil through the node meets "
+         "the group's target, instead of the bracketing search over well solves "
+         "(--network-solver=newton only).");
+    Parameters::Register<Parameters::NetworkAutochokeBracketSamples>
+        ("Samples the legacy autochoke search takes across its bracket before the root find; each "
+         "sample solves every well in the group. 300 is the historical value; a dozen finds the "
+         "same root for a fraction of the well solves.");
+    Parameters::Register<Parameters::NetworkComplementarity>
+        ("Close each production well's own limits -- rate, tubing, bhp -- with one complementarity "
+         "row in the simultaneous network solve instead of an active set, so nothing switches "
+         "(--network-solver=newton with --network-analytic-jacobian=true).");
+    Parameters::Register<Parameters::GasLiftNetworkResponse>
+        ("Answer the gas lift optimiser's trial evaluations from the simultaneous network "
+         "solve -- the well's rates with every node pressure responding to its lift gas -- "
+         "instead of a well solve at a fixed thp (--network-solver=newton only).");
+    Parameters::Register<Parameters::NetworkDumpFailures>
+        ("Path prefix for writing out each network system that fails to converge, for replay in "
+         "the standalone bench; empty disables it");
+    Parameters::Register<Parameters::NetworkPressureUpdateSecant>
+        ("Networks whose node pressures use the bracketing/secant update in the inner network "
+         "iterations instead of the damped update: injection, all or none");
     Parameters::Register<Parameters::NonlinearSolver>
         ("Choose nonlinear solver. Valid choices are newton or nldd.");
     Parameters::Register<Parameters::LocalSolveApproach>

opm/simulators/flow/BlackoilModelParameters.hpp

@@ -158,6 +158,15 @@ template<class Scalar>
 struct NetworkPressureUpdateDampingFactor { static constexpr Scalar value = 0.1; };
 template<class Scalar>
 struct NetworkMaxPressureUpdateInBars { static constexpr Scalar value = 5.0; };
+struct NetworkPressureUpdateSecant { static constexpr auto value = "injection"; };
+struct NetworkSolver { static constexpr auto value = "fixedpoint"; };
+struct NetworkAnalyticJacobian { static constexpr bool value = false; };
+struct NetworkGroupControl { static constexpr bool value = false; };
+struct NetworkAutochoke { static constexpr bool value = false; };
+struct NetworkAutochokeBracketSamples { static constexpr int value = 300; };
+struct NetworkComplementarity { static constexpr bool value = false; };
+struct GasLiftNetworkResponse { static constexpr bool value = false; };
+struct NetworkDumpFailures { static constexpr auto value = ""; };
 // Reservoir coupling: when false (default) the master exchanges node pressures
 // and slave rates with the slaves once per master inner network sub-iteration
 // (tight coupling).  When true, the exchange happens only once per master outer network
@@ -358,6 +367,28 @@ struct BlackoilModelParameters
     /// Maximum pressure update in the inner network pressure update iterations
     Scalar network_max_pressure_update_in_bars_;
 
+    /// Which networks use the bracketing/secant node-pressure update in the inner network
+    /// iterations instead of the damped one: "injection" (default), "all" or "none"
+    std::string network_pressure_update_secant_;
+
+    /// How the injection networks are solved: "fixedpoint" (default) relaxes the
+    /// node pressures against the wells; "newton" solves pressures and rates
+    /// simultaneously, falling back to the fixed point when it does not converge.
+    std::string network_solver_;
+
+    /// Assemble the network Jacobian from the VFP table derivatives.
+    bool network_analytic_jacobian_;
+
+    /// Let the network place the split of a group's injection total itself.
+    bool network_group_control_;
+
+    /// Path prefix for writing network systems that fail to converge; empty off.
+    std::string network_dump_failures_;
+    bool network_autochoke_ = false;
+    int network_autochoke_bracket_samples_ = 300;
+    bool network_complementarity_ = false;
+    bool gaslift_network_response_ = false;
+
     /// Reservoir coupling: use loose (per-outer-iteration) master/slave network
     /// coupling instead of the default tight (per-sub-iteration) coupling.
     bool rc_network_loose_coupling_;

opm/simulators/utils/UnsupportedFlowKeywords.cpp

@@ -228,7 +228,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"GINODE", {true, std::nullopt}},
         {"GLIFTLIM", {true, std::nullopt}},
         {"GNETDP", {true, std::nullopt}},
-        {"GNETINJE", {true, std::nullopt}},
         {"GNETPUMP", {true, std::nullopt}},
         {"GRADGRUP", {true, std::nullopt}},
         {"GRADRESV", {true, std::nullopt}},

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -1324,7 +1324,7 @@ groupAndNetworkData(const int reportStepIdx) const
     auto grp_nwrk_values = data::GroupAndNetworkValues{};
 
     this->assignGroupValues(reportStepIdx, grp_nwrk_values.groupData);
-    this->genNetwork_.assignNodeAndBranchValues(grp_nwrk_values.nodeData, grp_nwrk_values.branchData, grp_nwrk_values.convergedBranchData, reportStepIdx - 1); // Schedule state info at previous step
+    this->genNetwork_.assignNodeAndBranchValues(grp_nwrk_values, reportStepIdx - 1); // Schedule state info at previous step
 
     return grp_nwrk_values;
 }