修复test_mql事件后积分步长控制

PythonModels/components/amesim_pneumatic.py

@@ -82,12 +82,24 @@ class AmesimPneumaticVolume(DynamicComponent):
         gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS,
         p0: float = 101_325.0,
         T0: float = 293.15,
+        heat_transfer_coefficient: float = 0.0,
+        heat_transfer_area: float = 0.0,
+        external_temperature_k: float = 293.15,
     ) -> None:
         if volume <= 0.0:
             raise ValueError("volume must be positive.")
+        if heat_transfer_coefficient < 0.0:
+            raise ValueError("heat_transfer_coefficient must be non-negative.")
+        if heat_transfer_area < 0.0:
+            raise ValueError("heat_transfer_area must be non-negative.")
+        if external_temperature_k <= 0.0:
+            raise ValueError("external_temperature_k must be positive.")
         super().__init__(name=name)
         self.volume = volume
         self.gas = gas
+        self.heat_transfer_coefficient = heat_transfer_coefficient
+        self.heat_transfer_area = heat_transfer_area
+        self.external_temperature = external_temperature_k
         rho0 = gas.density(p0, T0)
         m0 = rho0 * volume
         U0 = m0 * gas.specific_internal_energy(T0)
@@ -103,8 +115,20 @@ class AmesimPneumaticVolume(DynamicComponent):
         gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS,
         p0: float = 101_325.0,
         T0: float = 293.15,
+        heat_transfer_coefficient: float = 0.0,
+        heat_transfer_area: float = 0.0,
+        external_temperature_k: float = 293.15,
     ) -> "AmesimPneumaticVolume":
-        return cls(name=name, volume=liters_to_m3(volume_liters), gas=gas, p0=p0, T0=T0)
+        return cls(
+            name=name,
+            volume=liters_to_m3(volume_liters),
+            gas=gas,
+            p0=p0,
+            T0=T0,
+            heat_transfer_coefficient=heat_transfer_coefficient,
+            heat_transfer_area=heat_transfer_area,
+            external_temperature_k=external_temperature_k,
+        )
 
     def get_state_vector(self) -> list[float]:
         return self.state.as_vector()
@@ -118,6 +142,14 @@ class AmesimPneumaticVolume(DynamicComponent):
     def volume_rate_m3_s(self) -> float:
         return 0.0
 
+    def thermal_energy_flow_w(self, temperature_k: float | None = None) -> float:
+        temperature = self.properties().T if temperature_k is None else temperature_k
+        return (
+            self.heat_transfer_coefficient
+            * self.heat_transfer_area
+            * (self.external_temperature - temperature)
+        )
+
     def gas_mass_g(self) -> float:
         return kg_to_g(self.state.m)
 
@@ -139,7 +171,10 @@ class AmesimPneumaticVolume(DynamicComponent):
         return ThermodynamicProperties(p=p, T=T, rho=rho, u=u, h=h)
 
     def derivatives(self, inlet_h: float, m_flow: float) -> VolumeState:
-        return VolumeState(m=m_flow, U=m_flow * inlet_h)
+        return VolumeState(
+            m=m_flow,
+            U=m_flow * inlet_h + self.thermal_energy_flow_w(),
+        )
 
     def derivatives_from_two_connections(
         self,
@@ -151,6 +186,7 @@ class AmesimPneumaticVolume(DynamicComponent):
         internal_h: float,
         volume_rate_m3_s: float | None = None,
     ) -> VolumeState:
+        properties = self.properties()
         inlet_h_a = self.connection_inlet_enthalpy(
             port_m_flow=port_a_m_flow,
             connected_h=connected_h_a,
@@ -166,7 +202,8 @@ class AmesimPneumaticVolume(DynamicComponent):
             U=(
                 port_a_m_flow * inlet_h_a
                 + port_b_m_flow * inlet_h_b
-                - self.properties().p * (
+                + self.thermal_energy_flow_w(properties.T)
+                - properties.p * (
                     self.volume_rate_m3_s()
                     if volume_rate_m3_s is None
                     else volume_rate_m3_s
@@ -186,6 +223,9 @@ class AmesimVariablePneumaticVolume(AmesimPneumaticVolume):
         p0: float = 101_325.0,
         T0: float = 293.15,
         external_volume: float = 0.0,
+        heat_transfer_coefficient: float = 0.0,
+        heat_transfer_area: float = 0.0,
+        external_temperature_k: float = 293.15,
     ) -> None:
         if dead_volume <= 0.0:
             raise ValueError("dead_volume must be positive.")
@@ -200,6 +240,9 @@ class AmesimVariablePneumaticVolume(AmesimPneumaticVolume):
             gas=gas,
             p0=p0,
             T0=T0,
+            heat_transfer_coefficient=heat_transfer_coefficient,
+            heat_transfer_area=heat_transfer_area,
+            external_temperature_k=external_temperature_k,
         )
 
     @classmethod
@@ -211,6 +254,9 @@ class AmesimVariablePneumaticVolume(AmesimPneumaticVolume):
         p0: float = 101_325.0,
         T0: float = 293.15,
         external_volume_liters: float = 0.0,
+        heat_transfer_coefficient: float = 0.0,
+        heat_transfer_area: float = 0.0,
+        external_temperature_k: float = 293.15,
     ) -> "AmesimVariablePneumaticVolume":
         return cls(
             name=name,
@@ -219,6 +265,9 @@ class AmesimVariablePneumaticVolume(AmesimPneumaticVolume):
             p0=p0,
             T0=T0,
             external_volume=liters_to_m3(external_volume_liters),
+            heat_transfer_coefficient=heat_transfer_coefficient,
+            heat_transfer_area=heat_transfer_area,
+            external_temperature_k=external_temperature_k,
         )
 
     def volume_rate_m3_s(self) -> float:

PythonModels/components/amesim_pneumatic_line.py

@@ -48,7 +48,7 @@ class _DarcyPipeResistanceMixin:
             if upper > 1.0e3:
                 raise ValueError("unable to bracket pneumatic pipe resistance flow")
         lower = 0.0
-        for _ in range(80):
+        for _ in range(48):
             middle = 0.5 * (lower + upper)
             if self._darcy_pressure_drop(
                 middle,
@@ -255,15 +255,19 @@ class AmesimPnl0001Pipe(_DarcyPipeResistanceMixin, DynamicComponent):
         connected_h_2: float,
     ) -> VolumeState:
         internal = self.properties()
-        inlet_h_1 = self.connection_inlet_enthalpy(
+        # PNL0001 is a fixed-volume distributed line store. Its transported
-            port_m_flow=port_1_m_flow,
+        # energy variable therefore follows specific internal energy, not the
-            connected_h=connected_h_1,
+        # chamber-style stagnation enthalpy contract. For this ideal gas,
-            internal_h=internal.h,
+        # h = gamma * u. Outflow always carries the local u.
+        inlet_u_1 = (
+            connected_h_1 / self.gas.gamma
+            if port_1_m_flow > 0.0
+            else internal.u
         )
-        inlet_h_2 = self.connection_inlet_enthalpy(
+        inlet_u_2 = (
-            port_m_flow=port_2_m_flow,
+            connected_h_2 / self.gas.gamma
-            connected_h=connected_h_2,
+            if port_2_m_flow > 0.0
-            internal_h=internal.h,
+            else internal.u
         )
         heat_flow = (
             self.heat_transfer_coefficient
@@ -272,7 +276,7 @@ class AmesimPnl0001Pipe(_DarcyPipeResistanceMixin, DynamicComponent):
         )
         return VolumeState(
             m=port_1_m_flow + port_2_m_flow,
-            U=port_1_m_flow * inlet_h_1 + port_2_m_flow * inlet_h_2 + heat_flow,
+            U=port_1_m_flow * inlet_u_1 + port_2_m_flow * inlet_u_2 + heat_flow,
         )
 
 

PythonModels/core/solver.py

@@ -10,8 +10,9 @@ class SolveIVPConfig:
     t_stop: float = 20.0
     method: str = "BDF"
     rtol: float = 1e-6
-    atol: float = 1e-8
+    atol: float = 1e-10
     max_step: float = 1e-3
+    first_step: float | None = None
 
 
 @dataclass(frozen=True)
@@ -91,12 +92,16 @@ def integrate_ode(
     except ImportError:
         return _runge_kutta_4(rhs, initial_state, config, t_eval)
 
-    return solve_ivp(
+    solve_options = {
-        fun=rhs,
+        "fun": rhs,
-        t_span=(config.t_start, config.t_stop),
+        "t_span": (config.t_start, config.t_stop),
-        y0=initial_state,
+        "y0": initial_state,
-        method=config.method,
+        "method": config.method,
-        rtol=config.rtol,
+        "rtol": config.rtol,
-        atol=config.atol,
+        "atol": config.atol,
-        t_eval=t_eval,
+        "max_step": config.max_step,
-    )
+        "t_eval": t_eval,
+    }
+    if config.first_step is not None:
+        solve_options["first_step"] = config.first_step
+    return solve_ivp(**solve_options)

PythonModels/scripts/run_test_mql_full_state_comparison.py

@@ -559,6 +559,8 @@ def format_test_mql_full_state_comparison_summary(
                 f"  - mass_derivative_kg_s={chamber_diagnostic.mass_derivative_kg_s}",
                 f"  - port_a_energy_flow_w={chamber_diagnostic.port_a_energy_flow_w}",
                 f"  - boundary_work_w={chamber_diagnostic.boundary_work_w}",
+                f"  - thermal_energy_flow_w="
+                f"{chamber_diagnostic.thermal_energy_flow_w}",
                 f"  - energy_derivative_w={chamber_diagnostic.energy_derivative_w}",
             ]
         )

PythonModels/systems/test_mql.py

@@ -3870,6 +3870,7 @@ class TestMqlVariableChamberRhsDiagnostic:
     port_a_energy_flow_w: float
     port_b_energy_flow_w: float
     boundary_work_w: float
+    thermal_energy_flow_w: float
     energy_derivative_w: float
 
 
@@ -4120,6 +4121,9 @@ class TestMqlFullStateClosure:
         port_a_energy_flow = port_a_mass_flow * port_a_inlet_h
         port_b_energy_flow = port_b_mass_flow * port_b_inlet_h
         boundary_work = -chamber_properties.p * chamber.volume_rate_m3_s()
+        thermal_energy_flow = chamber.thermal_energy_flow_w(
+            chamber_properties.T
+        )
         return TestMqlVariableChamberRhsDiagnostic(
             chamber_alias=chamber_alias,
             piston_alias=piston_alias,
@@ -4133,8 +4137,12 @@ class TestMqlFullStateClosure:
             port_a_energy_flow_w=port_a_energy_flow,
             port_b_energy_flow_w=port_b_energy_flow,
             boundary_work_w=boundary_work,
+            thermal_energy_flow_w=thermal_energy_flow,
             energy_derivative_w=(
-                port_a_energy_flow + port_b_energy_flow + boundary_work
+                port_a_energy_flow
+                + port_b_energy_flow
+                + boundary_work
+                + thermal_energy_flow
             ),
         )
 

PythonModels/systems/test_mql_pneumatic.py

@@ -229,6 +229,9 @@ def _build_chamber(
             gas=gas,
             p0=initial_pressure_pa,
             T0=_component_temperature(component),
+            heat_transfer_coefficient=component.parameter_value("kth"),
+            heat_transfer_area=component.parameter_value("sth"),
+            external_temperature_k=_component_temperature(component),
         )
     return AmesimPneumaticVolume.from_liters(
         name=component.alias,
@@ -236,6 +239,9 @@ def _build_chamber(
         gas=gas,
         p0=initial_pressure_pa,
         T0=_component_temperature(component),
+        heat_transfer_coefficient=component.parameter_value("kth"),
+        heat_transfer_area=component.parameter_value("sth"),
+        external_temperature_k=_component_temperature(component),
     )
 
 

tests/test_amesim_pneumatic_components.py

@@ -79,6 +79,31 @@ class AmesimPneumaticComponentsTest(unittest.TestCase):
         )
         self.assertAlmostEqual(volume.port_a.p, volume.port_b.p)
 
+    def test_volume_applies_amesim_heat_exchange_to_energy_derivative(self) -> None:
+        volume = AmesimPneumaticVolume.from_liters(
+            name="heated_chamber",
+            volume_liters=15.0,
+            p0=100000.0,
+            T0=300.0,
+            heat_transfer_coefficient=1500.0,
+            heat_transfer_area=0.7,
+            external_temperature_k=293.15,
+        )
+        props = volume.properties()
+
+        derivative = volume.derivatives_from_two_connections(
+            port_a_m_flow=0.0,
+            connected_h_a=props.h,
+            port_b_m_flow=0.0,
+            connected_h_b=props.h,
+            internal_h=props.h,
+        )
+
+        self.assertAlmostEqual(
+            derivative.U,
+            1500.0 * 0.7 * (293.15 - 300.0),
+        )
+
     def test_variable_volume_tracks_external_volume(self) -> None:
         volume = AmesimVariablePneumaticVolume.from_liters(
             name="pn_c1_8",

tests/test_amesim_pnl0001_pipe.py

@@ -63,7 +63,8 @@ class AmesimPnl0001PipeTests(unittest.TestCase):
         self.assertAlmostEqual(derivative.m, 0.1)
         self.assertAlmostEqual(
             derivative.U,
-            0.2 * (internal.h + 1000.0) - 0.1 * internal.h,
+            0.2 * (internal.h + 1000.0) / self.pipe.gas.gamma
+            - 0.1 * internal.u,
         )
 
 

tests/test_core_solver.py

@@ -0,0 +1,68 @@
+import sys
+import types
+import unittest
+from unittest.mock import patch
+
+from PythonModels.core.solver import SolveIVPConfig, integrate_ode
+
+
+class IntegrateOdeTests(unittest.TestCase):
+    def test_scipy_solver_receives_step_size_controls(self) -> None:
+        calls: list[dict[str, object]] = []
+
+        def fake_solve_ivp(**kwargs):
+            calls.append(kwargs)
+            return object()
+
+        scipy_module = types.ModuleType("scipy")
+        integrate_module = types.ModuleType("scipy.integrate")
+        integrate_module.solve_ivp = fake_solve_ivp
+        scipy_module.integrate = integrate_module
+
+        with patch.dict(
+            sys.modules,
+            {"scipy": scipy_module, "scipy.integrate": integrate_module},
+        ):
+            integrate_ode(
+                rhs=lambda _time, state: state,
+                initial_state=[1.0],
+                config=SolveIVPConfig(
+                    t_start=0.0,
+                    t_stop=1.0,
+                    max_step=1.0e-4,
+                    first_step=1.0e-8,
+                ),
+                t_eval=[0.0, 1.0],
+            )
+
+        self.assertEqual(calls[0]["max_step"], 1.0e-4)
+        self.assertEqual(calls[0]["first_step"], 1.0e-8)
+
+    def test_scipy_solver_omits_unset_first_step(self) -> None:
+        calls: list[dict[str, object]] = []
+
+        def fake_solve_ivp(**kwargs):
+            calls.append(kwargs)
+            return object()
+
+        scipy_module = types.ModuleType("scipy")
+        integrate_module = types.ModuleType("scipy.integrate")
+        integrate_module.solve_ivp = fake_solve_ivp
+        scipy_module.integrate = integrate_module
+
+        with patch.dict(
+            sys.modules,
+            {"scipy": scipy_module, "scipy.integrate": integrate_module},
+        ):
+            integrate_ode(
+                rhs=lambda _time, state: state,
+                initial_state=[1.0],
+                config=SolveIVPConfig(t_start=0.0, t_stop=1.0),
+            )
+
+        self.assertEqual(calls[0]["max_step"], 1.0e-3)
+        self.assertNotIn("first_step", calls[0])
+
+
+if __name__ == "__main__":
+    unittest.main()

tests/test_run_test_mql_full_state_comparison.py

@@ -89,7 +89,8 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase):
             chamber_diagnostic.energy_derivative_w,
             chamber_diagnostic.port_a_energy_flow_w
             + chamber_diagnostic.port_b_energy_flow_w
-            + chamber_diagnostic.boundary_work_w,
+            + chamber_diagnostic.boundary_work_w
+            + chamber_diagnostic.thermal_energy_flow_w,
         )
         self.assertTrue(csv_summary_exists)
         self.assertIn("Mode: Python 132 full-state closure comparison", summary)
@@ -112,6 +113,7 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase):
         self.assertIn("piston_alias=pn_brp2_8", summary)
         self.assertIn("mass_derivative_kg_s=", summary)
         self.assertIn("boundary_work_w=", summary)
+        self.assertIn("thermal_energy_flow_w=", summary)
         self.assertIn("final_python=", summary)
         self.assertIn("final_amesim=", summary)
         self.assertIn("python.press@pn_c1_8", csv_header)

tests/test_test_mql_pneumatic.py

@@ -55,6 +55,9 @@ class TestMqlPneumaticAssemblyTests(unittest.TestCase):
         self.assertAlmostEqual(fixed_chamber.volume, 0.057)
         self.assertAlmostEqual(variable_chamber.dead_volume, 0.015)
         self.assertAlmostEqual(variable_chamber.volume, 0.015)
+        self.assertAlmostEqual(variable_chamber.heat_transfer_coefficient, 1500.0)
+        self.assertAlmostEqual(variable_chamber.heat_transfer_area, 0.7)
+        self.assertAlmostEqual(variable_chamber.external_temperature, 293.15)
         self.assertIs(fixed_chamber.gas, HELIUM_PNEUMATIC_GAS)
         self.assertIs(variable_chamber.gas, HELIUM_PNEUMATIC_GAS)
         self.assertAlmostEqual(fixed_chamber.properties().T, 293.15)

tests/test_test_mql_pnl0001_segment.py

@@ -709,7 +709,8 @@ class TestMqlPn3P4NodeChamberSegmentTests(unittest.TestCase):
             diagnostic.energy_derivative_w,
             diagnostic.port_a_energy_flow_w
             + diagnostic.port_b_energy_flow_w
-            + diagnostic.boundary_work_w,
+            + diagnostic.boundary_work_w
+            + diagnostic.thermal_energy_flow_w,
         )
 
     def test_full_state_closure_applies_step_controls_to_variable_orifices(self) -> None:
@@ -731,6 +732,27 @@ class TestMqlPn3P4NodeChamberSegmentTests(unittest.TestCase):
             open_snapshot.p4_port3_remote_node_to_primary_line_flow,
             0.0,
         )
+        node_inlet_h = (
+            open_snapshot.p4_port3_remote_balance.port_2_enthalpy_flow_w
+            / open_snapshot.p4_port3_remote_node_to_primary_line_flow
+        )
+        self.assertAlmostEqual(
+            node_inlet_h,
+            open_snapshot.p4_port3_remote_orifice_line_port_2.h,
+        )
+        rhs = closure.rhs_at(0.04, state)
+        line = closure.pneumatic_closure.components.p4_port3_remote_primary_line
+        expected_line_derivative = line.derivatives_from_connections(
+            port_1_m_flow=(
+                open_snapshot.p4_port3_remote_chamber_to_line_flow
+            ),
+            connected_h_1=open_snapshot.p4_port3_remote_primary_chamber.h,
+            port_2_m_flow=(
+                open_snapshot.p4_port3_remote_node_to_primary_line_flow
+            ),
+            connected_h_2=node_inlet_h,
+        )
+        self.assertAlmostEqual(rhs[45], expected_line_derivative.U)
         self.assertAlmostEqual(
             self.system.pneumatic_assembly.variable_orifices[
                 "pn_morifice_1"