1753 lines
72 KiB
Python
1753 lines
72 KiB
Python
from __future__ import annotations
|
|
|
|
from collections.abc import Callable
|
|
from dataclasses import dataclass, replace
|
|
from math import floor, isfinite
|
|
import os
|
|
from sys import maxsize
|
|
from typing import Literal
|
|
|
|
from app.simulation.core.base import Component, DynamicComponent
|
|
from app.simulation.core.metadata import ResultVariableMetadata
|
|
from app.simulation.core.ports import PortState
|
|
from app.simulation.performance import performance_span, profile_phase
|
|
from app.simulation.property_cache import with_property_cache
|
|
from app.simulation.solvers.algebraic import PressureFlowSolver
|
|
from app.simulation.solvers.algebraic_blocks import StreamPressureBlockSolver
|
|
from app.simulation.solvers.jacobian import (
|
|
SparseJacobianCompatibilityError,
|
|
SparseSecantJacobian,
|
|
)
|
|
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 (
|
|
IntegrationCancelled,
|
|
ODESolution,
|
|
SolveIVPConfig,
|
|
integrate_ode,
|
|
)
|
|
from app.simulation.solvers.signal import SignalResolver
|
|
from app.simulation.solvers.stream import StreamResolver
|
|
from app.simulation.solvers.tangent import (
|
|
ThreePistonTangentCompilation,
|
|
ThreePistonTangentProvider,
|
|
compile_supported_piston_tangent_provider,
|
|
)
|
|
from app.simulation.solvers.thermofluid import (
|
|
ThermofluidClosureDiagnostics,
|
|
ThermofluidClosureError,
|
|
ThermofluidClosureFailure,
|
|
ThermofluidClosureSuccess,
|
|
ThermofluidTransactionPlan,
|
|
)
|
|
from app.simulation.systems.network import Endpoint, SimulationNetwork
|
|
|
|
|
|
SimulationProgressCallback = Callable[[float, str], None]
|
|
SimulationCancellationCheck = Callable[[], bool]
|
|
SimulationRunStatus = Literal["completed", "cancelled", "failed"]
|
|
ODE_JACOBIAN_MODE_ENVIRONMENT_VARIABLE = "SIMULATION_ODE_JACOBIAN_MODE"
|
|
|
|
|
|
def _requested_ode_jacobian_mode() -> Literal[
|
|
"optimized",
|
|
"hybrid",
|
|
"semi-analytic",
|
|
"scipy",
|
|
]:
|
|
value = os.getenv(
|
|
ODE_JACOBIAN_MODE_ENVIRONMENT_VARIABLE,
|
|
"scipy",
|
|
).strip().lower()
|
|
if value in {"optimized", "colored"}:
|
|
return "optimized"
|
|
if value in {"hybrid", "secant"}:
|
|
return "hybrid"
|
|
if value in {"semi-analytic", "semi_analytic", "analytic"}:
|
|
return "semi-analytic"
|
|
if value in {
|
|
"scipy",
|
|
"native",
|
|
"finite-difference",
|
|
"0",
|
|
"false",
|
|
"no",
|
|
"off",
|
|
}:
|
|
return "scipy"
|
|
raise ValueError(
|
|
f"{ODE_JACOBIAN_MODE_ENVIRONMENT_VARIABLE} must be "
|
|
"'optimized', 'hybrid', 'semi-analytic', or 'scipy'."
|
|
)
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class _ThermofluidClosurePlan:
|
|
"""Static execution data for one compiled network.
|
|
|
|
The first pressure-flow solve remains global. Later fixed-point passes only
|
|
need the physical islands whose constitutive equations read stream-derived
|
|
enthalpy. An unclassified custom stream component deliberately falls back
|
|
to the original global solve.
|
|
"""
|
|
|
|
physical_ports: tuple[PortState, ...]
|
|
global_component_group: tuple[str, ...]
|
|
secondary_pressure_solvers: tuple[PressureFlowSolver, ...]
|
|
secondary_component_groups: tuple[tuple[str, ...], ...]
|
|
uses_conservative_global_solver: bool
|
|
conservative_fallback_reason: str | None
|
|
secondary_block_solvers: tuple[StreamPressureBlockSolver, ...] = ()
|
|
|
|
|
|
@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 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,
|
|
) -> list[float]:
|
|
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.",
|
|
)
|
|
|
|
ratio = duration / step
|
|
if not isfinite(ratio):
|
|
raise SimulationSampleTimeError(
|
|
"SIMULATION_SAMPLE_COUNT_UNREPRESENTABLE",
|
|
"Simulation sample count cannot be represented by this runtime; "
|
|
"increase sampleStep.",
|
|
)
|
|
interval_count = floor(ratio)
|
|
# There is no product-level point cap. Still reject a collection that the
|
|
# Python runtime cannot index before multiplying by the potentially huge
|
|
# interval count or allocating the output grid.
|
|
if interval_count > maxsize - 2:
|
|
raise SimulationSampleTimeError(
|
|
"SIMULATION_SAMPLE_COUNT_UNREPRESENTABLE",
|
|
"Simulation sample count cannot be represented by this runtime; "
|
|
"increase sampleStep.",
|
|
)
|
|
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 > maxsize:
|
|
raise SimulationSampleTimeError(
|
|
"SIMULATION_SAMPLE_COUNT_UNREPRESENTABLE",
|
|
"Simulation sample count cannot be represented by this runtime; "
|
|
"increase sampleStep.",
|
|
)
|
|
|
|
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._thermofluid_closure_plan = self._build_thermofluid_closure_plan()
|
|
self._thermofluid_transaction_plan = ThermofluidTransactionPlan.compile(
|
|
network,
|
|
diagnostic_owners=(
|
|
self.signal_resolver,
|
|
self.pneumatic_volume_resolver,
|
|
self.stream_resolver,
|
|
self.pressure_flow_solver,
|
|
*self._thermofluid_closure_plan.secondary_pressure_solvers,
|
|
),
|
|
)
|
|
self._thermofluid_closure_diagnostics = ThermofluidClosureDiagnostics()
|
|
self.algebraic_solve_count = 0
|
|
self.algebraic_seeded_solve_count = 0
|
|
self.algebraic_nonlinear_solve_count = 0
|
|
self.algebraic_optimizer_evaluation_count = 0
|
|
self.algebraic_residual_evaluation_count = 0
|
|
self.algebraic_block_fallback_count = 0
|
|
self.algebraic_dense_fallback_count = 0
|
|
self.thermofluid_pressure_pass_count = 0
|
|
self.max_algebraic_residual = 0.0
|
|
self.max_algebraic_evaluations = 0
|
|
self.max_algebraic_residual_evaluations = 0
|
|
self._last_algebraic_diagnostics = None
|
|
self._last_algebraic_scope: tuple[str, ...] = ()
|
|
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
|
|
self._ode_tangent_provider: ThreePistonTangentProvider | None = None
|
|
|
|
def _request_causal_residual_audit(self) -> None:
|
|
"""Make topology or mode boundaries verify the next causal closure."""
|
|
|
|
self.pressure_flow_solver.request_causal_audit()
|
|
for block_solver in (
|
|
self._thermofluid_closure_plan.secondary_block_solvers
|
|
):
|
|
block_solver.request_causal_audit()
|
|
|
|
@staticmethod
|
|
def _overrides_stream_update(component: Component) -> bool:
|
|
component_type = type(component)
|
|
return (
|
|
component_type.update_stream_outflows
|
|
is not Component.update_stream_outflows
|
|
or component_type.update_flow_temperature_references
|
|
is not Component.update_flow_temperature_references
|
|
)
|
|
|
|
def _physical_component_groups(self) -> tuple[tuple[str, ...], ...]:
|
|
"""Return physical islands in component insertion order."""
|
|
|
|
physical_names = tuple(
|
|
component.name
|
|
for component in self.network.components.values()
|
|
if any(
|
|
definition.kind == "physical"
|
|
for definition in component.active_port_definitions
|
|
)
|
|
)
|
|
adjacency = {name: set() for name in physical_names}
|
|
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)
|
|
|
|
groups: list[tuple[str, ...]] = []
|
|
visited: set[str] = set()
|
|
for root in physical_names:
|
|
if root in visited:
|
|
continue
|
|
members = {root}
|
|
pending = [root]
|
|
visited.add(root)
|
|
while pending:
|
|
current = pending.pop()
|
|
for neighbor in adjacency[current]:
|
|
if neighbor in visited:
|
|
continue
|
|
visited.add(neighbor)
|
|
members.add(neighbor)
|
|
pending.append(neighbor)
|
|
groups.append(tuple(name for name in physical_names if name in members))
|
|
return tuple(groups)
|
|
|
|
def _network_for_component_group(
|
|
self,
|
|
component_names: tuple[str, ...],
|
|
all_physical_names: frozenset[str],
|
|
) -> SimulationNetwork:
|
|
if frozenset(component_names) == all_physical_names:
|
|
return self.network
|
|
|
|
selected = frozenset(component_names)
|
|
subnetwork = SimulationNetwork(
|
|
name=f"{self.network.name}:thermofluid:{len(component_names)}"
|
|
)
|
|
for component in self.network.components.values():
|
|
if component.name in selected:
|
|
subnetwork.add_component(component)
|
|
subnetwork.connections.extend(
|
|
connection
|
|
for connection in self.network.connections
|
|
if connection.kind == "physical"
|
|
and connection.endpoint_a.component in selected
|
|
and connection.endpoint_b.component in selected
|
|
)
|
|
return subnetwork
|
|
|
|
def _pressure_solver_for_component_group(
|
|
self,
|
|
component_names: tuple[str, ...],
|
|
all_physical_names: frozenset[str],
|
|
) -> PressureFlowSolver:
|
|
subnetwork = self._network_for_component_group(
|
|
component_names,
|
|
all_physical_names,
|
|
)
|
|
if subnetwork is self.network:
|
|
return self.pressure_flow_solver
|
|
return PressureFlowSolver(
|
|
subnetwork,
|
|
residual_tolerance=self.pressure_flow_solver.residual_tolerance,
|
|
max_evaluations=self.pressure_flow_solver.max_evaluations,
|
|
scope_kind="physicalIsland",
|
|
)
|
|
|
|
def _build_thermofluid_closure_plan(self) -> _ThermofluidClosurePlan:
|
|
physical_ports = tuple(
|
|
component.get_port(definition.name)
|
|
for component in self.network.components.values()
|
|
for definition in component.active_port_definitions
|
|
if definition.kind == "physical"
|
|
)
|
|
physical_groups = self._physical_component_groups()
|
|
all_physical_names = frozenset(
|
|
name for group in physical_groups for name in group
|
|
)
|
|
all_physical_order = tuple(
|
|
name for group in physical_groups for name in group
|
|
)
|
|
|
|
sensitive_names: set[str] = set()
|
|
has_unclassified_stream_component = False
|
|
has_invalid_dependency_declaration = False
|
|
for name in all_physical_names:
|
|
component = self.network.components[name]
|
|
# Only an exact-class declaration opts into pruning. A custom
|
|
# subclass cannot accidentally inherit a purity promise after
|
|
# changing its stream hook or constitutive equations.
|
|
declared = type(component).__dict__.get(
|
|
"PRESSURE_FLOW_DEPENDS_ON_STREAM"
|
|
)
|
|
if declared is True:
|
|
sensitive_names.add(name)
|
|
elif declared is False:
|
|
continue
|
|
elif declared is None and self._overrides_stream_update(component):
|
|
# Preserve the exact legacy behavior for custom components that
|
|
# receive stream values but have not declared equation purity.
|
|
has_unclassified_stream_component = True
|
|
elif declared is not None:
|
|
has_invalid_dependency_declaration = True
|
|
|
|
if has_invalid_dependency_declaration:
|
|
return _ThermofluidClosurePlan(
|
|
physical_ports=physical_ports,
|
|
global_component_group=all_physical_order,
|
|
secondary_pressure_solvers=(self.pressure_flow_solver,),
|
|
secondary_component_groups=(all_physical_order,),
|
|
uses_conservative_global_solver=True,
|
|
conservative_fallback_reason="invalidDependencyDeclaration",
|
|
)
|
|
|
|
if has_unclassified_stream_component:
|
|
return _ThermofluidClosurePlan(
|
|
physical_ports=physical_ports,
|
|
global_component_group=all_physical_order,
|
|
secondary_pressure_solvers=(self.pressure_flow_solver,),
|
|
secondary_component_groups=(all_physical_order,),
|
|
uses_conservative_global_solver=True,
|
|
conservative_fallback_reason="unclassifiedStreamComponent",
|
|
)
|
|
|
|
component_names = set(self.network.components)
|
|
compiled_equations = self.pressure_flow_solver.equation_templates
|
|
for equation in compiled_equations:
|
|
if equation.owner != "component":
|
|
continue
|
|
referenced_components = {
|
|
parts[0]
|
|
for variable in equation.variables
|
|
if len(parts := variable.rsplit(".", 2)) == 3
|
|
and parts[0] in component_names
|
|
}
|
|
if referenced_components - {equation.owner_id}:
|
|
# Catalog equations are component-local and connectors carry
|
|
# cross-component constraints. A custom residual may violate
|
|
# that convention, so retain the unsplit global problem.
|
|
return _ThermofluidClosurePlan(
|
|
physical_ports=physical_ports,
|
|
global_component_group=all_physical_order,
|
|
secondary_pressure_solvers=(self.pressure_flow_solver,),
|
|
secondary_component_groups=(all_physical_order,),
|
|
uses_conservative_global_solver=True,
|
|
conservative_fallback_reason="crossComponentEquation",
|
|
)
|
|
|
|
for group in physical_groups:
|
|
selected = frozenset(group)
|
|
connection_ids = {
|
|
connection.id
|
|
for connection in self.network.connections
|
|
if connection.kind == "physical"
|
|
and connection.endpoint_a.component in selected
|
|
and connection.endpoint_b.component in selected
|
|
}
|
|
unknown_count = sum(
|
|
unknown.component in selected
|
|
for unknown in self.pressure_flow_solver.unknowns
|
|
)
|
|
equation_count = sum(
|
|
(
|
|
equation.owner == "component"
|
|
and equation.owner_id in selected
|
|
)
|
|
or (
|
|
equation.owner == "connection"
|
|
and equation.owner_id in connection_ids
|
|
)
|
|
for equation in compiled_equations
|
|
)
|
|
if unknown_count != equation_count:
|
|
# The full network can be square even when two disconnected
|
|
# rectangular islands happen to cancel each other's equation
|
|
# count. Preserve the original global least-squares problem in
|
|
# that unusual case rather than changing its solution space.
|
|
return _ThermofluidClosurePlan(
|
|
physical_ports=physical_ports,
|
|
global_component_group=all_physical_order,
|
|
secondary_pressure_solvers=(self.pressure_flow_solver,),
|
|
secondary_component_groups=(all_physical_order,),
|
|
uses_conservative_global_solver=True,
|
|
conservative_fallback_reason="nonSquarePhysicalIsland",
|
|
)
|
|
|
|
coupled_groups = tuple(
|
|
group for group in physical_groups if sensitive_names.intersection(group)
|
|
)
|
|
secondary_pressure_solvers = tuple(
|
|
self._pressure_solver_for_component_group(group, all_physical_names)
|
|
for group in coupled_groups
|
|
)
|
|
secondary_block_solvers = tuple(
|
|
StreamPressureBlockSolver(
|
|
pressure_solver,
|
|
tuple(name for name in group if name in sensitive_names),
|
|
)
|
|
for pressure_solver, group in zip(
|
|
secondary_pressure_solvers,
|
|
coupled_groups,
|
|
)
|
|
)
|
|
return _ThermofluidClosurePlan(
|
|
physical_ports=physical_ports,
|
|
global_component_group=all_physical_order,
|
|
secondary_pressure_solvers=secondary_pressure_solvers,
|
|
secondary_component_groups=coupled_groups,
|
|
uses_conservative_global_solver=False,
|
|
conservative_fallback_reason=None,
|
|
secondary_block_solvers=secondary_block_solvers,
|
|
)
|
|
|
|
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)
|
|
)
|
|
|
|
@staticmethod
|
|
def _entry_has_pneumatic_state(entry: object) -> bool:
|
|
"""Return whether one reduced ODE entry owns pneumatic state.
|
|
|
|
Mechanical constraint groups are synthetic state owners. Every other
|
|
entry is a dynamic component, so its active port metadata is the
|
|
topology-level way to classify it without depending on model names.
|
|
"""
|
|
|
|
if isinstance(entry, MechanicalConstraintGroup):
|
|
return False
|
|
return any(
|
|
definition.kind == "physical" and definition.domain == "pneumatic"
|
|
for definition in entry.active_port_definitions
|
|
)
|
|
|
|
def _add_pneumatic_volume_state_dependencies(
|
|
self,
|
|
dependencies: list[set[int]],
|
|
entries: tuple[object, ...],
|
|
owner_by_component: dict[str, int],
|
|
) -> None:
|
|
"""Close the cross-domain dependency hidden by external volume ports.
|
|
|
|
A pneumatic-volume source such as a piston writes swept volume from
|
|
mechanical coordinates into a connected storage component before the
|
|
pressure-flow closure. The ordinary physical-path walk intentionally
|
|
stops at a storage state. Consequently, a second storage connected to
|
|
that chamber can depend on the piston even though the path crosses the
|
|
chamber state, and that derivative was previously omitted from the BDF
|
|
sparsity pattern.
|
|
|
|
Reuse the resolver's compiled output/connection plan to locate each
|
|
receiving storage. Mechanical states already found from that receiver,
|
|
the volume source's own ODE state (when it has one), and pneumatic states
|
|
whose local closure reaches the receiver form one conservative
|
|
cross-domain dependency set. Add it in both directions. If executable
|
|
custom/source metadata cannot bound those drivers, use a dense pattern.
|
|
"""
|
|
|
|
resolver = self.pneumatic_volume_resolver
|
|
pneumatic_entries = tuple(
|
|
self._entry_has_pneumatic_state(entry) for entry in entries
|
|
)
|
|
all_entry_indexes = set(range(len(entries)))
|
|
|
|
def use_conservative_dense_pattern() -> None:
|
|
for entry_dependencies in dependencies:
|
|
entry_dependencies.update(all_entry_indexes)
|
|
|
|
for component in resolver._output_components:
|
|
# ``pneumatic_volume_outputs`` is executable code rather than an
|
|
# equation-level dependency declaration. Catalog components with
|
|
# no directed signal input can be bounded by their own ODE state
|
|
# and the mechanical states already connected through topology.
|
|
# Custom/output components with an external signal driver keep the
|
|
# implicit integrator safe by disabling sparsity for this system.
|
|
if (
|
|
not type(component).__module__.startswith(
|
|
"app.simulation.components."
|
|
)
|
|
or any(
|
|
(
|
|
definition.kind == "signal"
|
|
and definition.nominal_role == "input"
|
|
)
|
|
or (
|
|
definition.kind == "physical"
|
|
and definition.domain
|
|
not in {"pneumatic", "mechanical"}
|
|
)
|
|
for definition in component.active_port_definitions
|
|
)
|
|
):
|
|
use_conservative_dense_pattern()
|
|
return
|
|
|
|
source_index = owner_by_component.get(component.name)
|
|
receiver_indexes: set[int] = set()
|
|
for definition in component.active_port_definitions:
|
|
if (
|
|
definition.kind != "physical"
|
|
or definition.domain != "pneumatic"
|
|
):
|
|
continue
|
|
binding = resolver._connected_endpoint.get(
|
|
Endpoint(component.name, definition.name)
|
|
)
|
|
if binding is None:
|
|
continue
|
|
receiver_index = owner_by_component.get(
|
|
binding.connected_endpoint.component
|
|
)
|
|
if receiver_index is not None and pneumatic_entries[receiver_index]:
|
|
receiver_indexes.add(receiver_index)
|
|
|
|
for receiver_index in receiver_indexes:
|
|
driver_indexes = {
|
|
entry_index
|
|
for entry_index in dependencies[receiver_index]
|
|
if isinstance(entries[entry_index], MechanicalConstraintGroup)
|
|
}
|
|
if source_index is not None:
|
|
driver_indexes.add(source_index)
|
|
if not driver_indexes:
|
|
use_conservative_dense_pattern()
|
|
return
|
|
coupled_pneumatic_indexes = {
|
|
entry_index
|
|
for entry_index, is_pneumatic in enumerate(pneumatic_entries)
|
|
if is_pneumatic
|
|
and (
|
|
entry_index == receiver_index
|
|
or receiver_index in dependencies[entry_index]
|
|
)
|
|
}
|
|
for pneumatic_index in coupled_pneumatic_indexes:
|
|
dependencies[pneumatic_index].update(driver_indexes)
|
|
for driver_index in driver_indexes:
|
|
dependencies[driver_index].update(
|
|
coupled_pneumatic_indexes
|
|
)
|
|
|
|
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)
|
|
|
|
self._add_pneumatic_volume_state_dependencies(
|
|
dependencies,
|
|
entries,
|
|
owner_by_component,
|
|
)
|
|
|
|
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,
|
|
}
|
|
|
|
def _exact_ode_jacobian_rows(self) -> dict[int, dict[int, float]]:
|
|
"""Return mode-independent kinematic rows safe to evaluate exactly."""
|
|
|
|
rows: dict[int, dict[int, float]] = {}
|
|
cursor = 0
|
|
for entry in self.mechanical_state_reducer.state_entries:
|
|
if isinstance(entry, MechanicalConstraintGroup):
|
|
# A discrete endstop can replace x' = v with x' = 0 for the
|
|
# active constrained mode. Keep those rows numerical; free
|
|
# mechanical groups always have d(x')/d(v) = 1.
|
|
if not entry.discrete_endstop_components:
|
|
rows[cursor + 1] = {cursor: 1.0}
|
|
cursor += 2
|
|
else:
|
|
cursor += entry.state_size
|
|
return rows
|
|
|
|
@profile_phase(
|
|
"simulation.closure",
|
|
minimum_mode="audit",
|
|
reset_property_shadow=True,
|
|
)
|
|
def _close_current_state(
|
|
self,
|
|
time: float,
|
|
*,
|
|
record_rhs_outcome: bool = False,
|
|
) -> dict[str, dict[str, float]]:
|
|
transaction = self._thermofluid_transaction_plan.capture()
|
|
last_algebraic_diagnostics = self._last_algebraic_diagnostics
|
|
last_algebraic_scope = self._last_algebraic_scope
|
|
try:
|
|
connected_h, success = self._close_current_state_unchecked(time)
|
|
except ThermofluidClosureError as exc:
|
|
transaction.restore()
|
|
self._last_algebraic_diagnostics = last_algebraic_diagnostics
|
|
self._last_algebraic_scope = last_algebraic_scope
|
|
self._request_causal_residual_audit()
|
|
if record_rhs_outcome:
|
|
failure = self._thermofluid_closure_diagnostics.record_failure(
|
|
exc.diagnostics
|
|
)
|
|
raise ThermofluidClosureError(failure) from None
|
|
raise
|
|
except BaseException:
|
|
transaction.restore()
|
|
self._last_algebraic_diagnostics = last_algebraic_diagnostics
|
|
self._last_algebraic_scope = last_algebraic_scope
|
|
self._request_causal_residual_audit()
|
|
raise
|
|
if record_rhs_outcome:
|
|
self._thermofluid_closure_diagnostics.record_success(success)
|
|
return connected_h
|
|
|
|
def _close_current_state_unchecked(
|
|
self,
|
|
time: float,
|
|
) -> tuple[
|
|
dict[str, dict[str, float]],
|
|
ThermofluidClosureSuccess,
|
|
]:
|
|
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()
|
|
initial_algebraic = self.pressure_flow_solver.solve(
|
|
effort_variables=("p",),
|
|
)
|
|
algebraic_diagnostics = [initial_algebraic]
|
|
pressure_flow_solve_count = 1
|
|
self.thermofluid_pressure_pass_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.
|
|
# The compiled closure plan keeps custom stream-aware components on the
|
|
# legacy global path. For catalog models, only stream-sensitive physical
|
|
# islands are revisited; independent islands keep the first solve.
|
|
closure_plan = self._thermofluid_closure_plan
|
|
self._last_algebraic_diagnostics = initial_algebraic
|
|
self._last_algebraic_scope = closure_plan.global_component_group
|
|
secondary_pressure_solvers = closure_plan.secondary_pressure_solvers
|
|
secondary_block_solvers = closure_plan.secondary_block_solvers
|
|
connected_h: dict[str, dict[str, float]] = {}
|
|
stream_diagnostics = []
|
|
coupling_deltas = []
|
|
max_coupling_iterations = 25
|
|
flow_relative_tolerance = 1.0e-12
|
|
for coupling_iteration in range(1, max_coupling_iterations + 1):
|
|
previous_flows = self._thermofluid_transaction_plan.flow_values()
|
|
stream, connected_h = self.stream_resolver.solve(
|
|
dynamic_ports_are_current=True,
|
|
)
|
|
stream_diagnostics.append(stream)
|
|
for component in self.dynamic_components:
|
|
component.update_stream_outflows(connected_h[component.name])
|
|
self.stream_resolver.refresh_flow_temperature_references()
|
|
if secondary_pressure_solvers:
|
|
self.thermofluid_pressure_pass_count += 1
|
|
block_scale_context = (
|
|
self.pressure_flow_solver.scale_context()
|
|
if any(
|
|
solver is not self.pressure_flow_solver
|
|
for solver in secondary_pressure_solvers
|
|
)
|
|
else None
|
|
)
|
|
if (
|
|
secondary_block_solvers
|
|
and not closure_plan.uses_conservative_global_solver
|
|
):
|
|
# Stream propagation only invalidates equations that explicitly
|
|
# consume the new enthalpy/temperature references. Re-solve the
|
|
# exact equation/unknown blocks containing those equations; the
|
|
# first global pass above remains the causalization boundary for
|
|
# mechanics, contact, and all stream-independent pneumatic blocks.
|
|
for block_solver in secondary_block_solvers:
|
|
block_result = block_solver.solve(
|
|
scale_context=block_scale_context,
|
|
)
|
|
# One public secondary closure is one logical solve. The
|
|
# block solver folds every local attempt and a possible
|
|
# accepted global fallback into this single diagnostic, so
|
|
# evaluations and blockFallbackUsed are counted exactly
|
|
# once here rather than once per internal equation block.
|
|
(algebraic,) = block_result.diagnostics
|
|
(scope,) = block_result.scopes
|
|
algebraic_diagnostics.append(algebraic)
|
|
self._last_algebraic_diagnostics = algebraic
|
|
self._last_algebraic_scope = scope
|
|
pressure_flow_solve_count += 1
|
|
else:
|
|
for pressure_solver, component_group in zip(
|
|
secondary_pressure_solvers,
|
|
closure_plan.secondary_component_groups,
|
|
):
|
|
algebraic = pressure_solver.solve(
|
|
effort_variables=(),
|
|
scale_context=(
|
|
block_scale_context
|
|
if pressure_solver is not self.pressure_flow_solver
|
|
else None
|
|
),
|
|
)
|
|
algebraic_diagnostics.append(algebraic)
|
|
self._last_algebraic_diagnostics = algebraic
|
|
self._last_algebraic_scope = component_group
|
|
pressure_flow_solve_count += 1
|
|
coupling_delta = self._thermofluid_transaction_plan.measure_flow_delta(
|
|
previous_flows,
|
|
iteration=coupling_iteration,
|
|
relative_tolerance=flow_relative_tolerance,
|
|
)
|
|
coupling_deltas.append(coupling_delta)
|
|
if (
|
|
not secondary_pressure_solvers
|
|
or coupling_delta.max_delta <= coupling_delta.tolerance
|
|
):
|
|
break
|
|
else:
|
|
raise ThermofluidClosureError(
|
|
ThermofluidClosureFailure.from_iterations(
|
|
time,
|
|
coupling_deltas,
|
|
)
|
|
)
|
|
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
|
|
seeded_count = sum(
|
|
item.jacobian_mode == "seeded" for item in algebraic_diagnostics
|
|
)
|
|
self.algebraic_seeded_solve_count += seeded_count
|
|
self.algebraic_nonlinear_solve_count += (
|
|
len(algebraic_diagnostics) - seeded_count
|
|
)
|
|
self.algebraic_optimizer_evaluation_count += sum(
|
|
item.evaluations for item in algebraic_diagnostics
|
|
)
|
|
self.algebraic_residual_evaluation_count += sum(
|
|
item.residual_evaluations for item in algebraic_diagnostics
|
|
)
|
|
self.algebraic_block_fallback_count += sum(
|
|
item.block_fallback_used for item in algebraic_diagnostics
|
|
)
|
|
self.algebraic_dense_fallback_count += sum(
|
|
item.dense_fallback_used for item in algebraic_diagnostics
|
|
)
|
|
self.max_algebraic_residual = max(
|
|
self.max_algebraic_residual,
|
|
*(item.max_scaled_residual for item in algebraic_diagnostics),
|
|
)
|
|
self.max_algebraic_evaluations = max(
|
|
self.max_algebraic_evaluations,
|
|
*(item.evaluations for item in algebraic_diagnostics),
|
|
)
|
|
self.max_algebraic_residual_evaluations = max(
|
|
self.max_algebraic_residual_evaluations,
|
|
*(item.residual_evaluations for item in algebraic_diagnostics),
|
|
)
|
|
self.max_stream_iterations = max(
|
|
self.max_stream_iterations,
|
|
*(item.iterations for item in stream_diagnostics),
|
|
)
|
|
return connected_h, ThermofluidClosureSuccess.from_iteration(
|
|
time,
|
|
coupling_deltas[-1],
|
|
)
|
|
|
|
@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,
|
|
record_rhs_outcome=True,
|
|
)
|
|
derivatives = self._state_derivatives(connected_h)
|
|
provider = self._ode_tangent_provider
|
|
if provider is not None:
|
|
provider.record_primal(_time, state_vector, connected_h)
|
|
return derivatives
|
|
|
|
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)
|
|
|
|
@with_property_cache
|
|
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
|
|
next_signal_audit_index = 0
|
|
|
|
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]:
|
|
nonlocal next_signal_audit_index
|
|
while (
|
|
next_signal_audit_index < len(signal_event_times)
|
|
and float(time) >= signal_event_times[next_signal_audit_index]
|
|
):
|
|
self._request_causal_residual_audit()
|
|
next_signal_audit_index += 1
|
|
if cancel_check is None:
|
|
report_solver_time(time)
|
|
return self.rhs(time, state_vector)
|
|
|
|
jacobian = None
|
|
jacobian_fallback_reason: str | None = None
|
|
tangent_compilation: ThreePistonTangentCompilation | None = None
|
|
selected_tangent_provider: ThreePistonTangentProvider | None = None
|
|
self._ode_tangent_provider = None
|
|
requested_jacobian_mode = (
|
|
_requested_ode_jacobian_mode()
|
|
if jac_sparsity is not None
|
|
else "scipy"
|
|
)
|
|
if (
|
|
jac_sparsity is not None
|
|
and requested_jacobian_mode
|
|
in {"optimized", "hybrid", "semi-analytic"}
|
|
):
|
|
state_count = int(jac_sparsity.shape[0])
|
|
if (
|
|
requested_jacobian_mode == "hybrid"
|
|
and not self.pressure_flow_solver.causal_fast_path_enabled
|
|
):
|
|
jacobian_fallback_reason = "causalAlgebraicExecutionUnavailable"
|
|
elif int(jac_sparsity.nnz) >= state_count * state_count:
|
|
jacobian_fallback_reason = "denseStateDependencyPattern"
|
|
else:
|
|
exact_columns = None
|
|
if requested_jacobian_mode == "semi-analytic":
|
|
tangent_compilation = (
|
|
compile_supported_piston_tangent_provider(self)
|
|
)
|
|
if tangent_compilation.eligible:
|
|
provider = tangent_compilation.provider
|
|
if provider is None:
|
|
raise RuntimeError(
|
|
"An eligible tangent compilation has no provider."
|
|
)
|
|
selected_tangent_provider = provider
|
|
exact_columns = (
|
|
tangent_compilation.columns,
|
|
provider,
|
|
)
|
|
else:
|
|
jacobian_fallback_reason = (
|
|
f"semiAnalytic:{tangent_compilation.reason}"
|
|
)
|
|
|
|
if (
|
|
requested_jacobian_mode != "semi-analytic"
|
|
or tangent_compilation is not None
|
|
and tangent_compilation.eligible
|
|
):
|
|
def evaluate_jacobian_rhs(time, state):
|
|
if cancel_check is not None and cancel_check():
|
|
raise IntegrationCancelled
|
|
return monitored_rhs(
|
|
time,
|
|
[float(value) for value in state],
|
|
)
|
|
|
|
try:
|
|
jacobian = SparseSecantJacobian(
|
|
evaluate_jacobian_rhs,
|
|
jac_sparsity,
|
|
integration_config.atol,
|
|
exact_rows=self._exact_ode_jacobian_rows(),
|
|
exact_columns=exact_columns,
|
|
max_consecutive_reuses=(
|
|
1
|
|
if requested_jacobian_mode == "hybrid"
|
|
else 0
|
|
),
|
|
)
|
|
self._ode_tangent_provider = (
|
|
selected_tangent_provider
|
|
)
|
|
except SparseJacobianCompatibilityError as exc:
|
|
jacobian_fallback_reason = (
|
|
f"scipyCompatibility:{type(exc).__name__}"
|
|
)
|
|
|
|
def handle_state_transition(*args):
|
|
transition = self.mechanical_state_reducer.state_transition(*args)
|
|
if transition is not None:
|
|
self._request_causal_residual_audit()
|
|
return transition
|
|
|
|
try:
|
|
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=(
|
|
handle_state_transition
|
|
if self.mechanical_state_reducer.has_state_events
|
|
else None
|
|
),
|
|
jac_sparsity=jac_sparsity,
|
|
jac=jacobian,
|
|
recoverable_trial_retries=True,
|
|
)
|
|
finally:
|
|
self._ode_tangent_provider = None
|
|
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,
|
|
}
|
|
]
|
|
if jacobian is not None:
|
|
direct_jacobian = jacobian.diagnostics()
|
|
solver_segment_diagnostics[0].update(
|
|
{
|
|
"jacobianEvaluationCount": int(
|
|
direct_jacobian["jacobianEvaluationCount"]
|
|
),
|
|
"jacobianFullBuildCount": int(
|
|
direct_jacobian["fullBuildCount"]
|
|
),
|
|
"jacobianSecantReuseCount": int(
|
|
direct_jacobian["secantReuseCount"]
|
|
),
|
|
"jacobianAuditFailureCount": int(
|
|
direct_jacobian["auditFailureCount"]
|
|
),
|
|
"finiteDifferenceRhsEvaluationCount": int(
|
|
direct_jacobian[
|
|
"finiteDifferenceRhsEvaluationCount"
|
|
]
|
|
),
|
|
"jacobianBaseRhsEvaluationCount": int(
|
|
direct_jacobian["baseRhsEvaluationCount"]
|
|
),
|
|
"jacobianJvAuditRhsEvaluationCount": int(
|
|
direct_jacobian["jvAuditEvaluationCount"]
|
|
),
|
|
"exactColumnBuildCount": int(
|
|
direct_jacobian["exactColumnBuildCount"]
|
|
),
|
|
"exactColumnFallbackCount": int(
|
|
direct_jacobian["exactColumnFallbackCount"]
|
|
),
|
|
"jacobianAssemblySeconds": float(
|
|
direct_jacobian["assemblySeconds"]
|
|
),
|
|
}
|
|
)
|
|
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_work_keys = (
|
|
"jacobianEvaluationCount",
|
|
"jacobianFullBuildCount",
|
|
"jacobianSecantReuseCount",
|
|
"jacobianAuditFailureCount",
|
|
"finiteDifferenceRhsEvaluationCount",
|
|
"jacobianBaseRhsEvaluationCount",
|
|
"jacobianJvAuditRhsEvaluationCount",
|
|
"exactColumnBuildCount",
|
|
"exactColumnFallbackCount",
|
|
"jacobianAssemblySeconds",
|
|
)
|
|
for key in jacobian_work_keys:
|
|
if any(key in segment for segment in solver_segment_diagnostics):
|
|
solver_totals[key] = sum(
|
|
segment.get(key, 0)
|
|
for segment in solver_segment_diagnostics
|
|
)
|
|
jacobian_diagnostics = (
|
|
self.jacobian_sparsity_diagnostics()
|
|
if integration_config.method in {"BDF", "Radau"}
|
|
else None
|
|
)
|
|
runtime_jacobian_diagnostics: dict[str, object] | None = None
|
|
if jacobian_diagnostics is not None:
|
|
color_group_count = int(jacobian_diagnostics["colorGroupCount"])
|
|
if jacobian is None:
|
|
for segment in solver_segment_diagnostics:
|
|
segment["finiteDifferenceRhsEstimate"] = (
|
|
int(segment["njev"]) * color_group_count
|
|
)
|
|
runtime_jacobian_diagnostics = {
|
|
"mode": "scipySparseFiniteDifference",
|
|
"fallbackReason": jacobian_fallback_reason,
|
|
"jacobianEvaluationCount": int(solver_totals["njev"]),
|
|
"fullBuildCount": int(solver_totals["njev"]),
|
|
"finiteDifferenceRhsEstimateIsExact": False,
|
|
}
|
|
else:
|
|
for segment in solver_segment_diagnostics:
|
|
segment["finiteDifferenceRhsEstimate"] = int(
|
|
segment.get("finiteDifferenceRhsEvaluationCount", 0)
|
|
) + int(
|
|
segment.get("jacobianJvAuditRhsEvaluationCount", 0)
|
|
)
|
|
runtime_jacobian_diagnostics = dict(jacobian.diagnostics())
|
|
runtime_jacobian_diagnostics.update(
|
|
{
|
|
"fallbackReason": jacobian_fallback_reason,
|
|
"finiteDifferenceRhsEstimateIsExact": True,
|
|
}
|
|
)
|
|
if tangent_compilation is not None:
|
|
runtime_jacobian_diagnostics["tangentCompilation"] = (
|
|
tangent_compilation.diagnostics()
|
|
)
|
|
if tangent_compilation.eligible and jacobian is not None:
|
|
runtime_jacobian_diagnostics["mode"] = (
|
|
"semiAnalyticExactColumns"
|
|
)
|
|
exact_builds = int(
|
|
runtime_jacobian_diagnostics[
|
|
"exactColumnBuildCount"
|
|
]
|
|
)
|
|
exact_fallbacks = int(
|
|
runtime_jacobian_diagnostics[
|
|
"exactColumnFallbackCount"
|
|
]
|
|
)
|
|
if exact_builds == 0 and exact_fallbacks == 0:
|
|
effective_mode = "notEvaluated"
|
|
elif exact_builds == 0:
|
|
effective_mode = "numericalFallbackOnly"
|
|
elif exact_fallbacks:
|
|
effective_mode = "mixedExactAndNumericalFallback"
|
|
else:
|
|
effective_mode = "exactColumns"
|
|
runtime_jacobian_diagnostics["effectiveMode"] = (
|
|
effective_mode
|
|
)
|
|
if exact_fallbacks:
|
|
runtime_jacobian_diagnostics[
|
|
"runtimeFallbackReason"
|
|
] = runtime_jacobian_diagnostics[
|
|
"lastExactColumnFallbackReason"
|
|
]
|
|
solver_totals["finiteDifferenceRhsEstimate"] = sum(
|
|
int(segment["finiteDifferenceRhsEstimate"])
|
|
for segment in solver_segment_diagnostics
|
|
)
|
|
if jacobian is None:
|
|
solver_totals["jacobianRhsEvaluationCountEstimate"] = (
|
|
int(solver_totals["finiteDifferenceRhsEstimate"])
|
|
+ int(solver_totals["njev"])
|
|
)
|
|
else:
|
|
solver_totals["jacobianRhsEvaluationCount"] = (
|
|
int(
|
|
solver_totals.get(
|
|
"finiteDifferenceRhsEvaluationCount",
|
|
0,
|
|
)
|
|
)
|
|
+ int(
|
|
solver_totals.get(
|
|
"jacobianBaseRhsEvaluationCount",
|
|
0,
|
|
)
|
|
)
|
|
+ int(
|
|
solver_totals.get(
|
|
"jacobianJvAuditRhsEvaluationCount",
|
|
0,
|
|
)
|
|
)
|
|
)
|
|
|
|
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,
|
|
"jacobian": runtime_jacobian_diagnostics,
|
|
"segmentCount": len(solver_segment_diagnostics),
|
|
"segments": solver_segment_diagnostics,
|
|
"totals": solver_totals,
|
|
},
|
|
"pressureFlow": {
|
|
"solveCount": self.algebraic_solve_count,
|
|
"seededSolveCount": self.algebraic_seeded_solve_count,
|
|
"nonlinearSolveCount": self.algebraic_nonlinear_solve_count,
|
|
"fastPathHitRate": (
|
|
self.algebraic_seeded_solve_count
|
|
/ self.algebraic_solve_count
|
|
if self.algebraic_solve_count
|
|
else 0.0
|
|
),
|
|
"optimizerEvaluationCount": (
|
|
self.algebraic_optimizer_evaluation_count
|
|
),
|
|
"residualEvaluationCount": (
|
|
self.algebraic_residual_evaluation_count
|
|
),
|
|
"blockFallbackCount": self.algebraic_block_fallback_count,
|
|
"denseFallbackCount": self.algebraic_dense_fallback_count,
|
|
"closurePassCount": self.thermofluid_pressure_pass_count,
|
|
"secondaryPhysicalIslandCount": len(
|
|
self._thermofluid_closure_plan.secondary_pressure_solvers
|
|
),
|
|
"secondaryBlockCount": sum(
|
|
len(solver.blocks)
|
|
for solver in self._thermofluid_closure_plan.secondary_block_solvers
|
|
if solver.available
|
|
),
|
|
"secondaryUnknownCount": sum(
|
|
len(block.unknowns)
|
|
for solver in self._thermofluid_closure_plan.secondary_block_solvers
|
|
if solver.available
|
|
for block in solver.blocks
|
|
),
|
|
"equationBlockFallbackReasons": [
|
|
solver.fallback_reason
|
|
for solver in self._thermofluid_closure_plan.secondary_block_solvers
|
|
if solver.fallback_reason is not None
|
|
],
|
|
"usesConservativeGlobalCoupling": (
|
|
self._thermofluid_closure_plan.uses_conservative_global_solver
|
|
),
|
|
"couplingPlanFallbackReason": (
|
|
self._thermofluid_closure_plan.conservative_fallback_reason
|
|
),
|
|
"maxScaledResidual": self.max_algebraic_residual,
|
|
"maxEvaluationsPerSolve": self.max_algebraic_evaluations,
|
|
"maxResidualEvaluationsPerSolve": (
|
|
self.max_algebraic_residual_evaluations
|
|
),
|
|
"lastScope": list(self._last_algebraic_scope),
|
|
"last": (
|
|
self._last_algebraic_diagnostics.as_dict()
|
|
if self._last_algebraic_diagnostics is not None
|
|
else None
|
|
),
|
|
"causalExecution": (
|
|
self.pressure_flow_solver.causal_execution_diagnostics()
|
|
),
|
|
"secondaryCausalExecution": [
|
|
solver.causal_execution_diagnostics()
|
|
for solver in self._thermofluid_closure_plan.secondary_block_solvers
|
|
],
|
|
},
|
|
"stream": {
|
|
"maxIterationsPerSolve": self.max_stream_iterations,
|
|
"maxThermofluidIterations": self.max_thermofluid_iterations,
|
|
"thermofluidClosure": {
|
|
**self._thermofluid_closure_diagnostics.as_dict(),
|
|
"transaction": (
|
|
self._thermofluid_transaction_plan.diagnostics()
|
|
),
|
|
},
|
|
"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,
|
|
)
|