对齐PNL0002上游温度并稳定热流体闭合

app/simulation/systems/generic.py

@@ -40,6 +40,10 @@ class SimulationPreparationError(ValueError):
         self.issues = issues
 
 
+class ThermofluidClosureError(RuntimeError):
+    """Raised when stream enthalpy and pressure-flow do not reach one fixed point."""
+
+
 @dataclass(frozen=True)
 class GenericSimulationResult:
     success: bool
@@ -260,6 +264,7 @@ class GenericFluidSystem:
         self.max_algebraic_residual = 0.0
         self.max_algebraic_evaluations = 0
         self.max_stream_iterations = 0
+        self.max_thermofluid_iterations = 0
         self.signal_propagation_count = 0
         self.pneumatic_volume_propagation_count = 0
 
@@ -280,20 +285,62 @@ class GenericFluidSystem:
         for component in self.dynamic_components:
             component.refresh_thermodynamic_ports()
         algebraic = self.pressure_flow_solver.solve()
+        pressure_flow_solve_count = 1
         pneumatic_volume = self.pneumatic_volume_resolver.solve()
         self.pneumatic_volume_propagation_count += pneumatic_volume.propagated
         if pneumatic_volume.propagated:
             for component in self.dynamic_components:
                 component.refresh_thermodynamic_ports()
             algebraic = self.pressure_flow_solver.solve()
-        stream, connected_h = self.stream_resolver.solve()
+            pressure_flow_solve_count += 1
+
         # Some constitutive flow laws recover their upstream temperature from
-        # the connected stream enthalpy.  Stream propagation updates that
-        # cache after the first pressure-flow pass, so refresh explicit flows
-        # once more before evaluating state derivatives and result variables.
-        algebraic = self.pressure_flow_solver.solve()
+        # connected stream enthalpy, while junction stream mixing itself depends
+        # on the resulting mass flows.  A single stream -> pressure-flow refresh
+        # leaves that two-way coupling to the next RHS call, making the ODE RHS
+        # depend on evaluation history and corrupting finite-difference
+        # Jacobians.  Close both layers to one fixed point inside this call.
+        physical_ports = tuple(
+            port
+            for component in self.network.components.values()
+            for definition in component.port_definitions
+            if definition.kind == "physical"
+            for port in (component.get_port(definition.name),)
+        )
+        connected_h: dict[str, dict[str, float]] = {}
+        max_coupling_iterations = 25
+        flow_relative_tolerance = 1.0e-12
+        for coupling_iteration in range(1, max_coupling_iterations + 1):
+            previous_flows = tuple(port.m_flow for port in physical_ports)
+            stream, connected_h = self.stream_resolver.solve()
+            for component in self.dynamic_components:
+                component.update_stream_outflows(connected_h[component.name])
+            algebraic = self.pressure_flow_solver.solve()
+            pressure_flow_solve_count += 1
+            current_flows = tuple(port.m_flow for port in physical_ports)
+            flow_scale = max(
+                [abs(value) for value in (*previous_flows, *current_flows)] + [1.0]
+            )
+            max_flow_delta = max(
+                (
+                    abs(current - previous)
+                    for previous, current in zip(previous_flows, current_flows)
+                ),
+                default=0.0,
+            )
+            if max_flow_delta <= flow_relative_tolerance * flow_scale:
+                break
+        else:
+            raise ThermofluidClosureError(
+                "Stream enthalpy and pressure-flow coupling did not converge "
+                f"after {max_coupling_iterations} iterations."
+            )
+        self.max_thermofluid_iterations = max(
+            self.max_thermofluid_iterations,
+            coupling_iteration,
+        )
         self.mechanical_state_reducer.update_constraint_accelerations()
-        self.algebraic_solve_count += 2 + int(bool(pneumatic_volume.propagated))
+        self.algebraic_solve_count += pressure_flow_solve_count
         self.max_algebraic_residual = max(
             self.max_algebraic_residual,
             algebraic.max_scaled_residual,
@@ -467,6 +514,7 @@ class GenericFluidSystem:
             },
             "stream": {
                 "maxIterationsPerSolve": self.max_stream_iterations,
+                "maxThermofluidIterations": self.max_thermofluid_iterations,
                 "last": (
                     self.stream_resolver.last_diagnostics.as_dict()
                     if self.stream_resolver.last_diagnostics is not None