Implement injection networks

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
@@ -1238,6 +1239,9 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/BlackoilWellModelNetwork_impl.hpp
   opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp
   opm/simulators/wells/BlackoilWellModelNetworkPressureComputation.hpp
+  opm/simulators/wells/NetworkAndersonAcceleration.hpp
+  opm/simulators/wells/NetworkNodePressureUpdater.hpp
+  opm/simulators/wells/NetworkSystem.hpp
   opm/simulators/wells/BlackoilWellModelNldd.hpp
   opm/simulators/wells/BlackoilWellModelNldd_impl.hpp
   opm/simulators/wells/BlackoilWellModelRescoup.hpp

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_pressure_update_acceleration_ = Parameters::Get<Parameters::NetworkPressureUpdateAcceleration>();
+    network_anderson_depth_ = Parameters::Get<Parameters::NetworkAndersonDepth>();
+    network_well_proxy_ = Parameters::Get<Parameters::NetworkWellProxy>();
+    network_well_proxy_max_iterations_ = Parameters::Get<Parameters::NetworkWellProxyMaxIterations>();
+    network_solver_ = Parameters::Get<Parameters::NetworkSolver>();
+    network_analytic_jacobian_ = Parameters::Get<Parameters::NetworkAnalyticJacobian>();
+    network_group_control_ = Parameters::Get<Parameters::NetworkGroupControl>();
+    network_dump_failures_ = Parameters::Get<Parameters::NetworkDumpFailures>();
     local_domains_ordering_ = domainOrderingMeasureFromString(Parameters::Get<Parameters::LocalDomainsOrderingMeasure>());
     write_partitions_ = Parameters::Get<Parameters::DebugEmitCellPartition>();
 
@@ -276,6 +285,32 @@ 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::NetworkPressureUpdateAcceleration>
+        ("Acceleration of the network node-pressure iteration, applied to the whole pressure "
+         "vector of a network instead of the per-node update: none or anderson");
+    Parameters::Register<Parameters::NetworkAndersonDepth>
+        ("Number of past iterates kept by Anderson acceleration of the network pressures");
+    Parameters::Register<Parameters::NetworkWellProxy>
+        ("Balance the injection networks against the wells' inflow-performance linearisation "
+         "(q = A - B*bhp) before re-solving the wells: none or ipr");
+    Parameters::Register<Parameters::NetworkWellProxyMaxIterations>
+        ("Iteration cap for the inflow-performance network balance");
+    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::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 NetworkPressureUpdateAcceleration { static constexpr auto value = "none"; };
+struct NetworkAndersonDepth { static constexpr int value = 4; };
+struct NetworkWellProxy { static constexpr auto value = "none"; };
+struct NetworkSolver { static constexpr auto value = "fixedpoint"; };
+struct NetworkAnalyticJacobian { static constexpr bool value = false; };
+struct NetworkGroupControl { static constexpr bool value = false; };
+struct NetworkDumpFailures { static constexpr auto value = ""; };
+struct NetworkWellProxyMaxIterations { static constexpr int value = 50; };
 // 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,38 @@ 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_;
+
+    /// Acceleration applied to the whole node-pressure vector of a network instead of the
+    /// per-node update: "none" (default) or "anderson"
+    std::string network_pressure_update_acceleration_;
+
+    /// Number of past iterates Anderson acceleration keeps
+    int network_anderson_depth_;
+
+    /// Balance the injection networks against the wells' well-index linearisation
+    /// before re-solving them: "none" (default) or "ipr"
+    std::string network_well_proxy_;
+
+    /// Iteration cap for that inner balance
+    int network_well_proxy_max_iterations_;
+
+    /// 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_;
+
     /// 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;
 }

opm/simulators/wells/BlackoilWellModelNetwork.hpp

@@ -28,6 +28,8 @@
 
 #include <opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp>
 
+#include <optional>
+
 #include <map>
 #include <string>
 #include <tuple>
@@ -66,6 +68,26 @@ class BlackoilWellModelNetwork :
     void doPreStepRebalance(DeferredLogger& deferred_logger);
 
 protected:
+    /// Balance the injection networks against the wells' well-index linearisation
+    /// (q = ipr_b*bhp - ipr_a, from the converged well Jacobian) instead of re-solving the
+    /// well equations, so a residual evaluation costs a handful of VFP lookups. Returns the
+    /// last imbalance. Experimental, off unless --network-well-proxy=ipr.
+    Scalar proxyBalance(const int episodeIdx,
+                        const double dt,
+                        const int max_iterations,
+                        const Scalar damping_factor,
+                        const Scalar max_pressure_update,
+                        const bool use_secant,
+                        const bool secant_production,
+                        DeferredLogger& deferred_logger);
+
+    /// Rate this injector would take at its current THP constraint, from the well-index
+    /// linearisation alone. nullopt when it admits no solution there.
+    std::optional<Scalar>
+    proxyInjectionRate(WellInterface<TypeTag>& well,
+                       const int phase_pos,
+                       DeferredLogger& deferred_logger) const;
+
     /// This function is to be used for well groups in an extended network that act as a subsea manifold
     /// The wells of such group should have a common THP and total phase rate(s) obeying (if possible)
     /// the well group constraint set by GCONPROD