from __future__ import annotations from dataclasses import dataclass, field from typing import Callable from PythonModels.components.cylinder import Cylinder from PythonModels.components.orifice import Orifice from PythonModels.components.pipe import Pipe from PythonModels.components.tank import Tank from PythonModels.components.tee import Tee from PythonModels.core.medium import IdealGasMedium, ThermodynamicProperties from PythonModels.core.state import VolumeState @dataclass(frozen=True) class BranchInletFlowDiagnostics: converged: bool iterations: int residual: float m_flow: float inlet_pressure: float @dataclass(frozen=True) class DownstreamPressureDiagnostics: converged: bool iterations: int residual: float pressure: float target_total_internal_energy: float @dataclass(frozen=True) class TestModelSolveDiagnostics: upper_branch_inlet: BranchInletFlowDiagnostics lower_branch_inlet: BranchInletFlowDiagnostics downstream_pressure_projection: DownstreamPressureDiagnostics | None @dataclass(frozen=True) class BranchClosureComponents: name: str orifice: Orifice pipe: Pipe @dataclass(frozen=True) class BranchClosureState: name: str pipe: ThermodynamicProperties inlet_flow: float outlet_flow: float inlet_h: float inlet_flow_diagnostics: BranchInletFlowDiagnostics @dataclass(frozen=True) class BranchSnapshot: name: str pipe: ThermodynamicProperties inlet_flow: float outlet_flow: float inlet_h: float inlet_flow_diagnostics: BranchInletFlowDiagnostics @dataclass(frozen=True) class TestModelSnapshot: cylinder: ThermodynamicProperties tank: ThermodynamicProperties tee_upstream_h: float tee_downstream_h: float branches: tuple[BranchSnapshot, ...] = field(default_factory=tuple) solve_diagnostics: TestModelSolveDiagnostics | None = None @property def pipe_upper(self) -> ThermodynamicProperties: return self.branches[0].pipe @property def pipe_lower(self) -> ThermodynamicProperties: return self.branches[1].pipe @property def branch_inlet_flows(self) -> tuple[float, ...]: return tuple(branch.inlet_flow for branch in self.branches) @property def branch_outlet_flows(self) -> tuple[float, ...]: return tuple(branch.outlet_flow for branch in self.branches) @dataclass(frozen=True) class InitializationDiagnostics: converged: bool iterations: int max_state_delta: float max_flow_delta: float max_enthalpy_delta: float downstream_pressure_spread: float state_vector: tuple[float, ...] @dataclass(frozen=True) class TestModelClosureComponents: cylinder: Cylinder upstream_tee: Tee upper_branch: BranchClosureComponents lower_branch: BranchClosureComponents downstream_tee: Tee tank: Tank def branches(self) -> tuple[BranchClosureComponents, BranchClosureComponents]: return (self.upper_branch, self.lower_branch) class TestModelClosure: """Owns Testmodel-specific closure, projection and port-writeback logic.""" def __init__( self, *, medium: IdealGasMedium, components: TestModelClosureComponents, initial_state_vector: Callable[[], list[float]], apply_state_vector: Callable[[list[float]], None], ) -> None: self.medium = medium self.components = components self._initial_state_vector = initial_state_vector self._apply_state_vector = apply_state_vector self.last_solve_diagnostics: TestModelSolveDiagnostics | None = None self.last_downstream_pressure_diagnostics: DownstreamPressureDiagnostics | None = None @staticmethod def _downstream_pressure_spread(snapshot: TestModelSnapshot) -> float: downstream_pressures = tuple(branch.pipe.p for branch in snapshot.branches) + ( snapshot.tank.p, ) return max(downstream_pressures) - min(downstream_pressures) @staticmethod def _initialization_flow_delta( previous_snapshot: TestModelSnapshot | None, current_snapshot: TestModelSnapshot, ) -> float: if previous_snapshot is None: return max(abs(branch.outlet_flow) for branch in current_snapshot.branches) return max( abs(curr - prev) for curr, prev in zip( current_snapshot.branch_outlet_flows, previous_snapshot.branch_outlet_flows, ) ) @staticmethod def _initialization_enthalpy_delta( previous_snapshot: TestModelSnapshot | None, current_snapshot: TestModelSnapshot, ) -> float: if previous_snapshot is None: return abs(current_snapshot.tee_downstream_h - current_snapshot.tank.h) return max( abs(current_snapshot.tee_upstream_h - previous_snapshot.tee_upstream_h), abs(current_snapshot.tee_downstream_h - previous_snapshot.tee_downstream_h), ) def consistent_initial_state_vector(self) -> list[float]: return list(self.initialize_consistent_state().state_vector) def initialize_consistent_state( self, max_iterations: int = 12, state_tolerance: float = 1e-9, flow_tolerance: float = 1e-9, enthalpy_tolerance: float = 1e-6, pressure_tolerance: float = 1e-6, strict_internal_solvers: bool = False, ) -> InitializationDiagnostics: raw_state = self._initial_state_vector() previous_snapshot: TestModelSnapshot | None = None diagnostics: InitializationDiagnostics | None = None for iteration in range(1, max_iterations + 1): state_before_projection = self._initial_state_vector() self.snapshot(state_before_projection, strict=strict_internal_solvers) self.project_downstream_pressure_constraints(strict=strict_internal_solvers) state_after_projection = self._initial_state_vector() snapshot_after_projection = self.snapshot( state_after_projection, strict=strict_internal_solvers, ) max_state_delta = max( abs(after - before) for before, after in zip(state_before_projection, state_after_projection) ) max_flow_delta = self._initialization_flow_delta( previous_snapshot, snapshot_after_projection, ) max_enthalpy_delta = self._initialization_enthalpy_delta( previous_snapshot, snapshot_after_projection, ) downstream_pressure_spread = self._downstream_pressure_spread( snapshot_after_projection, ) diagnostics = InitializationDiagnostics( converged=( max_state_delta <= state_tolerance and max_flow_delta <= flow_tolerance and max_enthalpy_delta <= enthalpy_tolerance and downstream_pressure_spread <= pressure_tolerance ), iterations=iteration, max_state_delta=max_state_delta, max_flow_delta=max_flow_delta, max_enthalpy_delta=max_enthalpy_delta, downstream_pressure_spread=downstream_pressure_spread, state_vector=tuple(state_after_projection), ) previous_snapshot = snapshot_after_projection if diagnostics.converged: self._apply_state_vector(raw_state) return diagnostics assert diagnostics is not None self._apply_state_vector(raw_state) return diagnostics def _solve_branch_inlet_flow( self, orifice: Orifice, pipe: Pipe, p_upstream: float, pipe_props: ThermodynamicProperties, *, strict: bool = False, ) -> tuple[float, BranchInletFlowDiagnostics]: m_flow = orifice.mass_flow(p_upstream, pipe_props.p) rho = max(pipe_props.rho, 1e-9) p_inlet = pipe.inlet_pressure(m_flow, rho, pipe_props.p) residual = abs(orifice.mass_flow(p_upstream, p_inlet) - m_flow) converged = False iterations = 0 for iteration in range(1, 9): p_inlet = pipe.inlet_pressure(m_flow, rho, pipe_props.p) next_m_flow = orifice.mass_flow(p_upstream, p_inlet) residual = abs(next_m_flow - m_flow) iterations = iteration if residual <= 1e-9 * max(1.0, abs(next_m_flow)): m_flow = next_m_flow converged = True break m_flow = next_m_flow diagnostics = BranchInletFlowDiagnostics( converged=converged, iterations=iterations, residual=residual, m_flow=m_flow, inlet_pressure=p_inlet, ) if strict and not diagnostics.converged: raise RuntimeError( f"Branch inlet flow solve did not converge for {pipe.name}: residual={residual:.6e}" ) return m_flow, diagnostics def _solve_downstream_branch_flows( self, cylinder: ThermodynamicProperties, tank: ThermodynamicProperties, branch_states: tuple[BranchClosureState, BranchClosureState], ) -> tuple[float, float]: return self._solve_downstream_branch_flows_from_state( inlet_h_upper=branch_states[0].inlet_h, inlet_h_lower=branch_states[1].inlet_h, pipe_upper_h=max(branch_states[0].pipe.h, 1e-9), pipe_lower_h=max(branch_states[1].pipe.h, 1e-9), tank_h=max(tank.h, 1e-9), q_in_upper=branch_states[0].inlet_flow, q_in_lower=branch_states[1].inlet_flow, ) def _project_volume_energy_to_pressure( self, component: Pipe | Tank, target_pressure: float, ) -> None: target_temperature = target_pressure * component.V / ( max(component.state.m, 1e-12) * self.medium.R_gas ) target_internal_energy = ( component.state.m * self.medium.specific_internal_energy(target_temperature) ) component.state = VolumeState(m=component.state.m, U=target_internal_energy) def _downstream_total_internal_energy_for_pressure( self, target_pressure: float, downstream_components: tuple[Pipe | Tank, ...], ) -> float: total_internal_energy = 0.0 for component in downstream_components: target_temperature = target_pressure * component.V / ( max(component.state.m, 1e-12) * self.medium.R_gas ) total_internal_energy += ( component.state.m * self.medium.specific_internal_energy(target_temperature) ) return total_internal_energy def _solve_downstream_common_pressure( self, downstream_components: tuple[Pipe | Tank, ...], target_total_internal_energy: float, *, strict: bool = False, ) -> tuple[float, DownstreamPressureDiagnostics]: lower_pressure = 1.0 upper_pressure = max(component.properties().p for component in downstream_components) upper_pressure = max(upper_pressure, 1e5) def residual(pressure: float) -> float: return ( self._downstream_total_internal_energy_for_pressure( pressure, downstream_components, ) - target_total_internal_energy ) upper_residual = residual(upper_pressure) iteration_count = 0 while upper_residual < 0.0: upper_pressure *= 2.0 upper_residual = residual(upper_pressure) final_pressure = 0.5 * (lower_pressure + upper_pressure) final_residual = residual(final_pressure) converged = False for iteration in range(1, 81): middle_pressure = 0.5 * (lower_pressure + upper_pressure) middle_residual = residual(middle_pressure) iteration_count = iteration final_pressure = middle_pressure final_residual = middle_residual if abs(middle_residual) <= 1e-12 * max(1.0, target_total_internal_energy): converged = True break if middle_residual > 0.0: upper_pressure = middle_pressure else: lower_pressure = middle_pressure diagnostics = DownstreamPressureDiagnostics( converged=converged, iterations=iteration_count, residual=final_residual, pressure=final_pressure, target_total_internal_energy=target_total_internal_energy, ) if strict and not diagnostics.converged: raise RuntimeError( "Downstream common-pressure solve did not converge: " f"residual={final_residual:.6e}" ) return final_pressure, diagnostics def project_downstream_pressure_constraints(self, *, strict: bool = False) -> None: downstream_components = ( self.components.upper_branch.pipe, self.components.lower_branch.pipe, self.components.tank, ) total_internal_energy = sum(component.state.U for component in downstream_components) common_pressure, diagnostics = self._solve_downstream_common_pressure( downstream_components, total_internal_energy, strict=strict, ) self.last_downstream_pressure_diagnostics = diagnostics for component in downstream_components: self._project_volume_energy_to_pressure(component, common_pressure) def _downstream_connection_enthalpy( self, q_out_upper: float, q_out_lower: float, pipe_upper_h: float, pipe_lower_h: float, tank_h: float, ) -> float: return self.components.downstream_tee.inlet_stream_enthalpy( q_out_lower, pipe_lower_h, q_out_upper, pipe_upper_h, fallback_h=tank_h, ) def _solve_downstream_branch_flows_from_state( self, *, inlet_h_upper: float, inlet_h_lower: float, pipe_upper_h: float, pipe_lower_h: float, tank_h: float, q_in_upper: float, q_in_lower: float, ) -> tuple[float, float]: return self.components.downstream_tee.solve_branch_outlet_flows_from_energy_balance( ratio_branch1=self.components.upper_branch.pipe.V / self.components.tank.V, ratio_branch2=self.components.lower_branch.pipe.V / self.components.tank.V, inlet_h_branch1=inlet_h_upper, inlet_h_branch2=inlet_h_lower, branch1_h=pipe_upper_h, branch2_h=pipe_lower_h, inlet_h=tank_h, q_in_branch1=q_in_upper, q_in_branch2=q_in_lower, ) def _evaluate_branch_states( self, cylinder: ThermodynamicProperties, ) -> tuple[BranchClosureState, BranchClosureState]: states: list[BranchClosureState] = [] for branch in self.components.branches(): pipe_properties = branch.pipe.properties() inlet_flow, inlet_flow_diagnostics = self._solve_branch_inlet_flow( branch.orifice, branch.pipe, cylinder.p, pipe_properties, ) inlet_h = branch.pipe.port_a_inlet_enthalpy( port_a_m_flow=inlet_flow, connected_h=cylinder.h, internal_h=pipe_properties.h, ) states.append( BranchClosureState( name=branch.name, pipe=pipe_properties, inlet_flow=inlet_flow, outlet_flow=0.0, inlet_h=inlet_h, inlet_flow_diagnostics=inlet_flow_diagnostics, ) ) return (states[0], states[1]) @staticmethod def _with_branch_outlet_flows( branch_states: tuple[BranchClosureState, BranchClosureState], outlet_flows: tuple[float, float], ) -> tuple[BranchClosureState, BranchClosureState]: return ( BranchClosureState( name=branch_states[0].name, pipe=branch_states[0].pipe, inlet_flow=branch_states[0].inlet_flow, outlet_flow=outlet_flows[0], inlet_h=branch_states[0].inlet_h, inlet_flow_diagnostics=branch_states[0].inlet_flow_diagnostics, ), BranchClosureState( name=branch_states[1].name, pipe=branch_states[1].pipe, inlet_flow=branch_states[1].inlet_flow, outlet_flow=outlet_flows[1], inlet_h=branch_states[1].inlet_h, inlet_flow_diagnostics=branch_states[1].inlet_flow_diagnostics, ), ) @staticmethod def _branch_snapshots( branch_states: tuple[BranchClosureState, BranchClosureState], ) -> tuple[BranchSnapshot, BranchSnapshot]: return ( BranchSnapshot( name=branch_states[0].name, pipe=branch_states[0].pipe, inlet_flow=branch_states[0].inlet_flow, outlet_flow=branch_states[0].outlet_flow, inlet_h=branch_states[0].inlet_h, inlet_flow_diagnostics=branch_states[0].inlet_flow_diagnostics, ), BranchSnapshot( name=branch_states[1].name, pipe=branch_states[1].pipe, inlet_flow=branch_states[1].inlet_flow, outlet_flow=branch_states[1].outlet_flow, inlet_h=branch_states[1].inlet_h, inlet_flow_diagnostics=branch_states[1].inlet_flow_diagnostics, ), ) def snapshot( self, state_vector: list[float] | None = None, *, strict: bool = False, ) -> TestModelSnapshot: if state_vector is not None: self._apply_state_vector(list(state_vector)) cylinder = self.components.cylinder.properties() tank = self.components.tank.properties() branch_states = self._evaluate_branch_states(cylinder) if strict: for branch_state in branch_states: if not branch_state.inlet_flow_diagnostics.converged: raise RuntimeError( "Branch inlet flow solve did not converge for " f"{branch_state.name}: residual=" f"{branch_state.inlet_flow_diagnostics.residual:.6e}" ) outlet_flows = self._solve_downstream_branch_flows(cylinder, tank, branch_states) branch_states = self._with_branch_outlet_flows(branch_states, outlet_flows) tee_upstream_h = self.components.upstream_tee.inlet_stream_enthalpy( -branch_states[0].inlet_flow, branch_states[0].pipe.h, -branch_states[1].inlet_flow, branch_states[1].pipe.h, fallback_h=cylinder.h, ) tee_downstream_h = self._downstream_connection_enthalpy( branch_states[0].outlet_flow, branch_states[1].outlet_flow, branch_states[0].pipe.h, branch_states[1].pipe.h, tank.h, ) self._write_port_states( cylinder, tank, branch_states, tee_upstream_h, tee_downstream_h, ) solve_diagnostics = TestModelSolveDiagnostics( upper_branch_inlet=branch_states[0].inlet_flow_diagnostics, lower_branch_inlet=branch_states[1].inlet_flow_diagnostics, downstream_pressure_projection=self.last_downstream_pressure_diagnostics, ) self.last_solve_diagnostics = solve_diagnostics branch_snapshots = self._branch_snapshots(branch_states) return TestModelSnapshot( cylinder=cylinder, tank=tank, tee_upstream_h=tee_upstream_h, tee_downstream_h=tee_downstream_h, branches=branch_snapshots, solve_diagnostics=solve_diagnostics, ) def _write_port_states( self, cylinder: ThermodynamicProperties, tank: ThermodynamicProperties, branch_states: tuple[BranchClosureState, BranchClosureState], tee_upstream_h: float, tee_downstream_h: float, ) -> None: cylinder_m_flow = -sum(branch_state.inlet_flow for branch_state in branch_states) tank_m_flow = sum(branch_state.outlet_flow for branch_state in branch_states) self.components.cylinder.port_b.m_flow = cylinder_m_flow self.components.upstream_tee.port_in.p = cylinder.p self.components.upstream_tee.port_out1.p = cylinder.p self.components.upstream_tee.port_out2.p = cylinder.p self.components.upstream_tee.port_in.m_flow = -cylinder_m_flow self.components.upstream_tee.port_in.h_outflow = tee_upstream_h self.components.upstream_tee.port_out1.h_outflow = cylinder.h self.components.upstream_tee.port_out2.h_outflow = cylinder.h self.components.upstream_tee.port_out1.m_flow = -branch_states[0].inlet_flow self.components.upstream_tee.port_out2.m_flow = -branch_states[1].inlet_flow for branch_components, branch_state in zip(self.components.branches(), branch_states): branch_components.orifice.port_a.p = cylinder.p branch_components.orifice.port_b.p = branch_components.pipe.inlet_pressure( branch_state.inlet_flow, max(branch_state.pipe.rho, 1e-9), branch_state.pipe.p, ) branch_components.orifice.port_a.m_flow = branch_state.inlet_flow branch_components.orifice.port_b.m_flow = -branch_state.inlet_flow branch_components.orifice.port_a.h_outflow = cylinder.h branch_components.orifice.port_b.h_outflow = branch_state.pipe.h branch_components.pipe.port_a.p = branch_components.orifice.port_b.p branch_components.pipe.port_a.m_flow = branch_state.inlet_flow branch_components.pipe.port_b.m_flow = -branch_state.outlet_flow branch_components.pipe.port_b.p = branch_state.pipe.p self.components.downstream_tee.port_in.p = tank.p self.components.downstream_tee.port_out1.p = tank.p self.components.downstream_tee.port_out2.p = tank.p self.components.downstream_tee.port_in.m_flow = -tank_m_flow self.components.downstream_tee.port_out1.m_flow = branch_states[1].outlet_flow self.components.downstream_tee.port_out2.m_flow = branch_states[0].outlet_flow self.components.downstream_tee.port_in.h_outflow = tee_downstream_h self.components.downstream_tee.port_out1.h_outflow = tank.h self.components.downstream_tee.port_out2.h_outflow = tank.h self.components.tank.port_a.m_flow = tank_m_flow def _branch_derivative_states( self, snapshot: TestModelSnapshot, ) -> tuple[VolumeState, VolumeState]: derivative_states: list[VolumeState] = [] for branch_components, branch_snapshot in zip(self.components.branches(), snapshot.branches): derivative_states.append( branch_components.pipe.derivatives_from_connections( port_a_m_flow=branch_snapshot.inlet_flow, connected_h_a=snapshot.cylinder.h, port_b_m_flow=-branch_snapshot.outlet_flow, connected_h_b=snapshot.tank.h, internal_h=branch_snapshot.pipe.h, ) ) return (derivative_states[0], derivative_states[1]) def rhs(self, state_vector: list[float]) -> list[float]: snapshot = self.snapshot(state_vector) cylinder_m_flow = -sum(branch.inlet_flow for branch in snapshot.branches) tank_m_flow = sum(branch.outlet_flow for branch in snapshot.branches) d_cylinder = self.components.cylinder.derivatives_from_connection( connected_h=snapshot.tee_upstream_h, port_m_flow=cylinder_m_flow, internal_h=snapshot.cylinder.h, ) branch_derivatives = self._branch_derivative_states(snapshot) d_tank = self.components.tank.derivatives_from_connection( connected_h=snapshot.tee_downstream_h, port_m_flow=tank_m_flow, internal_h=snapshot.tank.h, ) return [ d_cylinder.m, d_cylinder.U, branch_derivatives[0].m, branch_derivatives[0].U, branch_derivatives[1].m, branch_derivatives[1].U, d_tank.m, d_tank.U, ]