from __future__ import annotations from dataclasses import dataclass from math import floor, isfinite from PythonModels.core.algebraic import PressureFlowSolver from PythonModels.core.base import DynamicComponent from PythonModels.core.network import Endpoint, SimulationNetwork from PythonModels.core.solver import SolveIVPConfig, integrate_ode from PythonModels.core.stream import StreamResolver @dataclass(frozen=True) class SimulationPreparationIssue: code: str message: str def as_dict(self) -> dict[str, str]: return {"code": self.code, "message": self.message} class SimulationPreparationError(ValueError): def __init__(self, issues: tuple[SimulationPreparationIssue, ...]) -> None: super().__init__("The compiled model is not ready for simulation.") self.issues = issues @dataclass(frozen=True) class GenericSimulationResult: success: bool message: str series: dict[str, list[float]] final: dict[str, float] diagnostics: dict[str, object] def as_dict(self) -> dict[str, object]: return { "success": self.success, "message": self.message, "series": self.series, "final": self.final, "diagnostics": self.diagnostics, } class _UnionFind: def __init__(self, items: set[Endpoint]) -> None: self.parent = {item: item for item in items} def find(self, item: Endpoint) -> Endpoint: parent = self.parent[item] if parent != item: self.parent[item] = self.find(parent) return self.parent[item] def union(self, first: Endpoint, second: Endpoint) -> None: first_root = self.find(first) second_root = self.find(second) if first_root != second_root: self.parent[second_root] = first_root def _equation_port(component_name: str, variable: str) -> Endpoint | None: parts = variable.rsplit(".", 2) if len(parts) != 3: return None prefix, port_name, variable_name = parts if prefix != component_name or variable_name != "p": return None return Endpoint(component_name, port_name) def simulation_preparation_issues( network: SimulationNetwork, ) -> tuple[SimulationPreparationIssue, ...]: issues: list[SimulationPreparationIssue] = [] physical_endpoints = { Endpoint(component.name, definition.name) for component in network.components.values() for definition in component.port_definitions if definition.kind == "physical" } connected_endpoints = { endpoint for connection in network.connections if connection.kind == "physical" for endpoint in connection.endpoints } for endpoint in sorted(physical_endpoints - connected_endpoints, key=str): issues.append( SimulationPreparationIssue( "PORT_UNCONNECTED", f"Physical port {endpoint} must be connected before simulation.", ) ) if any( definition.kind == "signal" for component in network.components.values() for definition in component.port_definitions ): issues.append( SimulationPreparationIssue( "SIGNAL_PORT_UNSUPPORTED", "Signal-port simulation is not implemented in the current MVP solver.", ) ) structure = network.pressure_flow_structure_dict() if not structure["isSquare"]: issues.append( SimulationPreparationIssue( "PRESSURE_FLOW_SYSTEM_NOT_SQUARE", "Pressure-flow equation count does not match the unknown count: " f"{structure['equationCount']} equations for {structure['unknownCount']} unknowns.", ) ) dynamic_names = { component.name for component in network.components.values() if isinstance(component, DynamicComponent) } if not dynamic_names: issues.append( SimulationPreparationIssue( "DYNAMIC_STATE_MISSING", "Each simulated network requires at least one storage component.", ) ) adjacency = {name: set() for name in network.components} for connection in network.connections: first, second = connection.endpoints adjacency[first.component].add(second.component) adjacency[second.component].add(first.component) remaining = set(adjacency) while remaining: start = remaining.pop() group = {start} stack = [start] while stack: current = stack.pop() for neighbour in adjacency[current] - group: group.add(neighbour) remaining.discard(neighbour) stack.append(neighbour) if not (group & dynamic_names): issues.append( SimulationPreparationIssue( "ALGEBRAIC_ISLAND_HAS_NO_STORAGE", "A connected physical network has no pressure/enthalpy storage anchor: " + ", ".join(sorted(group)) + ".", ) ) if physical_endpoints: effort_groups = _UnionFind(physical_endpoints) for connection in network.connections: if connection.kind == "physical": effort_groups.union(*connection.endpoints) storage_ports: dict[Endpoint, str] = {} for component in network.components.values(): for equation in component.pressure_flow_equation_residuals(): pressure_ports = [ endpoint for variable in equation.variables if (endpoint := _equation_port(component.name, variable)) is not None ] if equation.relation == "equal" and len(pressure_ports) == 2: effort_groups.union(pressure_ports[0], pressure_ports[1]) if equation.relation == "state": for endpoint in pressure_ports: storage_ports[endpoint] = component.name storages_by_group: dict[Endpoint, set[str]] = {} for endpoint, component_name in storage_ports.items(): storages_by_group.setdefault(effort_groups.find(endpoint), set()).add( component_name ) for storage_names in storages_by_group.values(): if len(storage_names) > 1: issues.append( SimulationPreparationIssue( "IDEAL_STORAGE_COUPLING_UNSUPPORTED", "Storage components are connected without a resistance: " + ", ".join(sorted(storage_names)) + ". Insert an orifice or pipe between them.", ) ) return tuple(issues) def simulation_sample_times( config: SolveIVPConfig, step: float, *, max_points: int = 10001, ) -> list[float]: if step <= 0.0 or not isfinite(step): raise ValueError("Simulation sample step must be finite and greater than zero.") duration = config.t_stop - config.t_start if duration <= 0.0: raise ValueError("Simulation stop time must be greater than start time.") interval_count = int(floor(duration / step + 1e-12)) times = [config.t_start + index * step for index in range(interval_count + 1)] if times[-1] < config.t_stop - 1e-12: times.append(config.t_stop) else: times[-1] = config.t_stop if len(times) > max_points: raise ValueError( f"Simulation requests {len(times)} samples; the limit is {max_points}." ) return times class GenericFluidSystem: """Topology-driven, semi-explicit fluid simulation for registered components.""" def __init__(self, network: SimulationNetwork) -> None: issues = simulation_preparation_issues(network) if issues: raise SimulationPreparationError(issues) self.network = network self.dynamic_components = network.dynamic_components() self.pressure_flow_solver = PressureFlowSolver(network) self.stream_resolver = StreamResolver(network) self.algebraic_solve_count = 0 self.max_algebraic_residual = 0.0 self.max_algebraic_evaluations = 0 self.max_stream_iterations = 0 def initial_state_vector(self) -> list[float]: return self.network.initial_state_vector() def apply_state_vector(self, values: list[float]) -> None: self.network.apply_state_vector(values) def _close_current_state(self) -> dict[str, dict[str, float]]: for component in self.dynamic_components: component.refresh_thermodynamic_ports() algebraic = self.pressure_flow_solver.solve() stream, connected_h = self.stream_resolver.solve() self.algebraic_solve_count += 1 self.max_algebraic_residual = max( self.max_algebraic_residual, algebraic.max_scaled_residual, ) self.max_algebraic_evaluations = max( self.max_algebraic_evaluations, algebraic.evaluations, ) self.max_stream_iterations = max( self.max_stream_iterations, stream.iterations, ) return connected_h def consistent_initial_state_vector(self) -> list[float]: state = self.initial_state_vector() self.apply_state_vector(state) self._close_current_state() return state def rhs(self, _time: float, state_vector: list[float]) -> list[float]: self.apply_state_vector(state_vector) connected_h = self._close_current_state() derivatives: list[float] = [] for component in self.dynamic_components: derivatives.extend( component.state_derivative_from_ports(connected_h[component.name]) ) return derivatives def _append_current_state(self, series: dict[str, list[float]]) -> None: for component in self.dynamic_components: properties = component.refresh_thermodynamic_ports() state = component.get_state_vector() if len(state) >= 2: series.setdefault(f"{component.name}.m", []).append(float(state[0])) series.setdefault(f"{component.name}.U", []).append(float(state[1])) for name in ("p", "T", "rho", "u", "h"): if hasattr(properties, name): series.setdefault(f"{component.name}.{name}", []).append( float(getattr(properties, name)) ) for component in self.network.components.values(): for port_name, port in component.ports.items(): prefix = f"{component.name}.{port_name}" series.setdefault(f"{prefix}.p", []).append(float(port.p)) series.setdefault(f"{prefix}.m_flow", []).append(float(port.m_flow)) series.setdefault(f"{prefix}.h_outflow", []).append( float(port.h_outflow) ) def simulate( self, config: SolveIVPConfig, *, sample_step: float, ) -> GenericSimulationResult: t_eval = simulation_sample_times(config, sample_step) initial_state = self.consistent_initial_state_vector() solution = integrate_ode( rhs=self.rhs, initial_state=initial_state, config=config, t_eval=t_eval, ) times = [float(value) for value in solution.t] series: dict[str, list[float]] = {"time": times} for time_index in range(len(times)): state = [ float(solution.y[state_index][time_index]) for state_index in range(len(solution.y)) ] self.apply_state_vector(state) self._close_current_state() self._append_current_state(series) final = { key: values[-1] for key, values in series.items() if key != "time" and values } diagnostics = { "pressureFlow": { "solveCount": self.algebraic_solve_count, "maxScaledResidual": self.max_algebraic_residual, "maxEvaluationsPerSolve": self.max_algebraic_evaluations, "last": ( self.pressure_flow_solver.last_diagnostics.as_dict() if self.pressure_flow_solver.last_diagnostics is not None else None ), }, "stream": { "maxIterationsPerSolve": self.max_stream_iterations, "last": ( self.stream_resolver.last_diagnostics.as_dict() if self.stream_resolver.last_diagnostics is not None else None ), }, "stateCount": len(initial_state), "sampleCount": len(times), } return GenericSimulationResult( success=bool(solution.success), message=str(solution.message), series=series, final=final, diagnostics=diagnostics, )