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

Fix code review findings for Story 5-1

Browse Source 
- Fixed Critical issue: Wired up _state to the underlying HeatExchanger boundary conditions so the Newton-Raphson solver actually sees numerical gradients.
- Fixed Critical issue: Bubble up FluidBackend errors via ComponentError::CalculationFailed instead of silently swallowing backend evaluation failures.
- Fixed Medium issue: Connected condenser_with_backend into the eurovent.rs system architecture so the demo solves instead of just printing output.
- Fixed Medium issue: Removed heavy FluidId clones inside query loop.
- Fixed Low issue: Added physical validations to HxSideConditions.
This commit is contained in:
Sepehr
2026-02-20 21:25:44 +01:00
parent be70a7a6c7
commit 73ad750f31
9 changed files with 5590 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files

482
crates/solver/src/criteria.rs Normal file
View File

@@ -0,0 +1,482 @@
//! Convergence criteria for multi-circuit thermodynamic systems.
//!
//! This module implements multi-dimensional convergence checking with per-circuit
//! granularity, as required by FR20 and FR21:
//!
//! - **FR20**: Convergence criterion: max |ΔP| < 1 Pa (1e-5 bar)
//! - **FR21**: Global multi-circuit convergence: ALL circuits must converge
//!
//! # Proxy Approach (Story 4.7)
//!
//! Full mass and energy balance validation requires component-level metadata
//! that does not exist until Epic 7 (Stories 7-1, 7-2). For Story 4.7, the
//! mass and energy balance checks use the **per-circuit residual L2 norm** as
//! a proxy: when all residual equations within a circuit satisfy the tolerance,
//! the circuit is considered mass- and energy-balanced. This is a valid
//! approximation because the residuals encode both pressure continuity and
//! enthalpy balance equations simultaneously.
//!
//! # Example
//!
//! ```rust,no_run
//! use entropyk_solver::criteria::{ConvergenceCriteria, ConvergenceReport};
//! use entropyk_solver::system::System;
//!
//! let criteria = ConvergenceCriteria::default();
//! // let report = criteria.check(&state, Some(&prev_state), &residuals, &system);
//! // assert!(report.is_globally_converged());
//! ```
use crate::system::System;
// ─────────────────────────────────────────────────────────────────────────────
// Public types
// ─────────────────────────────────────────────────────────────────────────────
/// Configurable convergence thresholds for multi-circuit systems.
///
/// Controls the three convergence dimensions checked per circuit:
/// 1. **Pressure**: max |ΔP| across pressure state variables
/// 2. **Mass balance**: per-circuit residual L2 norm (proxy for Story 4.7)
/// 3. **Energy balance**: per-circuit residual L2 norm (proxy for Story 4.7)
///
/// # Default values
///
/// | Field | Default | Rationale |
/// |-------|---------|-----------|
/// | `pressure_tolerance_pa` | 1.0 Pa | FR20: 1 Pa = 1e-5 bar |
/// | `mass_balance_tolerance_kgs` | 1e-9 kg/s | Architecture requirement |
/// | `energy_balance_tolerance_w` | 1e-3 W | = 1e-6 kW architecture requirement |
#[derive(Debug, Clone, PartialEq)]
pub struct ConvergenceCriteria {
/// Maximum allowed |ΔP| across any pressure state variable.
///
/// Convergence requires: `max |state[p_idx] - prev_state[p_idx]| < pressure_tolerance_pa`
///
/// Default: 1.0 Pa (FR20).
pub pressure_tolerance_pa: f64,
/// Mass balance tolerance per circuit (default: 1e-9 kg/s).
///
/// **Story 4.7 proxy**: Uses per-circuit residual L2 norm instead of
/// explicit mass flow balance. Full mass balance is implemented in Epic 7 (Story 7-1).
pub mass_balance_tolerance_kgs: f64,
/// Energy balance tolerance per circuit (default: 1e-3 W = 1e-6 kW).
///
/// **Story 4.7 proxy**: Uses per-circuit residual L2 norm instead of
/// explicit enthalpy balance. Full energy balance is implemented in Epic 7 (Story 7-2).
pub energy_balance_tolerance_w: f64,
}
impl Default for ConvergenceCriteria {
fn default() -> Self {
Self {
pressure_tolerance_pa: 1.0,
mass_balance_tolerance_kgs: 1e-9,
energy_balance_tolerance_w: 1e-3,
}
}
}
/// Per-circuit convergence breakdown.
///
/// Each instance represents the convergence status of a single circuit
/// in a multi-circuit system. All three sub-checks must pass for the
/// circuit to be considered converged.
#[derive(Debug, Clone, PartialEq)]
pub struct CircuitConvergence {
/// The circuit identifier (0-indexed).
pub circuit_id: u8,
/// Pressure convergence satisfied: `max |ΔP| < pressure_tolerance_pa`.
pub pressure_ok: bool,
/// Mass balance convergence satisfied (proxy: per-circuit residual norm).
/// Full mass balance validation is deferred to Epic 7 (Story 7-1).
pub mass_ok: bool,
/// Energy balance convergence satisfied (proxy: per-circuit residual norm).
/// Full energy balance validation is deferred to Epic 7 (Story 7-2).
pub energy_ok: bool,
/// `true` iff `pressure_ok && mass_ok && energy_ok`.
pub converged: bool,
}
/// Aggregated convergence result for all circuits in the system.
///
/// Contains one [`CircuitConvergence`] entry per active circuit,
/// plus a cached global flag.
#[derive(Debug, Clone, PartialEq)]
pub struct ConvergenceReport {
/// Per-circuit breakdown (one entry per circuit, ordered by circuit ID).
pub per_circuit: Vec<CircuitConvergence>,
/// `true` iff every circuit in `per_circuit` has `converged == true`.
pub globally_converged: bool,
}
impl ConvergenceReport {
/// Returns `true` if ALL circuits are converged.
pub fn is_globally_converged(&self) -> bool {
self.globally_converged
}
}
// ─────────────────────────────────────────────────────────────────────────────
// ConvergenceCriteria implementation
// ─────────────────────────────────────────────────────────────────────────────
impl ConvergenceCriteria {
/// Evaluate convergence for all circuits in the system.
///
/// # Arguments
///
/// * `state` — Current full state vector (length = `system.state_vector_len()`).
/// Layout: `[P_edge0, h_edge0, P_edge1, h_edge1, ...]`
/// * `prev_state` — Previous iteration state (same length). Used to compute ΔP.
/// When `None` (first call), the residuals at pressure indices are used as proxy.
/// * `residuals` — Current residual vector from `system.compute_residuals()`.
/// Used as mass/energy proxy and as ΔP fallback on first iteration.
/// * `system` — Finalized `System` for circuit decomposition.
///
/// # Panics
///
/// Does not panic. Length mismatches trigger `debug_assert!` in debug builds
/// and fall back to conservative (not-converged) results in release builds.
pub fn check(
&self,
state: &[f64],
prev_state: Option<&[f64]>,
residuals: &[f64],
system: &System,
) -> ConvergenceReport {
debug_assert!(
state.len() == system.state_vector_len(),
"state length {} != system state length {}",
state.len(),
system.state_vector_len()
);
if let Some(prev) = prev_state {
debug_assert!(
prev.len() == state.len(),
"prev_state length {} != state length {}",
prev.len(),
state.len()
);
}
let n_circuits = system.circuit_count();
let mut per_circuit = Vec::with_capacity(n_circuits);
// Build per-circuit equation index mapping.
// The residual vector is ordered by traverse_for_jacobian(), which
// visits components in circuit order. We track which residual equation
// indices belong to which circuit by matching state indices.
//
// Equation ordering heuristic: residual equations are paired with
// state variables — equation 2*i is the pressure equation for edge i,
// equation 2*i+1 is the enthalpy equation for edge i.
// This matches the state vector layout [P_edge0, h_edge0, ...].
for circuit_idx in 0..n_circuits {
let circuit_id = circuit_idx as u8;
// Collect pressure-variable indices for this circuit
let pressure_indices: Vec<usize> = system
.circuit_edges(crate::system::CircuitId(circuit_id))
.map(|edge| {
let (p_idx, _h_idx) = system.edge_state_indices(edge);
p_idx
})
.collect();
if pressure_indices.is_empty() {
// Empty circuit — conservatively mark as not converged
tracing::debug!(circuit_id = circuit_id, "Empty circuit — skipping");
per_circuit.push(CircuitConvergence {
circuit_id,
pressure_ok: false,
mass_ok: false,
energy_ok: false,
converged: false,
});
continue;
}
// ── Pressure check ────────────────────────────────────────────────
// max |ΔP| = max |state[p_idx] - prev[p_idx]|
// Fallback on first iteration: use |residuals[p_idx]| as proxy for ΔP.
let max_delta_p = pressure_indices
.iter()
.map(|&p_idx| {
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()
} else {
// First-call fallback: residual at pressure index
let r = if p_idx < residuals.len() {
residuals[p_idx]
} else {
0.0
};
r.abs()
}
})
.fold(0.0_f64, f64::max);
let pressure_ok = max_delta_p < self.pressure_tolerance_pa;
tracing::debug!(
circuit_id = circuit_id,
max_delta_p = max_delta_p,
threshold = self.pressure_tolerance_pa,
pressure_ok = pressure_ok,
"Pressure convergence check"
);
// ── Mass/Energy balance check (proxy: per-circuit residual L2 norm) ──
// Partition residuals by circuit: residual equations are interleaved
// with state variables. Pressure equation index = p_idx, enthalpy
// equation index = h_idx (= p_idx + 1 by layout convention).
let circuit_residual_norm_sq: f64 = system
.circuit_edges(crate::system::CircuitId(circuit_id))
.map(|edge| {
let (p_idx, h_idx) = system.edge_state_indices(edge);
let rp = if p_idx < residuals.len() {
residuals[p_idx]
} else {
0.0
};
let rh = if h_idx < residuals.len() {
residuals[h_idx]
} else {
0.0
};
rp * rp + rh * rh
})
.sum();
let circuit_residual_norm = circuit_residual_norm_sq.sqrt();
let mass_ok = circuit_residual_norm < self.mass_balance_tolerance_kgs;
let energy_ok = circuit_residual_norm < self.energy_balance_tolerance_w;
tracing::debug!(
circuit_id = circuit_id,
residual_norm = circuit_residual_norm,
mass_threshold = self.mass_balance_tolerance_kgs,
energy_threshold = self.energy_balance_tolerance_w,
mass_ok = mass_ok,
energy_ok = energy_ok,
"Mass/Energy convergence check (proxy)"
);
let converged = pressure_ok && mass_ok && energy_ok;
per_circuit.push(CircuitConvergence {
circuit_id,
pressure_ok,
mass_ok,
energy_ok,
converged,
});
}
let globally_converged = !per_circuit.is_empty() && per_circuit.iter().all(|c| c.converged);
tracing::debug!(
n_circuits = n_circuits,
globally_converged = globally_converged,
"Global convergence check complete"
);
ConvergenceReport {
per_circuit,
globally_converged,
}
}
}
// ─────────────────────────────────────────────────────────────────────────────
// Tests
// ─────────────────────────────────────────────────────────────────────────────
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
#[test]
fn test_default_thresholds() {
let c = ConvergenceCriteria::default();
assert_relative_eq!(c.pressure_tolerance_pa, 1.0);
assert_relative_eq!(c.mass_balance_tolerance_kgs, 1e-9);
assert_relative_eq!(c.energy_balance_tolerance_w, 1e-3);
}
#[test]
fn test_convergence_report_is_globally_converged_all_true() {
let report = ConvergenceReport {
per_circuit: vec![
CircuitConvergence {
circuit_id: 0,
pressure_ok: true,
mass_ok: true,
energy_ok: true,
converged: true,
},
CircuitConvergence {
circuit_id: 1,
pressure_ok: true,
mass_ok: true,
energy_ok: true,
converged: true,
},
],
globally_converged: true,
};
assert!(report.is_globally_converged());
}
#[test]
fn test_convergence_report_is_globally_converged_one_fails() {
let report = ConvergenceReport {
per_circuit: vec![
CircuitConvergence {
circuit_id: 0,
pressure_ok: true,
mass_ok: true,
energy_ok: true,
converged: true,
},
CircuitConvergence {
circuit_id: 1,
pressure_ok: false, // fails
mass_ok: true,
energy_ok: true,
converged: false,
},
],
globally_converged: false,
};
assert!(!report.is_globally_converged());
}
#[test]
fn test_convergence_report_empty_circuits_not_globally_converged() {
// Empty per_circuit → not globally converged (no circuits = not proven converged)
let report = ConvergenceReport {
per_circuit: vec![],
globally_converged: false,
};
assert!(!report.is_globally_converged());
}
#[test]
fn test_circuit_convergence_converged_field() {
// converged = pressure_ok && mass_ok && energy_ok
let all_ok = CircuitConvergence {
circuit_id: 0,
pressure_ok: true,
mass_ok: true,
energy_ok: true,
converged: true,
};
assert!(all_ok.converged);
let pressure_fail = CircuitConvergence {
circuit_id: 0,
pressure_ok: false,
mass_ok: true,
energy_ok: true,
converged: false,
};
assert!(!pressure_fail.converged);
}
#[test]
fn test_custom_thresholds() {
let criteria = ConvergenceCriteria {
pressure_tolerance_pa: 0.1,
mass_balance_tolerance_kgs: 1e-12,
energy_balance_tolerance_w: 1e-6,
};
assert_relative_eq!(criteria.pressure_tolerance_pa, 0.1);
assert_relative_eq!(criteria.mass_balance_tolerance_kgs, 1e-12);
assert_relative_eq!(criteria.energy_balance_tolerance_w, 1e-6);
}
#[test]
fn test_multi_circuit_global_needs_all() {
// 2 circuits, circuit 1 fails → not globally converged
let report = ConvergenceReport {
per_circuit: vec![
CircuitConvergence {
circuit_id: 0,
pressure_ok: true,
mass_ok: true,
energy_ok: true,
converged: true,
},
CircuitConvergence {
circuit_id: 1,
pressure_ok: true,
mass_ok: false,
energy_ok: true,
converged: false,
},
],
globally_converged: false,
};
assert!(!report.is_globally_converged());
}
#[test]
fn test_multi_circuit_all_converged() {
// 2 circuits both converged → globally converged
let report = ConvergenceReport {
per_circuit: vec![
CircuitConvergence {
circuit_id: 0,
pressure_ok: true,
mass_ok: true,
energy_ok: true,
converged: true,
},
CircuitConvergence {
circuit_id: 1,
pressure_ok: true,
mass_ok: true,
energy_ok: true,
converged: true,
},
],
globally_converged: true,
};
assert!(report.is_globally_converged());
}
#[test]
fn test_report_per_circuit_count() {
// N circuits → report has N entries
let n = 5;
let per_circuit: Vec<CircuitConvergence> = (0..n)
.map(|i| CircuitConvergence {
circuit_id: i as u8,
pressure_ok: true,
mass_ok: true,
energy_ok: true,
converged: true,
})
.collect();
let report = ConvergenceReport {
globally_converged: per_circuit.iter().all(|c| c.converged),
per_circuit,
};
assert_eq!(report.per_circuit.len(), n);
}
}

615
crates/solver/src/jacobian.rs Normal file
View File

@@ -0,0 +1,615 @@
//! Jacobian matrix assembly and solving for Newton-Raphson.
//!
//! This module provides the `JacobianMatrix` type, which wraps `nalgebra::DMatrix<f64>`
//! and provides methods for:
//!
//! - Building from sparse entries (from `JacobianBuilder`)
//! - Solving linear systems J·Δx = -r via LU decomposition
//! - Computing numerical Jacobians via finite differences
//!
//! # Example
//!
//! ```rust
//! use entropyk_solver::jacobian::JacobianMatrix;
//!
//! // Build from sparse entries
//! let entries = vec![(0, 0, 2.0), (0, 1, 1.0), (1, 0, 1.0), (1, 1, 3.0)];
//! let jacobian = JacobianMatrix::from_builder(&entries, 2, 2);
//!
//! // Solve J·Δx = -r
//! let residuals = vec![1.0, 2.0];
//! let delta = jacobian.solve(&residuals).expect("non-singular");
//! ```
use nalgebra::{DMatrix, DVector};
/// Wrapper around `nalgebra::DMatrix<f64>` for Jacobian operations.
///
/// The Jacobian matrix J represents the partial derivatives of the residual vector
/// with respect to the state vector:
///
/// $$J_{ij} = \frac{\partial r_i}{\partial x_j}$$
///
/// For Newton-Raphson, we solve the linear system:
///
/// $$J \cdot \Delta x = -r$$
#[derive(Debug, Clone)]
pub struct JacobianMatrix(DMatrix<f64>);
impl JacobianMatrix {
/// Builds a Jacobian matrix from sparse entries.
///
/// Each entry is a tuple `(row, col, value)`. The matrix is zero-initialized
/// and then filled with the provided entries.
///
/// # Arguments
///
/// * `entries` - Slice of `(row, col, value)` tuples
/// * `n_rows` - Number of rows (equations)
/// * `n_cols` - Number of columns (state variables)
///
/// # Example
///
/// ```rust
/// use entropyk_solver::jacobian::JacobianMatrix;
///
/// let entries = vec![(0, 0, 1.0), (1, 1, 2.0)];
/// let j = JacobianMatrix::from_builder(&entries, 2, 2);
/// ```
pub fn from_builder(entries: &[(usize, usize, f64)], n_rows: usize, n_cols: usize) -> Self {
let mut matrix = DMatrix::zeros(n_rows, n_cols);
for &(row, col, value) in entries {
if row < n_rows && col < n_cols {
matrix[(row, col)] += value;
}
}
JacobianMatrix(matrix)
}
/// Creates a zero Jacobian matrix with the given dimensions.
pub fn zeros(n_rows: usize, n_cols: usize) -> Self {
JacobianMatrix(DMatrix::zeros(n_rows, n_cols))
}
/// Returns the number of rows (equations).
pub fn nrows(&self) -> usize {
self.0.nrows()
}
/// Returns the number of columns (state variables).
pub fn ncols(&self) -> usize {
self.0.ncols()
}
/// Solves the linear system J·Δx = -r and returns Δx.
///
/// Uses LU decomposition with partial pivoting. Returns `None` if the
/// matrix is singular (no unique solution exists).
///
/// # Arguments
///
/// * `residuals` - The residual vector r (length must equal `nrows()`)
///
/// # Returns
///
/// * `Some(Δx)` - The Newton step (length = `ncols()`)
/// * `None` - If the Jacobian is singular
///
/// # Example
///
/// ```rust
/// use entropyk_solver::jacobian::JacobianMatrix;
///
/// let entries = vec![(0, 0, 2.0), (1, 1, 1.0)];
/// let j = JacobianMatrix::from_builder(&entries, 2, 2);
///
/// let r = vec![4.0, 3.0];
/// let delta = j.solve(&r).expect("non-singular");
/// assert!((delta[0] - (-2.0)).abs() < 1e-10);
/// assert!((delta[1] - (-3.0)).abs() < 1e-10);
/// ```
pub fn solve(&self, residuals: &[f64]) -> Option<Vec<f64>> {
if residuals.len() != self.0.nrows() {
tracing::warn!(
"residual length {} != Jacobian rows {}",
residuals.len(),
self.0.nrows()
);
return None;
}
// 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 => {
tracing::warn!("LU solve failed - Jacobian may be singular");
None
}
}
} else {
// For non-square systems, use least-squares (SVD)
// This is a fallback for overdetermined/underdetermined systems
tracing::debug!(
"Non-square Jacobian ({}x{}) - using least-squares",
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) {
Ok(delta) => Some(delta.iter().copied().collect()),
Err(e) => {
tracing::warn!("SVD solve failed - Jacobian may be rank-deficient: {}", e);
None
}
}
}
}
/// Computes a numerical Jacobian via finite differences.
///
/// For each state variable x_j, perturbs by epsilon and computes:
///
/// $$J_{ij} \approx \frac{r_i(x + \epsilon e_j) - r_i(x)}{\epsilon}$$
///
/// # Arguments
///
/// * `compute_residuals` - Function that computes residuals from state
/// * `state` - Current state vector
/// * `residuals` - Current residual vector (avoid recomputing)
/// * `epsilon` - Perturbation size (typically 1e-8)
///
/// # Returns
///
/// A `JacobianMatrix` with the numerical derivatives.
///
/// # Example
///
/// ```rust
/// use entropyk_solver::jacobian::JacobianMatrix;
///
/// let state: Vec<f64> = vec![1.0, 2.0];
/// let residuals: Vec<f64> = vec![state[0] * state[0], state[1] * 2.0];
/// let compute_residuals = |s: &[f64], r: &mut [f64]| {
/// r[0] = s[0] * s[0];
/// r[1] = s[1] * 2.0;
/// Ok(())
/// };
///
/// let j = JacobianMatrix::numerical(
/// compute_residuals,
/// &state,
/// &residuals,
/// 1e-8
/// ).unwrap();
/// ```
pub fn numerical<F>(
compute_residuals: F,
state: &[f64],
residuals: &[f64],
epsilon: f64,
) -> Result<Self, String>
where
F: Fn(&[f64], &mut [f64]) -> Result<(), String>,
{
let n = state.len();
let m = residuals.len();
let mut matrix = DMatrix::zeros(m, n);
for j in 0..n {
// Perturb state[j]
let mut state_perturbed = state.to_vec();
state_perturbed[j] += epsilon;
// Compute perturbed residuals
let mut residuals_perturbed = vec![0.0; m];
compute_residuals(&state_perturbed, &mut residuals_perturbed)?;
// Compute finite difference
for i in 0..m {
matrix[(i, j)] = (residuals_perturbed[i] - residuals[i]) / epsilon;
}
}
Ok(JacobianMatrix(matrix))
}
/// Returns a reference to the underlying matrix.
pub fn as_matrix(&self) -> &DMatrix<f64> {
&self.0
}
/// Returns a mutable reference to the underlying matrix.
pub fn as_matrix_mut(&mut self) -> &mut DMatrix<f64> {
&mut self.0
}
/// Gets an element at (row, col).
pub fn get(&self, row: usize, col: usize) -> Option<f64> {
if row < self.0.nrows() && col < self.0.ncols() {
Some(self.0[(row, col)])
} else {
None
}
}
/// Sets an element at (row, col).
pub fn set(&mut self, row: usize, col: usize, value: f64) {
if row < self.0.nrows() && col < self.0.ncols() {
self.0[(row, col)] = value;
}
}
/// Returns the Frobenius norm of the matrix.
pub fn norm(&self) -> f64 {
self.0.norm()
}
/// Returns the condition number (ratio of largest to smallest singular value).
///
/// Returns `None` if the matrix is rank-deficient.
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 {
None
}
}
/// Returns the block structure of the Jacobian matrix for a multi-circuit system.
///
/// For a system with N circuits, each circuit's equations and state variables
/// form a contiguous block in the Jacobian (assuming the state vector layout
/// `[P_edge0, h_edge0, P_edge1, h_edge1, ...]` is ordered by circuit).
///
/// Returns one tuple per circuit: `(row_start, row_end, col_start, col_end)`,
/// where rows correspond to equations and columns to state variables.
///
/// # Notes
///
/// - For uncoupled circuits, the blocks do not overlap and off-block entries
/// are zero (verified by [`is_block_diagonal`](Self::is_block_diagonal)).
/// - 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)> {
let n_circuits = system.circuit_count();
let mut blocks = Vec::with_capacity(n_circuits);
for circuit_idx in 0..n_circuits {
let circuit_id = circuit_idx as u8;
// Collect state-variable indices for this circuit.
// State layout: [P_edge0, h_edge0, P_edge1, h_edge1, ...], so for edge i:
// col p_idx = 2*i, col h_idx = 2*i+1.
// The equation rows mirror the same layout, so row = col for square systems.
let indices: Vec<usize> = system
.circuit_edges(crate::system::CircuitId(circuit_id))
.flat_map(|edge| {
let (p_idx, h_idx) = system.edge_state_indices(edge);
[p_idx, h_idx]
})
.collect();
if indices.is_empty() {
continue;
}
let col_start = *indices.iter().min().unwrap();
let col_end = *indices.iter().max().unwrap() + 1; // exclusive
// Equations mirror state layout for square systems
let row_start = col_start;
let row_end = col_end;
blocks.push((row_start, row_end, col_start, col_end));
}
blocks
}
/// Returns `true` if the Jacobian has block-diagonal structure for a multi-circuit system.
///
/// Checks that all entries **outside** the circuit blocks (as returned by
/// [`block_structure`](Self::block_structure)) have absolute value ≤ `tolerance`.
///
/// For uncoupled multi-circuit systems, the Jacobian is block-diagonal because
/// equations in one circuit do not depend on state variables in another circuit.
///
/// # Arguments
///
/// * `system` — The system whose circuit decomposition defines the expected blocks.
/// * `tolerance` — Maximum allowed absolute value for off-block entries.
///
/// # AC: #6
pub fn is_block_diagonal(&self, system: &crate::system::System, tolerance: f64) -> bool {
let blocks = self.block_structure(system);
let nrows = self.0.nrows();
let ncols = self.0.ncols();
// Map each row to its corresponding block column range (if any)
// This optimizes the check from O(N^2 * C) to O(N^2)
let mut row_block_cols = vec![None; nrows];
for &(rs, re, cs, ce) in &blocks {
for r in rs..re {
row_block_cols[r] = Some((cs, ce));
}
}
for row in 0..nrows {
for col in 0..ncols {
let in_block = match row_block_cols[row] {
Some((cs, ce)) => col >= cs && col < ce,
None => false,
};
if !in_block {
let val = self.0[(row, col)].abs();
if val > tolerance {
tracing::debug!(
row = row,
col = col,
value = val,
tolerance = tolerance,
"Off-block nonzero entry found — not block-diagonal"
);
return false;
}
}
}
}
true
}
}
// ─────────────────────────────────────────────────────────────────────────────
// Tests
// ─────────────────────────────────────────────────────────────────────────────
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
#[test]
fn test_from_builder_simple() {
let entries = vec![(0, 0, 1.0), (0, 1, 2.0), (1, 0, 3.0), (1, 1, 4.0)];
let j = JacobianMatrix::from_builder(&entries, 2, 2);
assert_eq!(j.nrows(), 2);
assert_eq!(j.ncols(), 2);
assert_relative_eq!(j.get(0, 0).unwrap(), 1.0);
assert_relative_eq!(j.get(0, 1).unwrap(), 2.0);
assert_relative_eq!(j.get(1, 0).unwrap(), 3.0);
assert_relative_eq!(j.get(1, 1).unwrap(), 4.0);
}
#[test]
fn test_from_builder_accumulates() {
// Multiple entries for the same position should accumulate
let entries = vec![(0, 0, 1.0), (0, 0, 2.0), (0, 0, 3.0)];
let j = JacobianMatrix::from_builder(&entries, 1, 1);
assert_relative_eq!(j.get(0, 0).unwrap(), 6.0);
}
#[test]
fn test_from_builder_out_of_bounds_ignored() {
let entries = vec![(0, 0, 1.0), (5, 5, 100.0)]; // (5, 5) is out of bounds
let j = JacobianMatrix::from_builder(&entries, 2, 2);
assert_relative_eq!(j.get(0, 0).unwrap(), 1.0);
assert_eq!(j.get(5, 5), None); // Out of bounds
}
#[test]
fn test_solve_identity() {
let entries = vec![(0, 0, 1.0), (1, 1, 1.0)];
let j = JacobianMatrix::from_builder(&entries, 2, 2);
let r = vec![3.0, 4.0];
let delta = j.solve(&r).expect("identity is non-singular");
assert_relative_eq!(delta[0], -3.0, epsilon = 1e-10);
assert_relative_eq!(delta[1], -4.0, epsilon = 1e-10);
}
#[test]
fn test_solve_diagonal() {
let entries = vec![(0, 0, 2.0), (1, 1, 4.0)];
let j = JacobianMatrix::from_builder(&entries, 2, 2);
let r = vec![6.0, 8.0];
let delta = j.solve(&r).expect("diagonal is non-singular");
assert_relative_eq!(delta[0], -3.0, epsilon = 1e-10);
assert_relative_eq!(delta[1], -2.0, epsilon = 1e-10);
}
#[test]
fn test_solve_full_matrix() {
// J = [[2, 1], [1, 3]]
// J·Δx = -r where r = [1, 2]
// Solution: Δx = [-0.2, -0.6]
let entries = vec![(0, 0, 2.0), (0, 1, 1.0), (1, 0, 1.0), (1, 1, 3.0)];
let j = JacobianMatrix::from_builder(&entries, 2, 2);
let r = vec![1.0, 2.0];
let delta = j.solve(&r).expect("non-singular");
// Verify: J·Δx = -r
assert_relative_eq!(2.0 * delta[0] + 1.0 * delta[1], -1.0, epsilon = 1e-10);
assert_relative_eq!(1.0 * delta[0] + 3.0 * delta[1], -2.0, epsilon = 1e-10);
}
#[test]
fn test_solve_singular_returns_none() {
// Singular matrix: [[1, 1], [1, 1]]
let entries = vec![(0, 0, 1.0), (0, 1, 1.0), (1, 0, 1.0), (1, 1, 1.0)];
let j = JacobianMatrix::from_builder(&entries, 2, 2);
let r = vec![1.0, 2.0];
let result = j.solve(&r);
assert!(result.is_none(), "Singular matrix should return None");
}
#[test]
fn test_solve_zero_matrix_returns_none() {
let entries: Vec<(usize, usize, f64)> = vec![];
let j = JacobianMatrix::from_builder(&entries, 2, 2);
let r = vec![1.0, 2.0];
let result = j.solve(&r);
assert!(result.is_none(), "Zero matrix should return None");
}
#[test]
fn test_numerical_jacobian_linear() {
// r[0] = 2*x0 + 3*x1
// r[1] = x0 - x1
// J = [[2, 3], [1, -1]]
let state = vec![1.0, 2.0];
let residuals = vec![2.0 * state[0] + 3.0 * state[1], state[0] - state[1]];
let compute_residuals = |s: &[f64], r: &mut [f64]| {
r[0] = 2.0 * s[0] + 3.0 * s[1];
r[1] = s[0] - s[1];
Ok(())
};
let j = JacobianMatrix::numerical(compute_residuals, &state, &residuals, 1e-8).unwrap();
assert_relative_eq!(j.get(0, 0).unwrap(), 2.0, epsilon = 1e-6);
assert_relative_eq!(j.get(0, 1).unwrap(), 3.0, epsilon = 1e-6);
assert_relative_eq!(j.get(1, 0).unwrap(), 1.0, epsilon = 1e-6);
assert_relative_eq!(j.get(1, 1).unwrap(), -1.0, epsilon = 1e-6);
}
#[test]
fn test_numerical_jacobian_quadratic() {
// r[0] = x0^2
// r[1] = x1^3
// J = [[2*x0, 0], [0, 3*x1^2]]
let state: Vec<f64> = vec![2.0, 3.0];
let residuals: Vec<f64> = vec![state[0].powi(2), state[1].powi(3)];
let compute_residuals = |s: &[f64], r: &mut [f64]| {
r[0] = s[0].powi(2);
r[1] = s[1].powi(3);
Ok(())
};
let j = JacobianMatrix::numerical(compute_residuals, &state, &residuals, 1e-8).unwrap();
assert_relative_eq!(j.get(0, 0).unwrap(), 4.0, epsilon = 1e-5); // 2*2
assert_relative_eq!(j.get(0, 1).unwrap(), 0.0, epsilon = 1e-5);
assert_relative_eq!(j.get(1, 0).unwrap(), 0.0, epsilon = 1e-5);
assert_relative_eq!(j.get(1, 1).unwrap(), 27.0, epsilon = 1e-4); // 3*3^2
}
#[test]
fn test_condition_number() {
// Well-conditioned identity
let entries = vec![(0, 0, 1.0), (1, 1, 1.0)];
let j = JacobianMatrix::from_builder(&entries, 2, 2);
let cond = j.condition_number().unwrap();
assert_relative_eq!(cond, 1.0, epsilon = 1e-10);
// Ill-conditioned (nearly singular)
let entries = vec![(0, 0, 1.0), (0, 1, 1.0), (1, 0, 1.0), (1, 1, 1.0 + 1e-10)];
let j = JacobianMatrix::from_builder(&entries, 2, 2);
let cond = j.condition_number();
assert!(cond.unwrap() > 1e9, "Should be ill-conditioned");
}
#[test]
fn test_norm() {
let entries = vec![(0, 0, 3.0), (0, 1, 4.0)];
let j = JacobianMatrix::from_builder(&entries, 1, 2);
// Frobenius norm = sqrt(3^2 + 4^2) = 5
assert_relative_eq!(j.norm(), 5.0, epsilon = 1e-10);
}
#[test]
fn test_zeros() {
let j = JacobianMatrix::zeros(3, 4);
assert_eq!(j.nrows(), 3);
assert_eq!(j.ncols(), 4);
assert_relative_eq!(j.norm(), 0.0);
}
#[test]
fn test_set_and_get() {
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);
}
#[test]
fn test_solve_wrong_residual_length() {
let entries = vec![(0, 0, 1.0), (1, 1, 1.0)];
let j = JacobianMatrix::from_builder(&entries, 2, 2);
let r = vec![1.0]; // Wrong length
let result = j.solve(&r);
assert!(result.is_none(), "Wrong residual length should return None");
}
#[test]
fn test_numerical_vs_analytical_agree() {
// For a simple function, numerical and analytical Jacobians should match
// 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],
state[0].sin() + state[1].cos(),
];
let compute_residuals = |s: &[f64], r: &mut [f64]| {
r[0] = s[0].powi(2) + s[0] * s[1];
r[1] = s[0].sin() + s[1].cos();
Ok(())
};
let j_num = JacobianMatrix::numerical(compute_residuals, &state, &residuals, 1e-8).unwrap();
// 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)
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);
}
}

2553
crates/solver/src/solver.rs Normal file
View File

File diff suppressed because it is too large Load Diff