259 lines
8.9 KiB
Python
259 lines
8.9 KiB
Python
from __future__ import annotations
|
|
|
|
from dataclasses import dataclass
|
|
from math import sqrt
|
|
|
|
from app.simulation.core.ports import PortState, VariableRole
|
|
from app.simulation.systems.network import SimulationNetwork
|
|
|
|
|
|
class AlgebraicSolveError(RuntimeError):
|
|
def __init__(self, message: str, diagnostics: "AlgebraicSolveDiagnostics") -> None:
|
|
super().__init__(message)
|
|
self.diagnostics = diagnostics
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class AlgebraicUnknown:
|
|
component: str
|
|
port: str
|
|
variable: str
|
|
role: VariableRole
|
|
state: PortState
|
|
|
|
@property
|
|
def id(self) -> str:
|
|
return f"{self.component}.{self.port}.{self.variable}"
|
|
|
|
def read(self) -> float:
|
|
return float(getattr(self.state, self.variable))
|
|
|
|
def write(self, value: float) -> None:
|
|
setattr(self.state, self.variable, float(value))
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class AlgebraicSolveDiagnostics:
|
|
success: bool
|
|
message: str
|
|
evaluations: int
|
|
pressure_scale: float
|
|
flow_scale: float
|
|
max_scaled_residual: float
|
|
max_raw_residual: float
|
|
|
|
def as_dict(self) -> dict[str, object]:
|
|
return {
|
|
"success": self.success,
|
|
"message": self.message,
|
|
"evaluations": self.evaluations,
|
|
"pressureScale": self.pressure_scale,
|
|
"flowScale": self.flow_scale,
|
|
"maxScaledResidual": self.max_scaled_residual,
|
|
"maxRawResidual": self.max_raw_residual,
|
|
}
|
|
|
|
|
|
class PressureFlowSolver:
|
|
"""Solve the acausal pressure-flow subsystem for a compiled network."""
|
|
|
|
def __init__(
|
|
self,
|
|
network: SimulationNetwork,
|
|
*,
|
|
residual_tolerance: float = 1e-7,
|
|
max_evaluations: int = 500,
|
|
) -> None:
|
|
self.network = network
|
|
self.residual_tolerance = residual_tolerance
|
|
self.max_evaluations = max_evaluations
|
|
self.unknowns = self._build_unknowns()
|
|
self.last_diagnostics: AlgebraicSolveDiagnostics | None = None
|
|
|
|
def _build_unknowns(self) -> tuple[AlgebraicUnknown, ...]:
|
|
unknowns: list[AlgebraicUnknown] = []
|
|
for component in self.network.components.values():
|
|
for definition in component.port_definitions:
|
|
if definition.kind != "physical":
|
|
continue
|
|
state = component.get_port(definition.name)
|
|
for variable in definition.variables:
|
|
if variable.role not in {"effort", "flow"}:
|
|
continue
|
|
unknowns.append(
|
|
AlgebraicUnknown(
|
|
component=component.name,
|
|
port=definition.name,
|
|
variable=variable.name,
|
|
role=variable.role,
|
|
state=state,
|
|
)
|
|
)
|
|
return tuple(unknowns)
|
|
|
|
def _seed_equal_pressures(self) -> None:
|
|
for _ in range(max(2, len(self.network.connections))):
|
|
changed = False
|
|
for connection in self.network.connections:
|
|
if connection.kind != "physical":
|
|
continue
|
|
first = self.network.components[
|
|
connection.endpoint_a.component
|
|
].get_port(connection.endpoint_a.port)
|
|
second = self.network.components[
|
|
connection.endpoint_b.component
|
|
].get_port(connection.endpoint_b.port)
|
|
if first.p > 0.0 and second.p <= 0.0:
|
|
second.p = first.p
|
|
changed = True
|
|
elif second.p > 0.0 and first.p <= 0.0:
|
|
first.p = second.p
|
|
changed = True
|
|
|
|
for component in self.network.components.values():
|
|
equal_pressure_equations = [
|
|
equation
|
|
for equation in component.pressure_flow_equation_residuals()
|
|
if equation.relation == "equal" and equation.role == "effort"
|
|
]
|
|
for equation in equal_pressure_equations:
|
|
states = []
|
|
for variable in equation.variables:
|
|
_, port_name, variable_name = variable.rsplit(".", 2)
|
|
if variable_name == "p":
|
|
states.append(component.get_port(port_name))
|
|
if len(states) != 2:
|
|
continue
|
|
first, second = states
|
|
if first.p > 0.0 and second.p <= 0.0:
|
|
second.p = first.p
|
|
changed = True
|
|
elif second.p > 0.0 and first.p <= 0.0:
|
|
first.p = second.p
|
|
changed = True
|
|
if not changed:
|
|
break
|
|
|
|
def _scales(self) -> tuple[float, float]:
|
|
pressure_scale = max(
|
|
[
|
|
abs(unknown.read())
|
|
for unknown in self.unknowns
|
|
if unknown.role == "effort" and unknown.read() > 0.0
|
|
]
|
|
+ [1e5]
|
|
)
|
|
estimated_flows = [
|
|
abs(float(getattr(component, "K_eff"))) * sqrt(pressure_scale)
|
|
for component in self.network.components.values()
|
|
if hasattr(component, "K_eff")
|
|
]
|
|
flow_scale = max(
|
|
estimated_flows
|
|
+ [
|
|
abs(unknown.read())
|
|
for unknown in self.unknowns
|
|
if unknown.role == "flow"
|
|
]
|
|
+ [1e-3]
|
|
)
|
|
return pressure_scale, flow_scale
|
|
|
|
def solve(self) -> AlgebraicSolveDiagnostics:
|
|
try:
|
|
import numpy as np
|
|
from scipy.optimize import least_squares
|
|
except ImportError as exc:
|
|
raise RuntimeError(
|
|
"Topology-driven simulation requires SciPy; install requirements.txt."
|
|
) from exc
|
|
|
|
self._seed_equal_pressures()
|
|
pressure_scale, flow_scale = self._scales()
|
|
positive_pressures = [
|
|
unknown.read()
|
|
for unknown in self.unknowns
|
|
if unknown.role == "effort" and unknown.read() > 0.0
|
|
]
|
|
fallback_pressure = (
|
|
sum(positive_pressures) / len(positive_pressures)
|
|
if positive_pressures
|
|
else pressure_scale
|
|
)
|
|
|
|
def variable_scale(unknown: AlgebraicUnknown) -> float:
|
|
return pressure_scale if unknown.role == "effort" else flow_scale
|
|
|
|
x0 = np.asarray(
|
|
[
|
|
(
|
|
unknown.read()
|
|
if unknown.role != "effort" or unknown.read() > 0.0
|
|
else fallback_pressure
|
|
)
|
|
/ variable_scale(unknown)
|
|
for unknown in self.unknowns
|
|
],
|
|
dtype=float,
|
|
)
|
|
lower = np.asarray(
|
|
[
|
|
1.0 / pressure_scale if unknown.role == "effort" else -np.inf
|
|
for unknown in self.unknowns
|
|
]
|
|
)
|
|
upper = np.full(len(self.unknowns), np.inf)
|
|
|
|
def assign(values) -> None:
|
|
for unknown, value in zip(self.unknowns, values):
|
|
unknown.write(float(value) * variable_scale(unknown))
|
|
|
|
def scaled_residuals(values):
|
|
assign(values)
|
|
equations = self.network.pressure_flow_equation_residuals()
|
|
return np.asarray(
|
|
[
|
|
equation.value
|
|
/ (pressure_scale if equation.role == "effort" else flow_scale)
|
|
for equation in equations
|
|
],
|
|
dtype=float,
|
|
)
|
|
|
|
result = least_squares(
|
|
scaled_residuals,
|
|
x0,
|
|
bounds=(lower, upper),
|
|
x_scale="jac",
|
|
ftol=1e-10,
|
|
xtol=1e-10,
|
|
gtol=1e-10,
|
|
max_nfev=self.max_evaluations,
|
|
)
|
|
assign(result.x)
|
|
equations = self.network.pressure_flow_equation_residuals()
|
|
scaled = [
|
|
abs(
|
|
equation.value
|
|
/ (pressure_scale if equation.role == "effort" else flow_scale)
|
|
)
|
|
for equation in equations
|
|
]
|
|
success = bool(result.success) and max(scaled, default=0.0) <= self.residual_tolerance
|
|
diagnostics = AlgebraicSolveDiagnostics(
|
|
success=success,
|
|
message=str(result.message),
|
|
evaluations=int(result.nfev),
|
|
pressure_scale=pressure_scale,
|
|
flow_scale=flow_scale,
|
|
max_scaled_residual=max(scaled, default=0.0),
|
|
max_raw_residual=max((abs(item.value) for item in equations), default=0.0),
|
|
)
|
|
self.last_diagnostics = diagnostics
|
|
if not success:
|
|
raise AlgebraicSolveError(
|
|
"Pressure-flow equations did not converge to the requested tolerance.",
|
|
diagnostics,
|
|
)
|
|
return diagnostics
|