From 456c29b3b6ec529bbda67fb3d1e87440b935a316 Mon Sep 17 00:00:00 2001 From: huojiarong Date: Wed, 12 Aug 2026 11:57:42 +0000 Subject: [PATCH] =?UTF-8?q?=E9=AA=8C=E6=94=B6=E5=9B=9B=E8=B7=AF=E6=A8=A1?= =?UTF-8?q?=E5=9E=8B=E5=B9=B6=E4=BC=98=E5=8C=96=E6=8B=93=E6=89=91=E6=B1=82?= =?UTF-8?q?=E8=A7=A3=E6=80=A7=E8=83=BD?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- app/main.py | 8 +- .../components/amesim/flow/orifices.py | 2 + .../components/amesim/flow/pipes.py | 16 +- .../amesim/mechanical/translational.py | 11 +- .../components/amesim/media/mediums.py | 11 + app/simulation/core/equations.py | 2 +- app/simulation/reporting/pnl0002_replay.py | 263 ++++++++++ app/simulation/solvers/algebraic.py | 472 ++++++++++++++---- app/simulation/solvers/mechanical.py | 25 +- app/simulation/solvers/solver.py | 44 +- app/simulation/systems/generic.py | 119 ++++- tests/test_amesim_helium_medium.py | 37 ++ ...est_amesim_mechanical_public_components.py | 5 +- tests/test_amesim_mechanical_xml.py | 11 + .../test_amesim_pnl0002_pnl0003_component.py | 32 ++ tests/test_amesim_pnl00r_component.py | 2 +- tests/test_amesim_pnor001_component.py | 16 + tests/test_component_catalog.py | 8 +- tests/test_mechanical_solver_causalization.py | 30 ++ tests/test_pnl0002_replay.py | 108 ++++ 20 files changed, 1085 insertions(+), 137 deletions(-) create mode 100644 app/simulation/reporting/pnl0002_replay.py create mode 100644 tests/test_pnl0002_replay.py diff --git a/app/main.py b/app/main.py index 558e084..d1f23a9 100644 --- a/app/main.py +++ b/app/main.py @@ -689,10 +689,10 @@ def run_system_xml_simulation( t_stop=project.simulation.t_stop, method=project.simulation.method, # The pressure-flow closure is solved to a scaled 1e-7 - # residual. Asking the outer adaptive integrator for 1e-6 - # relative accuracy makes its finite-difference Jacobian chase - # algebraic solver noise after discontinuous signal events. - rtol=1.0e-5, + # residual. State-specific mechanical absolute tolerances now + # keep ideal-stop Jacobian perturbations stable, so the outer + # integrator can use its canonical 1e-6 relative accuracy. + rtol=1.0e-6, max_step=project.simulation.max_step, ), sample_step=project.simulation.step, diff --git a/app/simulation/components/amesim/flow/orifices.py b/app/simulation/components/amesim/flow/orifices.py index 2a9c56d..2d6e0d1 100644 --- a/app/simulation/components/amesim/flow/orifices.py +++ b/app/simulation/components/amesim/flow/orifices.py @@ -1,5 +1,6 @@ from __future__ import annotations +from functools import lru_cache from collections.abc import Mapping from math import isclose, log, sqrt, tanh @@ -277,6 +278,7 @@ class AmesimPnor001(AlgebraicComponent): ) ) + @lru_cache(maxsize=32768) def _one_way_flow_characteristics( self, *, diff --git a/app/simulation/components/amesim/flow/pipes.py b/app/simulation/components/amesim/flow/pipes.py index 8ebcc17..ce8c377 100644 --- a/app/simulation/components/amesim/flow/pipes.py +++ b/app/simulation/components/amesim/flow/pipes.py @@ -1,6 +1,7 @@ from __future__ import annotations from collections.abc import Mapping +from functools import lru_cache from math import isclose, log, log10, pi, sqrt, tanh from app.simulation.components.amesim.gases import ( @@ -47,7 +48,7 @@ class AmesimPnl00r(AlgebraicComponent): """ MODEL_TYPE = "amesim_pnl00r" - MODEL_VERSION = "0.2.0" + MODEL_VERSION = "0.3.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), @@ -217,8 +218,12 @@ class AmesimPnl00r(AlgebraicComponent): if self.rr <= 0.0: turbulent = smooth_turbulent else: + # Nikuradse's fully rough asymptote is the Re-independent limit + # of Colebrook. Haaland's rounded all-regime approximation is + # about 0.2003% high at the test_mql roughness values, enough to + # bias its long high-Re PNL0001 filling transient. fully_rough = 1.0 / ( - -1.8 * log10((self.rr / 3.7) ** 1.11) + -2.0 * log10(self.rr / 3.7) ) ** 2 roughness_reynolds = reynolds_number * self.rr roughness_weight = roughness_reynolds * roughness_reynolds / ( @@ -351,7 +356,7 @@ class AmesimPnl0001(ThermodynamicVolumeComponent): """AMESim PNL0001 C-R pneumatic pipe with compressibility and friction.""" MODEL_TYPE = "amesim_pnl0001" - MODEL_VERSION = "0.3.0" + MODEL_VERSION = "0.4.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), @@ -685,6 +690,7 @@ class AmesimPnl0001(ThermodynamicVolumeComponent): upper = middle return 0.5 * (lower + upper) + @lru_cache(maxsize=32768) def _one_way_pn2pipefr_mass_flow( self, *, @@ -890,7 +896,7 @@ class AmesimPnl0002(AmesimPnl0001): """AMESim PNL0002 R-C-R pneumatic pipe with one center compliance.""" MODEL_TYPE = "amesim_pnl0002" - MODEL_VERSION = "0.5.0" + MODEL_VERSION = "0.6.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), @@ -1123,7 +1129,7 @@ class AmesimPnl0003(DynamicComponent): state_size = 4 MODEL_TYPE = "amesim_pnl0003" - MODEL_VERSION = "0.3.0" + MODEL_VERSION = "0.4.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), diff --git a/app/simulation/components/amesim/mechanical/translational.py b/app/simulation/components/amesim/mechanical/translational.py index 5a52583..23a0c18 100644 --- a/app/simulation/components/amesim/mechanical/translational.py +++ b/app/simulation/components/amesim/mechanical/translational.py @@ -1082,7 +1082,14 @@ class AmesimLmechn1(AlgebraicComponent): @property def total_force(self) -> float: - return sum(self.get_port(port_name).f for port_name in self.active_ports) + # AMESim's ``tforce`` is the force transmitted by the summed branch + # ports (1..v1). Port 9 is the balancing/common port and is excluded + # from that reported value. + return sum(self.get_port(port_name).f for port_name in self.active_ports[:-1]) + + @property + def force_balance(self) -> float: + return self.total_force + self.port_9.f def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]: reference = self.port_9 @@ -1119,7 +1126,7 @@ class AmesimLmechn1(AlgebraicComponent): relation="sumToZero", variables=tuple(f"{self.name}.{port_name}.f" for port_name in self.active_ports), role="flow", - value=self.total_force, + value=self.force_balance, ) ) return tuple(residuals) diff --git a/app/simulation/components/amesim/media/mediums.py b/app/simulation/components/amesim/media/mediums.py index d9c54df..084def4 100644 --- a/app/simulation/components/amesim/media/mediums.py +++ b/app/simulation/components/amesim/media/mediums.py @@ -2,6 +2,7 @@ from __future__ import annotations from collections.abc import Callable from dataclasses import dataclass +from functools import lru_cache from typing import ClassVar from app.simulation.core.errors import RecoverableTrialStateError @@ -197,6 +198,7 @@ class AmesimHeliumPengRobinsonMedium(IdealGasMedium): h / self.R_gas - self.nasa_enthalpy_constant_K ) / self.nasa_cp_over_R + @lru_cache(maxsize=8192) def temperature_from_pressure_enthalpy(self, p: float, h: float) -> float: temperature = max(self.temperature_from_enthalpy(h), 2.2) for _iteration in range(16): @@ -220,12 +222,21 @@ class AmesimHeliumPengRobinsonMedium(IdealGasMedium): ) return self.temperature_from_internal_energy(U / m) + @lru_cache(maxsize=8192) def properties_from_mU( self, m: float, U: float, V: float, ) -> ThermodynamicProperties: + """Recover a real-gas state, reusing exact repeated evaluations. + + Implicit integration asks several component interfaces for the same + ``(m, U, V)`` state while closing one RHS evaluation and while building + finite-difference Jacobians. The calculation is pure and its result is + immutable, so an exact-key bounded cache avoids repeating the + Peng-Robinson temperature iteration without changing model semantics. + """ if m <= 0.0: raise RecoverableTrialStateError( "Mass must stay positive when recovering temperature." diff --git a/app/simulation/core/equations.py b/app/simulation/core/equations.py index 7c28b48..ece76f1 100644 --- a/app/simulation/core/equations.py +++ b/app/simulation/core/equations.py @@ -10,7 +10,7 @@ EquationOwner = Literal["connection", "component"] EquationRelation = Literal["equal", "sumToZero", "constitutive", "state"] -@dataclass(frozen=True) +@dataclass(frozen=True, slots=True) class EquationResidual: """One executable scalar equation in the pressure-flow subsystem.""" diff --git a/app/simulation/reporting/pnl0002_replay.py b/app/simulation/reporting/pnl0002_replay.py new file mode 100644 index 0000000..29806b4 --- /dev/null +++ b/app/simulation/reporting/pnl0002_replay.py @@ -0,0 +1,263 @@ +from __future__ import annotations + +import argparse +import json +from dataclasses import dataclass +from pathlib import Path +from typing import Mapping + +from app.simulation.components.amesim.flow.pipes import AmesimPnl0002 +from app.simulation.reporting.amesim_results import ( + AmesimResults, + load_test_mql_amesim_results, +) + + +@dataclass(frozen=True, slots=True) +class Pnl0002ReplayPaths: + """AMESim data paths needed to replay one PNL0002 resistance.""" + + center_pressure: str + center_temperature: str + port_1_pressure: str + port_1_temperature: str + port_1_mass_flow: str + port_2_pressure: str + port_2_temperature: str + port_2_mass_flow: str + reynolds: str + friction_factor: str + + +PNL83_REPLAY_PATHS = Pnl0002ReplayPaths( + center_pressure="pctr@pneumatic_83", + center_temperature="tctr@pneumatic_83", + port_1_pressure="press1@pnnode4_16", + port_1_temperature="temp1@pnnode4_16", + port_1_mass_flow="dm1@pneumatic_83", + port_2_pressure="press3@pnnode4_17", + port_2_temperature="temp3@pnnode4_17", + port_2_mass_flow="dm2@pneumatic_83", + reynolds="re@pneumatic_83", + friction_factor="ff@pneumatic_83", +) + + +def _required_series( + results: AmesimResults, + path: str, +) -> tuple[float, ...]: + try: + values = results.series(path) + except KeyError as exc: + raise ValueError(f"AMESim replay variable is not saved: {path}") from exc + if len(values) != len(results.times): + raise ValueError(f"AMESim replay variable has an invalid length: {path}") + return values + + +def _metric_summary( + rows: list[dict[str, float]], + key: str, +) -> dict[str, float]: + values = [float(row[key]) for row in rows] + max_index = max(range(len(values)), key=lambda index: abs(values[index])) + return { + "maxAbs": abs(values[max_index]), + "maxAbsTime": rows[max_index]["time"], + "finalSigned": values[-1], + } + + +def replay_pnl0002_amesim_states( + pipe: AmesimPnl0002, + results: AmesimResults, + paths: Pnl0002ReplayPaths, + *, + amesim_mass_flow_scale: float = -1.0e-3, +) -> dict[str, object]: + """Replay saved AMESim states through current PNL0002 flow functions. + + This is a calibration-only, no-integration calculation. It does not write + states into the pipe or alter the production simulation path. + """ + + series_by_field = { + field: _required_series(results, getattr(paths, field)) + for field in paths.__dataclass_fields__ + } + rows: list[dict[str, float]] = [] + for index, time_s in enumerate(results.times): + center_pressure = series_by_field["center_pressure"][index] + center_temperature = series_by_field["center_temperature"][index] + port_1_pressure = series_by_field["port_1_pressure"][index] + port_2_pressure = series_by_field["port_2_pressure"][index] + observed_flow_1 = ( + amesim_mass_flow_scale + * series_by_field["port_1_mass_flow"][index] + ) + observed_flow_2 = ( + amesim_mass_flow_scale + * series_by_field["port_2_mass_flow"][index] + ) + upstream_temperature_1 = ( + series_by_field["port_1_temperature"][index] + if observed_flow_1 >= 0.0 + else center_temperature + ) + upstream_temperature_2 = ( + series_by_field["port_2_temperature"][index] + if observed_flow_2 >= 0.0 + else center_temperature + ) + predicted_flow_1 = pipe.mass_flow( + port_1_pressure, + center_pressure, + upstream_temperature_1, + ) + predicted_flow_2 = pipe.mass_flow( + port_2_pressure, + center_pressure, + upstream_temperature_2, + ) + reynolds_1 = pipe.reynolds_number( + observed_flow_1, + upstream_temperature_1, + ) + reynolds_2 = pipe.reynolds_number( + observed_flow_2, + upstream_temperature_2, + ) + friction_1 = pipe.friction_factor(reynolds_1) + friction_2 = pipe.friction_factor(reynolds_2) + replay_reynolds = 0.5 * (reynolds_1 + reynolds_2) + replay_friction = 0.5 * (friction_1 + friction_2) + rows.append( + { + "time": float(time_s), + "centerPressure": center_pressure, + "centerTemperature": center_temperature, + "port1Pressure": port_1_pressure, + "port2Pressure": port_2_pressure, + "observedPort1MassFlow": observed_flow_1, + "observedPort2MassFlow": observed_flow_2, + "predictedPort1MassFlow": predicted_flow_1, + "predictedPort2MassFlow": predicted_flow_2, + "port1MassFlowError": predicted_flow_1 - observed_flow_1, + "port2MassFlowError": predicted_flow_2 - observed_flow_2, + "amesimReynolds": series_by_field["reynolds"][index], + "replayReynolds": replay_reynolds, + "reynoldsError": ( + replay_reynolds - series_by_field["reynolds"][index] + ), + "amesimFrictionFactor": series_by_field["friction_factor"][index], + "replayFrictionFactor": replay_friction, + "frictionFactorError": ( + replay_friction + - series_by_field["friction_factor"][index] + ), + } + ) + + metric_keys = ( + "port1MassFlowError", + "port2MassFlowError", + "reynoldsError", + "frictionFactorError", + ) + return { + "mode": "amesim-state-replay-no-integration", + "component": pipe.name, + "pointCount": len(rows), + "massFlowScale": amesim_mass_flow_scale, + "paths": { + field: getattr(paths, field) + for field in paths.__dataclass_fields__ + }, + "parameters": dict(pipe.parameter_values), + "summary": { + key: _metric_summary(rows, key) + for key in metric_keys + }, + "rows": rows, + } + + +def _compile_project_pipe( + project_path: Path, + component_name: str, +) -> AmesimPnl0002: + from app.main import ( + ReactFlowProjectPayload, + _compile_xml_document_or_422, + _validated_xml_document_or_422, + build_reactflow_system_xml, + validate_system_xml_document, + ) + + payload = ReactFlowProjectPayload.model_validate_json( + project_path.read_text(encoding="utf-8") + ) + xml_bytes = build_reactflow_system_xml(payload) + document = _validated_xml_document_or_422( + validate_system_xml_document(xml_bytes) + ) + _project, network = _compile_xml_document_or_422(document) + try: + component = network.components[component_name] + except KeyError as exc: + raise ValueError(f"Project component does not exist: {component_name}") from exc + if not isinstance(component, AmesimPnl0002): + raise ValueError(f"Project component is not PNL0002: {component_name}") + return component + + +def _paths_from_arguments(arguments: argparse.Namespace) -> Pnl0002ReplayPaths: + values: Mapping[str, str] = { + field: getattr(arguments, field) + for field in PNL83_REPLAY_PATHS.__dataclass_fields__ + } + return Pnl0002ReplayPaths(**values) + + +def main() -> None: + parser = argparse.ArgumentParser( + description="Replay saved AMESim p/T/m_flow through a PNL0002 model without integration." + ) + parser.add_argument("project", type=Path) + parser.add_argument("amesim_archive", type=Path) + parser.add_argument("output", type=Path) + parser.add_argument("--component", default="pneumatic_83") + parser.add_argument("--mass-flow-scale", type=float, default=-1.0e-3) + for field in PNL83_REPLAY_PATHS.__dataclass_fields__: + parser.add_argument( + "--" + field.replace("_", "-"), + dest=field, + default=getattr(PNL83_REPLAY_PATHS, field), + ) + arguments = parser.parse_args() + + pipe = _compile_project_pipe(arguments.project, arguments.component) + results = load_test_mql_amesim_results(arguments.amesim_archive) + report = replay_pnl0002_amesim_states( + pipe, + results, + _paths_from_arguments(arguments), + amesim_mass_flow_scale=arguments.mass_flow_scale, + ) + arguments.output.parent.mkdir(parents=True, exist_ok=True) + arguments.output.write_text( + json.dumps(report, ensure_ascii=False, indent=2) + "\n", + encoding="utf-8", + ) + print(json.dumps({key: report[key] for key in ( + "mode", + "component", + "pointCount", + "parameters", + "summary", + )}, ensure_ascii=False, indent=2)) + + +if __name__ == "__main__": + main() diff --git a/app/simulation/solvers/algebraic.py b/app/simulation/solvers/algebraic.py index 6130a86..bd65286 100644 --- a/app/simulation/solvers/algebraic.py +++ b/app/simulation/solvers/algebraic.py @@ -12,6 +12,7 @@ from app.simulation.components.amesim.flow.pipes import ( AmesimPnl0001, AmesimPnl0002, ) +from app.simulation.core.equations import EquationResidual from app.simulation.core.ports import PortState, VariableRole from app.simulation.systems.network import SimulationNetwork @@ -45,13 +46,45 @@ class AlgebraicUnknown: class ExplicitFlowAssignment: equation_id: str unknown: AlgebraicUnknown - evaluate: Callable[[], float] + evaluate: Callable[[], float] | None + component: object | None = None + equation_index: int | None = None + + +@dataclass(frozen=True) +class ExplicitFlowStage: + assignments: tuple[ExplicitFlowAssignment, ...] + @dataclass(frozen=True) class EffortAnchor: unknown: AlgebraicUnknown evaluate: Callable[[], float] + +@dataclass(frozen=True) +class ConnectionEquationEvaluation: + template: EquationResidual + evaluate: Callable[[], float] + + +@dataclass(frozen=True) +class PnorPnl0001SeriesBinding: + orifice: AmesimPnor001 + orifice_port: str + pipe: AmesimPnl0001 + pipe_port: str + + +@dataclass(frozen=True) +class ClosedResistancePressureBinding: + component: object + port_name: str + neighbor: object + neighbor_port: str + pressure_source_port: str | None + + @dataclass(frozen=True) class EffortEqualityGroup: variable: str @@ -105,10 +138,43 @@ class PressureFlowSolver: self.max_evaluations = max_evaluations self.unknowns = self._build_unknowns() self._unknowns_by_id = {unknown.id: unknown for unknown in self.unknowns} + self._unknowns_by_variable = { + variable: tuple( + unknown + for unknown in self.unknowns + if unknown.variable == variable + ) + for variable in ("p", "m_flow", "x", "v", "f") + } + self._component_equation_owners = tuple(network.components.values()) + self._estimated_flow_components = tuple( + component + for component in self._component_equation_owners + if hasattr(component, "K_eff") + ) + self._causal_contact_components = tuple( + component + for component in self._component_equation_owners + if getattr(component, "clear_causal_contact", None) is not None + ) + self._connection_equation_plan = tuple( + ConnectionEquationEvaluation( + template=equation, + evaluate=self._equation_value_reader(equation), + ) + for equation in network.connection_equation_residuals() + ) self._effort_groups = { variable: self._build_effort_equality_groups(variable) for variable in ("p", "x", "v") } + self._pnor_pnl0001_series_plan = ( + self._build_pnor_pnl0001_series_plan() + ) + self._closed_resistance_pressure_plan = ( + self._build_closed_resistance_pressure_plan() + ) + self._unilateral_contact_plan = self._build_unilateral_contact_plan() self._explicit_flow_plan = self._build_explicit_flow_plan() self.last_diagnostics: AlgebraicSolveDiagnostics | None = None @@ -143,7 +209,10 @@ class PressureFlowSolver: return None return component_name, port_name - def _seed_equal_efforts(self) -> None: + def _seed_equal_efforts( + self, + variables: tuple[str, ...] = ("p", "x", "v"), + ) -> None: """Lift state-owned efforts across their complete equality groups. Dynamic components refresh their own ports before each closure, while @@ -154,7 +223,19 @@ class PressureFlowSolver: before evaluating explicit flow laws. """ - for variable in ("p", "x", "v"): + self.propagate_equal_efforts(variables) + + def propagate_equal_efforts(self, variables: tuple[str, ...]) -> None: + """Propagate selected state-owned efforts without solving flows. + + Piston geometry needs current mechanical ``x``/``v`` before swept + volume propagation, but pressure and flow equations can wait until the + connected chamber has refreshed that volume. + """ + unknown = sorted(set(variables) - set(self._effort_groups)) + if unknown: + raise ValueError("Unsupported effort variables: " + ", ".join(unknown)) + for variable in variables: self._seed_equal_effort(variable) def _build_effort_equality_groups( @@ -538,10 +619,10 @@ class PressureFlowSolver: for binding in bindings: self._apply_unilateral_contact_binding(binding) - def _seed_unilateral_contacts( + def _build_unilateral_contact_plan( self, ) -> tuple[UnilateralContactBinding, ...]: - """Create local eliminations for contacts with one algebraic coordinate.""" + """Compile contacts that can eliminate one algebraic coordinate.""" position_groups = { unknown.id: group @@ -549,7 +630,6 @@ class PressureFlowSolver: for unknown in group.members } bindings: list[UnilateralContactBinding] = [] - bound_group_ids: set[int] = set() for component in self.network.components.values(): if component.model_type != "amesim_lstp00a": continue @@ -591,6 +671,18 @@ class PressureFlowSolver: # With both coordinates state-owned, penetration is a dynamic # result rather than an algebraic active-set choice. continue + bindings.append(binding) + + return tuple(bindings) + + def _seed_unilateral_contacts( + self, + ) -> tuple[UnilateralContactBinding, ...]: + """Apply compiled local contact eliminations for the current state.""" + + bindings: list[UnilateralContactBinding] = [] + bound_group_ids: set[int] = set() + for binding in self._unilateral_contact_plan: group_id = id(binding.algebraic_group) if group_id in bound_group_ids: # One relative contact law may eliminate a free coordinate. @@ -630,33 +722,87 @@ class PressureFlowSolver: component = self.network.components[equation.owner_id] equation_id = equation.id - - def read_component_equation() -> float: - for current in component.pressure_flow_equation_residuals(): - if current.id == equation_id: - return float(current.value) + equation_ids = tuple( + current.id + for current in component.pressure_flow_equation_residuals() + ) + try: + equation_index = equation_ids.index(equation_id) + except ValueError as exc: raise RuntimeError( f"Compiled algebraic equation disappeared at runtime: {equation_id}." - ) + ) from exc + + def read_component_equation() -> float: + current_equations = component.pressure_flow_equation_residuals() + if ( + equation_index >= len(current_equations) + or current_equations[equation_index].id != equation_id + ): + raise RuntimeError( + f"Compiled algebraic equation disappeared at runtime: {equation_id}." + ) + return float(current_equations[equation_index].value) return read_component_equation + def _pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]: + """Evaluate live values through a precompiled connector topology.""" + component_residuals = tuple( + residual + for component in self._component_equation_owners + for residual in component.pressure_flow_equation_residuals() + ) + connection_residuals = tuple( + EquationResidual( + id=item.template.id, + owner=item.template.owner, + owner_id=item.template.owner_id, + relation=item.template.relation, + variables=item.template.variables, + value=item.evaluate(), + role=item.template.role, + ) + for item in self._connection_equation_plan + ) + return component_residuals + connection_residuals - def _build_explicit_flow_plan(self) -> tuple[ExplicitFlowAssignment, ...]: - """Compile the legacy deterministic flow assignment order once.""" - assignments: list[ExplicitFlowAssignment] = [] + def _explicit_flow_assignment( + self, + equation, + unknown: AlgebraicUnknown, + ) -> ExplicitFlowAssignment: + if equation.owner == "connection": + return ExplicitFlowAssignment( + equation_id=equation.id, + unknown=unknown, + evaluate=self._equation_value_reader(equation), + ) + + component = self.network.components[equation.owner_id] + equations = component.pressure_flow_equation_residuals() + equation_ids = tuple(current.id for current in equations) + try: + equation_index = equation_ids.index(equation.id) + except ValueError as exc: + raise RuntimeError( + f"Compiled algebraic equation disappeared at runtime: {equation.id}." + ) from exc + return ExplicitFlowAssignment( + equation_id=equation.id, + unknown=unknown, + evaluate=None, + component=component, + equation_index=equation_index, + ) + + def _build_explicit_flow_plan(self) -> tuple[ExplicitFlowStage, ...]: + """Compile flow causalization into independent dependency stages.""" + + stages: list[ExplicitFlowStage] = [] seeded_ids: set[str] = set() - def append_assignment(equation, unknown: AlgebraicUnknown) -> None: - assignments.append( - ExplicitFlowAssignment( - equation_id=equation.id, - unknown=unknown, - evaluate=self._equation_value_reader(equation), - ) - ) - seeded_ids.add(unknown.id) - + initial_assignments: list[ExplicitFlowAssignment] = [] for component in self.network.components.values(): for equation in component.pressure_flow_equation_residuals(): if equation.relation != "constitutive" or equation.role != "flow": @@ -665,12 +811,19 @@ class PressureFlowSolver: if len(flow_unknowns) != 1: continue unknown = flow_unknowns[0] - if unknown.id not in seeded_ids: - append_assignment(equation, unknown) + if unknown.id in seeded_ids: + continue + initial_assignments.append( + self._explicit_flow_assignment(equation, unknown) + ) + seeded_ids.add(unknown.id) + if initial_assignments: + stages.append(ExplicitFlowStage(tuple(initial_assignments))) - equations = self.network.pressure_flow_equation_residuals() + equations = self._pressure_flow_equation_residuals() while True: - propagated = False + stage_assignments: list[ExplicitFlowAssignment] = [] + stage_unknown_ids: set[str] = set() for equation in equations: if equation.role != "flow" or equation.relation not in { "constitutive", @@ -689,12 +842,19 @@ class PressureFlowSolver: ) if len(unseeded) != 1: continue - append_assignment(equation, unseeded[0]) - propagated = True - break - if propagated: + unknown = unseeded[0] + if unknown.id in stage_unknown_ids: + continue + stage_assignments.append( + self._explicit_flow_assignment(equation, unknown) + ) + stage_unknown_ids.add(unknown.id) + if stage_assignments: + stages.append(ExplicitFlowStage(tuple(stage_assignments))) + seeded_ids.update(stage_unknown_ids) continue + fallback_assignment: ExplicitFlowAssignment | None = None for equation in equations: if equation.role != "flow" or equation.relation not in { "constitutive", @@ -711,40 +871,92 @@ class PressureFlowSolver: continue if len({unknown.variable for unknown in flow_unknowns}) != 1: continue - append_assignment(equation, unseeded[-1]) - propagated = True + unknown = unseeded[-1] + fallback_assignment = self._explicit_flow_assignment( + equation, + unknown, + ) + seeded_ids.add(unknown.id) break - if not propagated: + if fallback_assignment is None: break + stages.append(ExplicitFlowStage((fallback_assignment,))) - return tuple(assignments) + return tuple(stages) - def _solve_explicit_flow_unknowns(self) -> set[str]: - """Execute the precompiled explicit flow/force causalization plan.""" + @staticmethod + def _evaluate_explicit_flow_stage( + assignments: tuple[ExplicitFlowAssignment, ...], + ) -> dict[str, float]: + values: dict[str, float] = {} + assignments_by_component: dict[object, list[ExplicitFlowAssignment]] = {} + for assignment in assignments: + if assignment.component is None: + assert assignment.evaluate is not None + values[assignment.equation_id] = assignment.evaluate() + continue + assignments_by_component.setdefault(assignment.component, []).append( + assignment + ) + + for component, component_assignments in assignments_by_component.items(): + equations = component.pressure_flow_equation_residuals() + for assignment in component_assignments: + assert assignment.equation_index is not None + equation_index = assignment.equation_index + if ( + equation_index >= len(equations) + or equations[equation_index].id != assignment.equation_id + ): + raise RuntimeError( + "Compiled algebraic equation disappeared at runtime: " + f"{assignment.equation_id}." + ) + values[assignment.equation_id] = float( + equations[equation_index].value + ) + return values + + def _solve_explicit_flow_unknowns( + self, + variables: tuple[str, ...] = ("f", "m_flow"), + ) -> set[str]: + """Execute staged flow/force assignments without repeated equations.""" + + selected = frozenset(variables) for unknown in self.unknowns: - if unknown.variable in {"f", "m_flow"}: + if unknown.variable in selected: unknown.write(0.0) seeded_ids: set[str] = set() - for assignment in self._explicit_flow_plan: - target_value = assignment.unknown.read() - assignment.evaluate() - if not isfinite(target_value): - continue - assignment.unknown.write(target_value) - seeded_ids.add(assignment.unknown.id) + for stage in self._explicit_flow_plan: + assignments = tuple( + assignment + for assignment in stage.assignments + if assignment.unknown.variable in selected + ) + values = self._evaluate_explicit_flow_stage(assignments) + targets = tuple( + ( + assignment, + assignment.unknown.read() - values[assignment.equation_id], + ) + for assignment in assignments + ) + for assignment, target_value in targets: + if not isfinite(target_value): + continue + assignment.unknown.write(target_value) + seeded_ids.add(assignment.unknown.id) return seeded_ids - def _seed_closed_resistance_pressures(self) -> None: - """Seed a sealed resistance end at its zero-flow pressure. - - A PNPL01 fixes flow, not pressure. Starting a dead-ended Darcy branch - with the plug-side pressure at the medium reference can otherwise put - the nonlinear solver on the singular square-root part of the inverse - flow law. At zero flow, these AMESim pipe resistances have exactly zero - pressure drop, which gives a deterministic and physically exact seed. - """ + def _build_closed_resistance_pressure_plan( + self, + ) -> tuple[ClosedResistancePressureBinding, ...]: + """Compile sealed resistance ends whose zero-flow pressure is known.""" + bindings: list[ClosedResistancePressureBinding] = [] connected: dict[tuple[str, str], tuple[str, str]] = {} for connection in self.network.connections: if connection.kind != "physical" or connection.domain != "pneumatic": @@ -764,20 +976,50 @@ class PressureFlowSolver: if not isinstance(neighbor, AmesimPnpl01): continue if isinstance(component, AmesimPnl0002): - pressure = component.properties().p + pressure_source_port = None elif isinstance(component, AmesimPnl0001): if port_name != "port_1": continue - pressure = component.properties().p + pressure_source_port = None else: - other_port_name = "port_2" if port_name == "port_1" else "port_1" - pressure = component.get_port(other_port_name).p - component.get_port(port_name).p = pressure - neighbor.get_port(neighbor_key[1]).p = pressure + pressure_source_port = ( + "port_2" if port_name == "port_1" else "port_1" + ) + bindings.append( + ClosedResistancePressureBinding( + component=component, + port_name=port_name, + neighbor=neighbor, + neighbor_port=neighbor_key[1], + pressure_source_port=pressure_source_port, + ) + ) + return tuple(bindings) - def _seed_pnor_pnl0001_series_pressures(self) -> None: - """Causalize the pressure between a PNOR001 and PNL0001 R port.""" + def _seed_closed_resistance_pressures(self) -> None: + """Seed a sealed resistance end at its zero-flow pressure. + A PNPL01 fixes flow, not pressure. Starting a dead-ended Darcy branch + with the plug-side pressure at the medium reference can otherwise put + the nonlinear solver on the singular square-root part of the inverse + flow law. At zero flow, these AMESim pipe resistances have exactly zero + pressure drop, which gives a deterministic and physically exact seed. + """ + + for binding in self._closed_resistance_pressure_plan: + component = binding.component + pressure = ( + component.properties().p + if binding.pressure_source_port is None + else component.get_port(binding.pressure_source_port).p + ) + component.get_port(binding.port_name).p = pressure + binding.neighbor.get_port(binding.neighbor_port).p = pressure + + def _build_pnor_pnl0001_series_plan( + self, + ) -> tuple[PnorPnl0001SeriesBinding, ...]: + bindings: list[PnorPnl0001SeriesBinding] = [] for connection in self.network.connections: first_endpoint, second_endpoint = connection.endpoints first = self.network.components[first_endpoint.component] @@ -792,6 +1034,26 @@ class PressureFlowSolver: continue if isinstance(pipe, AmesimPnl0002) or pipe_port != "port_1": continue + bindings.append( + PnorPnl0001SeriesBinding( + orifice=orifice, + orifice_port=orifice_port, + pipe=pipe, + pipe_port=pipe_port, + ) + ) + return tuple(bindings) + + def _seed_pnor_pnl0001_series_pressures(self) -> None: + """Causalize the pressure between a PNOR001 and PNL0001 R port.""" + + from scipy.optimize import brentq + + for binding in self._pnor_pnl0001_series_plan: + orifice = binding.orifice + orifice_port = binding.orifice_port + pipe = binding.pipe + pipe_port = binding.pipe_port orifice_other = "port_2" if orifice_port == "port_1" else "port_1" pressure_a = orifice.get_port(orifice_other).p @@ -826,15 +1088,16 @@ class PressureFlowSolver: elif (lower_value < 0.0) == (upper_value < 0.0): continue else: - for _iteration in range(64): - middle = 0.5 * (lower + upper) - middle_value = mismatch(middle) - if (middle_value < 0.0) == (lower_value < 0.0): - lower = middle - lower_value = middle_value - else: - upper = middle - pressure = 0.5 * (lower + upper) + pressure = float( + brentq( + mismatch, + lower, + upper, + xtol=1.0e-6, + rtol=1.0e-12, + maxiter=32, + ) + ) orifice.get_port(orifice_port).p = pressure pipe.get_port(pipe_port).p = pressure @@ -842,22 +1105,20 @@ class PressureFlowSolver: pressure_scale = max( [ abs(unknown.read()) - for unknown in self.unknowns - if unknown.variable == "p" and unknown.read() > 0.0 + for unknown in self._unknowns_by_variable["p"] + if unknown.read() > 0.0 ] + [1e5] ) estimated_flows = [ abs(float(getattr(component, "K_eff"))) * sqrt(pressure_scale) - for component in self.network.components.values() - if hasattr(component, "K_eff") + for component in self._estimated_flow_components ] mass_flow_scale = max( estimated_flows + [ abs(unknown.read()) - for unknown in self.unknowns - if unknown.variable == "m_flow" + for unknown in self._unknowns_by_variable["m_flow"] ] + [1e-3] ) @@ -865,20 +1126,24 @@ class PressureFlowSolver: "p": pressure_scale, "m_flow": mass_flow_scale, "x": max( - [abs(unknown.read()) for unknown in self.unknowns if unknown.variable == "x"] + [abs(unknown.read()) for unknown in self._unknowns_by_variable["x"]] + [1.0] ), "v": max( - [abs(unknown.read()) for unknown in self.unknowns if unknown.variable == "v"] + [abs(unknown.read()) for unknown in self._unknowns_by_variable["v"]] + [1.0] ), "f": max( - [abs(unknown.read()) for unknown in self.unknowns if unknown.variable == "f"] + [abs(unknown.read()) for unknown in self._unknowns_by_variable["f"]] + [1.0] ), } - def solve(self) -> AlgebraicSolveDiagnostics: + def solve( + self, + *, + effort_variables: tuple[str, ...] = ("p", "x", "v"), + ) -> AlgebraicSolveDiagnostics: try: import numpy as np from scipy.optimize import least_squares @@ -887,18 +1152,16 @@ class PressureFlowSolver: "Topology-driven simulation requires SciPy; install requirements.txt." ) from exc - for component in self.network.components.values(): - clear_causal_contact = getattr(component, "clear_causal_contact", None) - if clear_causal_contact is not None: - clear_causal_contact() + for component in self._causal_contact_components: + component.clear_causal_contact() - self._seed_equal_efforts() + self._seed_equal_efforts(effort_variables) self._seed_closed_resistance_pressures() self._seed_pnor_pnl0001_series_pressures() self._solve_explicit_flow_unknowns() contact_bindings = self._seed_unilateral_contacts() if contact_bindings: - self._solve_explicit_flow_unknowns() + self._solve_explicit_flow_unknowns(("f",)) self._refresh_unilateral_contacts(contact_bindings) scales = self._scales() pressure_scale = scales["p"] @@ -911,21 +1174,10 @@ class PressureFlowSolver: ) for unknown in self.unknowns } - positive_pressures = [ - unknown.read() - for unknown in self.unknowns - if unknown.variable == "p" and unknown.read() > 0.0 - ] - fallback_pressure = ( - sum(positive_pressures) / len(positive_pressures) - if positive_pressures - else pressure_scale - ) - def variable_scale(unknown: AlgebraicUnknown) -> float: return unknown_scales[unknown.id] - seeded_equations = self.network.pressure_flow_equation_residuals() + seeded_equations = self._pressure_flow_equation_residuals() def initial_equation_scale(equation) -> float: variable_names = [ @@ -959,7 +1211,10 @@ class PressureFlowSolver: } def equation_scale(equation) -> float: - return equation_scales.get(equation.id, initial_equation_scale(equation)) + cached = equation_scales.get(equation.id) + if cached is not None: + return cached + return initial_equation_scale(equation) seeded_scaled = [ abs(equation.value / equation_scale(equation)) @@ -994,6 +1249,17 @@ class PressureFlowSolver: self.last_diagnostics = diagnostics return diagnostics + positive_pressures = [ + unknown.read() + for unknown in self._unknowns_by_variable["p"] + if unknown.read() > 0.0 + ] + fallback_pressure = ( + sum(positive_pressures) / len(positive_pressures) + if positive_pressures + else pressure_scale + ) + # A causal contact retains its small relative penetration around the # current absolute port coordinates. Keep that local coordinate during # nonlinear fallback: the contact law remains responsive to optimizer @@ -1027,7 +1293,7 @@ class PressureFlowSolver: def scaled_residuals(values): assign(values) self._refresh_unilateral_contacts(contact_bindings) - equations = self.network.pressure_flow_equation_residuals() + equations = self._pressure_flow_equation_residuals() return np.asarray( [ equation.value / equation_scale(equation) @@ -1048,7 +1314,7 @@ class PressureFlowSolver: ) assign(result.x) self._refresh_unilateral_contacts(contact_bindings) - equations = self.network.pressure_flow_equation_residuals() + equations = self._pressure_flow_equation_residuals() scaled = [ abs( equation.value / equation_scale(equation) diff --git a/app/simulation/solvers/mechanical.py b/app/simulation/solvers/mechanical.py index 404b50f..3827185 100644 --- a/app/simulation/solvers/mechanical.py +++ b/app/simulation/solvers/mechanical.py @@ -391,6 +391,27 @@ class MechanicalStateReducer: def has_state_events(self) -> bool: return any(group.discrete_endstop_components for group in self.groups) + def absolute_tolerances( + self, + default: float, + *, + mechanical: float = 1.0e-12, + ) -> list[float]: + """Return state-aligned tolerances with machine-scale mechanics. + + A scalar ``1e-8`` absolute tolerance makes SciPy perturb a zero-valued + endstop position across the much smaller unilateral boundary band while + constructing finite-difference Jacobians. Mechanical coordinates need + a tighter floor; thermodynamic states retain the caller's tolerance. + """ + values: list[float] = [] + for entry in self.state_entries: + if isinstance(entry, MechanicalConstraintGroup): + values.extend([min(default, mechanical)] * 2) + else: + values.extend([default] * entry.state_size) + return values + def reset_constraint_modes(self) -> None: for group in self.groups: group.reset_mode() @@ -518,7 +539,7 @@ class MechanicalStateReducer: candidates.append((previous_time, group, "lower", lower)) elif ( lower is not None - and previous_position > lower + and previous_position > lower + group._boundary_tolerance(lower) and current_position <= lower ): candidates.append( @@ -573,7 +594,7 @@ class MechanicalStateReducer: candidates.append((previous_time, group, "upper", upper)) elif ( upper is not None - and previous_position < upper + and previous_position < upper - group._boundary_tolerance(upper) and current_position >= upper ): candidates.append( diff --git a/app/simulation/solvers/solver.py b/app/simulation/solvers/solver.py index 389649e..5e23d4d 100644 --- a/app/simulation/solvers/solver.py +++ b/app/simulation/solvers/solver.py @@ -39,7 +39,7 @@ class SolveIVPConfig: t_stop: float = 20.0 method: str = "BDF" rtol: float = 1e-6 - atol: float = 1e-8 + atol: float | Sequence[float] = 1e-8 max_step: float = 1e-3 first_step: float | None = None @@ -88,6 +88,34 @@ def _append_or_replace_solution_sample( return _append_solution_sample(times, states, time, state) +def _project_nearby_pre_transition_sample( + times: list[float], + states: list[list[float]], + transition: StateTransition, + config: SolveIVPConfig, +) -> None: + """Resolve a sample/event ordering that is below solver time precision. + + An adaptive dense interpolant can place a discontinuous impact a few + nanoseconds after its analytically coincident output sample. Keep the + located event and restart time unchanged, but report that ambiguous sample + on the reset side of the discontinuity. + """ + if not times or not math.isfinite(config.max_step): + return + time_gap = float(transition.time) - times[-1] + tolerance = max( + 64.0 * math.ulp(max(abs(float(transition.time)), 1.0)), + min( + abs(float(config.max_step) * float(config.rtol)), + 1.0e-8, + ), + ) + if not 0.0 < time_gap <= tolerance: + return + for index, value in enumerate(transition.state): + states[index][-1] = float(value) + def _normalize_state_transition( transition: StateTransition, @@ -534,6 +562,7 @@ def _integrate_scipy_stepwise( accepted_step_callback: AcceptedStepCallback | None, breakpoints: Sequence[float] = (), state_transition_handler: StateTransitionHandler | None = None, + jac_sparsity=None, ) -> ODESolution: """Initial stepwise integration path for breakpoints and state resets. @@ -625,6 +654,8 @@ def _integrate_scipy_stepwise( "atol": config.atol, "max_step": segment_max_step, } + if jac_sparsity is not None and config.method in {"BDF", "Radau"}: + solver_options["jac_sparsity"] = jac_sparsity requested_first_step = ( 0.1 * segment_max_step if last_recoverable_error is not None @@ -666,6 +697,7 @@ def _integrate_scipy_stepwise( error = exc break + restart_at_transition = False restart_after_recoverable = False while solver.status == "running": @@ -806,6 +838,12 @@ def _integrate_scipy_stepwise( ) sample_index += 1 + _project_nearby_pre_transition_sample( + times, + states, + transition, + config, + ) last_accepted_time = transition.time last_accepted_state = list(transition.state) last_transition = transition @@ -923,6 +961,7 @@ def integrate_ode( accepted_step_callback: AcceptedStepCallback | None = None, breakpoints: Sequence[float] | None = None, state_transition_handler: StateTransitionHandler | None = None, + jac_sparsity=None, ): """Integrate an ODE, optionally restarting at equation discontinuities. @@ -992,6 +1031,7 @@ def integrate_ode( accepted_step_callback, normalized_breakpoints, state_transition_handler, + jac_sparsity, ) solve_options = { @@ -1006,4 +1046,6 @@ def integrate_ode( } if config.first_step is not None: solve_options["first_step"] = config.first_step + if jac_sparsity is not None and config.method in {"BDF", "Radau"}: + solve_options["jac_sparsity"] = jac_sparsity return solve_ivp(**solve_options) diff --git a/app/simulation/systems/generic.py b/app/simulation/systems/generic.py index 20fa330..a7fd842 100644 --- a/app/simulation/systems/generic.py +++ b/app/simulation/systems/generic.py @@ -1,14 +1,17 @@ from __future__ import annotations from collections.abc import Callable -from dataclasses import dataclass +from dataclasses import dataclass, replace from math import floor, isfinite from typing import Literal from app.simulation.core.base import DynamicComponent from app.simulation.core.metadata import ResultVariableMetadata from app.simulation.solvers.algebraic import PressureFlowSolver -from app.simulation.solvers.mechanical import MechanicalStateReducer +from app.simulation.solvers.mechanical import ( + MechanicalConstraintGroup, + MechanicalStateReducer, +) from app.simulation.solvers.pneumatic_storage import ( IdealPneumaticStorageReducer, ideal_storage_group_is_reducible, @@ -267,6 +270,7 @@ class GenericFluidSystem: self.max_thermofluid_iterations = 0 self.signal_propagation_count = 0 self.pneumatic_volume_propagation_count = 0 + self._jacobian_sparsity = None def initial_state_vector(self) -> list[float]: return self.pneumatic_storage_reducer.synchronize_state_vector( @@ -279,20 +283,92 @@ class GenericFluidSystem: self.pneumatic_storage_reducer.synchronize_state_vector(values) ) + def _build_jacobian_sparsity(self): + """Build a conservative state dependency graph for implicit solvers. + + Two state entries are coupled when a physical path connects them without + crossing a third storage state. This over-approximates the local + pressure-flow/mechanical closure while preserving branch sparsity. + """ + from scipy.sparse import lil_matrix + + entries = self.mechanical_state_reducer.state_entries + entry_components: list[set[str]] = [] + entry_sizes: list[int] = [] + for entry in entries: + if isinstance(entry, MechanicalConstraintGroup): + entry_components.append( + {component.name for component in entry.components} + ) + entry_sizes.append(2) + else: + entry_components.append({entry.name}) + entry_sizes.append(entry.state_size) + + owner_by_component = { + component_name: entry_index + for entry_index, component_names in enumerate(entry_components) + for component_name in component_names + } + adjacency = {name: set() for name in self.network.components} + for connection in self.network.connections: + if connection.kind != "physical": + continue + first, second = connection.endpoints + adjacency[first.component].add(second.component) + adjacency[second.component].add(first.component) + + dependencies: list[set[int]] = [] + for entry_index, component_names in enumerate(entry_components): + visited = set(component_names) + pending = list(component_names) + found = {entry_index} + while pending: + current = pending.pop() + for neighbour in adjacency[current] - visited: + visited.add(neighbour) + neighbour_entry = owner_by_component.get(neighbour) + if ( + neighbour_entry is not None + and neighbour_entry != entry_index + ): + found.add(neighbour_entry) + else: + pending.append(neighbour) + dependencies.append(found) + + offsets = [0] + for state_size in entry_sizes: + offsets.append(offsets[-1] + state_size) + sparsity = lil_matrix( + (offsets[-1], offsets[-1]), + dtype=bool, + ) + for row_entry, column_entries in enumerate(dependencies): + for column_entry in column_entries: + sparsity[ + offsets[row_entry] : offsets[row_entry + 1], + offsets[column_entry] : offsets[column_entry + 1], + ] = True + return sparsity.tocsr() + + def jacobian_sparsity(self): + if self._jacobian_sparsity is None: + self._jacobian_sparsity = self._build_jacobian_sparsity() + return self._jacobian_sparsity + def _close_current_state(self, time: float) -> dict[str, dict[str, float]]: signal = self.signal_resolver.solve(time) self.signal_propagation_count += signal.propagated - for component in self.dynamic_components: - component.refresh_thermodynamic_ports() - algebraic = self.pressure_flow_solver.solve() - pressure_flow_solve_count = 1 + self.pressure_flow_solver.propagate_equal_efforts(("x", "v")) 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() - pressure_flow_solve_count += 1 + for component in self.dynamic_components: + component.refresh_thermodynamic_ports() + algebraic = self.pressure_flow_solver.solve( + effort_variables=("p",), + ) + pressure_flow_solve_count = 1 # Some constitutive flow laws recover their upstream temperature from # connected stream enthalpy, while junction stream mixing itself depends @@ -321,7 +397,11 @@ class GenericFluidSystem: component.update_flow_temperature_references( temperature_reference_h[component.name] ) - algebraic = self.pressure_flow_solver.solve() + algebraic = self.pressure_flow_solver.solve( + effort_variables=( + ("p",) if pressure_flow_solve_count == 0 else () + ), + ) pressure_flow_solve_count += 1 current_flows = tuple(port.m_flow for port in physical_ports) flow_scale = max( @@ -415,6 +495,14 @@ class GenericFluidSystem: progress_callback(last_reported_progress, phase) report_progress(0.0, "initializing", force=True) + integration_config = config + if isinstance(config.atol, (int, float)): + integration_config = replace( + config, + atol=self.mechanical_state_reducer.absolute_tolerances( + float(config.atol) + ), + ) t_eval = simulation_sample_times(config, sample_step) signal_event_times = self.signal_resolver.event_times( config.t_start, @@ -443,7 +531,7 @@ class GenericFluidSystem: solution = integrate_ode( rhs=monitored_rhs, initial_state=initial_state, - config=config, + config=integration_config, t_eval=t_eval, cancel_check=cancel_check, accepted_step_callback=( @@ -455,6 +543,11 @@ class GenericFluidSystem: if self.mechanical_state_reducer.has_state_events else None ), + jac_sparsity=( + self.jacobian_sparsity() + if integration_config.method in {"BDF", "Radau"} + else None + ), ) if isinstance(solution, ODESolution): run_status: SimulationRunStatus = solution.status diff --git a/tests/test_amesim_helium_medium.py b/tests/test_amesim_helium_medium.py index 803eb9c..6920e53 100644 --- a/tests/test_amesim_helium_medium.py +++ b/tests/test_amesim_helium_medium.py @@ -93,6 +93,21 @@ class AmesimHeliumPengRobinsonMediumTests(unittest.TestCase): delta=1.0e-8, ) + medium.temperature_from_pressure_enthalpy.cache_clear() + first = medium.temperature_from_pressure_enthalpy( + pressure, + transport_enthalpy, + ) + after_first = medium.temperature_from_pressure_enthalpy.cache_info() + second = medium.temperature_from_pressure_enthalpy( + pressure, + transport_enthalpy, + ) + after_second = medium.temperature_from_pressure_enthalpy.cache_info() + self.assertEqual(second, first) + self.assertEqual(after_first.misses, 1) + self.assertEqual(after_second.hits, 1) + density = medium.density(pressure, temperature) volume = 0.01 properties = medium.properties_from_mU( @@ -120,6 +135,28 @@ class AmesimHeliumPengRobinsonMediumTests(unittest.TestCase): with self.assertRaises(RecoverableTrialStateError): medium.properties_from_mU(-1.0, 1.0, 1.0) + def test_exact_repeated_real_gas_state_recovery_is_cached(self) -> None: + medium = AmesimHeliumPengRobinsonMedium() + pressure = 15.3e6 + temperature = 293.15 + volume = 0.057 + mass = medium.density(pressure, temperature) * volume + internal_energy = mass * medium.specific_internal_energy_at_pressure( + pressure, + temperature, + ) + + first = medium.properties_from_mU(mass, internal_energy, volume) + second = medium.properties_from_mU(mass, internal_energy, volume) + changed = medium.properties_from_mU( + mass, + internal_energy, + volume * 1.01, + ) + + self.assertIs(second, first) + self.assertIsNot(changed, first) + def test_amesim_2404_real_gas_isentropic_factor_reference(self) -> None: medium = AmesimHeliumPengRobinsonMedium() diff --git a/tests/test_amesim_mechanical_public_components.py b/tests/test_amesim_mechanical_public_components.py index fc4253a..1593c3c 100644 --- a/tests/test_amesim_mechanical_public_components.py +++ b/tests/test_amesim_mechanical_public_components.py @@ -207,10 +207,11 @@ class AmesimMechanicalPublicComponentTests(unittest.TestCase): residuals = node.pressure_flow_equation_residuals() self.assertEqual(node.active_ports, ("port_1", "port_2", "port_9")) - self.assertAlmostEqual(node.total_force, 0.0) + self.assertAlmostEqual(node.total_force, 10.0) + self.assertAlmostEqual(node.force_balance, 0.0) self.assertEqual(len(residuals), 5) self.assertTrue(all(abs(residual.value) <= 1.0e-12 for residual in residuals)) - self.assertEqual(node.component_result_values(), {"tforce": 0.0}) + self.assertEqual(node.component_result_values(), {"tforce": 10.0}) def test_lmechn1_registry_rejects_fractional_integer_options(self) -> None: with self.assertRaisesRegex(ValueError, "v1 must be an integer"): diff --git a/tests/test_amesim_mechanical_xml.py b/tests/test_amesim_mechanical_xml.py index b601b61..35eceb1 100644 --- a/tests/test_amesim_mechanical_xml.py +++ b/tests/test_amesim_mechanical_xml.py @@ -10,6 +10,7 @@ from app.main import ( reactflow_project_storage_data, run_system_xml_simulation, ) +from app.simulation.systems.generic import GenericFluidSystem from app.system_xml import validate_system_xml_document from tests.test_amesim_pnvo001_signal_xml import signal_edge, signal_port from tests.test_generic_system_xml_simulation import component_node, physical_edge @@ -287,6 +288,16 @@ def force_node_mass_project() -> ReactFlowProjectPayload: class AmesimMechanicalXmlTests(unittest.TestCase): + def test_mechanical_states_use_machine_scale_absolute_tolerances(self) -> None: + system = GenericFluidSystem( + compile_reactflow_network(zero_force_mass_project()) + ) + + tolerances = system.mechanical_state_reducer.absolute_tolerances(1.0e-8) + + self.assertEqual(tolerances, [1.0e-12, 1.0e-12]) + self.assertEqual(len(tolerances), len(system.initial_state_vector())) + def test_sparse_legacy_reactflow_mecmas_defaults_are_canonicalized( self, diff --git a/tests/test_amesim_pnl0002_pnl0003_component.py b/tests/test_amesim_pnl0002_pnl0003_component.py index 3eee9d6..10e4c07 100644 --- a/tests/test_amesim_pnl0002_pnl0003_component.py +++ b/tests/test_amesim_pnl0002_pnl0003_component.py @@ -1,6 +1,7 @@ from __future__ import annotations import unittest +from math import log10 from app.simulation.components.amesim.flow.pipes import AmesimPnl0002, AmesimPnl0003 from app.simulation.core.medium import IdealGasMedium @@ -45,6 +46,20 @@ class AmesimPnl0002ComponentTests(unittest.TestCase): self.assertGreater(forward, 0.0) self.assertLess(reverse, 0.0) + def test_exact_repeated_pipe_flow_inversion_is_cached(self) -> None: + pipe = AmesimPnl0002("pnl_2", self.medium, p0=100000.0, T0=293.15) + pipe._one_way_pn2pipefr_mass_flow.cache_clear() + + first = pipe.mass_flow(15.3e6, 10.0e6, 293.15) + after_first = pipe._one_way_pn2pipefr_mass_flow.cache_info() + second = pipe.mass_flow(15.3e6, 10.0e6, 293.15) + after_second = pipe._one_way_pn2pipefr_mass_flow.cache_info() + + self.assertEqual(second, first) + self.assertEqual(after_first.misses, 1) + self.assertEqual(after_second.misses, 1) + self.assertEqual(after_second.hits, 1) + def test_external_upstream_flow_uses_connected_stream_temperature(self) -> None: pipe = AmesimPnl0002("pnl_2", self.medium, p0=100000.0, T0=350.0) center = pipe.properties() @@ -78,6 +93,23 @@ class AmesimPnl0002ComponentTests(unittest.TestCase): self.assertAlmostEqual(flow, expected) self.assertNotAlmostEqual(flow, center_temperature_flow) + def test_fully_rough_limit_matches_nikuradse_colebrook_asymptote(self) -> None: + relative_roughness = 0.045 / 20.0 + pipe = AmesimPnl0002( + "pnl_2", + self.medium, + diam=0.02, + le=2.0, + rr=relative_roughness, + ) + expected = 1.0 / (-2.0 * log10(relative_roughness / 3.7)) ** 2 + + self.assertAlmostEqual( + pipe.friction_factor(1.0e12), + expected, + delta=1.0e-10, + ) + def test_missing_stream_cache_preserves_center_temperature_fallback(self) -> None: pipe = AmesimPnl0002("pnl_2", self.medium, p0=100000.0, T0=350.0) center = pipe.properties() diff --git a/tests/test_amesim_pnl00r_component.py b/tests/test_amesim_pnl00r_component.py index a4e33b8..65437c4 100644 --- a/tests/test_amesim_pnl00r_component.py +++ b/tests/test_amesim_pnl00r_component.py @@ -83,7 +83,7 @@ class AmesimPnl00rComponentTests(unittest.TestCase): self.assertAlmostEqual( pipe_20mm.friction_factor(56_887.5547), 0.0215740061043, - delta=8.0e-5, + delta=1.0e-4, ) self.assertAlmostEqual( pipe_20mm.friction_factor(700_686.41), diff --git a/tests/test_amesim_pnor001_component.py b/tests/test_amesim_pnor001_component.py index 8766b96..1d65c71 100644 --- a/tests/test_amesim_pnor001_component.py +++ b/tests/test_amesim_pnor001_component.py @@ -52,6 +52,22 @@ class AmesimPnor001ComponentTests(unittest.TestCase): self.assertEqual(residuals["pressure_flow_relation"].value, 0.0) self.assertEqual(set(orifice.component_result_values()), {"cm", "gasvel"}) + def test_exact_repeated_flow_characteristics_are_cached(self) -> None: + orifice = AmesimPnor001("pnor_1", self.medium) + orifice._one_way_flow_characteristics.cache_clear() + + first = orifice._one_way_flow_characteristics( + upstream_pressure=500000.0, + downstream_pressure=100000.0, + upstream_temperature=293.15, + ) + after_first = orifice._one_way_flow_characteristics.cache_info() + second = orifice._one_way_flow_characteristics( + upstream_pressure=500000.0, + downstream_pressure=100000.0, + upstream_temperature=293.15, + ) + def test_mass_flow_is_bidirectional_by_pressure_order(self) -> None: orifice = AmesimPnor001("pnor_1", self.medium) diff --git a/tests/test_component_catalog.py b/tests/test_component_catalog.py index e9ab7b6..e126c07 100644 --- a/tests/test_component_catalog.py +++ b/tests/test_component_catalog.py @@ -464,9 +464,11 @@ class ComponentCatalogTests(unittest.TestCase): with self.subTest(model_type=model_type): component = self.components[model_type] model_parameters = parameters(model_type) - expected_version = ( - "0.5.0" if model_type == "amesim_pnl0002" else "0.3.0" - ) + expected_version = { + "amesim_pnl0001": "0.4.0", + "amesim_pnl0002": "0.6.0", + "amesim_pnl0003": "0.4.0", + }[model_type] self.assertEqual(component["modelVersion"], expected_version) self.assertEqual(model_parameters["mode"]["editor"], "choice") self.assertEqual(model_parameters["mode"]["options"], thermal_options) diff --git a/tests/test_mechanical_solver_causalization.py b/tests/test_mechanical_solver_causalization.py index fb2e0f0..73ec7d3 100644 --- a/tests/test_mechanical_solver_causalization.py +++ b/tests/test_mechanical_solver_causalization.py @@ -488,6 +488,36 @@ class MechanicalSolverCausalizationTests(unittest.TestCase): for value in result.series["mass.a"]: self.assertAlmostEqual(value, 0.0, places=12) + def test_ideal_stop_ignores_zero_velocity_roundoff_inside_boundary_band(self) -> None: + system, _mass = _single_mass_system( + 0.0, + stoptype=1.0, + x0=0.0, + xmin=0.0, + xmax=1.0, + ) + reducer = system.mechanical_state_reducer + previous_state = [0.0, 2.0e-32] + current_state = [0.0, -1.0e-32] + + def dense_state(time: float) -> list[float]: + fraction = float(time) + return [ + 0.0, + previous_state[1] + + fraction * (current_state[1] - previous_state[1]), + ] + + transition = reducer.state_transition( + 0.0, + previous_state, + 1.0, + current_state, + dense_state, + ) + + self.assertIsNone(transition) + def test_ideal_upper_stop_projects_a_high_speed_impact(self) -> None: system, _mass = _single_mass_system( 100.0, diff --git a/tests/test_pnl0002_replay.py b/tests/test_pnl0002_replay.py new file mode 100644 index 0000000..4995cda --- /dev/null +++ b/tests/test_pnl0002_replay.py @@ -0,0 +1,108 @@ +from __future__ import annotations + +import unittest + +from app.simulation.components.amesim.flow.pipes import AmesimPnl0002 +from app.simulation.core.medium import IdealGasMedium +from app.simulation.reporting.amesim_results import AmesimResults +from app.simulation.reporting.pnl0002_replay import ( + Pnl0002ReplayPaths, + replay_pnl0002_amesim_states, +) + + +class Pnl0002ReplayTests(unittest.TestCase): + def setUp(self) -> None: + self.pipe = AmesimPnl0002( + "pneumatic_83", + IdealGasMedium(), + diam=0.02, + le=2.0, + rr=0.00225, + ) + self.paths = Pnl0002ReplayPaths( + center_pressure="center_p", + center_temperature="center_T", + port_1_pressure="port_1_p", + port_1_temperature="port_1_T", + port_1_mass_flow="port_1_dm", + port_2_pressure="port_2_p", + port_2_temperature="port_2_T", + port_2_mass_flow="port_2_dm", + reynolds="re", + friction_factor="ff", + ) + + def _matching_results(self) -> AmesimResults: + center_pressure = 300000.0 + center_temperature = 320.0 + port_1_pressure = 500000.0 + port_1_temperature = 330.0 + port_2_pressure = 100000.0 + port_2_temperature = 310.0 + flow_1 = self.pipe.mass_flow( + port_1_pressure, + center_pressure, + port_1_temperature, + ) + flow_2 = self.pipe.mass_flow( + port_2_pressure, + center_pressure, + center_temperature, + ) + reynolds_1 = self.pipe.reynolds_number(flow_1, port_1_temperature) + reynolds_2 = self.pipe.reynolds_number(flow_2, center_temperature) + friction = 0.5 * ( + self.pipe.friction_factor(reynolds_1) + + self.pipe.friction_factor(reynolds_2) + ) + return AmesimResults( + times=(0.0,), + variables=(), + saved_variable_indices=(), + series_by_data_path={ + "center_p": (center_pressure,), + "center_T": (center_temperature,), + "port_1_p": (port_1_pressure,), + "port_1_T": (port_1_temperature,), + "port_1_dm": (flow_1 / -1.0e-3,), + "port_2_p": (port_2_pressure,), + "port_2_T": (port_2_temperature,), + "port_2_dm": (flow_2 / -1.0e-3,), + "re": (0.5 * (reynolds_1 + reynolds_2),), + "ff": (friction,), + }, + final_values_by_data_path={}, + ) + + def test_matching_saved_states_replay_without_integration(self) -> None: + initial_state = self.pipe.get_state_vector() + + report = replay_pnl0002_amesim_states( + self.pipe, + self._matching_results(), + self.paths, + ) + + self.assertEqual(report["mode"], "amesim-state-replay-no-integration") + self.assertEqual(report["pointCount"], 1) + self.assertEqual(self.pipe.get_state_vector(), initial_state) + row = report["rows"][0] + self.assertAlmostEqual(row["port1MassFlowError"], 0.0, places=15) + self.assertAlmostEqual(row["port2MassFlowError"], 0.0, places=15) + self.assertAlmostEqual(row["reynoldsError"], 0.0, places=12) + self.assertAlmostEqual(row["frictionFactorError"], 0.0, places=15) + + def test_missing_saved_variable_has_actionable_error(self) -> None: + results = self._matching_results() + results.series_by_data_path.pop("ff") + + with self.assertRaisesRegex( + ValueError, + "AMESim replay variable is not saved: ff", + ): + replay_pnl0002_amesim_states(self.pipe, results, self.paths) + + +if __name__ == "__main__": + unittest.main()