from __future__ import annotations from collections.abc import Callable 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.performance import performance_span, profile_phase from app.simulation.solvers.algebraic import PressureFlowSolver from app.simulation.solvers.mechanical import ( MechanicalConstraintGroup, MechanicalStateReducer, ) from app.simulation.solvers.pneumatic_storage import ( IdealPneumaticStorageReducer, ideal_storage_group_is_reducible, ) from app.simulation.solvers.pneumatic_volume import PneumaticVolumeResolver from app.simulation.solvers.solver import ODESolution, SolveIVPConfig, integrate_ode from app.simulation.solvers.signal import SignalResolver from app.simulation.solvers.stream import StreamResolver from app.simulation.systems.network import Endpoint, SimulationNetwork SimulationProgressCallback = Callable[[float, str], None] SimulationCancellationCheck = Callable[[], bool] SimulationRunStatus = Literal["completed", "cancelled", "failed"] @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 class ThermofluidClosureError(RuntimeError): """Raised when stream enthalpy and pressure-flow do not reach one fixed point.""" class SimulationSampleTimeError(ValueError): """Stable failure contract for an unsafe or unrepresentable sample grid.""" def __init__(self, code: str, message: str) -> None: super().__init__(message) self.code = code @dataclass(frozen=True) class GenericSimulationResult: success: bool status: SimulationRunStatus message: str simulated_until: float requested_stop_time: float variables: tuple[ResultVariableMetadata, ...] series: dict[str, list[float]] final: dict[str, float] diagnostics: dict[str, object] def as_dict(self) -> dict[str, object]: return { "success": self.success, "status": self.status, "partial": self.status != "completed", "message": self.message, "simulatedUntil": self.simulated_until, "requestedStopTime": self.requested_stop_time, "variables": [variable.as_dict() for variable in self.variables], "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, port_name) for component in network.components.values() for port_name in component.required_connection_ports } 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.", ) ) 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.", ) ) physical_component_names = { component.name for component in network.components.values() if any( definition.kind == "physical" for definition in component.active_port_definitions ) } adjacency = {name: set() for name in physical_component_names} for connection in network.connections: if connection.kind != "physical": continue 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, dict[Endpoint, str]] = {} for endpoint, component_name in storage_ports.items(): storages_by_group.setdefault(effort_groups.find(endpoint), {})[ endpoint ] = component_name for storage_endpoints in storages_by_group.values(): storage_names = set(storage_endpoints.values()) if len(storage_names) > 1: if ideal_storage_group_is_reducible( network, storage_endpoints, ): continue 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 max_points < 2: raise SimulationSampleTimeError( "SIMULATION_SAMPLE_LIMIT_INVALID", "Simulation sample limit must allow at least two points.", ) t_start = float(config.t_start) t_stop = float(config.t_stop) if not isfinite(t_start) or not isfinite(t_stop): raise SimulationSampleTimeError( "SIMULATION_VALUE_NOT_FINITE", "Simulation start and stop times must be finite.", ) if step <= 0.0 or not isfinite(step): raise SimulationSampleTimeError( "SIMULATION_SAMPLE_STEP_INVALID", "Simulation sample step must be finite and greater than zero.", ) duration = t_stop - t_start if not isfinite(duration): raise SimulationSampleTimeError( "SIMULATION_TIME_SPAN_NOT_FINITE", "Simulation time span must be finite.", ) if duration <= 0.0: raise SimulationSampleTimeError( "SIMULATION_TIME_RANGE_INVALID", "Simulation stop time must be greater than start time.", ) # Bound the grid before dividing by a potentially tiny step or allocating # the result list. This avoids both float-to-int overflow and an OOM-sized # ``range``/list when input comes from an external System XML document. maximum_interval_count = max_points - 1 if step < duration / maximum_interval_count: raise SimulationSampleTimeError( "SIMULATION_SAMPLE_COUNT_EXCEEDED", f"Simulation sample count exceeds the limit of {max_points}; " "increase sampleStep.", ) ratio = duration / step if not isfinite(ratio): raise SimulationSampleTimeError( "SIMULATION_SAMPLE_COUNT_EXCEEDED", f"Simulation sample count exceeds the limit of {max_points}; " "increase sampleStep.", ) interval_count = int(floor(ratio)) last_regular_time = t_start + interval_count * step append_stop = last_regular_time < t_stop requested_point_count = interval_count + 1 + int(append_stop) if requested_point_count > max_points: raise SimulationSampleTimeError( "SIMULATION_SAMPLE_COUNT_EXCEEDED", f"Simulation requests {requested_point_count} samples; " f"the limit is {max_points}.", ) times = [t_start] for index in range(1, interval_count + 1): candidate = t_start + index * step if not isfinite(candidate): raise SimulationSampleTimeError( "SIMULATION_SAMPLE_TIME_UNREPRESENTABLE", "Simulation sampleStep cannot be represented over the requested " "absolute time range.", ) if candidate >= t_stop: candidate = t_stop if candidate <= times[-1]: raise SimulationSampleTimeError( "SIMULATION_SAMPLE_TIME_UNREPRESENTABLE", "Simulation sampleStep is too small to advance floating-point " "time over the requested absolute time range.", ) times.append(candidate) if candidate == t_stop: break if times[-1] < t_stop: times.append(t_stop) if len(times) < 2 or any( current >= following for current, following in zip(times, times[1:]) ): raise SimulationSampleTimeError( "SIMULATION_SAMPLE_TIME_UNREPRESENTABLE", "Simulation sample times must contain at least two strictly " "increasing values.", ) return times class GenericFluidSystem: """Topology-driven, semi-explicit fluid simulation for registered components.""" @profile_phase("simulation.system_construction") 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.mechanical_state_reducer = MechanicalStateReducer( network, self.dynamic_components, ) self.pneumatic_storage_reducer = IdealPneumaticStorageReducer( network, self.mechanical_state_reducer, ) self.pressure_flow_solver = PressureFlowSolver(network) self.pneumatic_volume_resolver = PneumaticVolumeResolver(network) self.signal_resolver = SignalResolver(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 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( self.mechanical_state_reducer.initial_state_vector(), validate=True, ) def apply_state_vector(self, values: list[float]) -> None: self.mechanical_state_reducer.apply_state_vector( 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 jacobian_sparsity_diagnostics(self) -> dict[str, float | int]: from scipy.optimize._numdiff import group_columns sparsity = self.jacobian_sparsity() group_count = int(group_columns(sparsity).max(initial=-1)) + 1 state_count = int(sparsity.shape[0]) return { "nonzeroCount": int(sparsity.nnz), "density": ( float(sparsity.nnz) / float(state_count * state_count) if state_count else 0.0 ), "colorGroupCount": group_count, } @profile_phase( "simulation.closure", minimum_mode="audit", reset_property_shadow=True, ) def _close_current_state(self, time: float) -> dict[str, dict[str, float]]: signal = self.signal_resolver.solve(time) self.signal_propagation_count += signal.propagated self.pressure_flow_solver.propagate_equal_efforts(("x", "v")) pneumatic_volume = self.pneumatic_volume_resolver.solve() self.pneumatic_volume_propagation_count += pneumatic_volume.propagated self._refresh_dynamic_components() 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 # 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.active_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() temperature_reference_h = ( self.stream_resolver.connected_temperature_reference_enthalpies() ) for component in self.dynamic_components: component.update_stream_outflows(connected_h[component.name]) component.update_flow_temperature_references( temperature_reference_h[component.name] ) 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( [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 += pressure_flow_solve_count 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 @profile_phase("simulation.refresh", minimum_mode="audit") def _refresh_dynamic_components(self) -> None: for component in self.dynamic_components: component.refresh_thermodynamic_ports() @profile_phase("simulation.derivatives", minimum_mode="audit") def _state_derivatives( self, connected_h: dict[str, dict[str, float]], ) -> list[float]: return self.pneumatic_storage_reducer.coupled_derivatives( self.mechanical_state_reducer.state_derivatives(connected_h) ) def consistent_initial_state_vector(self, time: float = 0.0) -> list[float]: state = self.initial_state_vector() self.apply_state_vector(state) self._close_current_state(time) return state @profile_phase("simulation.rhs", minimum_mode="audit") def rhs(self, _time: float, state_vector: list[float]) -> list[float]: self.apply_state_vector(state_vector) connected_h = self._close_current_state(_time) return self._state_derivatives(connected_h) def _append_current_state(self, series: dict[str, list[float]]) -> None: for component in self.network.components.values(): for relative_key, value in component.result_values().items(): series.setdefault( f"{component.name}.{relative_key}", [] ).append(value) def simulate( self, config: SolveIVPConfig, *, sample_step: float, progress_callback: SimulationProgressCallback | None = None, cancel_check: SimulationCancellationCheck | None = None, ) -> GenericSimulationResult: last_reported_progress = -1.0 last_reported_phase = "" def report_progress( progress: float, phase: str, *, force: bool = False, ) -> None: nonlocal last_reported_phase, last_reported_progress if progress_callback is None: return bounded_progress = min(1.0, max(0.0, progress)) if ( force or phase != last_reported_phase or bounded_progress - last_reported_progress >= 0.0025 ): last_reported_phase = phase last_reported_progress = max( last_reported_progress, bounded_progress, ) progress_callback(last_reported_progress, phase) report_progress(0.0, "initializing", force=True) with performance_span("simulation.sample_initialization"): 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, config.t_stop, ) initial_state = self.consistent_initial_state_vector(config.t_start) jac_sparsity = ( self.jacobian_sparsity() if integration_config.method in {"BDF", "Radau"} else None ) report_progress(0.0, "integrating", force=True) duration = config.t_stop - config.t_start furthest_solver_time = config.t_start def report_solver_time(time: float) -> None: nonlocal furthest_solver_time furthest_solver_time = max(furthest_solver_time, float(time)) time_fraction = ( (furthest_solver_time - config.t_start) / duration if duration > 0.0 else 1.0 ) report_progress(time_fraction, "integrating") def monitored_rhs(time: float, state_vector: list[float]) -> list[float]: if cancel_check is None: report_solver_time(time) return self.rhs(time, state_vector) solution = integrate_ode( rhs=monitored_rhs, initial_state=initial_state, config=integration_config, t_eval=t_eval, cancel_check=cancel_check, accepted_step_callback=( report_solver_time if cancel_check is not None else None ), breakpoints=signal_event_times, state_transition_handler=( self.mechanical_state_reducer.state_transition if self.mechanical_state_reducer.has_state_events else None ), jac_sparsity=jac_sparsity, ) if isinstance(solution, ODESolution): run_status: SimulationRunStatus = solution.status integration_error = solution.error else: run_status = "completed" if bool(solution.success) else "failed" integration_error = None result_message = str(solution.message) postprocess_progress = ( 1.0 if run_status == "completed" else max(0.0, last_reported_progress) ) report_progress(postprocess_progress, "postprocessing", force=True) times = [float(value) for value in solution.t] if isinstance(solution, ODESolution): solver_segment_diagnostics = [ segment.as_dict() for segment in solution.solver_segments ] else: solver_segment_diagnostics = [ { "startTime": float(config.t_start), "requestedStopTime": float(config.t_stop), "simulatedUntil": times[-1] if times else float(config.t_start), "nfev": int(getattr(solution, "nfev", 0)), "njev": int(getattr(solution, "njev", 0)), "nlu": int(getattr(solution, "nlu", 0)), "acceptedStepCount": 0, "solverStartCount": 1, "stateTransitionCount": 0, "recoverableRetryCount": 0, } ] solver_total_keys = ( "nfev", "njev", "nlu", "acceptedStepCount", "solverStartCount", "stateTransitionCount", "recoverableRetryCount", ) solver_totals = { key: sum(int(segment[key]) for segment in solver_segment_diagnostics) for key in solver_total_keys } jacobian_diagnostics = ( self.jacobian_sparsity_diagnostics() if integration_config.method in {"BDF", "Radau"} else None ) if jacobian_diagnostics is not None: color_group_count = int(jacobian_diagnostics["colorGroupCount"]) for segment in solver_segment_diagnostics: segment["finiteDifferenceRhsEstimate"] = ( int(segment["njev"]) * color_group_count ) solver_totals["finiteDifferenceRhsEstimate"] = sum( int(segment["finiteDifferenceRhsEstimate"]) for segment in solver_segment_diagnostics ) with performance_span("simulation.postprocessing"): series: dict[str, list[float]] = {"time": []} postprocessing_error: Exception | None = None self.mechanical_state_reducer.reset_constraint_modes() for time_index in range(len(times)): if ( run_status == "completed" and cancel_check is not None and cancel_check() ): run_status = "cancelled" result_message = ( "Simulation was stopped while preparing partial results." ) break state = [ float(solution.y[state_index][time_index]) for state_index in range(len(solution.y)) ] try: self.apply_state_vector(state) self._close_current_state(times[time_index]) self._append_current_state(series) series["time"].append(times[time_index]) except Exception as exc: run_status = "failed" result_message = str(exc) postprocessing_error = exc break if len(series["time"]) < 2: if postprocessing_error is not None: raise postprocessing_error if integration_error is not None: raise integration_error with performance_span("simulation.result_assembly"): final = { key: values[-1] for key, values in series.items() if key != "time" and values } diagnostics = { "integration": { "method": integration_config.method, "jacobianSparsity": jacobian_diagnostics, "segmentCount": len(solver_segment_diagnostics), "segments": solver_segment_diagnostics, "totals": solver_totals, }, "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, "maxThermofluidIterations": self.max_thermofluid_iterations, "last": ( self.stream_resolver.last_diagnostics.as_dict() if self.stream_resolver.last_diagnostics is not None else None ), }, "signal": { "propagations": self.signal_propagation_count, "eventTimes": list(signal_event_times), "last": ( self.signal_resolver.last_diagnostics.as_dict() if self.signal_resolver.last_diagnostics is not None else None ), }, "pneumaticVolume": { "propagations": self.pneumatic_volume_propagation_count, "last": ( self.pneumatic_volume_resolver.last_diagnostics.as_dict() if self.pneumatic_volume_resolver.last_diagnostics is not None else None ), }, "stateCount": len(initial_state), "sampleCount": len(series["time"]), } variables = tuple( variable for variable in self.network.result_variable_metadata() if variable.key in series ) report_progress( 1.0 if run_status == "completed" else max(0.0, last_reported_progress), "complete" if run_status == "completed" else run_status, force=True, ) return GenericSimulationResult( success=run_status == "completed" and bool(solution.success), status=run_status, message=result_message, simulated_until=( float(series["time"][-1]) if series["time"] else float(config.t_start) ), requested_stop_time=float(config.t_stop), variables=variables, series=series, final=final, diagnostics=diagnostics, )