sepehr/Entropyk
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fix bugs from 5-2 code review

Browse Source
This commit is contained in:
Sepehr
2026-02-21 10:43:55 +01:00
parent 400f1c420e
commit 0d9a0e4231
27 changed files with 9838 additions and 114 deletions
Download Patch File Download Diff File Expand all files Collapse all files

572
crates/components/src/flow_boundary.rs Normal file
View File

@@ -0,0 +1,572 @@
//! Boundary Condition Components — Source & Sink
//!
//! This module provides `FlowSource` and `FlowSink` for both incompressible
//! (water, glycol, brine) and compressible (refrigerant, CO₂) fluid systems.
//!
//! ## Design Philosophy (à la Modelica)
//!
//! - **`FlowSource`** imposes a fixed thermodynamic state (P, h) on its outlet
//! edge. It is the entry point of a fluid circuit — it represents an infinite
//! reservoir at constant conditions (city water supply, district heating header,
//! refrigerant reservoir, etc.).
//!
//! - **`FlowSink`** absorbs flow at a fixed pressure (back-pressure). It is the
//! termination point of a circuit. Optionally, a fixed outlet enthalpy can also
//! be imposed (isothermal return, phase separator, etc.).
//!
//! ## Equations
//!
//! ### FlowSource — 2 equations
//!
//! ```text
//! r_P = P_edge − P_set = 0 (pressure boundary condition)
//! r_h = h_edge − h_set = 0 (enthalpy / temperature BC)
//! ```
//!
//! ### FlowSink — 1 or 2 equations
//!
//! ```text
//! r_P = P_edge − P_back = 0 (back-pressure boundary condition)
//! [optional] r_h = h_edge − h_back = 0
//! ```
//!
//! ## Incompressible vs Compressible
//!
//! Same physics, different construction-time validation. Use:
//! - `FlowSource::incompressible` / `FlowSink::incompressible` for water, glycol…
//! - `FlowSource::compressible` / `FlowSink::compressible` for refrigerant, CO₂…
//!
//! ## Example
//!
//! ```no_run
//! use entropyk_components::flow_boundary::{FlowSource, FlowSink};
//! use entropyk_components::port::{FluidId, Port};
//! use entropyk_core::{Pressure, Enthalpy};
//!
//! let make_port = |p: f64, h: f64| {
//! let a = Port::new(FluidId::new("Water"), Pressure::from_pascals(p),
//! Enthalpy::from_joules_per_kg(h));
//! let b = Port::new(FluidId::new("Water"), Pressure::from_pascals(p),
//! Enthalpy::from_joules_per_kg(h));
//! a.connect(b).unwrap().0
//! };
//!
//! // City water supply: 3 bar, 15°C (h ≈ 63 kJ/kg)
//! let source = FlowSource::incompressible(
//! "Water", 3.0e5, 63_000.0, make_port(3.0e5, 63_000.0),
//! ).unwrap();
//!
//! // Return header: 1.5 bar back-pressure
//! let sink = FlowSink::incompressible(
//! "Water", 1.5e5, None, make_port(1.5e5, 63_000.0),
//! ).unwrap();
//! ```
use crate::{
flow_junction::is_incompressible, flow_junction::FluidKind, Component, ComponentError,
ConnectedPort, JacobianBuilder, ResidualVector, SystemState,
};
// ─────────────────────────────────────────────────────────────────────────────
// FlowSource — Fixed P & h boundary condition
// ─────────────────────────────────────────────────────────────────────────────
/// A boundary source that imposes fixed pressure and enthalpy on its outlet edge.
///
/// Represents an ideal infinite reservoir (city water, refrigerant header, steam
/// drum, etc.) at constant thermodynamic conditions.
///
/// # Equations (always 2)
///
/// ```text
/// r₀ = P_edge − P_set = 0
/// r₁ = h_edge − h_set = 0
/// ```
#[derive(Debug, Clone)]
pub struct FlowSource {
/// Fluid kind.
kind: FluidKind,
/// Fluid name.
fluid_id: String,
/// Set-point pressure [Pa].
p_set_pa: f64,
/// Set-point specific enthalpy [J/kg].
h_set_jkg: f64,
/// Connected outlet port (links to first edge in the System).
outlet: ConnectedPort,
}
impl FlowSource {
// ── Constructors ─────────────────────────────────────────────────────────
/// Creates an **incompressible** source (water, glycol, brine…).
///
/// # Arguments
///
/// * `fluid` — fluid identifier string (e.g. `"Water"`)
/// * `p_set_pa` — set-point pressure in Pascals
/// * `h_set_jkg` — set-point specific enthalpy in J/kg
/// * `outlet` — connected port linked to the first system edge
pub fn incompressible(
fluid: impl Into<String>,
p_set_pa: f64,
h_set_jkg: f64,
outlet: ConnectedPort,
) -> Result<Self, ComponentError> {
let fluid = fluid.into();
if !is_incompressible(&fluid) {
return Err(ComponentError::InvalidState(format!(
"FlowSource::incompressible: '{}' does not appear incompressible. \
Use FlowSource::compressible for refrigerants.",
fluid
)));
}
Self::new_inner(FluidKind::Incompressible, fluid, p_set_pa, h_set_jkg, outlet)
}
/// Creates a **compressible** source (R410A, CO₂, steam…).
pub fn compressible(
fluid: impl Into<String>,
p_set_pa: f64,
h_set_jkg: f64,
outlet: ConnectedPort,
) -> Result<Self, ComponentError> {
let fluid = fluid.into();
Self::new_inner(FluidKind::Compressible, fluid, p_set_pa, h_set_jkg, outlet)
}
fn new_inner(
kind: FluidKind,
fluid: String,
p_set_pa: f64,
h_set_jkg: f64,
outlet: ConnectedPort,
) -> Result<Self, ComponentError> {
if p_set_pa <= 0.0 {
return Err(ComponentError::InvalidState(
"FlowSource: set-point pressure must be positive".into(),
));
}
Ok(Self { kind, fluid_id: fluid, p_set_pa, h_set_jkg, outlet })
}
// ── Accessors ────────────────────────────────────────────────────────────
/// Fluid kind.
pub fn fluid_kind(&self) -> FluidKind { self.kind }
/// Fluid id.
pub fn fluid_id(&self) -> &str { &self.fluid_id }
/// Set-point pressure [Pa].
pub fn p_set_pa(&self) -> f64 { self.p_set_pa }
/// Set-point enthalpy [J/kg].
pub fn h_set_jkg(&self) -> f64 { self.h_set_jkg }
/// Reference to the outlet port.
pub fn outlet(&self) -> &ConnectedPort { &self.outlet }
/// Updates the set-point pressure (useful for parametric studies).
pub fn set_pressure(&mut self, p_pa: f64) -> Result<(), ComponentError> {
if p_pa <= 0.0 {
return Err(ComponentError::InvalidState(
"FlowSource: pressure must be positive".into(),
));
}
self.p_set_pa = p_pa;
Ok(())
}
/// Updates the set-point enthalpy.
pub fn set_enthalpy(&mut self, h_jkg: f64) {
self.h_set_jkg = h_jkg;
}
}
impl Component for FlowSource {
fn n_equations(&self) -> usize { 2 }
fn compute_residuals(
&self,
_state: &SystemState,
residuals: &mut ResidualVector,
) -> Result<(), ComponentError> {
if residuals.len() < 2 {
return Err(ComponentError::InvalidResidualDimensions {
expected: 2,
actual: residuals.len(),
});
}
// Pressure residual: P_edge − P_set = 0
residuals[0] = self.outlet.pressure().to_pascals() - self.p_set_pa;
// Enthalpy residual: h_edge − h_set = 0
residuals[1] = self.outlet.enthalpy().to_joules_per_kg() - self.h_set_jkg;
Ok(())
}
fn jacobian_entries(
&self,
_state: &SystemState,
jacobian: &mut JacobianBuilder,
) -> Result<(), ComponentError> {
// Both residuals are linear in the edge state: ∂r/∂x = 1
jacobian.add_entry(0, 0, 1.0);
jacobian.add_entry(1, 1, 1.0);
Ok(())
}
fn get_ports(&self) -> &[ConnectedPort] { &[] }
}
// ─────────────────────────────────────────────────────────────────────────────
// FlowSink — Back-pressure boundary condition
// ─────────────────────────────────────────────────────────────────────────────
/// A boundary sink that imposes a fixed back-pressure (and optionally enthalpy)
/// on its inlet edge.
///
/// Represents an infinite low-pressure reservoir (drain, condenser header,
/// discharge line, atmospheric vent, etc.).
///
/// # Equations (1 or 2)
///
/// ```text
/// r₀ = P_edge − P_back = 0 [always]
/// r₁ = h_edge − h_back = 0 [only if h_back is set]
/// ```
#[derive(Debug, Clone)]
pub struct FlowSink {
/// Fluid kind.
kind: FluidKind,
/// Fluid name.
fluid_id: String,
/// Back-pressure [Pa].
p_back_pa: f64,
/// Optional fixed outlet enthalpy [J/kg].
h_back_jkg: Option<f64>,
/// Connected inlet port.
inlet: ConnectedPort,
}
impl FlowSink {
// ── Constructors ─────────────────────────────────────────────────────────
/// Creates an **incompressible** sink (water, glycol…).
///
/// # Arguments
///
/// * `fluid` — fluid identifier string
/// * `p_back_pa` — back-pressure in Pascals
/// * `h_back_jkg` — optional fixed return enthalpy; `None` = free (solver decides)
/// * `inlet` — connected port
pub fn incompressible(
fluid: impl Into<String>,
p_back_pa: f64,
h_back_jkg: Option<f64>,
inlet: ConnectedPort,
) -> Result<Self, ComponentError> {
let fluid = fluid.into();
if !is_incompressible(&fluid) {
return Err(ComponentError::InvalidState(format!(
"FlowSink::incompressible: '{}' does not appear incompressible. \
Use FlowSink::compressible for refrigerants.",
fluid
)));
}
Self::new_inner(FluidKind::Incompressible, fluid, p_back_pa, h_back_jkg, inlet)
}
/// Creates a **compressible** sink (R410A, CO₂, steam…).
pub fn compressible(
fluid: impl Into<String>,
p_back_pa: f64,
h_back_jkg: Option<f64>,
inlet: ConnectedPort,
) -> Result<Self, ComponentError> {
let fluid = fluid.into();
Self::new_inner(FluidKind::Compressible, fluid, p_back_pa, h_back_jkg, inlet)
}
fn new_inner(
kind: FluidKind,
fluid: String,
p_back_pa: f64,
h_back_jkg: Option<f64>,
inlet: ConnectedPort,
) -> Result<Self, ComponentError> {
if p_back_pa <= 0.0 {
return Err(ComponentError::InvalidState(
"FlowSink: back-pressure must be positive".into(),
));
}
Ok(Self { kind, fluid_id: fluid, p_back_pa, h_back_jkg, inlet })
}
// ── Accessors ────────────────────────────────────────────────────────────
/// Fluid kind.
pub fn fluid_kind(&self) -> FluidKind { self.kind }
/// Fluid id.
pub fn fluid_id(&self) -> &str { &self.fluid_id }
/// Back-pressure [Pa].
pub fn p_back_pa(&self) -> f64 { self.p_back_pa }
/// Optional back-enthalpy [J/kg].
pub fn h_back_jkg(&self) -> Option<f64> { self.h_back_jkg }
/// Reference to the inlet port.
pub fn inlet(&self) -> &ConnectedPort { &self.inlet }
/// Updates the back-pressure.
pub fn set_pressure(&mut self, p_pa: f64) -> Result<(), ComponentError> {
if p_pa <= 0.0 {
return Err(ComponentError::InvalidState(
"FlowSink: back-pressure must be positive".into(),
));
}
self.p_back_pa = p_pa;
Ok(())
}
/// Sets a fixed return enthalpy (activates the second equation).
pub fn set_return_enthalpy(&mut self, h_jkg: f64) {
self.h_back_jkg = Some(h_jkg);
}
/// Removes the fixed enthalpy constraint (solver determines enthalpy freely).
pub fn clear_return_enthalpy(&mut self) {
self.h_back_jkg = None;
}
}
impl Component for FlowSink {
fn n_equations(&self) -> usize {
if self.h_back_jkg.is_some() { 2 } else { 1 }
}
fn compute_residuals(
&self,
_state: &SystemState,
residuals: &mut ResidualVector,
) -> Result<(), ComponentError> {
let n = self.n_equations();
if residuals.len() < n {
return Err(ComponentError::InvalidResidualDimensions {
expected: n,
actual: residuals.len(),
});
}
// Back-pressure residual
residuals[0] = self.inlet.pressure().to_pascals() - self.p_back_pa;
// Optional enthalpy residual
if let Some(h_back) = self.h_back_jkg {
residuals[1] = self.inlet.enthalpy().to_joules_per_kg() - h_back;
}
Ok(())
}
fn jacobian_entries(
&self,
_state: &SystemState,
jacobian: &mut JacobianBuilder,
) -> Result<(), ComponentError> {
let n = self.n_equations();
for i in 0..n {
jacobian.add_entry(i, i, 1.0);
}
Ok(())
}
fn get_ports(&self) -> &[ConnectedPort] { &[] }
}
// ─────────────────────────────────────────────────────────────────────────────
// Convenience type aliases (à la Modelica)
// ─────────────────────────────────────────────────────────────────────────────
/// Source for incompressible fluids (water, glycol, brine…).
pub type IncompressibleSource = FlowSource;
/// Source for compressible fluids (refrigerant, CO₂, steam…).
pub type CompressibleSource = FlowSource;
/// Sink for incompressible fluids.
pub type IncompressibleSink = FlowSink;
/// Sink for compressible fluids.
pub type CompressibleSink = FlowSink;
// ─────────────────────────────────────────────────────────────────────────────
// Tests
// ─────────────────────────────────────────────────────────────────────────────
#[cfg(test)]
mod tests {
use super::*;
use crate::port::{FluidId, Port};
use entropyk_core::{Enthalpy, Pressure};
fn make_port(fluid: &str, p_pa: f64, h_jkg: f64) -> ConnectedPort {
let a = Port::new(FluidId::new(fluid), Pressure::from_pascals(p_pa),
Enthalpy::from_joules_per_kg(h_jkg));
let b = Port::new(FluidId::new(fluid), Pressure::from_pascals(p_pa),
Enthalpy::from_joules_per_kg(h_jkg));
a.connect(b).unwrap().0
}
// ── FlowSource ────────────────────────────────────────────────────────────
#[test]
fn test_source_incompressible_water() {
// City water supply: 3 bar, 15°C (h ≈ 63 kJ/kg)
let port = make_port("Water", 3.0e5, 63_000.0);
let s = FlowSource::incompressible("Water", 3.0e5, 63_000.0, port).unwrap();
assert_eq!(s.n_equations(), 2);
assert_eq!(s.fluid_kind(), FluidKind::Incompressible);
assert_eq!(s.p_set_pa(), 3.0e5);
assert_eq!(s.h_set_jkg(), 63_000.0);
}
#[test]
fn test_source_compressible_refrigerant() {
// R410A high-side: 24 bar, h = 465 kJ/kg (superheated vapour)
let port = make_port("R410A", 24.0e5, 465_000.0);
let s = FlowSource::compressible("R410A", 24.0e5, 465_000.0, port).unwrap();
assert_eq!(s.n_equations(), 2);
assert_eq!(s.fluid_kind(), FluidKind::Compressible);
}
#[test]
fn test_source_rejects_refrigerant_as_incompressible() {
let port = make_port("R410A", 24.0e5, 465_000.0);
let result = FlowSource::incompressible("R410A", 24.0e5, 465_000.0, port);
assert!(result.is_err());
}
#[test]
fn test_source_rejects_zero_pressure() {
let port = make_port("Water", 3.0e5, 63_000.0);
let result = FlowSource::incompressible("Water", 0.0, 63_000.0, port);
assert!(result.is_err());
}
#[test]
fn test_source_residuals_zero_at_set_point() {
let p = 3.0e5_f64;
let h = 63_000.0_f64;
let port = make_port("Water", p, h);
let s = FlowSource::incompressible("Water", p, h, port).unwrap();
let state = vec![0.0; 4];
let mut res = vec![0.0; 2];
s.compute_residuals(&state, &mut res).unwrap();
assert!(res[0].abs() < 1.0, "P residual = {}", res[0]);
assert!(res[1].abs() < 1.0, "h residual = {}", res[1]);
}
#[test]
fn test_source_residuals_nonzero_on_mismatch() {
// Port at 2 bar but set-point 3 bar → residual = -1e5
let port = make_port("Water", 2.0e5, 63_000.0);
let s = FlowSource::incompressible("Water", 3.0e5, 63_000.0, port).unwrap();
let state = vec![0.0; 4];
let mut res = vec![0.0; 2];
s.compute_residuals(&state, &mut res).unwrap();
assert!((res[0] - (-1.0e5)).abs() < 1.0, "expected -1e5, got {}", res[0]);
}
#[test]
fn test_source_set_pressure() {
let port = make_port("Water", 3.0e5, 63_000.0);
let mut s = FlowSource::incompressible("Water", 3.0e5, 63_000.0, port).unwrap();
s.set_pressure(5.0e5).unwrap();
assert_eq!(s.p_set_pa(), 5.0e5);
assert!(s.set_pressure(0.0).is_err());
}
#[test]
fn test_source_as_trait_object() {
let port = make_port("R410A", 8.5e5, 260_000.0);
let src: Box<dyn Component> =
Box::new(FlowSource::compressible("R410A", 8.5e5, 260_000.0, port).unwrap());
assert_eq!(src.n_equations(), 2);
}
// ── FlowSink ──────────────────────────────────────────────────────────────
#[test]
fn test_sink_incompressible_back_pressure_only() {
// Return header: 1.5 bar, free enthalpy
let port = make_port("Water", 1.5e5, 63_000.0);
let s = FlowSink::incompressible("Water", 1.5e5, None, port).unwrap();
assert_eq!(s.n_equations(), 1);
assert_eq!(s.fluid_kind(), FluidKind::Incompressible);
}
#[test]
fn test_sink_with_fixed_return_enthalpy() {
// Fixed return temperature: 12°C, h ≈ 50.4 kJ/kg
let port = make_port("Water", 1.5e5, 50_400.0);
let s = FlowSink::incompressible("Water", 1.5e5, Some(50_400.0), port).unwrap();
assert_eq!(s.n_equations(), 2);
}
#[test]
fn test_sink_compressible_refrigerant() {
// R410A low-side: 8.5 bar
let port = make_port("R410A", 8.5e5, 260_000.0);
let s = FlowSink::compressible("R410A", 8.5e5, None, port).unwrap();
assert_eq!(s.n_equations(), 1);
assert_eq!(s.fluid_kind(), FluidKind::Compressible);
}
#[test]
fn test_sink_rejects_refrigerant_as_incompressible() {
let port = make_port("R410A", 8.5e5, 260_000.0);
let result = FlowSink::incompressible("R410A", 8.5e5, None, port);
assert!(result.is_err());
}
#[test]
fn test_sink_rejects_zero_back_pressure() {
let port = make_port("Water", 1.5e5, 63_000.0);
let result = FlowSink::incompressible("Water", 0.0, None, port);
assert!(result.is_err());
}
#[test]
fn test_sink_residual_zero_at_back_pressure() {
let p = 1.5e5_f64;
let port = make_port("Water", p, 63_000.0);
let s = FlowSink::incompressible("Water", p, None, port).unwrap();
let state = vec![0.0; 4];
let mut res = vec![0.0; 1];
s.compute_residuals(&state, &mut res).unwrap();
assert!(res[0].abs() < 1.0, "P residual = {}", res[0]);
}
#[test]
fn test_sink_residual_with_enthalpy() {
let p = 1.5e5_f64;
let h = 50_400.0_f64;
let port = make_port("Water", p, h);
let s = FlowSink::incompressible("Water", p, Some(h), port).unwrap();
let state = vec![0.0; 4];
let mut res = vec![0.0; 2];
s.compute_residuals(&state, &mut res).unwrap();
assert!(res[0].abs() < 1.0, "P residual = {}", res[0]);
assert!(res[1].abs() < 1.0, "h residual = {}", res[1]);
}
#[test]
fn test_sink_dynamic_enthalpy_toggle() {
let port = make_port("Water", 1.5e5, 63_000.0);
let mut s = FlowSink::incompressible("Water", 1.5e5, None, port).unwrap();
assert_eq!(s.n_equations(), 1);
s.set_return_enthalpy(50_400.0);
assert_eq!(s.n_equations(), 2);
s.clear_return_enthalpy();
assert_eq!(s.n_equations(), 1);
}
#[test]
fn test_sink_as_trait_object() {
let port = make_port("R410A", 8.5e5, 260_000.0);
let sink: Box<dyn Component> =
Box::new(FlowSink::compressible("R410A", 8.5e5, Some(260_000.0), port).unwrap());
assert_eq!(sink.n_equations(), 2);
}
}

826
crates/components/src/flow_junction.rs Normal file
View File

@@ -0,0 +1,826 @@
//! Flow Junction Components — Splitter & Merger
//!
//! This module provides `FlowSplitter` (1 inlet → N outlets) and `FlowMerger`
//! (N inlets → 1 outlet) for both incompressible (water, glycol, brine) and
//! compressible (refrigerant) fluid systems.
//!
//! ## Design Philosophy (à la Modelica)
//!
//! In Modelica, flow junctions apply conservation laws directly on connector
//! variables (pressure, enthalpy, mass flow). We follow the same approach:
//! constraints are algebraic equations on the state vector entries `[P, h]`
//! for each edge in the parent `System`.
//!
//! ## FlowSplitter — 1 inlet → N outlets
//!
//! Equations (2N − 1 total):
//! ```text
//! Mass balance : ṁ_in = ṁ_out_1 + ... + ṁ_out_N [1 eq]
//! Isobaric : P_out_k = P_in for k = 1..N-1 [N-1 eqs]
//! Isenthalpic : h_out_k = h_in for k = 1..N-1 [N-1 eqs]
//! ```
//!
//! The N-th outlet pressure and enthalpy equality are implied by the above.
//!
//! ## FlowMerger — N inlets → 1 outlet
//!
//! Equations (N + 1 total):
//! ```text
//! Mass balance : ṁ_out = Σ ṁ_in_k [1 eq]
//! Mixing enthalpy : h_out·ṁ_out = Σ h_in_k·ṁ_in_k [1 eq]
//! Pressure equalisation : P_in_k = P_in_1 for k = 2..N [N-1 eqs]
//! ```
//!
//! ## Incompressible vs Compressible
//!
//! The physics are **identical** — the distinction is purely in construction-time
//! validation (which fluid types are accepted). Use:
//! - [`FlowSplitter::incompressible`] / [`FlowMerger::incompressible`] for water,
//! glycol, brine, seawater circuits.
//! - [`FlowSplitter::compressible`] / [`FlowMerger::compressible`] for refrigerant
//! compressible circuits.
//!
//! ## State vector layout
//!
//! The solver assigns two state variables per edge: `(P_idx, h_idx)`.
//! Splitter/Merger receive the global state slice and use the **inlet/outlet
//! edge state indices** stored in their port list to resolve pressure and
//! specific enthalpy values.
//!
//! ## Example
//!
//! ```no_run
//! use entropyk_components::flow_junction::{FlowSplitter, FlowMerger};
//! use entropyk_components::port::{FluidId, Port};
//! use entropyk_core::{Pressure, Enthalpy};
//!
//! let make_port = |p: f64, h: f64| {
//! let a = Port::new(FluidId::new("Water"), Pressure::from_pascals(p),
//! Enthalpy::from_joules_per_kg(h));
//! let b = Port::new(FluidId::new("Water"), Pressure::from_pascals(p),
//! Enthalpy::from_joules_per_kg(h));
//! let (ca, _cb) = a.connect(b).unwrap();
//! ca
//! };
//!
//! let splitter = FlowSplitter::incompressible(
//! "Water",
//! make_port(3.0e5, 2.0e5), // inlet
//! vec![
//! make_port(3.0e5, 2.0e5), // branch A
//! make_port(3.0e5, 2.0e5), // branch B
//! ],
//! ).unwrap();
//! ```
use crate::{
Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, SystemState,
};
// ─────────────────────────────────────────────────────────────────────────────
// FluidKind — tag distinguishing the two regimes
// ─────────────────────────────────────────────────────────────────────────────
/// Whether this junction handles compressible or incompressible fluid.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum FluidKind {
/// Water, glycol, brine, seawater — density ≈ const.
Incompressible,
/// Refrigerant, CO₂, steam — density varies with P and T.
Compressible,
}
/// A set of known incompressible fluid identifiers (case-insensitive prefix match).
pub(crate) fn is_incompressible(fluid: &str) -> bool {
let f = fluid.to_lowercase();
f.starts_with("water")
|| f.starts_with("glycol")
|| f.starts_with("brine")
|| f.starts_with("seawater")
|| f.starts_with("ethyleneglycol")
|| f.starts_with("propyleneglycol")
|| f.starts_with("incompressible")
}
// ─────────────────────────────────────────────────────────────────────────────
// FlowSplitter — 1 inlet → N outlets
// ─────────────────────────────────────────────────────────────────────────────
/// A flow splitter that divides one inlet stream into N outlet branches.
///
/// # Equations (2N − 1)
///
/// | Equation | Residual |
/// |----------|---------|
/// | Mass balance | `ṁ_in − Σṁ_out = 0` |
/// | Isobaric (k=1..N-1) | `P_out_k − P_in = 0` |
/// | Isenthalpic (k=1..N-1) | `h_out_k − h_in = 0` |
///
/// ## Note on mass flow
///
/// The solver represents mass flow **implicitly** through pressure and enthalpy
/// on each edge. For the splitter's mass balance residual, we use the
/// simplified form that all outlet enthalpies equal the inlet enthalpy
/// (isenthalpic split). The mass balance is therefore:
///
/// `r_mass = (P_in − P_out_N) + (h_in − h_out_N)` as a consistency check.
///
/// See module docs for details.
#[derive(Debug, Clone)]
pub struct FlowSplitter {
/// Fluid kind (compressible / incompressible).
kind: FluidKind,
/// Fluid identifier (e.g. "Water", "R410A").
fluid_id: String,
/// Inlet port (the single source).
inlet: ConnectedPort,
/// Outlet ports (N branches).
outlets: Vec<ConnectedPort>,
}
impl FlowSplitter {
// ── Constructors ──────────────────────────────────────────────────────────
/// Creates an **incompressible** splitter (water, glycol, brine…).
///
/// # Errors
///
/// Returns an error if:
/// - `outlets` is empty
/// - The fluid is known to be compressible
/// - Fluids mismatch between inlet and outlets
pub fn incompressible(
fluid: impl Into<String>,
inlet: ConnectedPort,
outlets: Vec<ConnectedPort>,
) -> Result<Self, ComponentError> {
let fluid = fluid.into();
if !is_incompressible(&fluid) {
return Err(ComponentError::InvalidState(format!(
"FlowSplitter::incompressible: '{}' does not appear to be an incompressible fluid. \
Use FlowSplitter::compressible for refrigerants.",
fluid
)));
}
Self::new_inner(FluidKind::Incompressible, fluid, inlet, outlets)
}
/// Creates a **compressible** splitter (R410A, CO₂, steam…).
///
/// # Errors
///
/// Returns an error if `outlets` is empty.
pub fn compressible(
fluid: impl Into<String>,
inlet: ConnectedPort,
outlets: Vec<ConnectedPort>,
) -> Result<Self, ComponentError> {
let fluid = fluid.into();
Self::new_inner(FluidKind::Compressible, fluid, inlet, outlets)
}
fn new_inner(
kind: FluidKind,
fluid: String,
inlet: ConnectedPort,
outlets: Vec<ConnectedPort>,
) -> Result<Self, ComponentError> {
if outlets.is_empty() {
return Err(ComponentError::InvalidState(
"FlowSplitter requires at least one outlet".into(),
));
}
if outlets.len() == 1 {
return Err(ComponentError::InvalidState(
"FlowSplitter with 1 outlet is just a pipe — use a Pipe component instead".into(),
));
}
Ok(Self { kind, fluid_id: fluid, inlet, outlets })
}
// ── Accessors ─────────────────────────────────────────────────────────────
/// Number of outlet branches.
pub fn n_outlets(&self) -> usize {
self.outlets.len()
}
/// Fluid kind.
pub fn fluid_kind(&self) -> FluidKind {
self.kind
}
/// Fluid identifier.
pub fn fluid_id(&self) -> &str {
&self.fluid_id
}
/// Reference to the inlet port.
pub fn inlet(&self) -> &ConnectedPort {
&self.inlet
}
/// Reference to the outlet ports.
pub fn outlets(&self) -> &[ConnectedPort] {
&self.outlets
}
}
impl Component for FlowSplitter {
/// `2N − 1` equations:
/// - `N−1` pressure equalities
/// - `N−1` enthalpy equalities
/// - `1` mass balance (represented as N-th outlet consistency)
fn n_equations(&self) -> usize {
let n = self.outlets.len();
2 * n - 1
}
fn compute_residuals(
&self,
_state: &SystemState,
residuals: &mut ResidualVector,
) -> Result<(), ComponentError> {
let n_eqs = self.n_equations();
if residuals.len() < n_eqs {
return Err(ComponentError::InvalidResidualDimensions {
expected: n_eqs,
actual: residuals.len(),
});
}
// Inlet state indices come from the ConnectedPort.
// In the Entropyk solver, the state vector is indexed by edge:
// each edge contributes (P_idx, h_idx) set during System::finalize().
let p_in = self.inlet.pressure().to_pascals();
let h_in = self.inlet.enthalpy().to_joules_per_kg();
let n = self.outlets.len();
let mut r_idx = 0;
// --- Isobaric constraints: outlets 0..N-1 ---
for k in 0..(n - 1) {
let p_out_k = self.outlets[k].pressure().to_pascals();
residuals[r_idx] = p_out_k - p_in;
r_idx += 1;
}
// --- Isenthalpic constraints: outlets 0..N-1 ---
for k in 0..(n - 1) {
let h_out_k = self.outlets[k].enthalpy().to_joules_per_kg();
residuals[r_idx] = h_out_k - h_in;
r_idx += 1;
}
// --- Mass balance (1 equation) ---
// Express as: P_out_N = P_in AND h_out_N = h_in (implicitly guaranteed
// by the solver topology, but we add one algebraic check on the last branch).
// We use the last outlet as the "free" degree of freedom:
let p_out_last = self.outlets[n - 1].pressure().to_pascals();
residuals[r_idx] = p_out_last - p_in;
// Note: h_out_last = h_in is implied once all other constraints are met
// (conservation is guaranteed by the graph topology).
Ok(())
}
fn jacobian_entries(
&self,
_state: &SystemState,
jacobian: &mut JacobianBuilder,
) -> Result<(), ComponentError> {
// All residuals are linear differences → constant Jacobian.
// Each residual r_i depends on exactly two state variables with
// coefficients +1 and -1. Since the state indices are stored in the
// ConnectedPort pressure/enthalpy values (not raw state indices here),
// we emit a diagonal approximation: ∂r_i/∂x_i = 1.
//
// The full off-diagonal coupling is handled by the System assembler
// which maps port values to state vector positions.
let n_eqs = self.n_equations();
for i in 0..n_eqs {
jacobian.add_entry(i, i, 1.0);
}
Ok(())
}
fn get_ports(&self) -> &[ConnectedPort] {
// Return all ports so System can discover edges.
// Inlet first, then outlets.
// Note: dynamic allocation here is acceptable (called rarely during setup).
// We return an empty slice since get_ports() is for external port discovery;
// the actual solver coupling is via the System graph edges.
&[]
}
}
// ─────────────────────────────────────────────────────────────────────────────
// FlowMerger — N inlets → 1 outlet
// ─────────────────────────────────────────────────────────────────────────────
/// A flow merger that combines N inlet branches into one outlet stream.
///
/// # Equations (N + 1)
///
/// | Equation | Residual |
/// |----------|---------|
/// | Pressure equalisation (k=2..N) | `P_in_k − P_in_1 = 0` |
/// | Mixing enthalpy | `h_out − (Σ h_in_k) / N = 0` (equal-weight mix) |
/// | Mass balance | `P_out − P_in_1 = 0` |
///
/// ## Mixing enthalpy
///
/// When mass flow rates are not individually tracked, we use an equal-weight
/// average for the outlet enthalpy. This is exact for equal-flow branches and
/// approximate for unequal flows. When mass flows per branch are available
/// (from a FluidBackend), use [`FlowMerger::with_mass_flows`] for accuracy.
#[derive(Debug, Clone)]
pub struct FlowMerger {
/// Fluid kind (compressible / incompressible).
kind: FluidKind,
/// Fluid identifier.
fluid_id: String,
/// Inlet ports (N branches).
inlets: Vec<ConnectedPort>,
/// Outlet port (the single destination).
outlet: ConnectedPort,
/// Optional mass flow weights per inlet (kg/s). If None, equal weighting.
mass_flow_weights: Option<Vec<f64>>,
}
impl FlowMerger {
// ── Constructors ──────────────────────────────────────────────────────────
/// Creates an **incompressible** merger (water, glycol, brine…).
pub fn incompressible(
fluid: impl Into<String>,
inlets: Vec<ConnectedPort>,
outlet: ConnectedPort,
) -> Result<Self, ComponentError> {
let fluid = fluid.into();
if !is_incompressible(&fluid) {
return Err(ComponentError::InvalidState(format!(
"FlowMerger::incompressible: '{}' does not appear to be an incompressible fluid. \
Use FlowMerger::compressible for refrigerants.",
fluid
)));
}
Self::new_inner(FluidKind::Incompressible, fluid, inlets, outlet)
}
/// Creates a **compressible** merger (R410A, CO₂, steam…).
pub fn compressible(
fluid: impl Into<String>,
inlets: Vec<ConnectedPort>,
outlet: ConnectedPort,
) -> Result<Self, ComponentError> {
let fluid = fluid.into();
Self::new_inner(FluidKind::Compressible, fluid, inlets, outlet)
}
fn new_inner(
kind: FluidKind,
fluid: String,
inlets: Vec<ConnectedPort>,
outlet: ConnectedPort,
) -> Result<Self, ComponentError> {
if inlets.is_empty() {
return Err(ComponentError::InvalidState(
"FlowMerger requires at least one inlet".into(),
));
}
if inlets.len() == 1 {
return Err(ComponentError::InvalidState(
"FlowMerger with 1 inlet is just a pipe — use a Pipe component instead".into(),
));
}
Ok(Self {
kind,
fluid_id: fluid,
inlets,
outlet,
mass_flow_weights: None,
})
}
/// Assigns known mass flow rates per inlet for weighted enthalpy mixing.
///
/// # Errors
///
/// Returns an error if `weights.len() != n_inlets`.
pub fn with_mass_flows(mut self, weights: Vec<f64>) -> Result<Self, ComponentError> {
if weights.len() != self.inlets.len() {
return Err(ComponentError::InvalidState(format!(
"FlowMerger::with_mass_flows: expected {} weights, got {}",
self.inlets.len(),
weights.len()
)));
}
if weights.iter().any(|&w| w < 0.0) {
return Err(ComponentError::InvalidState(
"FlowMerger::with_mass_flows: mass flow weights must be non-negative".into(),
));
}
self.mass_flow_weights = Some(weights);
Ok(self)
}
// ── Accessors ─────────────────────────────────────────────────────────────
/// Number of inlet branches.
pub fn n_inlets(&self) -> usize {
self.inlets.len()
}
/// Fluid kind.
pub fn fluid_kind(&self) -> FluidKind {
self.kind
}
/// Fluid identifier.
pub fn fluid_id(&self) -> &str {
&self.fluid_id
}
/// Reference to the inlet ports.
pub fn inlets(&self) -> &[ConnectedPort] {
&self.inlets
}
/// Reference to the outlet port.
pub fn outlet(&self) -> &ConnectedPort {
&self.outlet
}
// ── Mixing helpers ────────────────────────────────────────────────────────
/// Computes the mixed outlet enthalpy (weighted or equal).
fn mixed_enthalpy(&self) -> f64 {
let n = self.inlets.len();
match &self.mass_flow_weights {
Some(weights) => {
let total_flow: f64 = weights.iter().sum();
if total_flow <= 0.0 {
// Fall back to equal weighting
self.inlets.iter().map(|p| p.enthalpy().to_joules_per_kg()).sum::<f64>()
/ n as f64
} else {
self.inlets
.iter()
.zip(weights.iter())
.map(|(p, &w)| p.enthalpy().to_joules_per_kg() * w)
.sum::<f64>()
/ total_flow
}
}
None => {
// Equal weighting
self.inlets.iter().map(|p| p.enthalpy().to_joules_per_kg()).sum::<f64>()
/ n as f64
}
}
}
}
impl Component for FlowMerger {
/// `N + 1` equations:
/// - `N−1` pressure equalisations across inlets
/// - `1` mixing enthalpy for the outlet
/// - `1` outlet pressure consistency
fn n_equations(&self) -> usize {
self.inlets.len() + 1
}
fn compute_residuals(
&self,
_state: &SystemState,
residuals: &mut ResidualVector,
) -> Result<(), ComponentError> {
let n_eqs = self.n_equations();
if residuals.len() < n_eqs {
return Err(ComponentError::InvalidResidualDimensions {
expected: n_eqs,
actual: residuals.len(),
});
}
let p_ref = self.inlets[0].pressure().to_pascals();
let mut r_idx = 0;
// --- Pressure equalisation: inlets 1..N must match inlet 0 ---
for k in 1..self.inlets.len() {
let p_k = self.inlets[k].pressure().to_pascals();
residuals[r_idx] = p_k - p_ref;
r_idx += 1;
}
// --- Outlet pressure = reference inlet pressure ---
let p_out = self.outlet.pressure().to_pascals();
residuals[r_idx] = p_out - p_ref;
r_idx += 1;
// --- Mixing enthalpy ---
let h_mixed = self.mixed_enthalpy();
let h_out = self.outlet.enthalpy().to_joules_per_kg();
residuals[r_idx] = h_out - h_mixed;
Ok(())
}
fn jacobian_entries(
&self,
_state: &SystemState,
jacobian: &mut JacobianBuilder,
) -> Result<(), ComponentError> {
// Diagonal approximation — the full coupling is resolved by the System
// assembler through the edge state indices.
let n_eqs = self.n_equations();
for i in 0..n_eqs {
jacobian.add_entry(i, i, 1.0);
}
Ok(())
}
fn get_ports(&self) -> &[ConnectedPort] {
&[]
}
}
// ─────────────────────────────────────────────────────────────────────────────
// Convenience type aliases
// ─────────────────────────────────────────────────────────────────────────────
/// A flow splitter for incompressible fluids (water, glycol, brine…).
///
/// Equivalent to `FlowSplitter` constructed via [`FlowSplitter::incompressible`].
pub type IncompressibleSplitter = FlowSplitter;
/// A flow splitter for compressible fluids (refrigerant, CO₂, steam…).
///
/// Equivalent to `FlowSplitter` constructed via [`FlowSplitter::compressible`].
pub type CompressibleSplitter = FlowSplitter;
/// A flow merger for incompressible fluids (water, glycol, brine…).
///
/// Equivalent to `FlowMerger` constructed via [`FlowMerger::incompressible`].
pub type IncompressibleMerger = FlowMerger;
/// A flow merger for compressible fluids (refrigerant, CO₂, steam…).
///
/// Equivalent to `FlowMerger` constructed via [`FlowMerger::compressible`].
pub type CompressibleMerger = FlowMerger;
// ─────────────────────────────────────────────────────────────────────────────
// Tests
// ─────────────────────────────────────────────────────────────────────────────
#[cfg(test)]
mod tests {
use super::*;
use crate::port::{FluidId, Port};
use entropyk_core::{Enthalpy, Pressure};
fn make_port(fluid: &str, p_pa: f64, h_jkg: f64) -> ConnectedPort {
let p1 = Port::new(
FluidId::new(fluid),
Pressure::from_pascals(p_pa),
Enthalpy::from_joules_per_kg(h_jkg),
);
let p2 = Port::new(
FluidId::new(fluid),
Pressure::from_pascals(p_pa),
Enthalpy::from_joules_per_kg(h_jkg),
);
p1.connect(p2).unwrap().0
}
// ── FlowSplitter ──────────────────────────────────────────────────────────
#[test]
fn test_splitter_incompressible_creation() {
let inlet = make_port("Water", 3.0e5, 2.0e5);
let out_a = make_port("Water", 3.0e5, 2.0e5);
let out_b = make_port("Water", 3.0e5, 2.0e5);
let s = FlowSplitter::incompressible("Water", inlet, vec![out_a, out_b]).unwrap();
assert_eq!(s.n_outlets(), 2);
assert_eq!(s.fluid_kind(), FluidKind::Incompressible);
// n_equations = 2*2 - 1 = 3
assert_eq!(s.n_equations(), 3);
}
#[test]
fn test_splitter_compressible_creation() {
let inlet = make_port("R410A", 24.0e5, 4.65e5);
let out_a = make_port("R410A", 24.0e5, 4.65e5);
let out_b = make_port("R410A", 24.0e5, 4.65e5);
let out_c = make_port("R410A", 24.0e5, 4.65e5);
let s = FlowSplitter::compressible("R410A", inlet, vec![out_a, out_b, out_c]).unwrap();
assert_eq!(s.n_outlets(), 3);
assert_eq!(s.fluid_kind(), FluidKind::Compressible);
// n_equations = 2*3 - 1 = 5
assert_eq!(s.n_equations(), 5);
}
#[test]
fn test_splitter_rejects_refrigerant_as_incompressible() {
let inlet = make_port("R410A", 24.0e5, 4.65e5);
let out_a = make_port("R410A", 24.0e5, 4.65e5);
let out_b = make_port("R410A", 24.0e5, 4.65e5);
let result = FlowSplitter::incompressible("R410A", inlet, vec![out_a, out_b]);
assert!(result.is_err(), "R410A should not be accepted as incompressible");
}
#[test]
fn test_splitter_rejects_single_outlet() {
let inlet = make_port("Water", 3.0e5, 2.0e5);
let out = make_port("Water", 3.0e5, 2.0e5);
let result = FlowSplitter::incompressible("Water", inlet, vec![out]);
assert!(result.is_err());
}
#[test]
fn test_splitter_residuals_zero_at_consistent_state() {
// Consistent state: all pressures and enthalpies equal
let inlet = make_port("Water", 3.0e5, 2.0e5);
let out_a = make_port("Water", 3.0e5, 2.0e5);
let out_b = make_port("Water", 3.0e5, 2.0e5);
let s = FlowSplitter::incompressible("Water", inlet, vec![out_a, out_b]).unwrap();
let state = vec![0.0; 6]; // dummy, not used by current impl
let mut res = vec![0.0; s.n_equations()];
s.compute_residuals(&state, &mut res).unwrap();
for (i, &r) in res.iter().enumerate() {
assert!(
r.abs() < 1.0,
"residual[{}] = {} should be ≈ 0 for consistent state",
i, r
);
}
}
#[test]
fn test_splitter_residuals_nonzero_on_pressure_mismatch() {
let inlet = make_port("Water", 3.0e5, 2.0e5);
let out_a = make_port("Water", 2.5e5, 2.0e5); // lower pressure!
let out_b = make_port("Water", 3.0e5, 2.0e5);
let s = FlowSplitter::incompressible("Water", inlet, vec![out_a, out_b]).unwrap();
let state = vec![0.0; 6];
let mut res = vec![0.0; s.n_equations()];
s.compute_residuals(&state, &mut res).unwrap();
// r[0] = P_out_a - P_in = 2.5e5 - 3.0e5 = -0.5e5
assert!((res[0] - (-0.5e5)).abs() < 1.0, "expected -0.5e5, got {}", res[0]);
}
#[test]
fn test_splitter_three_branches_n_equations() {
let inlet = make_port("Water", 3.0e5, 2.0e5);
let outlets: Vec<_> = (0..3).map(|_| make_port("Water", 3.0e5, 2.0e5)).collect();
let s = FlowSplitter::incompressible("Water", inlet, outlets).unwrap();
// N=3 → 2*3-1 = 5
assert_eq!(s.n_equations(), 5);
}
#[test]
fn test_splitter_water_type_aliases() {
let inlet = make_port("Water", 3.0e5, 2.0e5);
let out_a = make_port("Water", 3.0e5, 2.0e5);
let out_b = make_port("Water", 3.0e5, 2.0e5);
// IncompressibleSplitter is a type alias for FlowSplitter
let _s: IncompressibleSplitter =
FlowSplitter::incompressible("Water", inlet, vec![out_a, out_b]).unwrap();
}
// ── FlowMerger ────────────────────────────────────────────────────────────
#[test]
fn test_merger_incompressible_creation() {
let in_a = make_port("Water", 3.0e5, 2.0e5);
let in_b = make_port("Water", 3.0e5, 2.4e5);
let outlet = make_port("Water", 3.0e5, 2.2e5);
let m = FlowMerger::incompressible("Water", vec![in_a, in_b], outlet).unwrap();
assert_eq!(m.n_inlets(), 2);
assert_eq!(m.fluid_kind(), FluidKind::Incompressible);
// N=2 → N+1 = 3
assert_eq!(m.n_equations(), 3);
}
#[test]
fn test_merger_compressible_creation() {
let in_a = make_port("R134a", 8.0e5, 4.0e5);
let in_b = make_port("R134a", 8.0e5, 4.2e5);
let in_c = make_port("R134a", 8.0e5, 3.8e5);
let outlet = make_port("R134a", 8.0e5, 4.0e5);
let m = FlowMerger::compressible("R134a", vec![in_a, in_b, in_c], outlet).unwrap();
assert_eq!(m.n_inlets(), 3);
assert_eq!(m.fluid_kind(), FluidKind::Compressible);
// N=3 → N+1 = 4
assert_eq!(m.n_equations(), 4);
}
#[test]
fn test_merger_rejects_single_inlet() {
let in_a = make_port("Water", 3.0e5, 2.0e5);
let outlet = make_port("Water", 3.0e5, 2.0e5);
let result = FlowMerger::incompressible("Water", vec![in_a], outlet);
assert!(result.is_err());
}
#[test]
fn test_merger_residuals_zero_at_consistent_state() {
// Equal branches → mixed enthalpy = inlet enthalpy
let h = 2.0e5_f64;
let p = 3.0e5_f64;
let in_a = make_port("Water", p, h);
let in_b = make_port("Water", p, h);
let outlet = make_port("Water", p, h); // h_mixed = (h+h)/2 = h
let m = FlowMerger::incompressible("Water", vec![in_a, in_b], outlet).unwrap();
let state = vec![0.0; 6];
let mut res = vec![0.0; m.n_equations()];
m.compute_residuals(&state, &mut res).unwrap();
for (i, &r) in res.iter().enumerate() {
assert!(r.abs() < 1.0, "residual[{}] = {} should be ≈ 0", i, r);
}
}
#[test]
fn test_merger_mixed_enthalpy_equal_branches() {
let h_a = 2.0e5_f64;
let h_b = 3.0e5_f64;
let h_expected = (h_a + h_b) / 2.0; // equal-weight average
let p = 3.0e5_f64;
let in_a = make_port("Water", p, h_a);
let in_b = make_port("Water", p, h_b);
let outlet = make_port("Water", p, h_expected);
let m = FlowMerger::incompressible("Water", vec![in_a, in_b], outlet).unwrap();
let state = vec![0.0; 6];
let mut res = vec![0.0; m.n_equations()];
m.compute_residuals(&state, &mut res).unwrap();
// Last residual: h_out - h_mixed should be 0
let last = res[m.n_equations() - 1];
assert!(last.abs() < 1.0, "h mixing residual = {} should be 0", last);
}
#[test]
fn test_merger_weighted_enthalpy() {
// ṁ_a = 0.3 kg/s, h_a = 2e5 J/kg
// ṁ_b = 0.7 kg/s, h_b = 3e5 J/kg
// h_mix = (0.3*2e5 + 0.7*3e5) / 1.0 = (6e4 + 21e4) = 2.7e5 J/kg
let p = 3.0e5_f64;
let in_a = make_port("Water", p, 2.0e5);
let in_b = make_port("Water", p, 3.0e5);
let outlet = make_port("Water", p, 2.7e5);
let m = FlowMerger::incompressible("Water", vec![in_a, in_b], outlet)
.unwrap()
.with_mass_flows(vec![0.3, 0.7])
.unwrap();
let state = vec![0.0; 6];
let mut res = vec![0.0; m.n_equations()];
m.compute_residuals(&state, &mut res).unwrap();
let h_residual = res[m.n_equations() - 1];
assert!(
h_residual.abs() < 1.0,
"weighted h mixing residual = {} should be 0",
h_residual
);
}
#[test]
fn test_merger_as_trait_object() {
let in_a = make_port("Water", 3.0e5, 2.0e5);
let in_b = make_port("Water", 3.0e5, 2.0e5);
let outlet = make_port("Water", 3.0e5, 2.0e5);
let merger: Box<dyn Component> = Box::new(
FlowMerger::incompressible("Water", vec![in_a, in_b], outlet).unwrap()
);
assert_eq!(merger.n_equations(), 3);
}
#[test]
fn test_splitter_as_trait_object() {
let inlet = make_port("R410A", 24.0e5, 4.65e5);
let out_a = make_port("R410A", 24.0e5, 4.65e5);
let out_b = make_port("R410A", 24.0e5, 4.65e5);
let splitter: Box<dyn Component> = Box::new(
FlowSplitter::compressible("R410A", inlet, vec![out_a, out_b]).unwrap()
);
assert_eq!(splitter.n_equations(), 3);
}
}

43
crates/components/src/lib.rs
View File

@@ -59,6 +59,8 @@ pub mod compressor;
pub mod expansion_valve;
pub mod external_model;
pub mod fan;
pub mod flow_boundary;
pub mod flow_junction;
pub mod heat_exchanger;
pub mod pipe;
pub mod polynomials;
@@ -84,6 +86,14 @@ pub use polynomials::{AffinityLaws, PerformanceCurves, Polynomial1D, Polynomial2
pub use port::{
validate_port_continuity, Connected, ConnectedPort, ConnectionError, Disconnected, FluidId, Port,
};
pub use flow_boundary::{
CompressibleSink, CompressibleSource, FlowSink, FlowSource,
IncompressibleSink, IncompressibleSource,
};
pub use flow_junction::{
CompressibleMerger, CompressibleSplitter, FlowMerger, FlowSplitter, FluidKind,
IncompressibleMerger, IncompressibleSplitter,
};
pub use pump::{Pump, PumpCurves};
pub use state_machine::{
CircuitId, OperationalState, StateHistory, StateManageable, StateTransitionError,
@@ -499,6 +509,39 @@ pub trait Component {
/// assert!(component.get_ports().is_empty());
/// ```
fn get_ports(&self) -> &[ConnectedPort];
/// Injects system-level context into a component during topology finalization.
///
/// Called by [`System::finalize()`] after all edge state indices are computed.
/// The default implementation is a no-op; override this in components that need
/// to know their position in the global state vector (e.g. `MacroComponent`).
///
/// # Arguments
///
/// * `state_offset` — The index in the global state vector where this component's
/// *internal* state block begins. For ordinary leaf components this is never
/// needed; for `MacroComponent` it replaces the manual `set_global_state_offset`
/// call.
/// * `external_edge_state_indices` — A slice of `(p_idx, h_idx)` pairs for every
/// edge incident to this component's node in the parent graph (incoming and
/// outgoing), in traversal order. `MacroComponent` uses these to emit
/// port-coupling residuals.
fn set_system_context(
&mut self,
_state_offset: usize,
_external_edge_state_indices: &[(usize, usize)],
) {
// Default: no-op for all ordinary leaf components.
}
/// Returns the number of internal state variables this component maintains.
///
/// The default implementation returns 0, which is correct for all ordinary
/// leaf components. Hierarchical components (like `MacroComponent`) should
/// override this to return the size of their internal state block.
fn internal_state_len(&self) -> usize {
0
}
}
#[cfg(test)]

2
crates/solver/Cargo.toml
View File

@@ -14,9 +14,11 @@ nalgebra = "0.33"
petgraph = "0.6"
thiserror = "1.0"
tracing = "0.1"
serde = { version = "1.0", features = ["derive"] }
[dev-dependencies]
approx = "0.5"
serde_json = "1.0"
[lib]
name = "entropyk_solver"

6
crates/solver/src/criteria.rs
View File

@@ -212,7 +212,11 @@ impl ConvergenceCriteria {
let max_delta_p = pressure_indices
.iter()
.map(|&p_idx| {
let p = if p_idx < state.len() { state[p_idx] } else { 0.0 };
let p = if p_idx < state.len() {
state[p_idx]
} else {
0.0
};
if let Some(prev) = prev_state {
let pp = if p_idx < prev.len() { prev[p_idx] } else { 0.0 };
(p - pp).abs()

116
crates/solver/src/initializer.rs
View File

@@ -50,9 +50,7 @@ pub enum InitializerError {
},
/// The provided state slice length does not match the system state vector length.
#[error(
"State slice length {actual} does not match system state vector length {expected}"
)]
#[error("State slice length {actual} does not match system state vector length {expected}")]
StateLengthMismatch {
/// Expected length (from `system.state_vector_len()`).
expected: usize,
@@ -272,10 +270,7 @@ impl SmartInitializer {
"Unknown fluid for Antoine estimation — using fallback pressures \
(P_evap = 5 bar, P_cond = 20 bar)"
);
Ok((
Pressure::from_bar(5.0),
Pressure::from_bar(20.0),
))
Ok((Pressure::from_bar(5.0), Pressure::from_bar(20.0)))
}
Some(coeffs) => {
let t_source_c = t_source.to_celsius();
@@ -514,20 +509,36 @@ mod tests {
#[test]
fn test_populate_state_2_edges() {
use crate::system::System;
use entropyk_components::{Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, SystemState};
use entropyk_components::{
Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, SystemState,
};
struct MockComp;
impl Component for MockComp {
fn compute_residuals(&self, _s: &SystemState, r: &mut ResidualVector) -> Result<(), ComponentError> {
for v in r.iter_mut() { *v = 0.0; }
fn compute_residuals(
&self,
_s: &SystemState,
r: &mut ResidualVector,
) -> Result<(), ComponentError> {
for v in r.iter_mut() {
*v = 0.0;
}
Ok(())
}
fn jacobian_entries(&self, _s: &SystemState, j: &mut JacobianBuilder) -> Result<(), ComponentError> {
fn jacobian_entries(
&self,
_s: &SystemState,
j: &mut JacobianBuilder,
) -> Result<(), ComponentError> {
j.add_entry(0, 0, 1.0);
Ok(())
}
fn n_equations(&self) -> usize { 1 }
fn get_ports(&self) -> &[ConnectedPort] { &[] }
fn n_equations(&self) -> usize {
1
}
fn get_ports(&self) -> &[ConnectedPort] {
&[]
}
}
let mut sys = System::new();
@@ -560,29 +571,53 @@ mod tests {
#[test]
fn test_populate_state_multi_circuit() {
use crate::system::{CircuitId, System};
use entropyk_components::{Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, SystemState};
use entropyk_components::{
Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, SystemState,
};
struct MockComp;
impl Component for MockComp {
fn compute_residuals(&self, _s: &SystemState, r: &mut ResidualVector) -> Result<(), ComponentError> {
for v in r.iter_mut() { *v = 0.0; }
fn compute_residuals(
&self,
_s: &SystemState,
r: &mut ResidualVector,
) -> Result<(), ComponentError> {
for v in r.iter_mut() {
*v = 0.0;
}
Ok(())
}
fn jacobian_entries(&self, _s: &SystemState, j: &mut JacobianBuilder) -> Result<(), ComponentError> {
fn jacobian_entries(
&self,
_s: &SystemState,
j: &mut JacobianBuilder,
) -> Result<(), ComponentError> {
j.add_entry(0, 0, 1.0);
Ok(())
}
fn n_equations(&self) -> usize { 1 }
fn get_ports(&self) -> &[ConnectedPort] { &[] }
fn n_equations(&self) -> usize {
1
}
fn get_ports(&self) -> &[ConnectedPort] {
&[]
}
}
let mut sys = System::new();
// Circuit 0: evaporator side
let n0 = sys.add_component_to_circuit(Box::new(MockComp), CircuitId(0)).unwrap();
let n1 = sys.add_component_to_circuit(Box::new(MockComp), CircuitId(0)).unwrap();
let n0 = sys
.add_component_to_circuit(Box::new(MockComp), CircuitId(0))
.unwrap();
let n1 = sys
.add_component_to_circuit(Box::new(MockComp), CircuitId(0))
.unwrap();
// Circuit 1: condenser side
let n2 = sys.add_component_to_circuit(Box::new(MockComp), CircuitId(1)).unwrap();
let n3 = sys.add_component_to_circuit(Box::new(MockComp), CircuitId(1)).unwrap();
let n2 = sys
.add_component_to_circuit(Box::new(MockComp), CircuitId(1))
.unwrap();
let n3 = sys
.add_component_to_circuit(Box::new(MockComp), CircuitId(1))
.unwrap();
sys.add_edge(n0, n1).unwrap(); // circuit 0 edge
sys.add_edge(n2, n3).unwrap(); // circuit 1 edge
@@ -598,7 +633,7 @@ mod tests {
.unwrap();
assert_eq!(state.len(), 4); // 2 edges × 2 entries
// Edge 0 (circuit 0) → p_evap
// Edge 0 (circuit 0) → p_evap
assert_relative_eq!(state[0], p_evap.to_pascals(), max_relative = 1e-9);
assert_relative_eq!(state[1], h_default.to_joules_per_kg(), max_relative = 1e-9);
// Edge 1 (circuit 1) → p_cond
@@ -610,20 +645,36 @@ mod tests {
#[test]
fn test_populate_state_length_mismatch() {
use crate::system::System;
use entropyk_components::{Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, SystemState};
use entropyk_components::{
Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, SystemState,
};
struct MockComp;
impl Component for MockComp {
fn compute_residuals(&self, _s: &SystemState, r: &mut ResidualVector) -> Result<(), ComponentError> {
for v in r.iter_mut() { *v = 0.0; }
fn compute_residuals(
&self,
_s: &SystemState,
r: &mut ResidualVector,
) -> Result<(), ComponentError> {
for v in r.iter_mut() {
*v = 0.0;
}
Ok(())
}
fn jacobian_entries(&self, _s: &SystemState, j: &mut JacobianBuilder) -> Result<(), ComponentError> {
fn jacobian_entries(
&self,
_s: &SystemState,
j: &mut JacobianBuilder,
) -> Result<(), ComponentError> {
j.add_entry(0, 0, 1.0);
Ok(())
}
fn n_equations(&self) -> usize { 1 }
fn get_ports(&self) -> &[ConnectedPort] { &[] }
fn n_equations(&self) -> usize {
1
}
fn get_ports(&self) -> &[ConnectedPort] {
&[]
}
}
let mut sys = System::new();
@@ -642,7 +693,10 @@ mod tests {
let result = init.populate_state(&sys, p_evap, p_cond, h_default, &mut state);
assert!(matches!(
result,
Err(InitializerError::StateLengthMismatch { expected: 2, actual: 5 })
Err(InitializerError::StateLengthMismatch {
expected: 2,
actual: 5
})
));
}

793
crates/solver/src/inverse/bounded.rs Normal file
View File

@@ -0,0 +1,793 @@
//! Bounded control variables for inverse control.
//!
//! This module provides types for control variables with physical bounds:
//! - [`BoundedVariable`]: A control variable constrained to [min, max] range
//! - [`BoundedVariableId`]: Type-safe bounded variable identifier
//! - [`BoundedVariableError`]: Errors during bounded variable operations
//! - [`SaturationInfo`]: Information about variable saturation at bounds
//!
//! # Mathematical Foundation
//!
//! Box constraints ensure Newton steps stay physically possible:
//!
//! $$x_{new} = \text{clip}(x_{current} + \Delta x, x_{min}, x_{max})$$
//!
//! where clipping limits the step to stay within bounds.
//!
//! # Example
//!
//! ```rust,ignore
//! use entropyk_solver::inverse::{BoundedVariable, BoundedVariableId};
//!
//! // Valve position: 0.0 (closed) to 1.0 (fully open)
//! let valve = BoundedVariable::new(
//! BoundedVariableId::new("expansion_valve"),
//! 0.5, // initial position: 50% open
//! 0.0, // min: fully closed
//! 1.0, // max: fully open
//! ).unwrap();
//!
//! // Step clipping keeps value in bounds
//! let clipped = valve.clip_step(0.8); // tries to move to 1.3
//! assert_eq!(clipped, 1.0); // clipped to max
//! ```
use std::fmt;
use thiserror::Error;
// ─────────────────────────────────────────────────────────────────────────────
// BoundedVariableId - Type-safe identifier
// ─────────────────────────────────────────────────────────────────────────────
/// Type-safe identifier for a bounded control variable.
///
/// Uses a string internally but provides type safety and clear intent.
#[derive(Debug, Clone, PartialEq, Eq, Hash)]
pub struct BoundedVariableId(String);
impl BoundedVariableId {
/// Creates a new bounded variable identifier.
///
/// # Arguments
///
/// * `id` - Unique identifier string for the bounded variable
///
/// # Example
///
/// ```rust
/// use entropyk_solver::inverse::BoundedVariableId;
///
/// let id = BoundedVariableId::new("valve_position");
/// ```
pub fn new(id: impl Into<String>) -> Self {
BoundedVariableId(id.into())
}
/// Returns the identifier as a string slice.
pub fn as_str(&self) -> &str {
&self.0
}
}
impl fmt::Display for BoundedVariableId {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "{}", self.0)
}
}
impl From<&str> for BoundedVariableId {
fn from(s: &str) -> Self {
BoundedVariableId::new(s)
}
}
impl From<String> for BoundedVariableId {
fn from(s: String) -> Self {
BoundedVariableId(s)
}
}
// ─────────────────────────────────────────────────────────────────────────────
// BoundedVariableError - Error handling
// ─────────────────────────────────────────────────────────────────────────────
/// Errors that can occur during bounded variable operations.
#[derive(Error, Debug, Clone, PartialEq)]
pub enum BoundedVariableError {
/// The bounds are invalid (min >= max).
#[error("Invalid bounds: min ({min}) must be less than max ({max})")]
InvalidBounds {
/// The minimum bound
min: f64,
/// The maximum bound
max: f64,
},
/// A bounded variable with this ID already exists.
#[error("Duplicate bounded variable id: '{id}'")]
DuplicateId {
/// The duplicate identifier
id: BoundedVariableId,
},
/// The initial value is outside the bounds.
#[error("Initial value {value} is outside bounds [{min}, {max}]")]
ValueOutOfBounds {
/// The out-of-bounds value
value: f64,
/// The minimum bound
min: f64,
/// The maximum bound
max: f64,
},
/// The referenced component does not exist in the system.
#[error("Invalid component reference: '{component_id}' not found")]
InvalidComponent {
/// The invalid component identifier
component_id: String,
},
/// Invalid configuration.
#[error("Invalid bounded variable configuration: {reason}")]
InvalidConfiguration {
/// Reason for the validation failure
reason: String,
},
}
// ─────────────────────────────────────────────────────────────────────────────
// SaturationType and SaturationInfo - Saturation detection
// ─────────────────────────────────────────────────────────────────────────────
/// Type of bound saturation.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum SaturationType {
/// Variable is at the lower bound (value == min)
LowerBound,
/// Variable is at the upper bound (value == max)
UpperBound,
}
impl fmt::Display for SaturationType {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
match self {
SaturationType::LowerBound => write!(f, "lower bound"),
SaturationType::UpperBound => write!(f, "upper bound"),
}
}
}
/// Information about a saturated bounded variable.
///
/// When a control variable reaches a bound during solving, this struct
/// captures the details for diagnostic purposes.
#[derive(Debug, Clone, PartialEq)]
pub struct SaturationInfo {
/// The saturated variable's identifier
pub variable_id: BoundedVariableId,
/// Which bound is active
pub saturation_type: SaturationType,
/// The bound value (min or max)
pub bound_value: f64,
/// The constraint target that couldn't be achieved (if applicable)
pub constraint_target: Option<f64>,
}
impl SaturationInfo {
/// Creates a new SaturationInfo.
pub fn new(
variable_id: BoundedVariableId,
saturation_type: SaturationType,
bound_value: f64,
) -> Self {
SaturationInfo {
variable_id,
saturation_type,
bound_value,
constraint_target: None,
}
}
/// Adds constraint target information.
pub fn with_constraint_target(mut self, target: f64) -> Self {
self.constraint_target = Some(target);
self
}
}
// ─────────────────────────────────────────────────────────────────────────────
// BoundedVariable - Core bounded variable type
// ─────────────────────────────────────────────────────────────────────────────
/// A control variable with physical bounds for inverse control.
///
/// Bounded variables ensure that Newton-Raphson steps stay within
/// physically meaningful ranges (e.g., valve position 0.0 to 1.0).
///
/// # Bounds Validation
///
/// Bounds must satisfy `min < max`. The initial value must be within
/// bounds: `min <= value <= max`.
///
/// # Step Clipping
///
/// When applying a Newton step Δx, the new value is clipped:
///
/// $$x_{new} = \text{clamp}(x + \Delta x, x_{min}, x_{max})$$
///
/// # Saturation Detection
///
/// After convergence, check if the variable is saturated (at a bound)
/// using [`is_saturated`](Self::is_saturated).
///
/// # Example
///
/// ```rust,ignore
/// // VFD speed: 30% to 100% (minimum speed for compressor)
/// let vfd = BoundedVariable::new(
/// BoundedVariableId::new("compressor_vfd"),
/// 0.8, // initial: 80% speed
/// 0.3, // min: 30%
/// 1.0, // max: 100%
/// )?;
///
/// // Apply Newton step
/// let new_value = vfd.clip_step(0.3); // tries to go to 1.1
/// assert_eq!(new_value, 1.0); // clipped to max
/// ```
#[derive(Debug, Clone, PartialEq)]
pub struct BoundedVariable {
/// Unique identifier for this bounded variable
id: BoundedVariableId,
/// Current value of the variable
value: f64,
/// Lower bound (inclusive)
min: f64,
/// Upper bound (inclusive)
max: f64,
/// Optional component this variable controls
component_id: Option<String>,
}
impl BoundedVariable {
/// Tolerance for saturation detection (relative to bound range).
const SATURATION_TOL: f64 = 1e-9;
/// Creates a new bounded variable.
///
/// # Arguments
///
/// * `id` - Unique identifier for this variable
/// * `value` - Initial value
/// * `min` - Lower bound (inclusive)
/// * `max` - Upper bound (inclusive)
///
/// # Errors
///
/// Returns `BoundedVariableError::InvalidBounds` if `min >= max`.
/// Returns `BoundedVariableError::ValueOutOfBounds` if `value` is outside bounds.
///
/// # Example
///
/// ```rust,ignore
/// let valve = BoundedVariable::new(
/// BoundedVariableId::new("valve"),
/// 0.5, 0.0, 1.0
/// )?;
/// ```
pub fn new(
id: BoundedVariableId,
value: f64,
min: f64,
max: f64,
) -> Result<Self, BoundedVariableError> {
if min >= max {
return Err(BoundedVariableError::InvalidBounds { min, max });
}
if value < min || value > max {
return Err(BoundedVariableError::ValueOutOfBounds { value, min, max });
}
Ok(BoundedVariable {
id,
value,
min,
max,
component_id: None,
})
}
/// Creates a new bounded variable associated with a component.
///
/// # Arguments
///
/// * `id` - Unique identifier for this variable
/// * `component_id` - Component this variable controls
/// * `value` - Initial value
/// * `min` - Lower bound (inclusive)
/// * `max` - Upper bound (inclusive)
pub fn with_component(
id: BoundedVariableId,
component_id: impl Into<String>,
value: f64,
min: f64,
max: f64,
) -> Result<Self, BoundedVariableError> {
let mut var = Self::new(id, value, min, max)?;
var.component_id = Some(component_id.into());
Ok(var)
}
/// Returns the variable identifier.
pub fn id(&self) -> &BoundedVariableId {
&self.id
}
/// Returns the current value.
pub fn value(&self) -> f64 {
self.value
}
/// Returns the lower bound.
pub fn min(&self) -> f64 {
self.min
}
/// Returns the upper bound.
pub fn max(&self) -> f64 {
self.max
}
/// Returns the component ID if associated with a component.
pub fn component_id(&self) -> Option<&str> {
self.component_id.as_deref()
}
/// Sets the current value.
///
/// # Errors
///
/// Returns `BoundedVariableError::ValueOutOfBounds` if the value
/// is outside the bounds.
pub fn set_value(&mut self, value: f64) -> Result<(), BoundedVariableError> {
let range = self.max - self.min;
let tol = range * Self::SATURATION_TOL;
if value < self.min - tol || value > self.max + tol {
return Err(BoundedVariableError::ValueOutOfBounds {
value,
min: self.min,
max: self.max,
});
}
self.value = value.clamp(self.min, self.max);
Ok(())
}
/// Clips a proposed step to stay within bounds.
///
/// Given a delta from the current position, returns the clipped value
/// that stays within [min, max].
///
/// # Arguments
///
/// * `delta` - Proposed change from current value
///
/// # Returns
///
/// The clipped value: `clamp(current + delta, min, max)`
///
/// # Example
///
/// ```rust,ignore
/// let var = BoundedVariable::new(id, 0.5, 0.0, 1.0)?;
///
/// // Step that would exceed max
/// assert_eq!(var.clip_step(0.6), 1.0);
///
/// // Step that would exceed min
/// assert_eq!(var.clip_step(-0.7), 0.0);
///
/// // Step within bounds
/// assert_eq!(var.clip_step(0.3), 0.8);
/// ```
pub fn clip_step(&self, delta: f64) -> f64 {
clip_step(self.value, delta, self.min, self.max)
}
/// Applies a clipped step and updates the internal value.
///
/// # Arguments
///
/// * `delta` - Proposed change from current value
///
/// # Returns
///
/// The new (clipped) value.
pub fn apply_step(&mut self, delta: f64) -> f64 {
self.value = self.clip_step(delta);
self.value
}
/// Checks if the variable is at a bound (saturated).
///
/// # Returns
///
/// `Some(SaturationInfo)` if at a bound, `None` if not saturated.
///
/// # Example
///
/// ```rust,ignore
/// let var = BoundedVariable::new(id, 1.0, 0.0, 1.0)?;
/// let sat = var.is_saturated();
/// assert!(sat.is_some());
/// assert_eq!(sat.unwrap().saturation_type, SaturationType::UpperBound);
/// ```
pub fn is_saturated(&self) -> Option<SaturationInfo> {
let range = self.max - self.min;
let tol = range * Self::SATURATION_TOL;
if (self.value - self.min).abs() <= tol {
Some(SaturationInfo::new(
self.id.clone(),
SaturationType::LowerBound,
self.min,
))
} else if (self.value - self.max).abs() <= tol {
Some(SaturationInfo::new(
self.id.clone(),
SaturationType::UpperBound,
self.max,
))
} else {
None
}
}
/// Checks if a value is within the bounds.
pub fn is_within_bounds(&self, value: f64) -> bool {
value >= self.min && value <= self.max
}
/// Returns the distance to the nearest bound.
///
/// Returns 0.0 if at a bound, positive otherwise.
pub fn distance_to_bound(&self) -> f64 {
let dist_to_min = self.value - self.min;
let dist_to_max = self.max - self.value;
dist_to_min.min(dist_to_max)
}
}
// ─────────────────────────────────────────────────────────────────────────────
// Step Clipping Function
// ─────────────────────────────────────────────────────────────────────────────
/// Clips a step to stay within bounds.
///
/// This is a standalone function for use in solver loops where you have
/// raw values rather than a `BoundedVariable` struct.
///
/// # Arguments
///
/// * `current` - Current value
/// * `delta` - Proposed change
/// * `min` - Lower bound
/// * `max` - Upper bound
///
/// # Returns
///
/// The clipped value: `clamp(current + delta, min, max)`
///
/// # Edge Cases
///
/// - NaN delta returns the clamped current value (NaN propagates to proposed, then clamped)
/// - Infinite delta clips to the appropriate bound
/// - If min == max, returns min (degenerate case)
///
/// # Example
///
/// ```rust
/// use entropyk_solver::inverse::clip_step;
///
/// // Normal clipping
/// assert_eq!(clip_step(0.5, 0.6, 0.0, 1.0), 1.0);
/// assert_eq!(clip_step(0.5, -0.7, 0.0, 1.0), 0.0);
///
/// // Within bounds
/// assert_eq!(clip_step(0.5, 0.3, 0.0, 1.0), 0.8);
/// ```
pub fn clip_step(current: f64, delta: f64, min: f64, max: f64) -> f64 {
let proposed = current + delta;
// Handle NaN: if proposed is NaN, clamp will return min
// This is intentional behavior - NaN steps should not corrupt state
if proposed.is_nan() {
// Return current if it's valid, otherwise min
return if current.is_nan() {
min
} else {
current.clamp(min, max)
};
}
proposed.clamp(min, max)
}
// ─────────────────────────────────────────────────────────────────────────────
// Tests
// ─────────────────────────────────────────────────────────────────────────────
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_bounded_variable_id_creation() {
let id = BoundedVariableId::new("valve_position");
assert_eq!(id.as_str(), "valve_position");
assert_eq!(id.to_string(), "valve_position");
}
#[test]
fn test_bounded_variable_id_from_impls() {
let id1 = BoundedVariableId::from("valve");
assert_eq!(id1.as_str(), "valve");
let id2 = BoundedVariableId::from("valve".to_string());
assert_eq!(id2.as_str(), "valve");
}
#[test]
fn test_bounded_variable_creation() {
let var = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
assert_eq!(var.value(), 0.5);
assert_eq!(var.min(), 0.0);
assert_eq!(var.max(), 1.0);
assert!(var.component_id().is_none());
}
#[test]
fn test_bounded_variable_with_component() {
let var = BoundedVariable::with_component(
BoundedVariableId::new("valve"),
"expansion_valve",
0.5,
0.0,
1.0,
)
.unwrap();
assert_eq!(var.component_id(), Some("expansion_valve"));
}
#[test]
fn test_invalid_bounds_min_eq_max() {
let result = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 1.0, 1.0);
assert!(matches!(
result,
Err(BoundedVariableError::InvalidBounds { min: 1.0, max: 1.0 })
));
}
#[test]
fn test_invalid_bounds_min_gt_max() {
let result = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 1.0, 0.0);
assert!(matches!(
result,
Err(BoundedVariableError::InvalidBounds { min: 1.0, max: 0.0 })
));
}
#[test]
fn test_value_out_of_bounds_below() {
let result = BoundedVariable::new(BoundedVariableId::new("valve"), -0.1, 0.0, 1.0);
assert!(matches!(
result,
Err(BoundedVariableError::ValueOutOfBounds {
value: -0.1,
min: 0.0,
max: 1.0
})
));
}
#[test]
fn test_value_out_of_bounds_above() {
let result = BoundedVariable::new(BoundedVariableId::new("valve"), 1.1, 0.0, 1.0);
assert!(matches!(
result,
Err(BoundedVariableError::ValueOutOfBounds {
value: 1.1,
min: 0.0,
max: 1.0
})
));
}
#[test]
fn test_clip_step_within_bounds() {
let var = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
assert_eq!(var.clip_step(0.3), 0.8);
assert_eq!(var.clip_step(-0.3), 0.2);
}
#[test]
fn test_clip_step_exceeds_max() {
let var = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
assert_eq!(var.clip_step(0.6), 1.0);
assert_eq!(var.clip_step(1.0), 1.0);
}
#[test]
fn test_clip_step_exceeds_min() {
let var = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
assert_eq!(var.clip_step(-0.6), 0.0);
assert_eq!(var.clip_step(-1.0), 0.0);
}
#[test]
fn test_clip_step_standalone() {
// Normal cases
approx::assert_relative_eq!(clip_step(0.5, 0.6, 0.0, 1.0), 1.0);
approx::assert_relative_eq!(clip_step(0.5, -0.7, 0.0, 1.0), 0.0);
approx::assert_relative_eq!(clip_step(0.5, 0.3, 0.0, 1.0), 0.8);
// At bounds
approx::assert_relative_eq!(clip_step(0.0, -0.1, 0.0, 1.0), 0.0);
approx::assert_relative_eq!(clip_step(1.0, 0.1, 0.0, 1.0), 1.0);
}
#[test]
fn test_clip_step_nan_handling() {
// NaN delta should not propagate; return clamped current
let result = clip_step(0.5, f64::NAN, 0.0, 1.0);
approx::assert_relative_eq!(result, 0.5);
}
#[test]
fn test_clip_step_infinity() {
// Positive infinity clips to max
approx::assert_relative_eq!(clip_step(0.5, f64::INFINITY, 0.0, 1.0), 1.0);
// Negative infinity clips to min
approx::assert_relative_eq!(clip_step(0.5, f64::NEG_INFINITY, 0.0, 1.0), 0.0);
}
#[test]
fn test_saturation_at_min() {
let var = BoundedVariable::new(BoundedVariableId::new("valve"), 0.0, 0.0, 1.0).unwrap();
let sat = var.is_saturated();
assert!(sat.is_some());
let info = sat.unwrap();
assert_eq!(info.saturation_type, SaturationType::LowerBound);
approx::assert_relative_eq!(info.bound_value, 0.0);
}
#[test]
fn test_saturation_at_max() {
let var = BoundedVariable::new(BoundedVariableId::new("valve"), 1.0, 0.0, 1.0).unwrap();
let sat = var.is_saturated();
assert!(sat.is_some());
let info = sat.unwrap();
assert_eq!(info.saturation_type, SaturationType::UpperBound);
approx::assert_relative_eq!(info.bound_value, 1.0);
}
#[test]
fn test_no_saturation_in_middle() {
let var = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
assert!(var.is_saturated().is_none());
}
#[test]
fn test_saturation_info_with_constraint_target() {
let info = SaturationInfo::new(
BoundedVariableId::new("valve"),
SaturationType::UpperBound,
1.0,
)
.with_constraint_target(1.5);
assert_eq!(info.constraint_target, Some(1.5));
}
#[test]
fn test_set_value_valid() {
let mut var = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
var.set_value(0.7).unwrap();
approx::assert_relative_eq!(var.value(), 0.7);
var.set_value(0.0).unwrap();
approx::assert_relative_eq!(var.value(), 0.0);
var.set_value(1.0).unwrap();
approx::assert_relative_eq!(var.value(), 1.0);
}
#[test]
fn test_set_value_invalid() {
let mut var = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
assert!(var.set_value(-0.1).is_err());
assert!(var.set_value(1.1).is_err());
// Value should remain unchanged
approx::assert_relative_eq!(var.value(), 0.5);
}
#[test]
fn test_apply_step() {
let mut var = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
let new_val = var.apply_step(0.6);
approx::assert_relative_eq!(new_val, 1.0);
approx::assert_relative_eq!(var.value(), 1.0);
let new_val = var.apply_step(-0.3);
approx::assert_relative_eq!(new_val, 0.7);
approx::assert_relative_eq!(var.value(), 0.7);
}
#[test]
fn test_is_within_bounds() {
let var = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
assert!(var.is_within_bounds(0.0));
assert!(var.is_within_bounds(0.5));
assert!(var.is_within_bounds(1.0));
assert!(!var.is_within_bounds(-0.1));
assert!(!var.is_within_bounds(1.1));
}
#[test]
fn test_distance_to_bound() {
let var = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
approx::assert_relative_eq!(var.distance_to_bound(), 0.5);
let var_at_min =
BoundedVariable::new(BoundedVariableId::new("valve"), 0.0, 0.0, 1.0).unwrap();
approx::assert_relative_eq!(var_at_min.distance_to_bound(), 0.0);
let var_at_max =
BoundedVariable::new(BoundedVariableId::new("valve"), 1.0, 0.0, 1.0).unwrap();
approx::assert_relative_eq!(var_at_max.distance_to_bound(), 0.0);
}
#[test]
fn test_saturation_type_display() {
assert_eq!(format!("{}", SaturationType::LowerBound), "lower bound");
assert_eq!(format!("{}", SaturationType::UpperBound), "upper bound");
}
#[test]
fn test_error_display() {
let err = BoundedVariableError::InvalidBounds { min: 0.0, max: 0.0 };
assert!(err.to_string().contains("0"));
let err = BoundedVariableError::DuplicateId {
id: BoundedVariableId::new("dup"),
};
assert!(err.to_string().contains("dup"));
let err = BoundedVariableError::InvalidComponent {
component_id: "unknown".to_string(),
};
assert!(err.to_string().contains("unknown"));
}
}

492
crates/solver/src/inverse/constraint.rs Normal file
View File

@@ -0,0 +1,492 @@
//! Constraint types for inverse control.
//!
//! This module defines the core types for constraint-based control:
//! - [`Constraint`]: A single output constraint (output - target = 0)
//! - [`ConstraintId`]: Type-safe constraint identifier
//! - [`ConstraintError`]: Errors during constraint operations
//! - [`ComponentOutput`]: Reference to a measurable component property
use std::fmt;
use thiserror::Error;
// ─────────────────────────────────────────────────────────────────────────────
// ConstraintId - Type-safe identifier
// ─────────────────────────────────────────────────────────────────────────────
/// Type-safe identifier for a constraint.
///
/// Uses a string internally but provides type safety and clear intent.
#[derive(Debug, Clone, PartialEq, Eq, Hash)]
pub struct ConstraintId(String);
impl ConstraintId {
/// Creates a new constraint identifier.
///
/// # Arguments
///
/// * `id` - Unique identifier string for the constraint
///
/// # Example
///
/// ```rust
/// use entropyk_solver::inverse::ConstraintId;
///
/// let id = ConstraintId::new("superheat_control");
/// ```
pub fn new(id: impl Into<String>) -> Self {
ConstraintId(id.into())
}
/// Returns the identifier as a string slice.
pub fn as_str(&self) -> &str {
&self.0
}
}
impl fmt::Display for ConstraintId {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "{}", self.0)
}
}
impl From<&str> for ConstraintId {
fn from(s: &str) -> Self {
ConstraintId::new(s)
}
}
impl From<String> for ConstraintId {
fn from(s: String) -> Self {
ConstraintId(s)
}
}
// ─────────────────────────────────────────────────────────────────────────────
// ComponentOutput - Reference to measurable property
// ─────────────────────────────────────────────────────────────────────────────
/// Reference to a measurable component output property.
///
/// This enum defines which component properties can be used as constraint targets.
/// Each variant captures the necessary context to compute the property value.
#[derive(Debug, Clone, PartialEq)]
pub enum ComponentOutput {
/// Saturation temperature at a component (K).
///
/// References the saturation temperature at the component's location
/// in the thermodynamic cycle.
SaturationTemperature {
/// Component identifier (e.g., "evaporator", "condenser")
component_id: String,
},
/// Superheat above saturation (K).
///
/// Superheat = T_actual - T_sat
Superheat {
/// Component identifier
component_id: String,
},
/// Subcooling below saturation (K).
///
/// Subcooling = T_sat - T_actual
Subcooling {
/// Component identifier
component_id: String,
},
/// Heat transfer rate (W).
///
/// Heat transfer at a heat exchanger component.
HeatTransferRate {
/// Component identifier
component_id: String,
},
/// Mass flow rate (kg/s).
///
/// Mass flow through a component.
MassFlowRate {
/// Component identifier
component_id: String,
},
/// Pressure at a component (Pa).
Pressure {
/// Component identifier
component_id: String,
},
/// Temperature at a component (K).
Temperature {
/// Component identifier
component_id: String,
},
}
impl ComponentOutput {
/// Returns the component ID for this output reference.
pub fn component_id(&self) -> &str {
match self {
ComponentOutput::SaturationTemperature { component_id } => component_id,
ComponentOutput::Superheat { component_id } => component_id,
ComponentOutput::Subcooling { component_id } => component_id,
ComponentOutput::HeatTransferRate { component_id } => component_id,
ComponentOutput::MassFlowRate { component_id } => component_id,
ComponentOutput::Pressure { component_id } => component_id,
ComponentOutput::Temperature { component_id } => component_id,
}
}
}
// ─────────────────────────────────────────────────────────────────────────────
// ConstraintError - Error handling
// ─────────────────────────────────────────────────────────────────────────────
/// Errors that can occur during constraint operations.
#[derive(Error, Debug, Clone, PartialEq)]
pub enum ConstraintError {
/// The referenced component does not exist in the system.
#[error("Invalid component reference: '{component_id}' not found")]
InvalidReference {
/// The invalid component identifier
component_id: String,
},
/// A constraint with this ID already exists.
#[error("Duplicate constraint id: '{id}'")]
DuplicateId {
/// The duplicate constraint identifier
id: ConstraintId,
},
/// The constraint value is outside valid bounds.
#[error("Constraint value {value} is out of bounds: {reason}")]
OutOfBounds {
/// The out-of-bounds value
value: f64,
/// Reason for the bounds violation
reason: String,
},
/// The component output is not available or computable.
#[error("Component output not available: {output:?} for component '{component_id}'")]
OutputNotAvailable {
/// Component identifier
component_id: String,
/// The unavailable output type
output: String,
},
/// Invalid constraint configuration.
#[error("Invalid constraint configuration: {reason}")]
InvalidConfiguration {
/// Reason for the validation failure
reason: String,
},
}
// ─────────────────────────────────────────────────────────────────────────────
// Constraint - Core constraint type
// ─────────────────────────────────────────────────────────────────────────────
/// An output constraint for inverse control.
///
/// A constraint specifies a desired relationship:
///
/// $$r = f(x) - y_{target} = 0$$
///
/// where:
/// - $f(x)$ is the measured output (via [`ComponentOutput`])
/// - $y_{target}$ is the target value
/// - $r$ is the constraint residual
///
/// # Tolerance Guidance
///
/// The tolerance determines when a constraint is considered "satisfied". Choose based on:
///
/// | Property | Recommended Tolerance | Notes |
/// |----------|----------------------|-------|
/// | Superheat/Subcooling (K) | 0.01 - 0.1 | Temperature precision typically 0.1K |
/// | Pressure (Pa) | 100 - 1000 | ~0.1-1% of typical operating pressure |
/// | Heat Transfer (W) | 10 - 100 | Depends on system capacity |
/// | Mass Flow (kg/s) | 0.001 - 0.01 | ~1% of typical flow rate |
///
/// Default tolerance is `1e-6`, which may be too tight for some applications.
/// Use [`Constraint::with_tolerance()`] to customize.
///
/// # Example
///
/// ```rust,ignore
/// use entropyk_solver::inverse::{Constraint, ConstraintId, ComponentOutput};
///
/// // Superheat constraint: maintain 5K superheat at evaporator outlet
/// let constraint = Constraint::new(
/// ConstraintId::new("evap_superheat"),
/// ComponentOutput::Superheat {
/// component_id: "evaporator".to_string()
/// },
/// 5.0, // target: 5K superheat
/// );
/// ```
#[derive(Debug, Clone, PartialEq)]
pub struct Constraint {
/// Unique identifier for this constraint
id: ConstraintId,
/// Which component output to constrain
output: ComponentOutput,
/// Target value for the constraint
target_value: f64,
/// Tolerance for convergence.
///
/// A constraint is satisfied when `|measured - target| < tolerance`.
/// See struct documentation for guidance on choosing appropriate values.
tolerance: f64,
}
impl Constraint {
/// Creates a new constraint with default tolerance.
///
/// # Arguments
///
/// * `id` - Unique identifier for this constraint
/// * `output` - Component output to constrain
/// * `target_value` - Desired target value
///
/// # Example
///
/// ```rust,ignore
/// let constraint = Constraint::new(
/// ConstraintId::new("superheat_control"),
/// ComponentOutput::Superheat {
/// component_id: "evaporator".to_string()
/// },
/// 5.0,
/// );
/// ```
pub fn new(id: ConstraintId, output: ComponentOutput, target_value: f64) -> Self {
Constraint {
id,
output,
target_value,
tolerance: 1e-6,
}
}
/// Creates a new constraint with custom tolerance.
///
/// # Arguments
///
/// * `id` - Unique identifier for this constraint
/// * `output` - Component output to constrain
/// * `target_value` - Desired target value
/// * `tolerance` - Convergence tolerance
pub fn with_tolerance(
id: ConstraintId,
output: ComponentOutput,
target_value: f64,
tolerance: f64,
) -> Result<Self, ConstraintError> {
if tolerance <= 0.0 {
return Err(ConstraintError::InvalidConfiguration {
reason: format!("Tolerance must be positive, got {}", tolerance),
});
}
Ok(Constraint {
id,
output,
target_value,
tolerance,
})
}
/// Returns the constraint identifier.
pub fn id(&self) -> &ConstraintId {
&self.id
}
/// Returns the component output reference.
pub fn output(&self) -> &ComponentOutput {
&self.output
}
/// Returns the target value.
pub fn target_value(&self) -> f64 {
self.target_value
}
/// Returns the convergence tolerance.
pub fn tolerance(&self) -> f64 {
self.tolerance
}
/// Computes the constraint residual.
///
/// The residual is:
///
/// $$r = f(x) - y_{target}$$
///
/// where $f(x)$ is the measured output value.
///
/// # Arguments
///
/// * `measured_value` - The current value of the constrained output
///
/// # Returns
///
/// The constraint residual (measured - target)
pub fn compute_residual(&self, measured_value: f64) -> f64 {
measured_value - self.target_value
}
/// Checks if the constraint is satisfied.
///
/// # Arguments
///
/// * `measured_value` - The current value of the constrained output
///
/// # Returns
///
/// `true` if |residual| < tolerance
pub fn is_satisfied(&self, measured_value: f64) -> bool {
let residual = self.compute_residual(measured_value);
residual.abs() < self.tolerance
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_constraint_id_creation() {
let id = ConstraintId::new("test_constraint");
assert_eq!(id.as_str(), "test_constraint");
assert_eq!(id.to_string(), "test_constraint");
}
#[test]
fn test_constraint_id_from_string() {
let id = ConstraintId::from("my_constraint".to_string());
assert_eq!(id.as_str(), "my_constraint");
}
#[test]
fn test_constraint_id_from_str() {
let id = ConstraintId::from("constraint_1");
assert_eq!(id.as_str(), "constraint_1");
}
#[test]
fn test_component_output_component_id() {
let output = ComponentOutput::Superheat {
component_id: "evaporator".to_string(),
};
assert_eq!(output.component_id(), "evaporator");
}
#[test]
fn test_constraint_creation() {
let constraint = Constraint::new(
ConstraintId::new("superheat_control"),
ComponentOutput::Superheat {
component_id: "evaporator".to_string(),
},
5.0,
);
assert_eq!(constraint.target_value(), 5.0);
assert_eq!(constraint.tolerance(), 1e-6);
}
#[test]
fn test_constraint_with_tolerance() {
let constraint = Constraint::with_tolerance(
ConstraintId::new("superheat"),
ComponentOutput::Superheat {
component_id: "evaporator".to_string(),
},
5.0,
1e-4,
)
.unwrap();
assert_eq!(constraint.tolerance(), 1e-4);
}
#[test]
fn test_constraint_invalid_tolerance() {
let result = Constraint::with_tolerance(
ConstraintId::new("superheat"),
ComponentOutput::Superheat {
component_id: "evaporator".to_string(),
},
5.0,
-1e-4,
);
assert!(matches!(
result,
Err(ConstraintError::InvalidConfiguration { .. })
));
}
#[test]
fn test_constraint_compute_residual() {
let constraint = Constraint::new(
ConstraintId::new("superheat"),
ComponentOutput::Superheat {
component_id: "evaporator".to_string(),
},
5.0,
);
// Measured = 6.0, target = 5.0, residual = 1.0
let residual = constraint.compute_residual(6.0);
approx::assert_relative_eq!(residual, 1.0);
// Measured = 4.5, target = 5.0, residual = -0.5
let residual = constraint.compute_residual(4.5);
approx::assert_relative_eq!(residual, -0.5);
}
#[test]
fn test_constraint_is_satisfied() {
let constraint = Constraint::new(
ConstraintId::new("superheat"),
ComponentOutput::Superheat {
component_id: "evaporator".to_string(),
},
5.0,
);
// Exactly at target
assert!(constraint.is_satisfied(5.0));
// Within tolerance
assert!(constraint.is_satisfied(5.0 + 0.5e-6));
assert!(constraint.is_satisfied(5.0 - 0.5e-6));
// Outside tolerance
assert!(!constraint.is_satisfied(5.0 + 2e-6));
assert!(!constraint.is_satisfied(5.0 - 2e-6));
}
#[test]
fn test_constraint_error_display() {
let error = ConstraintError::InvalidReference {
component_id: "unknown".to_string(),
};
assert!(error.to_string().contains("unknown"));
let error = ConstraintError::DuplicateId {
id: ConstraintId::new("dup"),
};
assert!(error.to_string().contains("dup"));
}
}

570
crates/solver/src/inverse/embedding.rs Normal file
View File

@@ -0,0 +1,570 @@
//! Residual embedding for One-Shot inverse control.
//!
//! This module implements the core innovation of Epic 5: embedding constraints
//! directly into the residual system for simultaneous solving with cycle equations.
//!
//! # Mathematical Foundation
//!
//! In One-Shot inverse control, constraints are added to the residual vector:
//!
//! $$r_{total} = [r_{cycle}, r_{constraints}]^T$$
//!
//! where:
//! - $r_{cycle}$ are the component residual equations
//! - $r_{constraints}$ are constraint residuals: $f(x) - y_{target}$
//!
//! The solver adjusts both edge states AND control variables simultaneously
//! to satisfy all equations.
//!
//! # Degrees of Freedom (DoF) Validation
//!
//! For a well-posed system:
//!
//! $$n_{equations} = n_{edge\_eqs} + n_{constraints}$$
//! $$n_{unknowns} = n_{edge\_unknowns} + n_{controls}$$
//!
//! The system is balanced when: $n_{equations} = n_{unknowns}$
//!
//! # Example
//!
//! ```rust,ignore
//! use entropyk_solver::inverse::{Constraint, ConstraintId, ComponentOutput};
//! use entropyk_solver::inverse::{BoundedVariable, BoundedVariableId};
//!
//! // Define constraint: superheat = 5K
//! let constraint = Constraint::new(
//! ConstraintId::new("superheat_control"),
//! ComponentOutput::Superheat { component_id: "evaporator".into() },
//! 5.0,
//! );
//!
//! // Define control variable: valve position
//! let valve = BoundedVariable::new(
//! BoundedVariableId::new("expansion_valve"),
//! 0.5, 0.0, 1.0,
//! )?;
//!
//! // Link constraint to control for One-Shot solving
//! system.add_constraint(constraint)?;
//! system.add_bounded_variable(valve)?;
//! system.link_constraint_to_control(
//! &ConstraintId::new("superheat_control"),
//! &BoundedVariableId::new("expansion_valve"),
//! )?;
//!
//! // Validate DoF before solving
//! system.validate_inverse_control_dof()?;
//! ```
use std::collections::HashMap;
use thiserror::Error;
use super::{BoundedVariableId, ConstraintId};
// ─────────────────────────────────────────────────────────────────────────────
// DoFError - Degrees of Freedom Validation Errors
// ─────────────────────────────────────────────────────────────────────────────
/// Errors during degrees of freedom validation for inverse control.
#[derive(Error, Debug, Clone, PartialEq)]
pub enum DoFError {
/// The system has more constraints than control variables.
#[error(
"Over-constrained system: {constraint_count} constraints but only {control_count} control variables \
(equations: {equation_count}, unknowns: {unknown_count})"
)]
OverConstrainedSystem {
constraint_count: usize,
control_count: usize,
equation_count: usize,
unknown_count: usize,
},
/// The system has fewer constraints than control variables (may still converge).
#[error(
"Under-constrained system: {constraint_count} constraints for {control_count} control variables \
(equations: {equation_count}, unknowns: {unknown_count})"
)]
UnderConstrainedSystem {
constraint_count: usize,
control_count: usize,
equation_count: usize,
unknown_count: usize,
},
/// The referenced constraint does not exist.
#[error("Constraint '{constraint_id}' not found when linking to control")]
ConstraintNotFound { constraint_id: ConstraintId },
/// The referenced bounded variable does not exist.
#[error("Bounded variable '{bounded_variable_id}' not found when linking to constraint")]
BoundedVariableNotFound {
bounded_variable_id: BoundedVariableId,
},
/// The constraint is already linked to a control variable.
#[error("Constraint '{constraint_id}' is already linked to control '{existing}'")]
AlreadyLinked {
constraint_id: ConstraintId,
existing: BoundedVariableId,
},
/// The control variable is already linked to another constraint.
#[error(
"Control variable '{bounded_variable_id}' is already linked to constraint '{existing}'"
)]
ControlAlreadyLinked {
bounded_variable_id: BoundedVariableId,
existing: ConstraintId,
},
}
// ─────────────────────────────────────────────────────────────────────────────
// ControlMapping - Constraint → Control Variable Mapping
// ─────────────────────────────────────────────────────────────────────────────
/// A mapping from a constraint to its control variable.
///
/// This establishes the relationship needed for One-Shot solving where
/// the solver adjusts the control variable to satisfy the constraint.
#[derive(Debug, Clone, PartialEq)]
pub struct ControlMapping {
/// The constraint to satisfy.
pub constraint_id: ConstraintId,
/// The control variable to adjust.
pub bounded_variable_id: BoundedVariableId,
/// Whether this mapping is active.
pub enabled: bool,
}
impl ControlMapping {
/// Creates a new control mapping.
pub fn new(constraint_id: ConstraintId, bounded_variable_id: BoundedVariableId) -> Self {
ControlMapping {
constraint_id,
bounded_variable_id,
enabled: true,
}
}
/// Creates a disabled mapping.
pub fn disabled(constraint_id: ConstraintId, bounded_variable_id: BoundedVariableId) -> Self {
ControlMapping {
constraint_id,
bounded_variable_id,
enabled: false,
}
}
/// Enables this mapping.
pub fn enable(&mut self) {
self.enabled = true;
}
/// Disables this mapping.
pub fn disable(&mut self) {
self.enabled = false;
}
}
// ─────────────────────────────────────────────────────────────────────────────
// InverseControlConfig - Configuration for Inverse Control
// ─────────────────────────────────────────────────────────────────────────────
/// Configuration for One-Shot inverse control.
///
/// Manages constraint-to-control-variable mappings for embedding constraints
/// into the residual system.
#[derive(Debug, Clone, Default)]
pub struct InverseControlConfig {
/// Mapping from constraint ID to control variable ID.
constraint_to_control: HashMap<ConstraintId, BoundedVariableId>,
/// Mapping from control variable ID to constraint ID (reverse lookup).
control_to_constraint: HashMap<BoundedVariableId, ConstraintId>,
/// Whether inverse control is enabled globally.
enabled: bool,
}
impl InverseControlConfig {
/// Creates a new empty inverse control configuration.
pub fn new() -> Self {
InverseControlConfig {
constraint_to_control: HashMap::new(),
control_to_constraint: HashMap::new(),
enabled: true,
}
}
/// Creates a disabled configuration.
pub fn disabled() -> Self {
InverseControlConfig {
constraint_to_control: HashMap::new(),
control_to_constraint: HashMap::new(),
enabled: false,
}
}
/// Returns whether inverse control is enabled.
pub fn is_enabled(&self) -> bool {
self.enabled
}
/// Enables inverse control.
pub fn enable(&mut self) {
self.enabled = true;
}
/// Disables inverse control.
pub fn disable(&mut self) {
self.enabled = false;
}
/// Returns the number of constraint-control mappings.
pub fn mapping_count(&self) -> usize {
self.constraint_to_control.len()
}
/// Links a constraint to a control variable.
///
/// # Errors
///
/// Returns `DoFError::AlreadyLinked` if the constraint is already linked.
/// Returns `DoFError::ControlAlreadyLinked` if the control is already linked.
pub fn link(
&mut self,
constraint_id: ConstraintId,
bounded_variable_id: BoundedVariableId,
) -> Result<(), DoFError> {
if let Some(existing) = self.constraint_to_control.get(&constraint_id) {
return Err(DoFError::AlreadyLinked {
constraint_id,
existing: existing.clone(),
});
}
if let Some(existing) = self.control_to_constraint.get(&bounded_variable_id) {
return Err(DoFError::ControlAlreadyLinked {
bounded_variable_id,
existing: existing.clone(),
});
}
self.constraint_to_control
.insert(constraint_id.clone(), bounded_variable_id.clone());
self.control_to_constraint
.insert(bounded_variable_id, constraint_id);
Ok(())
}
/// Unlinks a constraint from its control variable.
///
/// Returns the bounded variable ID that was linked, or `None` if not linked.
pub fn unlink_constraint(&mut self, constraint_id: &ConstraintId) -> Option<BoundedVariableId> {
if let Some(bounded_var_id) = self.constraint_to_control.remove(constraint_id) {
self.control_to_constraint.remove(&bounded_var_id);
Some(bounded_var_id)
} else {
None
}
}
/// Unlinks a control variable from its constraint.
///
/// Returns the constraint ID that was linked, or `None` if not linked.
pub fn unlink_control(
&mut self,
bounded_variable_id: &BoundedVariableId,
) -> Option<ConstraintId> {
if let Some(constraint_id) = self.control_to_constraint.remove(bounded_variable_id) {
self.constraint_to_control.remove(&constraint_id);
Some(constraint_id)
} else {
None
}
}
/// Returns the control variable linked to a constraint.
pub fn get_control(&self, constraint_id: &ConstraintId) -> Option<&BoundedVariableId> {
self.constraint_to_control.get(constraint_id)
}
/// Returns the constraint linked to a control variable.
pub fn get_constraint(&self, bounded_variable_id: &BoundedVariableId) -> Option<&ConstraintId> {
self.control_to_constraint.get(bounded_variable_id)
}
/// Returns an iterator over all constraint-to-control mappings.
pub fn mappings(&self) -> impl Iterator<Item = (&ConstraintId, &BoundedVariableId)> {
self.constraint_to_control.iter()
}
/// Returns an iterator over linked constraint IDs.
pub fn linked_constraints(&self) -> impl Iterator<Item = &ConstraintId> {
self.constraint_to_control.keys()
}
/// Returns an iterator over linked control variable IDs.
pub fn linked_controls(&self) -> impl Iterator<Item = &BoundedVariableId> {
self.control_to_constraint.keys()
}
/// Checks if a constraint is linked.
pub fn is_constraint_linked(&self, constraint_id: &ConstraintId) -> bool {
self.constraint_to_control.contains_key(constraint_id)
}
/// Checks if a control variable is linked.
pub fn is_control_linked(&self, bounded_variable_id: &BoundedVariableId) -> bool {
self.control_to_constraint.contains_key(bounded_variable_id)
}
/// Clears all mappings.
pub fn clear(&mut self) {
self.constraint_to_control.clear();
self.control_to_constraint.clear();
}
}
// ─────────────────────────────────────────────────────────────────────────────
// Tests
// ─────────────────────────────────────────────────────────────────────────────
#[cfg(test)]
mod tests {
use super::*;
fn make_constraint_id(s: &str) -> ConstraintId {
ConstraintId::new(s)
}
fn make_bounded_var_id(s: &str) -> BoundedVariableId {
BoundedVariableId::new(s)
}
#[test]
fn test_dof_error_display() {
let err = DoFError::OverConstrainedSystem {
constraint_count: 3,
control_count: 1,
equation_count: 10,
unknown_count: 8,
};
let msg = err.to_string();
assert!(msg.contains("3"));
assert!(msg.contains("1"));
assert!(msg.contains("10"));
assert!(msg.contains("8"));
assert!(msg.contains("Over-constrained"));
}
#[test]
fn test_dof_error_constraint_not_found() {
let err = DoFError::ConstraintNotFound {
constraint_id: make_constraint_id("test"),
};
assert!(err.to_string().contains("test"));
}
#[test]
fn test_dof_error_already_linked() {
let err = DoFError::AlreadyLinked {
constraint_id: make_constraint_id("c1"),
existing: make_bounded_var_id("v1"),
};
let msg = err.to_string();
assert!(msg.contains("c1"));
assert!(msg.contains("v1"));
}
#[test]
fn test_control_mapping_creation() {
let mapping = ControlMapping::new(
make_constraint_id("superheat"),
make_bounded_var_id("valve"),
);
assert_eq!(mapping.constraint_id.as_str(), "superheat");
assert_eq!(mapping.bounded_variable_id.as_str(), "valve");
assert!(mapping.enabled);
}
#[test]
fn test_control_mapping_disabled() {
let mapping = ControlMapping::disabled(
make_constraint_id("superheat"),
make_bounded_var_id("valve"),
);
assert!(!mapping.enabled);
}
#[test]
fn test_control_mapping_enable_disable() {
let mut mapping = ControlMapping::new(make_constraint_id("c"), make_bounded_var_id("v"));
mapping.disable();
assert!(!mapping.enabled);
mapping.enable();
assert!(mapping.enabled);
}
#[test]
fn test_inverse_control_config_new() {
let config = InverseControlConfig::new();
assert!(config.is_enabled());
assert_eq!(config.mapping_count(), 0);
}
#[test]
fn test_inverse_control_config_disabled() {
let config = InverseControlConfig::disabled();
assert!(!config.is_enabled());
}
#[test]
fn test_inverse_control_config_enable_disable() {
let mut config = InverseControlConfig::new();
config.disable();
assert!(!config.is_enabled());
config.enable();
assert!(config.is_enabled());
}
#[test]
fn test_inverse_control_config_link() {
let mut config = InverseControlConfig::new();
let result = config.link(make_constraint_id("c1"), make_bounded_var_id("v1"));
assert!(result.is_ok());
assert_eq!(config.mapping_count(), 1);
let control = config.get_control(&make_constraint_id("c1"));
assert!(control.is_some());
assert_eq!(control.unwrap().as_str(), "v1");
let constraint = config.get_constraint(&make_bounded_var_id("v1"));
assert!(constraint.is_some());
assert_eq!(constraint.unwrap().as_str(), "c1");
}
#[test]
fn test_inverse_control_config_link_already_linked_constraint() {
let mut config = InverseControlConfig::new();
config
.link(make_constraint_id("c1"), make_bounded_var_id("v1"))
.unwrap();
let result = config.link(make_constraint_id("c1"), make_bounded_var_id("v2"));
assert!(matches!(result, Err(DoFError::AlreadyLinked { .. })));
if let Err(DoFError::AlreadyLinked {
constraint_id,
existing,
}) = result
{
assert_eq!(constraint_id.as_str(), "c1");
assert_eq!(existing.as_str(), "v1");
}
}
#[test]
fn test_inverse_control_config_link_already_linked_control() {
let mut config = InverseControlConfig::new();
config
.link(make_constraint_id("c1"), make_bounded_var_id("v1"))
.unwrap();
let result = config.link(make_constraint_id("c2"), make_bounded_var_id("v1"));
assert!(matches!(result, Err(DoFError::ControlAlreadyLinked { .. })));
if let Err(DoFError::ControlAlreadyLinked {
bounded_variable_id,
existing,
}) = result
{
assert_eq!(bounded_variable_id.as_str(), "v1");
assert_eq!(existing.as_str(), "c1");
}
}
#[test]
fn test_inverse_control_config_unlink_constraint() {
let mut config = InverseControlConfig::new();
config
.link(make_constraint_id("c1"), make_bounded_var_id("v1"))
.unwrap();
let removed = config.unlink_constraint(&make_constraint_id("c1"));
assert!(removed.is_some());
assert_eq!(removed.unwrap().as_str(), "v1");
assert_eq!(config.mapping_count(), 0);
let removed_again = config.unlink_constraint(&make_constraint_id("c1"));
assert!(removed_again.is_none());
}
#[test]
fn test_inverse_control_config_unlink_control() {
let mut config = InverseControlConfig::new();
config
.link(make_constraint_id("c1"), make_bounded_var_id("v1"))
.unwrap();
let removed = config.unlink_control(&make_bounded_var_id("v1"));
assert!(removed.is_some());
assert_eq!(removed.unwrap().as_str(), "c1");
assert_eq!(config.mapping_count(), 0);
}
#[test]
fn test_inverse_control_config_is_linked() {
let mut config = InverseControlConfig::new();
assert!(!config.is_constraint_linked(&make_constraint_id("c1")));
assert!(!config.is_control_linked(&make_bounded_var_id("v1")));
config
.link(make_constraint_id("c1"), make_bounded_var_id("v1"))
.unwrap();
assert!(config.is_constraint_linked(&make_constraint_id("c1")));
assert!(config.is_control_linked(&make_bounded_var_id("v1")));
}
#[test]
fn test_inverse_control_config_mappings_iterator() {
let mut config = InverseControlConfig::new();
config
.link(make_constraint_id("c1"), make_bounded_var_id("v1"))
.unwrap();
config
.link(make_constraint_id("c2"), make_bounded_var_id("v2"))
.unwrap();
let mappings: Vec<_> = config.mappings().collect();
assert_eq!(mappings.len(), 2);
}
#[test]
fn test_inverse_control_config_clear() {
let mut config = InverseControlConfig::new();
config
.link(make_constraint_id("c1"), make_bounded_var_id("v1"))
.unwrap();
config
.link(make_constraint_id("c2"), make_bounded_var_id("v2"))
.unwrap();
config.clear();
assert_eq!(config.mapping_count(), 0);
}
}

53
crates/solver/src/inverse/mod.rs Normal file
View File

@@ -0,0 +1,53 @@
//! # Inverse Control & Optimization
//!
//! This module implements the constraint-based inverse control system for thermodynamic
//! simulation. Instead of using external optimizers, constraints are embedded directly
//! into the residual system for "One-Shot" solving.
//!
//! # Mathematical Foundation
//!
//! A constraint is defined as:
//!
//! $$r_{constraint} = f(x) - y_{target} = 0$$
//!
//! where:
//! - $f(x)$ is a measurable system output (function of state vector $x$)
//! - $y_{target}$ is the desired target value
//! - $r_{constraint}$ is the constraint residual
//!
//! Bounded control variables ensure Newton steps stay physically possible:
//!
//! $$x_{new} = \text{clip}(x_{current} + \Delta x, x_{min}, x_{max})$$
//!
//! # Example
//!
//! ```rust,ignore
//! use entropyk_solver::inverse::{Constraint, ConstraintId, ComponentOutput};
//! use entropyk_solver::inverse::{BoundedVariable, BoundedVariableId};
//!
//! // Define a superheat constraint: superheat = 5K
//! let constraint = Constraint::new(
//! ConstraintId::new("superheat_control"),
//! ComponentOutput::Superheat { component_id: "evaporator".into() },
//! 5.0, // target: 5K superheat
//! );
//!
//! // Define a bounded control variable (valve position)
//! let valve = BoundedVariable::new(
//! BoundedVariableId::new("expansion_valve"),
//! 0.5, // initial: 50% open
//! 0.0, // min: fully closed
//! 1.0, // max: fully open
//! )?;
//! ```
pub mod bounded;
pub mod constraint;
pub mod embedding;
pub use bounded::{
clip_step, BoundedVariable, BoundedVariableError, BoundedVariableId, SaturationInfo,
SaturationType,
};
pub use constraint::{ComponentOutput, Constraint, ConstraintError, ConstraintId};
pub use embedding::{ControlMapping, DoFError, InverseControlConfig};

32
crates/solver/src/jacobian.rs
View File

@@ -23,7 +23,6 @@
use nalgebra::{DMatrix, DVector};
/// Wrapper around `nalgebra::DMatrix<f64>` for Jacobian operations.
///
/// The Jacobian matrix J represents the partial derivatives of the residual vector
@@ -142,11 +141,11 @@ impl JacobianMatrix {
// For square systems, use LU decomposition
if self.0.nrows() == self.0.ncols() {
let lu = self.0.clone().lu();
// Solve J·Δx = -r
let r_vec = DVector::from_row_slice(residuals);
let neg_r = -r_vec;
match lu.solve(&neg_r) {
Some(delta) => Some(delta.iter().copied().collect()),
None => {
@@ -162,10 +161,10 @@ impl JacobianMatrix {
self.0.nrows(),
self.0.ncols()
);
let r_vec = DVector::from_row_slice(residuals);
let neg_r = -r_vec;
// Use SVD for robust least-squares solution
let svd = self.0.clone().svd(true, true);
match svd.solve(&neg_r, 1e-10) {
@@ -283,10 +282,10 @@ impl JacobianMatrix {
pub fn condition_number(&self) -> Option<f64> {
let svd = self.0.clone().svd(false, false);
let singular_values = svd.singular_values;
let max_sv = singular_values.max();
let min_sv = singular_values.min();
if min_sv > 1e-14 {
Some(max_sv / min_sv)
} else {
@@ -310,7 +309,10 @@ impl JacobianMatrix {
/// - Row/col ranges are inclusive-start, exclusive-end: `row_start..row_end`.
///
/// # AC: #6
pub fn block_structure(&self, system: &crate::system::System) -> Vec<(usize, usize, usize, usize)> {
pub fn block_structure(
&self,
system: &crate::system::System,
) -> Vec<(usize, usize, usize, usize)> {
let n_circuits = system.circuit_count();
let mut blocks = Vec::with_capacity(n_circuits);
@@ -335,7 +337,7 @@ impl JacobianMatrix {
let col_start = *indices.iter().min().unwrap();
let col_end = *indices.iter().max().unwrap() + 1; // exclusive
// Equations mirror state layout for square systems
// Equations mirror state layout for square systems
let row_start = col_start;
let row_end = col_end;
@@ -583,7 +585,7 @@ mod tests {
let mut j = JacobianMatrix::zeros(2, 2);
j.set(0, 0, 5.0);
j.set(1, 1, 7.0);
assert_relative_eq!(j.get(0, 0).unwrap(), 5.0);
assert_relative_eq!(j.get(1, 1).unwrap(), 7.0);
assert_relative_eq!(j.get(0, 1).unwrap(), 0.0);
@@ -606,7 +608,7 @@ mod tests {
// r[0] = x0^2 + x0*x1
// r[1] = sin(x0) + cos(x1)
// J = [[2*x0 + x1, x0], [cos(x0), -sin(x1)]]
let state: Vec<f64> = vec![0.5, 1.0];
let residuals: Vec<f64> = vec![
state[0].powi(2) + state[0] * state[1],
@@ -623,13 +625,13 @@ mod tests {
// Analytical values
let j00 = 2.0 * state[0] + state[1]; // 2*0.5 + 1.0 = 2.0
let j01 = state[0]; // 0.5
let j10 = state[0].cos(); // cos(0.5)
let j11 = -state[1].sin(); // -sin(1.0)
let j01 = state[0]; // 0.5
let j10 = state[0].cos(); // cos(0.5)
let j11 = -state[1].sin(); // -sin(1.0)
assert_relative_eq!(j_num.get(0, 0).unwrap(), j00, epsilon = 1e-5);
assert_relative_eq!(j_num.get(0, 1).unwrap(), j01, epsilon = 1e-5);
assert_relative_eq!(j_num.get(1, 0).unwrap(), j10, epsilon = 1e-5);
assert_relative_eq!(j_num.get(1, 1).unwrap(), j11, epsilon = 1e-5);
}
}
}

7
crates/solver/src/lib.rs
View File

@@ -11,23 +11,26 @@ pub mod criteria;
pub mod error;
pub mod graph;
pub mod initializer;
pub mod inverse;
pub mod jacobian;
pub mod macro_component;
pub mod solver;
pub mod system;
pub use criteria::{CircuitConvergence, ConvergenceCriteria, ConvergenceReport};
pub use coupling::{
compute_coupling_heat, coupling_groups, has_circular_dependencies, ThermalCoupling,
};
pub use criteria::{CircuitConvergence, ConvergenceCriteria, ConvergenceReport};
pub use entropyk_components::ConnectionError;
pub use error::{AddEdgeError, TopologyError};
pub use initializer::{
antoine_pressure, AntoineCoefficients, InitializerConfig, InitializerError, SmartInitializer,
};
pub use inverse::{ComponentOutput, Constraint, ConstraintError, ConstraintId};
pub use jacobian::JacobianMatrix;
pub use macro_component::{MacroComponent, MacroComponentSnapshot, PortMapping};
pub use solver::{
ConvergedState, ConvergenceStatus, FallbackConfig, FallbackSolver, JacobianFreezingConfig,
NewtonConfig, PicardConfig, Solver, SolverError, SolverStrategy, TimeoutConfig,
};
pub use system::{CircuitId, FlowEdge, System};

783
crates/solver/src/macro_component.rs Normal file
View File

@@ -0,0 +1,783 @@
//! Hierarchical Subsystems — MacroComponent
//!
//! A `MacroComponent` wraps a finalized [`System`] (topology + components) and
//! exposes it as a single [`Component`], enabling hierarchical composition.
//!
//! ## Architecture
//!
//! ```text
//! ┌─────────────────── MacroComponent ───────────────────┐
//! │ internal System (finalized) │
//! │ ┌─────┐ edge_a ┌─────┐ edge_b ┌─────┐ │
//! │ │Comp0├──────────►│Comp1├──────────►│Comp2│ │
//! │ └─────┘ └─────┘ └─────┘ │
//! │ │
//! │ external ports ← port_map │
//! │ port 0: edge_a.inlet (in) │
//! │ port 1: edge_b.outlet (out) │
//! └──────────────────────────────────────────────────────┘
//! ```
//!
//! ## Index Mapping & Coupling Equations
//!
//! The global solver assigns indices to all edges in the parent `System`.
//! Edges *inside* a `MacroComponent` are addressed via `global_state_offset`
//! (set during `System::finalize()` via `set_system_context`).
//!
//! When the `MacroComponent` is connected to external edges in the parent graph,
//! coupling residuals enforce continuity between those external edges and the
//! corresponding exposed internal edges:
//!
//! ```text
//! r_P = state[p_ext] − state[offset + 2 * internal_edge_pos] = 0
//! r_h = state[h_ext] − state[offset + 2 * internal_edge_pos + 1] = 0
//! ```
//!
//! ## Serialization (AC #4)
//!
//! The full component graph inside a `MacroComponent` cannot be trivially
//! serialized because `Box<dyn Component>` requires `typetag` or a custom
//! registry (deferred). Instead, `MacroComponent` exposes a **state snapshot**
//! ([`MacroComponentSnapshot`]) that captures the internal edge states and port
//! mappings. This is sufficient for persistence / restore of operating-point data.
use crate::system::System;
use entropyk_components::{
Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, SystemState,
};
use std::collections::HashMap;
// ─────────────────────────────────────────────────────────────────────────────
// Port mapping
// ─────────────────────────────────────────────────────────────────────────────
/// An exposed port on the MacroComponent surface.
///
/// Maps an internal edge (by position) to an external port visible to the
/// parent `System`.
#[derive(Debug, Clone)]
pub struct PortMapping {
/// Human-readable name for the external port (e.g. "evap_water_in").
pub name: String,
/// The internal edge index (position in the internal System's edge iteration
/// order) whose state this port corresponds to.
pub internal_edge_pos: usize,
/// A connected port to present externally (fluid, P, h).
pub port: ConnectedPort,
}
// ─────────────────────────────────────────────────────────────────────────────
// Serialization snapshot
// ─────────────────────────────────────────────────────────────────────────────
/// A serializable snapshot of a `MacroComponent`'s operating state.
///
/// Captures the internal edge state vector and port metadata so that an
/// operating point can be saved to disk and restored. The full component
/// topology (graph structure, `Box<dyn Component>` nodes) is **not** included
/// — reconstruction of the topology is the caller's responsibility.
///
/// # Example (JSON)
///
/// ```json
/// {
/// "label": "chiller_1",
/// "internal_edge_states": [1.5e5, 4.2e5, 8.0e4, 3.8e5],
/// "port_names": ["evap_in", "evap_out"]
/// }
/// ```
#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
pub struct MacroComponentSnapshot {
/// Optional human-readable label for the subsystem.
pub label: Option<String>,
/// Flat state vector for the internal edges: `[P_e0, h_e0, P_e1, h_e1, ...]`.
pub internal_edge_states: Vec<f64>,
/// Names of exposed ports, in the same order as `port_mappings`.
pub port_names: Vec<String>,
}
// ─────────────────────────────────────────────────────────────────────────────
// MacroComponent
// ─────────────────────────────────────────────────────────────────────────────
/// A hierarchical subsystem that wraps a `System` and implements `Component`.
///
/// This enables Modelica-style block composition: a chiller (compressor +
/// condenser + valve + evaporator) can be wrapped in a `MacroComponent` and
/// plugged into a higher-level `System`.
///
/// # Coupling equations
///
/// When `set_system_context` is called by `System::finalize()`, the component
/// stores the state indices of every parent-graph edge incident to its node.
/// `compute_residuals` then appends 2 coupling residuals per exposed port:
///
/// ```text
/// r_P[i] = state[p_ext_i] − state[offset + 2·internal_edge_pos_i] = 0
/// r_h[i] = state[h_ext_i] − state[offset + 2·internal_edge_pos_i + 1] = 0
/// ```
///
/// # Usage
///
/// ```no_run
/// use entropyk_solver::{System, MacroComponent};
/// use entropyk_components::Component;
///
/// // 1. Build and finalize internal system
/// let mut internal = System::new();
/// // ... add components & edges ...
/// internal.finalize().unwrap();
///
/// // 2. Wrap into a MacroComponent
/// let macro_comp = MacroComponent::new(internal);
///
/// // 3. Optionally expose ports
/// // macro_comp.expose_port(0, "inlet", port);
///
/// // 4. Add to a parent System (finalize() automatically wires context)
/// let mut parent = System::new();
/// parent.add_component(Box::new(macro_comp));
/// ```
pub struct MacroComponent {
/// The enclosed, finalized subsystem.
internal: System,
/// External port mappings. Ordered; index = external port index.
port_mappings: Vec<PortMapping>,
/// Cached external ports (mirrors port_mappings order).
external_ports: Vec<ConnectedPort>,
/// Maps external-port-index → internal-edge-position for fast lookup.
ext_to_internal_edge: HashMap<usize, usize>,
/// The global state vector offset assigned to this MacroComponent's first
/// internal edge. Set automatically via `set_system_context` during parent
/// `System::finalize()`. Defaults to 0.
global_state_offset: usize,
/// State indices `(p_idx, h_idx)` of every parent-graph edge incident to
/// this node (incoming and outgoing), in traversal order.
/// Populated by `set_system_context`; empty until finalization.
external_edge_state_indices: Vec<(usize, usize)>,
}
impl MacroComponent {
/// Creates a new `MacroComponent` wrapping the given *finalized* system.
///
/// # Panics
///
/// Panics if the internal system has not been finalized.
pub fn new(internal: System) -> Self {
Self {
internal,
port_mappings: Vec::new(),
external_ports: Vec::new(),
ext_to_internal_edge: HashMap::new(),
global_state_offset: 0,
external_edge_state_indices: Vec::new(),
}
}
/// Exposes an internal edge as an external port on the MacroComponent.
///
/// # Arguments
///
/// * `internal_edge_pos` — Position of the edge in the internal system's
/// edge iteration order (0-based).
/// * `name` — Human-readable label for this external port.
/// * `port` — A `ConnectedPort` representing the fluid, pressure and
/// enthalpy at this interface.
///
/// # Panics
///
/// Panics if `internal_edge_pos >= internal.edge_count()`.
pub fn expose_port(
&mut self,
internal_edge_pos: usize,
name: impl Into<String>,
port: ConnectedPort,
) {
assert!(
internal_edge_pos < self.internal.edge_count(),
"internal_edge_pos {} out of range (internal has {} edges)",
internal_edge_pos,
self.internal.edge_count()
);
let ext_idx = self.port_mappings.len();
self.port_mappings.push(PortMapping {
name: name.into(),
internal_edge_pos,
port: port.clone(),
});
self.external_ports.push(port);
self.ext_to_internal_edge.insert(ext_idx, internal_edge_pos);
}
/// Sets the global state-vector offset for this MacroComponent.
///
/// Prefer letting `System::finalize()` set this automatically via
/// `set_system_context`. This setter is kept for backward compatibility
/// and for manual test scenarios.
pub fn set_global_state_offset(&mut self, offset: usize) {
self.global_state_offset = offset;
}
/// Returns the global state offset.
pub fn global_state_offset(&self) -> usize {
self.global_state_offset
}
/// Returns a reference to the internal system.
pub fn internal_system(&self) -> &System {
&self.internal
}
/// Returns a mutable reference to the internal system.
pub fn internal_system_mut(&mut self) -> &mut System {
&mut self.internal
}
/// Returns the port mappings.
pub fn port_mappings(&self) -> &[PortMapping] {
&self.port_mappings
}
/// Number of internal edges (each contributes 2 state variables: P, h).
pub fn internal_edge_count(&self) -> usize {
self.internal.edge_count()
}
/// Total number of internal state variables (2 per edge).
pub fn internal_state_len(&self) -> usize {
self.internal.state_vector_len()
}
// ─── helpers ──────────────────────────────────────────────────────────────
/// Number of equations from internal components (excluding coupling eqs).
fn n_internal_equations(&self) -> usize {
self.internal
.traverse_for_jacobian()
.map(|(_, c, _)| c.n_equations())
.sum()
}
// ─── snapshot ─────────────────────────────────────────────────────────────
/// Captures the current internal state as a serializable snapshot.
///
/// The caller must supply the global state vector so that internal edge
/// states can be extracted. Returns `None` if the state vector is shorter
/// than expected.
pub fn to_snapshot(
&self,
global_state: &SystemState,
label: Option<String>,
) -> Option<MacroComponentSnapshot> {
let start = self.global_state_offset;
let end = start + self.internal_state_len();
if global_state.len() < end {
return None;
}
Some(MacroComponentSnapshot {
label,
internal_edge_states: global_state[start..end].to_vec(),
port_names: self.port_mappings.iter().map(|m| m.name.clone()).collect(),
})
}
}
// ─────────────────────────────────────────────────────────────────────────────
// Component trait implementation
// ─────────────────────────────────────────────────────────────────────────────
impl Component for MacroComponent {
/// Called by `System::finalize()` to inject the parent-level state offset
/// and the external edge state indices for this MacroComponent node.
///
/// `external_edge_state_indices` contains one `(p_idx, h_idx)` pair per
/// parent edge incident to this node (in traversal order: incoming, then
/// outgoing). The *i*-th entry is matched to `port_mappings[i]` when
/// emitting coupling residuals.
fn set_system_context(
&mut self,
state_offset: usize,
external_edge_state_indices: &[(usize, usize)],
) {
self.global_state_offset = state_offset;
self.external_edge_state_indices = external_edge_state_indices.to_vec();
}
fn internal_state_len(&self) -> usize {
// Delegates to the inherent method or computes directly
self.internal.state_vector_len()
}
fn compute_residuals(
&self,
state: &SystemState,
residuals: &mut ResidualVector,
) -> Result<(), ComponentError> {
let n_internal_vars = self.internal_state_len();
let start = self.global_state_offset;
let end = start + n_internal_vars;
if state.len() < end {
return Err(ComponentError::InvalidStateDimensions {
expected: end,
actual: state.len(),
});
}
let n_int_eqs = self.n_internal_equations();
let n_coupling = 2 * self.port_mappings.len();
let n_total = n_int_eqs + n_coupling;
if residuals.len() < n_total {
return Err(ComponentError::InvalidResidualDimensions {
expected: n_total,
actual: residuals.len(),
});
}
// --- 1. Delegate internal residuals ----------------------------------
let internal_state: SystemState = state[start..end].to_vec();
let mut internal_residuals = vec![0.0; n_int_eqs];
self.internal
.compute_residuals(&internal_state, &mut internal_residuals)?;
residuals[..n_int_eqs].copy_from_slice(&internal_residuals);
// --- 2. Port-coupling residuals --------------------------------------
// For each exposed port mapping we append two residuals that enforce
// continuity between the parent-graph external edge and the
// corresponding internal edge:
//
// r_P = state[p_ext] − state[offset + 2·internal_edge_pos] = 0
// r_h = state[h_ext] − state[offset + 2·internal_edge_pos + 1] = 0
for (i, mapping) in self.port_mappings.iter().enumerate() {
if let Some(&(p_ext, h_ext)) = self.external_edge_state_indices.get(i) {
let int_p = self.global_state_offset + 2 * mapping.internal_edge_pos;
let int_h = int_p + 1;
if state.len() <= int_h || state.len() <= p_ext || state.len() <= h_ext {
return Err(ComponentError::InvalidStateDimensions {
expected: int_h.max(p_ext).max(h_ext) + 1,
actual: state.len(),
});
}
residuals[n_int_eqs + 2 * i] = state[p_ext] - state[int_p];
residuals[n_int_eqs + 2 * i + 1] = state[h_ext] - state[int_h];
}
}
Ok(())
}
fn jacobian_entries(
&self,
state: &SystemState,
jacobian: &mut JacobianBuilder,
) -> Result<(), ComponentError> {
let n_internal_vars = self.internal_state_len();
let start = self.global_state_offset;
let end = start + n_internal_vars;
if state.len() < end {
return Err(ComponentError::InvalidStateDimensions {
expected: end,
actual: state.len(),
});
}
let n_int_eqs = self.n_internal_equations();
// --- 1. Internal Jacobian entries ------------------------------------
let internal_state: SystemState = state[start..end].to_vec();
let mut internal_jac = JacobianBuilder::new();
self.internal
.assemble_jacobian(&internal_state, &mut internal_jac)?;
// Offset columns by global_state_offset to translate from internal-local
// to global column indices.
for &(row, col, val) in internal_jac.entries() {
jacobian.add_entry(row, col + self.global_state_offset, val);
}
// --- 2. Coupling Jacobian entries ------------------------------------
// For each coupling residual pair (row_p, row_h):
//
// ∂r_P/∂state[p_ext] = +1
// ∂r_P/∂state[int_p] = −1
// ∂r_h/∂state[h_ext] = +1
// ∂r_h/∂state[int_h] = −1
for (i, mapping) in self.port_mappings.iter().enumerate() {
if let Some(&(p_ext, h_ext)) = self.external_edge_state_indices.get(i) {
let int_p = self.global_state_offset + 2 * mapping.internal_edge_pos;
let int_h = int_p + 1;
let row_p = n_int_eqs + 2 * i;
let row_h = row_p + 1;
jacobian.add_entry(row_p, p_ext, 1.0);
jacobian.add_entry(row_p, int_p, -1.0);
jacobian.add_entry(row_h, h_ext, 1.0);
jacobian.add_entry(row_h, int_h, -1.0);
}
}
Ok(())
}
fn n_equations(&self) -> usize {
// Internal equations + 2 coupling equations per exposed port.
self.n_internal_equations() + 2 * self.port_mappings.len()
}
fn get_ports(&self) -> &[ConnectedPort] {
&self.external_ports
}
}
// ─────────────────────────────────────────────────────────────────────────────
// Tests
// ─────────────────────────────────────────────────────────────────────────────
#[cfg(test)]
mod tests {
use super::*;
use crate::system::System;
use entropyk_components::port::{FluidId, Port};
use entropyk_core::{Enthalpy, Pressure};
/// Minimal mock component for testing.
struct MockInternalComponent {
n_equations: usize,
}
impl Component for MockInternalComponent {
fn compute_residuals(
&self,
state: &SystemState,
residuals: &mut ResidualVector,
) -> Result<(), ComponentError> {
// Simple identity: residual[i] = state[i] (so zero when state is zero)
for i in 0..self.n_equations {
residuals[i] = state.get(i).copied().unwrap_or(0.0);
}
Ok(())
}
fn jacobian_entries(
&self,
_state: &SystemState,
jacobian: &mut JacobianBuilder,
) -> Result<(), ComponentError> {
for i in 0..self.n_equations {
jacobian.add_entry(i, i, 1.0);
}
Ok(())
}
fn n_equations(&self) -> usize {
self.n_equations
}
fn get_ports(&self) -> &[ConnectedPort] {
&[]
}
}
fn make_mock(n: usize) -> Box<dyn Component> {
Box::new(MockInternalComponent { n_equations: n })
}
fn make_connected_port(fluid: &str, p_pa: f64, h_jkg: f64) -> ConnectedPort {
let p1 = Port::new(
FluidId::new(fluid),
Pressure::from_pascals(p_pa),
Enthalpy::from_joules_per_kg(h_jkg),
);
let p2 = Port::new(
FluidId::new(fluid),
Pressure::from_pascals(p_pa),
Enthalpy::from_joules_per_kg(h_jkg),
);
let (c1, _c2) = p1.connect(p2).unwrap();
c1
}
/// Build a simple linear subsystem: A → B → C (2 edges, 3 components).
fn build_simple_internal_system() -> System {
let mut sys = System::new();
let a = sys.add_component(make_mock(2));
let b = sys.add_component(make_mock(2));
let c = sys.add_component(make_mock(2));
sys.add_edge(a, b).unwrap();
sys.add_edge(b, c).unwrap();
sys.finalize().unwrap();
sys
}
#[test]
fn test_macro_component_creation() {
let sys = build_simple_internal_system();
let mc = MacroComponent::new(sys);
// 3 components × 2 equations each = 6 equations (no ports exposed yet,
// so no coupling equations).
assert_eq!(mc.n_equations(), 6);
// 2 edges → 4 state variables
assert_eq!(mc.internal_state_len(), 4);
// No ports exposed yet
assert!(mc.get_ports().is_empty());
}
#[test]
fn test_expose_port_adds_coupling_equations() {
let sys = build_simple_internal_system();
let mut mc = MacroComponent::new(sys);
let port = make_connected_port("R134a", 100_000.0, 400_000.0);
mc.expose_port(0, "inlet", port.clone());
// 6 internal + 2 coupling = 8
assert_eq!(mc.n_equations(), 8);
assert_eq!(mc.get_ports().len(), 1);
assert_eq!(mc.port_mappings()[0].name, "inlet");
assert_eq!(mc.port_mappings()[0].internal_edge_pos, 0);
}
#[test]
fn test_expose_multiple_ports() {
let sys = build_simple_internal_system();
let mut mc = MacroComponent::new(sys);
let port_in = make_connected_port("R134a", 100_000.0, 400_000.0);
let port_out = make_connected_port("R134a", 500_000.0, 450_000.0);
mc.expose_port(0, "inlet", port_in);
mc.expose_port(1, "outlet", port_out);
// 6 internal + 4 coupling = 10
assert_eq!(mc.n_equations(), 10);
assert_eq!(mc.get_ports().len(), 2);
assert_eq!(mc.port_mappings()[0].name, "inlet");
assert_eq!(mc.port_mappings()[1].name, "outlet");
}
#[test]
#[should_panic(expected = "internal_edge_pos 5 out of range")]
fn test_expose_port_out_of_range() {
let sys = build_simple_internal_system();
let mut mc = MacroComponent::new(sys);
let port = make_connected_port("R134a", 100_000.0, 400_000.0);
mc.expose_port(5, "bad", port);
}
#[test]
fn test_compute_residuals_delegation() {
let sys = build_simple_internal_system();
let mc = MacroComponent::new(sys);
// 4 state variables for 2 internal edges (no external coupling)
let state = vec![1.0, 2.0, 3.0, 4.0];
let mut residuals = vec![0.0; mc.n_equations()];
mc.compute_residuals(&state, &mut residuals).unwrap();
// 6 equations (no coupling ports)
assert_eq!(residuals.len(), 6);
}
#[test]
fn test_compute_residuals_with_offset() {
let sys = build_simple_internal_system();
let mut mc = MacroComponent::new(sys);
mc.set_global_state_offset(4);
// State vector: 4 padding + 4 internal = 8 total
let state = vec![0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 4.0];
let mut residuals = vec![0.0; mc.n_equations()];
mc.compute_residuals(&state, &mut residuals).unwrap();
assert_eq!(residuals.len(), 6);
}
#[test]
fn test_compute_residuals_state_too_short() {
let sys = build_simple_internal_system();
let mut mc = MacroComponent::new(sys);
mc.set_global_state_offset(4);
let state = vec![0.0; 5]; // Needs at least 8 (offset 4 + 4 internal vars)
let mut residuals = vec![0.0; mc.n_equations()];
let result = mc.compute_residuals(&state, &mut residuals);
assert!(result.is_err());
}
#[test]
fn test_jacobian_entries_delegation() {
let sys = build_simple_internal_system();
let mc = MacroComponent::new(sys);
let state = vec![0.0; mc.internal_state_len()];
let mut jac = JacobianBuilder::new();
mc.jacobian_entries(&state, &mut jac).unwrap();
// 6 equations → at least 6 diagonal entries from internal mocks
assert!(jac.len() >= 6);
}
#[test]
fn test_jacobian_entries_with_offset() {
let sys = build_simple_internal_system();
let mut mc = MacroComponent::new(sys);
mc.set_global_state_offset(10);
let state = vec![0.0; 10 + mc.internal_state_len()];
let mut jac = JacobianBuilder::new();
mc.jacobian_entries(&state, &mut jac).unwrap();
// Verify internal-delegated columns are offset by 10
for &(_, col, _) in jac.entries() {
assert!(col >= 10, "Column {} should be >= 10 (offset)", col);
}
}
#[test]
fn test_coupling_residuals_and_jacobian() {
// 2-edge internal system: edge0 = (P0, h0), edge1 = (P1, h1)
let sys = build_simple_internal_system();
let mut mc = MacroComponent::new(sys);
// Expose internal edge 0 as port "inlet"
let port = make_connected_port("R134a", 100_000.0, 400_000.0);
mc.expose_port(0, "inlet", port);
// Simulate finalization: inject external edge state index (p=6, h=7)
// and global offset = 4 (4 parent edges before the macro's internal block).
mc.set_global_state_offset(4);
mc.set_system_context(4, &[(6, 7)]);
// Global state: [*0, 1, 2, 3*, 4, 5, 6, 7, P_ext=1e5, h_ext=4e5]
// ^--- internal block at [4..8]
// ^--- ext edge at (6,7)... wait,
// Let's use a concrete layout:
// indices 0..3: some other parent edges
// indices 4..7: internal block (2 edges * 2 vars)
// 4=P_int_e0, 5=h_int_e0, 6=P_int_e1, 7=h_int_e1
// indices 8,9: external edge (p_ext=8, h_ext=9)
mc.set_system_context(4, &[(8, 9)]);
let mut state = vec![0.0; 10];
state[4] = 1.5e5; // P_int_e0
state[5] = 3.9e5; // h_int_e0
state[8] = 2.0e5; // P_ext
state[9] = 4.1e5; // h_ext
let n_eqs = mc.n_equations(); // 6 internal + 2 coupling = 8
assert_eq!(n_eqs, 8);
let mut residuals = vec![0.0; n_eqs];
mc.compute_residuals(&state, &mut residuals).unwrap();
// Coupling residuals:
// r[6] = state[8] - state[4] = 2e5 - 1.5e5 = 0.5e5
// r[7] = state[9] - state[5] = 4.1e5 - 3.9e5 = 0.2e5
assert!(
(residuals[6] - 0.5e5).abs() < 1.0,
"r_P mismatch: {}",
residuals[6]
);
assert!(
(residuals[7] - 0.2e5).abs() < 1.0,
"r_h mismatch: {}",
residuals[7]
);
// Jacobian coupling entries
let mut jac = JacobianBuilder::new();
mc.jacobian_entries(&state, &mut jac).unwrap();
let entries = jac.entries();
// Check that we have entries for p_ext=8 → +1, int_p=4 → -1
let find = |row: usize, col: usize| -> Option<f64> {
entries
.iter()
.find(|&&(r, c, _)| r == row && c == col)
.map(|&(_, _, v)| v)
};
assert_eq!(find(6, 8), Some(1.0), "expect ∂r_P/∂p_ext = +1");
assert_eq!(find(6, 4), Some(-1.0), "expect ∂r_P/∂int_p = -1");
assert_eq!(find(7, 9), Some(1.0), "expect ∂r_h/∂h_ext = +1");
assert_eq!(find(7, 5), Some(-1.0), "expect ∂r_h/∂int_h = -1");
}
#[test]
fn test_n_equations_empty_system() {
let mut sys = System::new();
let a = sys.add_component(make_mock(0));
let b = sys.add_component(make_mock(0));
sys.add_edge(a, b).unwrap();
sys.finalize().unwrap();
let mc = MacroComponent::new(sys);
assert_eq!(mc.n_equations(), 0);
}
#[test]
fn test_macro_component_as_trait_object() {
let sys = build_simple_internal_system();
let mc = MacroComponent::new(sys);
// Verify it can be used as Box<dyn Component>
let component: Box<dyn Component> = Box::new(mc);
assert_eq!(component.n_equations(), 6);
}
#[test]
fn test_macro_component_in_parent_system() {
// Build an internal subsystem
let internal = build_simple_internal_system();
let mc = MacroComponent::new(internal);
// Place it in a parent system alongside another component
let mut parent = System::new();
let mc_node = parent.add_component(Box::new(mc));
let other = parent.add_component(make_mock(2));
parent.add_edge(mc_node, other).unwrap();
parent.finalize().unwrap();
// Parent should have 2 nodes and 1 edge
assert_eq!(parent.node_count(), 2);
assert_eq!(parent.edge_count(), 1);
}
// ── Serialization snapshot ─────────────────────────────────────────────────
#[test]
fn test_snapshot_round_trip() {
let sys = build_simple_internal_system();
let mut mc = MacroComponent::new(sys);
mc.set_global_state_offset(0);
let port = make_connected_port("R134a", 1e5, 4e5);
mc.expose_port(0, "inlet", port);
// Fake global state with known values in the internal block [0..4]
let global_state = vec![1.5e5, 3.9e5, 8.0e4, 4.2e5];
let snap = mc
.to_snapshot(&global_state, Some("chiller_1".into()))
.unwrap();
assert_eq!(snap.label.as_deref(), Some("chiller_1"));
assert_eq!(snap.internal_edge_states.len(), 4);
assert_eq!(snap.port_names, vec!["inlet"]);
// Round-trip through JSON
let json = serde_json::to_string(&snap).unwrap();
let restored: MacroComponentSnapshot = serde_json::from_str(&json).unwrap();
assert_eq!(restored.internal_edge_states, snap.internal_edge_states);
}
}

78
crates/solver/src/solver.rs
View File

@@ -112,6 +112,9 @@ pub enum ConvergenceStatus {
/// The solver stopped due to timeout but returned the best-known state.
/// (Used by Story 4.5 — Time-Budgeted Solving.)
TimedOutWithBestState,
/// The solver converged, but one or more control variables saturated at bounds.
/// (Used by Story 5.2 — Bounded Control Variables)
ControlSaturation,
}
/// The result of a successful (or best-effort) solve.
@@ -177,7 +180,10 @@ impl ConvergedState {
/// Returns `true` if the solver fully converged (not just best-effort).
pub fn is_converged(&self) -> bool {
self.status == ConvergenceStatus::Converged
matches!(
self.status,
ConvergenceStatus::Converged | ConvergenceStatus::ControlSaturation
)
}
}
@@ -689,11 +695,13 @@ impl Solver for NewtonConfig {
);
// Get system dimensions
let n_state = system.state_vector_len();
let n_state = system.full_state_vector_len();
let n_equations: usize = system
.traverse_for_jacobian()
.map(|(_, c, _)| c.n_equations())
.sum();
.sum::<usize>()
+ system.constraints().count()
+ system.coupling_residual_count();
// Validate system
if n_state == 0 || n_equations == 0 {
@@ -759,6 +767,12 @@ impl Solver for NewtonConfig {
// Check if already converged
if current_norm < self.tolerance {
let status = if !system.saturated_variables().is_empty() {
ConvergenceStatus::ControlSaturation
} else {
ConvergenceStatus::Converged
};
// Criteria check with no prev_state (first call)
if let Some(ref criteria) = self.convergence_criteria {
let report = criteria.check(&state, None, &residuals, system);
@@ -769,11 +783,7 @@ impl Solver for NewtonConfig {
"System already converged at initial state (criteria)"
);
return Ok(ConvergedState::with_report(
state,
0,
current_norm,
ConvergenceStatus::Converged,
report,
state, 0, current_norm, status, report,
));
}
} else {
@@ -783,10 +793,7 @@ impl Solver for NewtonConfig {
"System already converged at initial state"
);
return Ok(ConvergedState::new(
state,
0,
current_norm,
ConvergenceStatus::Converged,
state, 0, current_norm, status,
));
}
}
@@ -806,7 +813,7 @@ impl Solver for NewtonConfig {
best_residual = best_residual,
"Solver timed out"
);
// Story 4.5 - AC: #2, #6: Return best state or error based on config
return self.handle_timeout(best_state, best_residual, iteration - 1, timeout);
}
@@ -847,8 +854,9 @@ impl Solver for NewtonConfig {
result.map(|_| ()).map_err(|e| format!("{:?}", e))
};
// Rather than creating a new matrix, compute it and assign
let jm = JacobianMatrix::numerical(compute_residuals_fn, &state, &residuals, 1e-8)
.map_err(|e| SolverError::InvalidSystem {
let jm =
JacobianMatrix::numerical(compute_residuals_fn, &state, &residuals, 1e-8)
.map_err(|e| SolverError::InvalidSystem {
message: format!("Failed to compute numerical Jacobian: {}", e),
})?;
// Deep copy elements to existing matrix (DMatrix::copy_from does not reallocate)
@@ -866,10 +874,7 @@ impl Solver for NewtonConfig {
frozen_count = 0;
force_recompute = false;
tracing::debug!(
iteration = iteration,
"Fresh Jacobian computed"
);
tracing::debug!(iteration = iteration, "Fresh Jacobian computed");
} else {
// Reuse the frozen Jacobian (Story 4.8 — AC: #2)
frozen_count += 1;
@@ -969,19 +974,21 @@ impl Solver for NewtonConfig {
// Check convergence (AC: #1, Story 4.7 — criteria-aware)
let converged = if let Some(ref criteria) = self.convergence_criteria {
let report = criteria.check(&state, Some(&prev_iteration_state), &residuals, system);
let report =
criteria.check(&state, Some(&prev_iteration_state), &residuals, system);
if report.is_globally_converged() {
let status = if !system.saturated_variables().is_empty() {
ConvergenceStatus::ControlSaturation
} else {
ConvergenceStatus::Converged
};
tracing::info!(
iterations = iteration,
final_residual = current_norm,
"Newton-Raphson converged (criteria)"
);
return Ok(ConvergedState::with_report(
state,
iteration,
current_norm,
ConvergenceStatus::Converged,
report,
state, iteration, current_norm, status, report,
));
}
false
@@ -990,16 +997,18 @@ impl Solver for NewtonConfig {
};
if converged {
let status = if !system.saturated_variables().is_empty() {
ConvergenceStatus::ControlSaturation
} else {
ConvergenceStatus::Converged
};
tracing::info!(
iterations = iteration,
final_residual = current_norm,
"Newton-Raphson converged"
);
return Ok(ConvergedState::new(
state,
iteration,
current_norm,
ConvergenceStatus::Converged,
state, iteration, current_norm, status,
));
}
@@ -1364,7 +1373,7 @@ impl Solver for PicardConfig {
best_residual = best_residual,
"Solver timed out"
);
// Story 4.5 - AC: #2, #6: Return best state or error based on config
return self.handle_timeout(best_state, best_residual, iteration - 1, timeout);
}
@@ -1403,7 +1412,8 @@ impl Solver for PicardConfig {
// Check convergence (AC: #1, Story 4.7 — criteria-aware)
let converged = if let Some(ref criteria) = self.convergence_criteria {
let report = criteria.check(&state, Some(&prev_iteration_state), &residuals, system);
let report =
criteria.check(&state, Some(&prev_iteration_state), &residuals, system);
if report.is_globally_converged() {
tracing::info!(
iterations = iteration,
@@ -1686,7 +1696,7 @@ impl FallbackSolver {
Ok(converged) => {
// Update best state tracking (Story 4.5 - AC: #4)
state.update_best_state(&converged.state, converged.final_residual);
tracing::info!(
solver = match state.current_solver {
CurrentSolver::Newton => "NewtonRaphson",
@@ -1701,8 +1711,8 @@ impl FallbackSolver {
}
Err(SolverError::Timeout { timeout_ms }) => {
// Story 4.5 - AC: #4: Return best state on timeout if available
if let (Some(best_state), Some(best_residual)) =
(state.best_state.clone(), state.best_residual)
if let (Some(best_state), Some(best_residual)) =
(state.best_state.clone(), state.best_residual)
{
tracing::info!(
best_residual = best_residual,

1644
crates/solver/src/system.rs
View File

File diff suppressed because it is too large Load Diff

402
crates/solver/tests/macro_component_integration.rs Normal file
View File

@@ -0,0 +1,402 @@
//! Integration tests for MacroComponent (Story 3.6).
//!
//! Tests cover:
//! - AC #1: MacroComponent implements Component trait
//! - AC #2: External ports correctly mapped to internal edges
//! - AC #3: Residuals and Jacobian delegated with proper coupling equations
//! - AC #4: Serialization snapshot round-trip
use entropyk_components::{
Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, SystemState,
};
use entropyk_solver::{MacroComponent, MacroComponentSnapshot, System};
// ─────────────────────────────────────────────────────────────────────────────
// Test helpers
// ─────────────────────────────────────────────────────────────────────────────
/// A simple zero-residual pass-through mock component.
struct PassThrough {
n_eq: usize,
}
impl Component for PassThrough {
fn compute_residuals(
&self,
_state: &SystemState,
residuals: &mut ResidualVector,
) -> Result<(), ComponentError> {
for r in residuals.iter_mut().take(self.n_eq) {
*r = 0.0;
}
Ok(())
}
fn jacobian_entries(
&self,
_state: &SystemState,
jacobian: &mut JacobianBuilder,
) -> Result<(), ComponentError> {
for i in 0..self.n_eq {
jacobian.add_entry(i, i, 1.0);
}
Ok(())
}
fn n_equations(&self) -> usize {
self.n_eq
}
fn get_ports(&self) -> &[ConnectedPort] {
&[]
}
}
fn pass(n: usize) -> Box<dyn Component> {
Box::new(PassThrough { n_eq: n })
}
fn make_port(fluid: &str, p: f64, h: f64) -> ConnectedPort {
use entropyk_components::port::{FluidId, Port};
use entropyk_core::{Enthalpy, Pressure};
let p1 = Port::new(FluidId::new(fluid), Pressure::from_pascals(p), Enthalpy::from_joules_per_kg(h));
let p2 = Port::new(FluidId::new(fluid), Pressure::from_pascals(p), Enthalpy::from_joules_per_kg(h));
p1.connect(p2).unwrap().0
}
/// Build a 4-component refrigerant cycle: A→B→C→D→A (4 edges).
fn build_4_component_cycle() -> System {
let mut sys = System::new();
let a = sys.add_component(pass(2)); // compressor
let b = sys.add_component(pass(2)); // condenser
let c = sys.add_component(pass(2)); // valve
let d = sys.add_component(pass(2)); // evaporator
sys.add_edge(a, b).unwrap();
sys.add_edge(b, c).unwrap();
sys.add_edge(c, d).unwrap();
sys.add_edge(d, a).unwrap();
sys.finalize().unwrap();
sys
}
// ─────────────────────────────────────────────────────────────────────────────
// AC #1 & #2 — MacroComponent wraps 4-component cycle correctly
// ─────────────────────────────────────────────────────────────────────────────
#[test]
fn test_4_component_cycle_macro_creation() {
let internal = build_4_component_cycle();
let mc = MacroComponent::new(internal);
// 4 components × 2 eqs = 8 internal equations, 0 exposed ports
assert_eq!(mc.n_equations(), 8,
"should have 8 internal equations with no exposed ports");
// 4 edges × 2 vars = 8 internal state vars
assert_eq!(mc.internal_state_len(), 8);
assert!(mc.get_ports().is_empty());
}
#[test]
fn test_4_component_cycle_expose_two_ports() {
let internal = build_4_component_cycle();
let mut mc = MacroComponent::new(internal);
// Expose edge 0 as "refrig_in" and edge 2 as "refrig_out"
mc.expose_port(0, "refrig_in", make_port("R134a", 1e5, 4e5));
mc.expose_port(2, "refrig_out", make_port("R134a", 5e5, 4.5e5));
// 8 internal + 4 coupling (2 per port) = 12 equations
assert_eq!(mc.n_equations(), 12,
"should have 12 equations with 2 exposed ports");
assert_eq!(mc.get_ports().len(), 2);
assert_eq!(mc.port_mappings()[0].name, "refrig_in");
assert_eq!(mc.port_mappings()[1].name, "refrig_out");
}
#[test]
fn test_4_component_cycle_in_parent_system() {
// Wrap cycle in MacroComponent and place in a parent system
let internal = build_4_component_cycle();
let mc = MacroComponent::new(internal);
let mut parent = System::new();
let _mc_node = parent.add_component(Box::new(mc));
// Single-node system (no edges) would fail validation,
// so we add a second node and an edge.
let other = parent.add_component(pass(1));
// For finalize to succeed, all nodes must have at least one edge
// (system topology requires connected nodes).
// We skip finalize here since the topology is valid (2 nodes, 1 edge).
// Actually the validation requires an edge:
parent.add_edge(_mc_node, other).unwrap();
let result = parent.finalize();
assert!(result.is_ok(), "parent finalize should succeed: {:?}", result.err());
// Parent has 2 nodes, 1 edge
assert_eq!(parent.node_count(), 2);
assert_eq!(parent.edge_count(), 1);
// Parent state vector: 1 edge × 2 = 2 state vars
assert_eq!(parent.state_vector_len(), 2);
}
// ─────────────────────────────────────────────────────────────────────────────
// AC #3 — Residuals and Jacobian delegated with coupling equations
// ─────────────────────────────────────────────────────────────────────────────
#[test]
fn test_coupling_residuals_are_zero_at_consistent_state() {
// Build cycle, expose 1 port, inject consistent external state
let internal = build_4_component_cycle();
let mut mc = MacroComponent::new(internal);
mc.expose_port(0, "refrig_in", make_port("R134a", 1e5, 4e5));
// Internal block starts at offset 2 (2 parent-edge state vars before it).
// External edge for port 0 is at (p=0, h=1).
mc.set_global_state_offset(2);
mc.set_system_context(2, &[(0, 1)]);
// State layout: [P_ext=1e5, h_ext=4e5, P_int_e0=1e5, h_int_e0=4e5, ...]
// indices: 0 1 2 3
let mut state = vec![0.0; 2 + 8]; // 2 parent + 8 internal
state[0] = 1.0e5; // P_ext
state[1] = 4.0e5; // h_ext
state[2] = 1.0e5; // P_int_e0 (consistent with port)
state[3] = 4.0e5; // h_int_e0
let n_eqs = mc.n_equations(); // 8 + 2 = 10
let mut residuals = vec![0.0; n_eqs];
mc.compute_residuals(&state, &mut residuals).unwrap();
// Coupling residuals at indices 8, 9 should be zero (consistent state)
assert!(
residuals[8].abs() < 1e-10,
"P coupling residual should be 0, got {}",
residuals[8]
);
assert!(
residuals[9].abs() < 1e-10,
"h coupling residual should be 0, got {}",
residuals[9]
);
}
#[test]
fn test_coupling_residuals_nonzero_at_inconsistent_state() {
let internal = build_4_component_cycle();
let mut mc = MacroComponent::new(internal);
mc.expose_port(0, "refrig_in", make_port("R134a", 1e5, 4e5));
mc.set_global_state_offset(2);
mc.set_system_context(2, &[(0, 1)]);
let mut state = vec![0.0; 10];
state[0] = 2.0e5; // P_ext (different from internal)
state[1] = 5.0e5; // h_ext
state[2] = 1.0e5; // P_int_e0
state[3] = 4.0e5; // h_int_e0
let n_eqs = mc.n_equations();
let mut residuals = vec![0.0; n_eqs];
mc.compute_residuals(&state, &mut residuals).unwrap();
// Coupling: r[8] = P_ext - P_int = 2e5 - 1e5 = 1e5
assert!(
(residuals[8] - 1.0e5).abs() < 1.0,
"P coupling residual mismatch: {}",
residuals[8]
);
assert!(
(residuals[9] - 1.0e5).abs() < 1.0,
"h coupling residual mismatch: {}",
residuals[9]
);
}
#[test]
fn test_jacobian_coupling_entries_correct() {
let internal = build_4_component_cycle();
let mut mc = MacroComponent::new(internal);
mc.expose_port(0, "refrig_in", make_port("R134a", 1e5, 4e5));
// external edge: (p_ext=0, h_ext=1), internal starts at offset=2
mc.set_global_state_offset(2);
mc.set_system_context(2, &[(0, 1)]);
let state = vec![0.0; 10];
let mut jac = JacobianBuilder::new();
mc.jacobian_entries(&state, &mut jac).unwrap();
let entries = jac.entries();
let find = |row: usize, col: usize| -> Option<f64> {
entries.iter().find(|&&(r, c, _)| r == row && c == col).map(|&(_, _, v)| v)
};
// Coupling rows 8 (P) and 9 (h)
assert_eq!(find(8, 0), Some(1.0), "∂r_P/∂p_ext should be +1");
assert_eq!(find(8, 2), Some(-1.0), "∂r_P/∂int_p should be -1");
assert_eq!(find(9, 1), Some(1.0), "∂r_h/∂h_ext should be +1");
assert_eq!(find(9, 3), Some(-1.0), "∂r_h/∂int_h should be -1");
}
// ─────────────────────────────────────────────────────────────────────────────
// AC #4 — Serialization snapshot
// ─────────────────────────────────────────────────────────────────────────────
#[test]
fn test_macro_component_snapshot_serialization() {
let internal = build_4_component_cycle();
let mut mc = MacroComponent::new(internal);
mc.expose_port(0, "refrig_in", make_port("R134a", 1e5, 4e5));
mc.expose_port(2, "refrig_out", make_port("R134a", 5e5, 4.5e5));
mc.set_global_state_offset(0);
// Simulate a converged global state (8 internal vars, all nonzero)
let global_state: Vec<f64> = (0..8).map(|i| (i as f64 + 1.0) * 1e4).collect();
let snap = mc
.to_snapshot(&global_state, Some("chiller_A".into()))
.expect("snapshot should succeed");
assert_eq!(snap.label.as_deref(), Some("chiller_A"));
assert_eq!(snap.internal_edge_states.len(), 8);
assert_eq!(snap.port_names, vec!["refrig_in", "refrig_out"]);
// JSON round-trip
let json = serde_json::to_string_pretty(&snap).expect("must serialize");
let restored: MacroComponentSnapshot =
serde_json::from_str(&json).expect("must deserialize");
assert_eq!(restored.label, snap.label);
assert_eq!(restored.internal_edge_states, snap.internal_edge_states);
assert_eq!(restored.port_names, snap.port_names);
}
#[test]
fn test_snapshot_fails_on_short_state() {
let internal = build_4_component_cycle();
let mut mc = MacroComponent::new(internal);
mc.set_global_state_offset(0);
// Only 4 values, but internal needs 8
let short_state = vec![0.0; 4];
let snap = mc.to_snapshot(&short_state, None);
assert!(snap.is_none(), "should return None for short state vector");
}
// ─────────────────────────────────────────────────────────────────────────────
// Two MacroComponent chillers in parallel
// ─────────────────────────────────────────────────────────────────────────────
#[test]
fn test_two_macro_chillers_in_parallel_topology() {
// Build two identical 4-component chiller MacroComponents.
let chiller_a = {
let internal = build_4_component_cycle();
let mut mc = MacroComponent::new(internal);
mc.expose_port(0, "in_a", make_port("R134a", 1e5, 4e5));
mc.expose_port(2, "out_a", make_port("R134a", 5e5, 4.5e5));
mc
};
let chiller_b = {
let internal = build_4_component_cycle();
let mut mc = MacroComponent::new(internal);
mc.expose_port(0, "in_b", make_port("R134a", 1e5, 4e5));
mc.expose_port(2, "out_b", make_port("R134a", 5e5, 4.5e5));
mc
};
// Place both into a parent system with a splitter and merger mock.
let mut parent = System::new();
let ca = parent.add_component(Box::new(chiller_a));
let cb = parent.add_component(Box::new(chiller_b));
// Simple pass-through splitter & merger
let splitter = parent.add_component(pass(1));
let merger = parent.add_component(pass(1));
// Topology: splitter → chiller_a → merger
// → chiller_b → merger
parent.add_edge(splitter, ca).unwrap();
parent.add_edge(splitter, cb).unwrap();
parent.add_edge(ca, merger).unwrap();
parent.add_edge(cb, merger).unwrap();
let result = parent.finalize();
assert!(result.is_ok(), "parallel chiller topology should finalize cleanly: {:?}", result.err());
// 4 parent edges × 2 = 8 state variables in the parent
// 2 chillers × 8 internal variables = 16 internal variables
// Total state vector length = 24
assert_eq!(parent.state_vector_len(), 24);
// 4 nodes
assert_eq!(parent.node_count(), 4);
// 4 edges
assert_eq!(parent.edge_count(), 4);
// Total equations:
// chiller_a: 8 internal + 4 coupling (2 ports) = 12
// chiller_b: 8 internal + 4 coupling (2 ports) = 12
// splitter: 1
// merger: 1
// total: 26
let total_eqs: usize = parent
.traverse_for_jacobian()
.map(|(_, c, _)| c.n_equations())
.sum();
assert_eq!(total_eqs, 26, "total equation count mismatch: {}", total_eqs);
}
#[test]
fn test_two_macro_chillers_residuals_are_computable() {
let chiller_a = {
let internal = build_4_component_cycle();
let mut mc = MacroComponent::new(internal);
mc.expose_port(0, "in_a", make_port("R134a", 1e5, 4e5));
mc.expose_port(2, "out_a", make_port("R134a", 5e5, 4.5e5));
mc
};
let chiller_b = {
let internal = build_4_component_cycle();
let mut mc = MacroComponent::new(internal);
mc.expose_port(0, "in_b", make_port("R134a", 1e5, 4e5));
mc.expose_port(2, "out_b", make_port("R134a", 5e5, 4.5e5));
mc
};
// Each chiller has 8 internal state variables (4 edges × 2)
let internal_state_len_each = chiller_a.internal_state_len(); // = 8
let mut parent = System::new();
let ca = parent.add_component(Box::new(chiller_a));
let cb = parent.add_component(Box::new(chiller_b));
let splitter = parent.add_component(pass(1));
let merger = parent.add_component(pass(1));
parent.add_edge(splitter, ca).unwrap();
parent.add_edge(splitter, cb).unwrap();
parent.add_edge(ca, merger).unwrap();
parent.add_edge(cb, merger).unwrap();
parent.finalize().unwrap();
// The parent's own state vector covers its 4 edges (8 vars).
// Each MacroComponent's internal state block starts at offsets assigned cumulatively
// by System::finalize().
// chiller_a offset = 8
// chiller_b offset = 16
// Total state len = 8 parent + 8 chiller_a + 8 chiller_b = 24 total.
let full_state_len = parent.state_vector_len();
assert_eq!(full_state_len, 24);
let state = vec![0.0; full_state_len];
let total_eqs: usize = parent
.traverse_for_jacobian()
.map(|(_, c, _)| c.n_equations())
.sum();
let mut residuals = vec![0.0; total_eqs];
let result = parent.compute_residuals(&state, &mut residuals);
assert!(
result.is_ok(),
"residual computation should not error on zero state: {:?}",
result.err()
);
}