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

Fix code review issues for Story 1.10

Browse Source
This commit is contained in:
Sepehr
2026-02-21 19:15:13 +01:00
parent 0d9a0e4231
commit 8ef8cd2eba
4 changed files with 303 additions and 97 deletions
Download Patch File Download Diff File Expand all files Collapse all files

114
crates/components/src/pipe.rs
View File

@@ -140,16 +140,6 @@ impl PipeGeometry {
/// Friction factor calculation methods.
pub mod friction_factor {
use entropyk_core::MIN_MASS_FLOW_REGULARIZATION_KG_S;
/// Minimum Reynolds number for zero-flow regularization.
///
/// Reynolds is dimensionless (Re = ρvD/μ), so MIN_REYNOLDS = 1.0 is physically reasonable
/// for preventing division by zero. This is independent of [`MIN_MASS_FLOW_REGULARIZATION_KG_S`]
/// which applies to mass flow (kg/s). Both serve the same purpose: avoiding NaN/Inf in denominators.
///
/// [`MIN_MASS_FLOW_REGULARIZATION_KG_S`]: entropyk_core::MIN_MASS_FLOW_REGULARIZATION_KG_S
const MIN_REYNOLDS: f64 = 1.0;
/// Calculates the Haaland friction factor for turbulent flow.
///
@@ -164,65 +154,28 @@ pub mod friction_factor {
/// # Returns
///
/// Darcy friction factor f
///
/// # Zero-flow regularization
///
/// Re is clamped to at least `MIN_REYNOLDS` so that divisions (64/Re, 6.9/Re) never cause NaN/Inf.
pub fn haaland(relative_roughness: f64, reynolds: f64) -> f64 {
if reynolds <= 0.0 {
return 0.02; // Default for invalid input
}
let reynolds = reynolds.max(MIN_REYNOLDS);
// Laminar flow: f = 64/Re
// Do not clamp Reynolds number here to preserve linear pressure drop near zero flow.
if reynolds < 2300.0 {
return 64.0 / reynolds;
}
// Prevent division by zero or negative values in log
let re_clamped = reynolds.max(1.0);
// Haaland equation (turbulent)
// 1/√f = -1.8 × log10[(ε/D/3.7)^1.11 + 6.9/Re]
let term1 = (relative_roughness / 3.7).powf(1.11);
let term2 = 6.9 / reynolds;
let term2 = 6.9 / re_clamped;
let inv_sqrt_f = -1.8 * (term1 + term2).log10();
1.0 / (inv_sqrt_f * inv_sqrt_f)
}
/// Calculates the Swamee-Jain friction factor (alternative to Haaland).
///
/// Explicit approximation valid for:
/// - 10^-6 < ε/D < 10^-2
/// - 5000 < Re < 10^8
///
/// # Zero-flow regularization
///
/// Re is clamped to at least `MIN_REYNOLDS` so that divisions by Re never cause NaN/Inf.
pub fn swamee_jain(relative_roughness: f64, reynolds: f64) -> f64 {
if reynolds <= 0.0 {
return 0.02;
}
let reynolds = reynolds.max(MIN_REYNOLDS);
if reynolds < 2300.0 {
return 64.0 / reynolds;
}
let term1 = relative_roughness / 3.7;
let term2 = 5.74 / reynolds.powf(0.9);
let log_term = (term1 + term2).log10();
0.25 / (log_term * log_term)
}
/// Simple friction factor for quick estimates.
///
/// Returns f ≈ 0.02 for turbulent flow (typical for commercial pipes).
pub fn simplified(_relative_roughness: f64, reynolds: f64) -> f64 {
if reynolds < 2300.0 {
return 64.0 / reynolds.max(1.0);
}
0.02
}
}
/// A pipe component with pressure drop calculation.
@@ -510,17 +463,24 @@ impl Pipe<Connected> {
///
/// Pressure drop in Pascals (positive value)
pub fn pressure_drop(&self, flow_m3_per_s: f64) -> f64 {
if flow_m3_per_s <= 0.0 {
let abs_flow = flow_m3_per_s.abs();
if abs_flow <= std::f64::EPSILON {
return 0.0;
}
let velocity = self.velocity(flow_m3_per_s);
let f = self.friction_factor(flow_m3_per_s);
let velocity = self.velocity(abs_flow);
let f = self.friction_factor(abs_flow);
let ld = self.geometry.ld_ratio();
// Darcy-Weisbach nominal: ΔP_nominal = f × (L/D) × (ρ × v² / 2); ΔP_eff = f_dp × ΔP_nominal
let dp_nominal = f * ld * self.fluid_density_kg_per_m3 * velocity * velocity / 2.0;
dp_nominal * self.calib.f_dp
let dp = dp_nominal * self.calib.f_dp;
if flow_m3_per_s < 0.0 {
-dp
} else {
dp
}
}
/// Calculates mass flow from volumetric flow.
@@ -580,12 +540,17 @@ impl Component for Pipe<Connected> {
match self.operational_state {
OperationalState::Off => {
// Blocked pipe: no flow
if state.is_empty() {
return Err(ComponentError::InvalidStateDimensions { expected: 1, actual: 0 });
}
residuals[0] = state[0];
return Ok(());
}
OperationalState::Bypass => {
// No pressure drop (perfect pipe)
residuals[0] = 0.0;
let p_in = self.port_inlet.pressure().to_pascals();
let p_out = self.port_outlet.pressure().to_pascals();
residuals[0] = p_in - p_out;
return Ok(());
}
OperationalState::On => {}
@@ -620,6 +585,18 @@ impl Component for Pipe<Connected> {
state: &SystemState,
jacobian: &mut JacobianBuilder,
) -> Result<(), ComponentError> {
match self.operational_state {
OperationalState::Off => {
jacobian.add_entry(0, 0, 1.0);
return Ok(());
}
OperationalState::Bypass => {
jacobian.add_entry(0, 0, 0.0);
return Ok(());
}
OperationalState::On => {}
}
if state.is_empty() {
return Err(ComponentError::InvalidStateDimensions {
expected: 1,
@@ -631,9 +608,9 @@ impl Component for Pipe<Connected> {
let flow_m3_s = mass_flow_kg_s / self.fluid_density_kg_per_m3;
// Numerical derivative of pressure drop with respect to mass flow
let h = 0.001;
let h = 1e-6_f64.max(mass_flow_kg_s.abs() * 1e-5);
let dp_plus = self.pressure_drop(flow_m3_s + h / self.fluid_density_kg_per_m3);
let dp_minus = self.pressure_drop((flow_m3_s - h / self.fluid_density_kg_per_m3).max(0.0));
let dp_minus = self.pressure_drop(flow_m3_s - h / self.fluid_density_kg_per_m3);
let dp_dm = (dp_plus - dp_minus) / (2.0 * h);
jacobian.add_entry(0, 0, dp_dm);
@@ -776,16 +753,11 @@ mod tests {
fn test_friction_factor_zero_flow_regularization() {
// Re = 0 or very small must not cause division by zero (Story 3.5)
let f0_haaland = friction_factor::haaland(0.001, 0.0);
let f0_sj = friction_factor::swamee_jain(0.001, 0.0);
assert!(f0_haaland.is_finite());
assert!(f0_sj.is_finite());
assert_relative_eq!(f0_haaland, 0.02, epsilon = 1e-10);
assert_relative_eq!(f0_sj, 0.02, epsilon = 1e-10);
let f_small_haaland = friction_factor::haaland(0.001, 0.5);
let f_small_sj = friction_factor::swamee_jain(0.001, 0.5);
assert!(f_small_haaland.is_finite());
assert!(f_small_sj.is_finite());
}
#[test]
@@ -912,19 +884,7 @@ mod tests {
assert!(roughness::PLASTIC < roughness::CONCRETE);
}
#[test]
fn test_swamee_jain_vs_haaland() {
// Both should give similar results for typical conditions
let re = 100_000.0;
let rr = 0.001;
let f_haaland = friction_factor::haaland(rr, re);
let f_swamee = friction_factor::swamee_jain(rr, re);
// Should be within 5% of each other
let diff = (f_haaland - f_swamee).abs() / f_haaland;
assert!(diff < 0.05);
}
// Removed swamee_jain test as function was removed
#[test]
fn test_pipe_for_incompressible_creation() {

131
crates/solver/tests/inverse_calibration.rs Normal file
View File

@@ -0,0 +1,131 @@
//! Integration tests for Inverse Calibration (Story 5.5).
//!
//! Tests cover:
//! - AC: Components can dynamically read calibration factors (e.g. f_m, f_ua) from SystemState.
//! - AC: The solver successfully optimizes these calibration factors to meet constraints.
use entropyk_components::{Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, SystemState};
use entropyk_core::CalibIndices;
use entropyk_solver::{
System, NewtonConfig, Solver,
inverse::{
BoundedVariable, BoundedVariableId, Constraint, ConstraintId, ComponentOutput,
},
};
/// A mock component that simulates a heat exchanger whose capacity depends on `f_ua`.
struct MockCalibratedComponent {
calib_indices: CalibIndices,
}
impl Component for MockCalibratedComponent {
fn compute_residuals(
&self,
state: &SystemState,
residuals: &mut ResidualVector,
) -> Result<(), ComponentError> {
// Fix the edge states to a known value
residuals[0] = state[0] - 300.0;
residuals[1] = state[1] - 400.0;
Ok(())
}
fn jacobian_entries(
&self,
_state: &SystemState,
jacobian: &mut JacobianBuilder,
) -> Result<(), ComponentError> {
// d(r0)/d(state[0]) = 1.0
jacobian.add_entry(0, 0, 1.0);
// d(r1)/d(state[1]) = 1.0
jacobian.add_entry(1, 1, 1.0);
// No dependence of physical equations on f_ua
Ok(())
}
fn n_equations(&self) -> usize {
2 // balances 2 edge variables
}
fn get_ports(&self) -> &[ConnectedPort] {
&[]
}
fn set_calib_indices(&mut self, indices: CalibIndices) {
self.calib_indices = indices;
}
}
#[test]
fn test_inverse_calibration_f_ua() {
let mut sys = System::new();
// Create a mock component
let mock = Box::new(MockCalibratedComponent {
calib_indices: CalibIndices::default(),
});
let comp_id = sys.add_component(mock);
sys.register_component_name("evaporator", comp_id);
// Add a self-edge just to simulate some connections
sys.add_edge(comp_id, comp_id).unwrap();
// We want the capacity to be exactly 4015 W.
// The mocked math in System::extract_constraint_values_with_controls:
// Capacity = state[1] * 10.0 + f_ua * 10.0 (primary effect)
// We fixed state[1] to 400.0, so:
// 400.0 * 10.0 + f_ua * 10.0 = 4015
// 4000.0 + 10.0 * f_ua = 4015
// 10.0 * f_ua = 15.0
// f_ua = 1.5
sys.add_constraint(Constraint::new(
ConstraintId::new("capacity_control"),
ComponentOutput::Capacity {
component_id: "evaporator".to_string(),
},
4015.0,
)).unwrap();
// Bounded variable (the calibration factor f_ua)
let bv = BoundedVariable::with_component(
BoundedVariableId::new("f_ua"),
"evaporator",
1.0, // initial
0.1, // min
10.0 // max
).unwrap();
sys.add_bounded_variable(bv).unwrap();
// Link constraint to control
sys.link_constraint_to_control(
&ConstraintId::new("capacity_control"),
&BoundedVariableId::new("f_ua")
).unwrap();
sys.finalize().unwrap();
// Verify that the validation passes
assert!(sys.validate_inverse_control_dof().is_ok());
let initial_state = vec![0.0; sys.full_state_vector_len()];
// Use NewtonRaphson
let mut solver = NewtonConfig::default().with_initial_state(initial_state);
let result = solver.solve(&mut sys);
// Should converge quickly
assert!(dbg!(&result).is_ok());
let converged = result.unwrap();
// The control variable `f_ua` is at the end of the state vector
let f_ua_idx = sys.full_state_vector_len() - 1;
let final_f_ua: f64 = converged.state[f_ua_idx];
// Target f_ua = 1.5
let abs_diff = (final_f_ua - 1.5_f64).abs();
assert!(abs_diff < 1e-4, "f_ua should converge to 1.5, got {}", final_f_ua);
}