diff --git a/apps/web/src/lib/componentMeta.test.ts b/apps/web/src/lib/componentMeta.test.ts index cb76c77..774791b 100644 --- a/apps/web/src/lib/componentMeta.test.ts +++ b/apps/web/src/lib/componentMeta.test.ts @@ -206,8 +206,14 @@ describe("COMPONENTS catalogue integrity", () => { "channel_spacing_mm", "chevron_angle_deg", "correlation", + "dp_correlation", ]), ); + const dp = COMPONENT_BY_TYPE[type].params.find((p) => p.key === "dp_correlation"); + expect(dp?.options?.map((o) => o.value)).toEqual([ + "SimplifiedChannel", + "Martin1996", + ]); } }); diff --git a/apps/web/src/lib/componentMeta.ts b/apps/web/src/lib/componentMeta.ts index 3dc0e98..72b6688 100644 --- a/apps/web/src/lib/componentMeta.ts +++ b/apps/web/src/lib/componentMeta.ts @@ -1076,6 +1076,19 @@ export const COMPONENTS: ComponentMeta[] = [ { value: "Shah2021", label: "Shah 2021" }, ], }, + { + key: "dp_correlation", + label: "Pressure-drop correlation", + kind: "string", + default: "SimplifiedChannel", + section: "Pressure drop", + description: + "Channel friction on both hot and cold sides (independent of HTC). Default SimplifiedChannel; Martin1996 uses chevron angle.", + options: [ + { value: "SimplifiedChannel", label: "Simplified channel (default)" }, + { value: "Martin1996", label: "Martin 1996 (chevron)" }, + ], + }, { key: "ua", label: "UA override", kind: "number", unit: "W/K", section: "Calibration", min: 0.0, advanced: true }, { key: "calib_sst_c", @@ -1163,6 +1176,19 @@ export const COMPONENTS: ComponentMeta[] = [ { value: "Shah2021", label: "Shah 2021" }, ], }, + { + key: "dp_correlation", + label: "Pressure-drop correlation", + kind: "string", + default: "SimplifiedChannel", + section: "Pressure drop", + description: + "Channel friction on both hot and cold sides (independent of HTC). Default SimplifiedChannel; Martin1996 uses chevron angle.", + options: [ + { value: "SimplifiedChannel", label: "Simplified channel (default)" }, + { value: "Martin1996", label: "Martin 1996 (chevron)" }, + ], + }, { key: "ua", label: "UA override", kind: "number", unit: "W/K", section: "Calibration", min: 0.0, advanced: true }, { key: "calib_sdt_c", diff --git a/apps/web/src/lib/configBuilder.ts b/apps/web/src/lib/configBuilder.ts index d2f31df..32b18ac 100644 --- a/apps/web/src/lib/configBuilder.ts +++ b/apps/web/src/lib/configBuilder.ts @@ -141,6 +141,7 @@ const PARAM_ALIASES: Record = { rh: ["relative_humidity"], p_set_bar: ["pressure_bar"], p_back_bar: ["pressure_bar"], + dp_correlation: ["DpCorrelation", "dpCorrelation"], }; function getParam( diff --git a/crates/cli/examples/bphx_evaporator_condenser.json b/crates/cli/examples/bphx_evaporator_condenser.json index a907bca..f428145 100644 --- a/crates/cli/examples/bphx_evaporator_condenser.json +++ b/crates/cli/examples/bphx_evaporator_condenser.json @@ -33,6 +33,7 @@ "target_subcooling_k": 5.0, "emergent_pressure": true, "correlation": "Longo2004", + "dp_correlation": "SimplifiedChannel", "ua": 2500.0 }, { @@ -52,6 +53,7 @@ "target_superheat_k": 5.0, "emergent_pressure": true, "correlation": "Longo2004", + "dp_correlation": "SimplifiedChannel", "ua": 2000.0 }, { diff --git a/crates/cli/src/batch.rs b/crates/cli/src/batch.rs index 4bba9e5..2a201bb 100644 --- a/crates/cli/src/batch.rs +++ b/crates/cli/src/batch.rs @@ -315,7 +315,7 @@ fn process_single_file( initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), } } } @@ -425,7 +425,7 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, SimulationResult { input: "test2.json".to_string(), @@ -439,7 +439,7 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, ]; @@ -483,7 +483,7 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, SimulationResult { input: "test2.json".to_string(), @@ -497,7 +497,7 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, ]; @@ -529,7 +529,7 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, SimulationResult { input: "test2.json".to_string(), @@ -543,7 +543,7 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, ]; @@ -571,7 +571,7 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, SimulationResult { input: "test2.json".to_string(), @@ -585,7 +585,7 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 1000, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, ]; @@ -643,7 +643,7 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 100, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }], }; @@ -667,7 +667,7 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, SimulationResult { input: "ok2.json".to_string(), @@ -681,7 +681,7 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 60, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, SimulationResult { input: "fail.json".to_string(), @@ -695,7 +695,7 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, ]; diff --git a/crates/cli/src/run.rs b/crates/cli/src/run.rs index 34d0f66..cb6cddc 100644 --- a/crates/cli/src/run.rs +++ b/crates/cli/src/run.rs @@ -76,6 +76,13 @@ pub struct SimulationResult { /// solve via the `ENTROPYK_INIT_STATE_JSON` environment variable). #[serde(default, skip_serializing_if = "Option::is_none")] pub raw_state_vector: Option>, + /// Named solver-computed unknowns (free actuators + calibration factors + /// like `z_ua`, `opening`, `z_dp`) read out of `raw_state_vector` and + /// labeled with their stable id, owning component, and `[min, max]` + /// bounds. Empty for systems with no bounded control variables; skipped + /// on serialization when empty so older consumers keep working. + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub solved_variables: Vec, } /// Compact DoF summary for CLI / web API consumers. @@ -365,7 +372,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } @@ -400,7 +407,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } }; @@ -710,7 +717,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } }, @@ -730,7 +737,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } } @@ -775,7 +782,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } }, @@ -795,7 +802,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } None => 1, @@ -816,7 +823,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } }, @@ -836,7 +843,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } None => 0, @@ -860,7 +867,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } } @@ -880,7 +887,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } } @@ -950,7 +957,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } } @@ -976,7 +983,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; let find_component = |name: &str| { config @@ -1054,7 +1061,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } system.add_saturated_controller(controller); @@ -1072,7 +1079,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } } @@ -1246,7 +1253,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } system.add_free_actuator(var_id); @@ -1268,7 +1275,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } } @@ -1289,7 +1296,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } @@ -1320,7 +1327,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: Some(summary), elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } summary @@ -1517,7 +1524,7 @@ fn execute_simulation( initialization_diagnostics: None, dof: Some(dof_summary), elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; } @@ -1746,6 +1753,7 @@ fn execute_simulation( dof: Some(dof_summary.clone()), elapsed_ms, raw_state_vector: Some(converged.state.clone()), + solved_variables: entropyk::extract_solved_variables(&system, &converged.state), } } Err(e) => { @@ -1791,7 +1799,7 @@ fn execute_simulation( initialization_diagnostics: Some(initialization_diagnostics), dof: Some(dof_summary), elapsed_ms, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), } } } @@ -2939,6 +2947,32 @@ fn bphx_correlation_from_params( Ok(Some(correlation)) } +/// Parse BPHX frictional `dp_correlation` / `DpCorrelation` (independent of HTC). +/// +/// Default when absent: `SimplifiedChannel`. Unknown ids are rejected (no silent isobaric). +fn bphx_dp_correlation_from_params( + params: &std::collections::HashMap, + component_name: &str, +) -> CliResult { + use entropyk_components::heat_exchanger::BphxDpCorrelation; + + let Some(value) = params + .get("dp_correlation") + .or_else(|| params.get("DpCorrelation")) + .or_else(|| params.get("dpCorrelation")) + else { + return Ok(BphxDpCorrelation::default()); + }; + let Some(raw) = value.as_str() else { + return Err(CliError::Config(format!( + "{component_name}: dp_correlation must be a string (got {value})" + ))); + }; + BphxDpCorrelation::parse(raw).map_err(|msg| { + CliError::Config(format!("{component_name}: {msg}")) + }) +} + /// Extract calibration Z-factors for BphxEvaporator/BphxCondenser from JSON params. /// /// Errors if `ua_nominal == 0` and an explicit `ua` override is provided (geometry is @@ -5095,7 +5129,11 @@ fn create_component( .with_target_superheat(target_superheat) .with_refrigerant(refrigerant) .with_secondary_fluid(secondary_fluid) - .with_fluid_backend(Arc::clone(&backend)); + .with_fluid_backend(Arc::clone(&backend)) + .with_dp_correlation(bphx_dp_correlation_from_params( + params, + "BphxEvaporator", + )?); if let Some(correlation) = bphx_correlation_from_params(params, "BphxEvaporator")? { @@ -5151,7 +5189,11 @@ fn create_component( .with_refrigerant(refrigerant) .with_secondary_fluid(secondary_fluid) .with_fluid_backend(Arc::clone(&backend)) - .with_target_subcooling(target_subcooling); + .with_target_subcooling(target_subcooling) + .with_dp_correlation(bphx_dp_correlation_from_params( + params, + "BphxCondenser", + )?); if let Some(correlation) = bphx_correlation_from_params(params, "BphxCondenser")? { @@ -5630,7 +5672,7 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; let json = serde_json::to_string_pretty(&result).unwrap(); @@ -5787,6 +5829,56 @@ mod tests { ); } + #[test] + fn test_bphx_dp_correlation_defaults_and_aliases() { + use entropyk_components::heat_exchanger::BphxDpCorrelation; + use serde_json::json; + + assert_eq!( + bphx_dp_correlation_from_params(&std::collections::HashMap::new(), "BphxEvaporator") + .unwrap(), + BphxDpCorrelation::SimplifiedChannel + ); + let params = std::collections::HashMap::from([( + "DpCorrelation".to_string(), + json!("Martin1996"), + )]); + assert_eq!( + bphx_dp_correlation_from_params(¶ms, "BphxCondenser").unwrap(), + BphxDpCorrelation::Martin1996 + ); + } + + #[test] + fn test_bphx_dp_correlation_rejects_unknown() { + use serde_json::json; + + let params = std::collections::HashMap::from([( + "dp_correlation".to_string(), + json!("Amalfi2016"), + )]); + let error = bphx_dp_correlation_from_params(¶ms, "BphxEvaporator").unwrap_err(); + assert!( + error.to_string().contains("unsupported dp_correlation"), + "{error}" + ); + } + + #[test] + fn test_bphx_dp_correlation_rejects_non_string() { + use serde_json::json; + + let params = std::collections::HashMap::from([( + "dp_correlation".to_string(), + json!(42), + )]); + let error = bphx_dp_correlation_from_params(¶ms, "BphxEvaporator").unwrap_err(); + assert!( + error.to_string().contains("must be a string"), + "{error}" + ); + } + fn comp_from_json(v: serde_json::Value) -> crate::config::ComponentConfig { serde_json::from_value(v).expect("valid component JSON") } diff --git a/crates/cli/tests/batch_execution.rs b/crates/cli/tests/batch_execution.rs index 0f2e882..d03fd5c 100644 --- a/crates/cli/tests/batch_execution.rs +++ b/crates/cli/tests/batch_execution.rs @@ -92,7 +92,7 @@ fn test_simulation_result_statuses() { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, SimulationResult { input: "fail.json".to_string(), @@ -106,7 +106,7 @@ fn test_simulation_result_statuses() { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, SimulationResult { input: "timeout.json".to_string(), @@ -120,7 +120,7 @@ fn test_simulation_result_statuses() { initialization_diagnostics: None, dof: None, elapsed_ms: 1000, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, ]; @@ -163,7 +163,7 @@ fn test_batch_aggregator_csv_output() { initialization_diagnostics: None, dof: None, elapsed_ms: 150, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, SimulationResult { input: "scenario2.json".to_string(), @@ -183,7 +183,7 @@ fn test_batch_aggregator_csv_output() { initialization_diagnostics: None, dof: None, elapsed_ms: 200, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, SimulationResult { input: "scenario3.json".to_string(), @@ -197,7 +197,7 @@ fn test_batch_aggregator_csv_output() { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, ]; @@ -229,7 +229,7 @@ fn test_batch_aggregator_json_summary() { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, SimulationResult { input: "test2.json".to_string(), @@ -243,7 +243,7 @@ fn test_batch_aggregator_json_summary() { initialization_diagnostics: None, dof: None, elapsed_ms: 75, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, SimulationResult { input: "test3.json".to_string(), @@ -257,7 +257,7 @@ fn test_batch_aggregator_json_summary() { initialization_diagnostics: None, dof: None, elapsed_ms: 5000, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }, ]; @@ -322,7 +322,7 @@ fn test_batch_summary_csv_with_convergence() { initialization_diagnostics: None, dof: None, elapsed_ms: 300, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }]; let summary = BatchSummary { diff --git a/crates/cli/tests/single_run.rs b/crates/cli/tests/single_run.rs index 8ba0445..4dce070 100644 --- a/crates/cli/tests/single_run.rs +++ b/crates/cli/tests/single_run.rs @@ -29,7 +29,7 @@ fn test_simulation_result_serialization() { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; let json = serde_json::to_string_pretty(&result).unwrap(); @@ -70,7 +70,7 @@ fn test_error_result_serialization() { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; let json = serde_json::to_string(&result).unwrap(); @@ -96,7 +96,7 @@ fn test_error_result_serializes_failure_diagnostics() { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; let json = serde_json::to_string(&result).unwrap(); @@ -2303,7 +2303,7 @@ fn test_structural_failure_serializes_without_diagnostics() { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; let json = serde_json::to_string(&result).unwrap(); @@ -2348,7 +2348,7 @@ fn test_success_result_does_not_include_failure_diagnostics() { initialization_diagnostics: None, dof: None, elapsed_ms: 120, - raw_state_vector: None, + raw_state_vector: None, solved_variables: Vec::new(), }; let json = serde_json::to_string(&result).unwrap(); diff --git a/crates/components/src/heat_exchanger/bphx_condenser.rs b/crates/components/src/heat_exchanger/bphx_condenser.rs index 23b411c..03b6cc7 100644 --- a/crates/components/src/heat_exchanger/bphx_condenser.rs +++ b/crates/components/src/heat_exchanger/bphx_condenser.rs @@ -24,6 +24,7 @@ //! ``` use super::bphx_correlation::{BphxCorrelation, CorrelationEvaluation}; +use super::bphx_dp::BphxDpCorrelation; use super::bphx_exchanger::BphxExchanger; use super::bphx_geometry::{BphxGeometry, BphxType}; use super::correlation_registry::CorrelationSelectionError; @@ -150,6 +151,17 @@ impl BphxCondenser { self } + /// Sets the channel frictional pressure-drop correlation (independent of HTC). + pub fn with_dp_correlation(mut self, correlation: BphxDpCorrelation) -> Self { + self.inner = self.inner.with_dp_correlation(correlation); + self + } + + /// Returns the configured pressure-drop correlation. + pub fn dp_correlation(&self) -> BphxDpCorrelation { + self.inner.dp_correlation() + } + /// Sets the target subcooling in Kelvin. /// /// # Panics diff --git a/crates/components/src/heat_exchanger/bphx_correlation.rs b/crates/components/src/heat_exchanger/bphx_correlation.rs index d4014c9..994115d 100644 --- a/crates/components/src/heat_exchanger/bphx_correlation.rs +++ b/crates/components/src/heat_exchanger/bphx_correlation.rs @@ -344,7 +344,11 @@ impl BphxCorrelation { CorrelationId::Cooper1984 | CorrelationId::Mostinski1963 | CorrelationId::Friedel1979 - | CorrelationId::MullerSteinhagenHeck1986 => None, + | CorrelationId::MullerSteinhagenHeck1986 + | CorrelationId::SimplifiedChannel + | CorrelationId::Martin1996 + | CorrelationId::AmalfiThome2016 + | CorrelationId::TaoInfanteFerreiraSurvey => None, } } diff --git a/crates/components/src/heat_exchanger/bphx_dp.rs b/crates/components/src/heat_exchanger/bphx_dp.rs new file mode 100644 index 0000000..93ccf6f --- /dev/null +++ b/crates/components/src/heat_exchanger/bphx_dp.rs @@ -0,0 +1,422 @@ +//! BPHX channel pressure-drop correlations (friction), separate from HTC. +//! +//! Pressure-drop form (Fanning-like `f`, matching existing SimplifiedChannel): +//! `ΔP = z_dp · 2 · f · L · G² / (ρ · d_h)`. +//! +//! Martin (1996/1999) returns a Fanning friction factor from the reciprocal- +//! square-root chevron model (same intermediate `f` as `fluids` / +//! `ht.friction_plate_Martin_1999` before the ×4 Darcy conversion). + +use super::bphx_geometry::BphxGeometry; +use crate::ComponentError; +use std::f64::consts::PI; + +/// Reference dynamic viscosity [Pa·s] used for Re in channel DP models. +/// Matches the historical SimplifiedChannel implementation (Newton-friendly constant). +pub const BPHX_DP_MU_REF_PA_S: f64 = 0.0002; + +/// Selectable BPHX frictional pressure-drop correlation. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] +pub enum BphxDpCorrelation { + /// Existing Blasius/laminar C¹ blend (default). + #[default] + SimplifiedChannel, + /// Holger Martin chevron-plate single-phase friction (1996 / 1999 form). + Martin1996, +} + +impl BphxDpCorrelation { + /// Stable serialization name. + pub const fn as_str(self) -> &'static str { + match self { + Self::SimplifiedChannel => "SimplifiedChannel", + Self::Martin1996 => "Martin1996", + } + } + + /// Parse CLI/UI identifiers. Rejects unknown strings (no silent fallback). + pub fn parse(name: &str) -> Result { + let key = name.trim().to_ascii_lowercase().replace(['-', '_', ' '], ""); + match key.as_str() { + "simplifiedchannel" | "simplified" | "channel" | "default" => { + Ok(Self::SimplifiedChannel) + } + "martin1996" | "martin1999" | "martin" => Ok(Self::Martin1996), + _ => Err(format!( + "unsupported dp_correlation '{name}'. Use 'SimplifiedChannel' or 'Martin1996'." + )), + } + } +} + +/// Evaluated channel pressure drop and ∂ΔP/∂G at fixed ρ and z_dp. +#[derive(Debug, Clone, Copy)] +pub struct BphxDpEval { + /// Pressure drop [Pa] (always ≥ 0 for the channel friction model). + pub delta_p_pa: f64, + /// ∂ΔP/∂G [Pa / (kg/(m²·s))]. + pub d_delta_p_d_g: f64, +} + +/// Fanning friction factor and ∂f/∂Re for the selected correlation. +fn fanning_friction( + corr: BphxDpCorrelation, + re: f64, + chevron_angle_deg: f64, +) -> Result<(f64, f64), ComponentError> { + match corr { + BphxDpCorrelation::SimplifiedChannel => Ok(simplified_channel_fanning(re)), + BphxDpCorrelation::Martin1996 => { + validate_martin_chevron(chevron_angle_deg)?; + Ok(martin1996_fanning(re, chevron_angle_deg)) + } + } +} + +fn validate_martin_chevron(angle_deg: f64) -> Result<(), ComponentError> { + if !angle_deg.is_finite() || !(0.0..90.0).contains(&angle_deg) { + return Err(ComponentError::InvalidState(format!( + "Martin1996 requires chevron_angle in (0, 90) degrees (got {angle_deg})" + ))); + } + Ok(()) +} + +/// SimplifiedChannel: laminar 64/Re, Blasius 0.079 Re^−0.25, C¹ blend on [2300, 4000]. +/// +/// Note: laminar uses the Darcy 64/Re coefficient inside the Fanning ΔP form — +/// historical Entropyk convention; kept identical for continuity. +fn simplified_channel_fanning(re: f64) -> (f64, f64) { + const RE_LAMINAR: f64 = 2300.0; + const RE_TURBULENT: f64 = 4000.0; + let re_safe = re.max(1e-9); + let f_lam = 64.0 / re_safe; + let df_lam = -64.0 / (re_safe * re_safe); + + if re <= RE_LAMINAR { + return (f_lam, df_lam); + } + + let f_turb = 0.079 * re_safe.powf(-0.25); + let df_turb = 0.079 * (-0.25) * re_safe.powf(-1.25); + + if re >= RE_TURBULENT { + return (f_turb, df_turb); + } + + let f = entropyk_core::smoothing::cubic_blend( + f_lam, + f_turb, + re, + RE_LAMINAR, + RE_TURBULENT, + ); + // d/dRe of blend(f_lam(Re), f_turb(Re), Re): product rule on cubic Hermite. + let df_blend_dx = entropyk_core::smoothing::cubic_blend_derivative( + f_lam, + f_turb, + re, + RE_LAMINAR, + RE_TURBULENT, + ); + // cubic_blend treats endpoints as constants w.r.t x; add chain-rule terms for + // f_lam(Re) and f_turb(Re) via the complementary weights. + let t = ((re - RE_LAMINAR) / (RE_TURBULENT - RE_LAMINAR)).clamp(0.0, 1.0); + let w_b = t * t * (3.0 - 2.0 * t); // weight on b (= f_turb) + let w_a = 1.0 - w_b; + let df = df_blend_dx + w_a * df_lam + w_b * df_turb; + (f, df) +} + +/// Martin 1996/1999 Fanning friction factor vs Re and chevron angle φ [deg]. +/// +/// ```text +/// 1/√f = cosφ / √(0.045 tanφ + 0.09 sinφ + f0/cosφ) + (1−cosφ) / √(3.8 f1) +/// ``` +/// with piecewise `f0`, `f1` (Shah/Martin 1999 appendix form). +/// Laminar/turbulent `f0`/`f1` branches are C¹-blended on Re ∈ [1800, 2200] +/// (same Newton rationale as SimplifiedChannel’s [2300, 4000] blend). +fn martin1996_fanning(re: f64, chevron_angle_deg: f64) -> (f64, f64) { + const RE_LO: f64 = 1800.0; + const RE_HI: f64 = 2200.0; + let re_safe = re.max(1e-9); + + let lam = |re_v: f64| -> (f64, f64, f64, f64) { + let f0 = 16.0 / re_v; + let df0 = -16.0 / (re_v * re_v); + let f1 = 149.0 / re_v + 0.9625; + let df1 = -149.0 / (re_v * re_v); + (f0, df0, f1, df1) + }; + let turb = |re_v: f64| -> (f64, f64, f64, f64) { + let u = 1.56 * re_v.ln() - 3.0; + let u_safe = if u.abs() < 1e-12 { + if u >= 0.0 { + 1e-12 + } else { + -1e-12 + } + } else { + u + }; + let f0 = u_safe.powi(-2); + let df0 = -2.0 * u_safe.powi(-3) * (1.56 / re_v); + let f1 = 9.75 * re_v.powf(-0.289); + let df1 = 9.75 * (-0.289) * re_v.powf(-1.289); + (f0, df0, f1, df1) + }; + + let (f0, df0, f1, df1) = if re_safe <= RE_LO { + lam(re_safe) + } else if re_safe >= RE_HI { + turb(re_safe) + } else { + let (f0_l, df0_l, f1_l, df1_l) = lam(re_safe); + let (f0_t, df0_t, f1_t, df1_t) = turb(re_safe); + let f0 = entropyk_core::smoothing::cubic_blend(f0_l, f0_t, re_safe, RE_LO, RE_HI); + let f1 = entropyk_core::smoothing::cubic_blend(f1_l, f1_t, re_safe, RE_LO, RE_HI); + let t = ((re_safe - RE_LO) / (RE_HI - RE_LO)).clamp(0.0, 1.0); + let w_b = t * t * (3.0 - 2.0 * t); + let w_a = 1.0 - w_b; + let df0_blend = entropyk_core::smoothing::cubic_blend_derivative( + f0_l, f0_t, re_safe, RE_LO, RE_HI, + ); + let df1_blend = entropyk_core::smoothing::cubic_blend_derivative( + f1_l, f1_t, re_safe, RE_LO, RE_HI, + ); + let df0 = df0_blend + w_a * df0_l + w_b * df0_t; + let df1 = df1_blend + w_a * df1_l + w_b * df1_t; + (f0, df0, f1, df1) + }; + + let phi = chevron_angle_deg * PI / 180.0; + let cos_p = phi.cos(); + let sin_p = phi.sin(); + let tan_p = phi.tan(); + // Guard grazing angles so /cosφ stays finite inside (0, 90). + let cos_safe = cos_p.signum() * cos_p.abs().max(1e-12); + + let a = 0.045 * tan_p + 0.09 * sin_p + f0 / cos_safe; + let b = 3.8 * f1.max(1e-30); + let a_safe = a.max(1e-30); + + let term0 = cos_p / a_safe.sqrt(); + let term1 = (1.0 - cos_p) / b.sqrt(); + let rhs = term0 + term1; + let rhs_safe = rhs.max(1e-30); + let f = rhs_safe.powi(-2); + + let da_dre = df0 / cos_safe; + let db_dre = 3.8 * df1; + let dterm0 = cos_p * (-0.5) * a_safe.powf(-1.5) * da_dre; + let dterm1 = (1.0 - cos_p) * (-0.5) * b.powf(-1.5) * db_dre; + let drhs = dterm0 + dterm1; + let df = -2.0 * rhs_safe.powi(-3) * drhs; + (f, df) +} + +/// Channel frictional ΔP and ∂ΔP/∂G for the given correlation. +/// +/// `mass_flux` may be signed; Re uses `|G|`, ΔP uses `G²` (existing convention). +pub fn evaluate_channel_pressure_drop( + corr: BphxDpCorrelation, + geometry: &BphxGeometry, + mass_flux: f64, + rho: f64, + z_dp: f64, +) -> Result { + if !mass_flux.is_finite() { + return Ok(BphxDpEval { + delta_p_pa: 0.0, + d_delta_p_d_g: 0.0, + }); + } + if !rho.is_finite() || rho < 1e-10 || geometry.dh < 1e-10 { + return Ok(BphxDpEval { + delta_p_pa: 0.0, + d_delta_p_d_g: 0.0, + }); + } + if !geometry.plate_length.is_finite() || geometry.plate_length <= 0.0 { + return Ok(BphxDpEval { + delta_p_pa: 0.0, + d_delta_p_d_g: 0.0, + }); + } + if !z_dp.is_finite() || z_dp < 0.0 { + return Err(ComponentError::InvalidState(format!( + "BPHX z_dp must be finite and >= 0 (got {z_dp})" + ))); + } + + let g_abs = mass_flux.abs(); + let re = g_abs * geometry.dh / BPHX_DP_MU_REF_PA_S; + let (f, df_dre) = fanning_friction(corr, re, geometry.chevron_angle)?; + + // ΔP = z_dp · 2 f L G² / (ρ dh) + let k = z_dp * 2.0 * geometry.plate_length / (rho * geometry.dh); + let g2 = mass_flux * mass_flux; + let delta_p_pa = k * f * g2; + + // dRe/dG = sign(G) · dh/μ (0 at stagnation) + let d_re_d_g = if g_abs < 1e-30 { + 0.0 + } else { + mass_flux.signum() * geometry.dh / BPHX_DP_MU_REF_PA_S + }; + // d(G²)/dG = 2 G + let d_delta_p_d_g = k * (df_dre * d_re_d_g * g2 + f * 2.0 * mass_flux); + + Ok(BphxDpEval { + delta_p_pa, + d_delta_p_d_g, + }) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::heat_exchanger::BphxGeometry; + + fn geo() -> BphxGeometry { + BphxGeometry::from_dh_area(0.003, 0.5, 20) + } + + #[test] + fn parse_accepts_aliases_rejects_unknown() { + assert_eq!( + BphxDpCorrelation::parse("SimplifiedChannel").unwrap(), + BphxDpCorrelation::SimplifiedChannel + ); + assert_eq!( + BphxDpCorrelation::parse("martin-1996").unwrap(), + BphxDpCorrelation::Martin1996 + ); + assert!(BphxDpCorrelation::parse("Amalfi2016").is_err()); + assert!(BphxDpCorrelation::parse("friedel").is_err()); + } + + #[test] + fn simplified_matches_historical_scale() { + let g = 30.0; + let rho = 1100.0; + let eval = evaluate_channel_pressure_drop( + BphxDpCorrelation::SimplifiedChannel, + &geo(), + g, + rho, + 1.0, + ) + .unwrap(); + assert!(eval.delta_p_pa > 0.0); + assert!(eval.delta_p_pa.is_finite()); + } + + #[test] + fn martin_differs_from_simplified_same_state() { + let g = 80.0; + let rho = 1000.0; + let simp = evaluate_channel_pressure_drop( + BphxDpCorrelation::SimplifiedChannel, + &geo(), + g, + rho, + 1.0, + ) + .unwrap(); + let martin = evaluate_channel_pressure_drop( + BphxDpCorrelation::Martin1996, + &geo(), + g, + rho, + 1.0, + ) + .unwrap(); + assert!( + (simp.delta_p_pa - martin.delta_p_pa).abs() > 1.0, + "Martin and Simplified should differ: {} vs {}", + simp.delta_p_pa, + martin.delta_p_pa + ); + } + + #[test] + fn z_dp_scales_linearly() { + let g = 40.0; + let rho = 1100.0; + let a = evaluate_channel_pressure_drop( + BphxDpCorrelation::SimplifiedChannel, + &geo(), + g, + rho, + 1.0, + ) + .unwrap(); + let b = evaluate_channel_pressure_drop( + BphxDpCorrelation::SimplifiedChannel, + &geo(), + g, + rho, + 0.5, + ) + .unwrap(); + assert!((b.delta_p_pa - 0.5 * a.delta_p_pa).abs() < 1e-9 * a.delta_p_pa.max(1.0)); + } + + #[test] + fn zero_flow_smooth() { + let eval = evaluate_channel_pressure_drop( + BphxDpCorrelation::SimplifiedChannel, + &geo(), + 0.0, + 1100.0, + 1.0, + ) + .unwrap(); + assert_eq!(eval.delta_p_pa, 0.0); + assert!(eval.d_delta_p_d_g.is_finite()); + } + + #[test] + fn martin_rejects_invalid_chevron() { + let mut g = geo(); + g.chevron_angle = f64::NAN; + let err = evaluate_channel_pressure_drop( + BphxDpCorrelation::Martin1996, + &g, + 30.0, + 1000.0, + 1.0, + ) + .unwrap_err(); + assert!(format!("{err}").contains("chevron")); + } + + #[test] + fn analytic_d_dg_matches_central_difference() { + let geo = geo(); + let g0 = 55.0; + let rho = 1050.0; + for corr in [ + BphxDpCorrelation::SimplifiedChannel, + BphxDpCorrelation::Martin1996, + ] { + let eval = evaluate_channel_pressure_drop(corr, &geo, g0, rho, 1.0).unwrap(); + let h = 1e-4; + let plus = + evaluate_channel_pressure_drop(corr, &geo, g0 + h, rho, 1.0).unwrap(); + let minus = + evaluate_channel_pressure_drop(corr, &geo, g0 - h, rho, 1.0).unwrap(); + let fd = (plus.delta_p_pa - minus.delta_p_pa) / (2.0 * h); + let rel = (eval.d_delta_p_d_g - fd).abs() / fd.abs().max(1.0); + assert!( + rel < 1e-4, + "{:?}: analytic {} vs fd {} (rel {})", + corr, + eval.d_delta_p_d_g, + fd, + rel + ); + } + } +} diff --git a/crates/components/src/heat_exchanger/bphx_evaporator.rs b/crates/components/src/heat_exchanger/bphx_evaporator.rs index 51ec976..dd3de2e 100644 --- a/crates/components/src/heat_exchanger/bphx_evaporator.rs +++ b/crates/components/src/heat_exchanger/bphx_evaporator.rs @@ -24,6 +24,7 @@ //! ``` use super::bphx_correlation::{BphxCorrelation, CorrelationEvaluation}; +use super::bphx_dp::BphxDpCorrelation; use super::bphx_exchanger::BphxExchanger; use super::bphx_geometry::{BphxGeometry, BphxType}; use super::correlation_registry::CorrelationSelectionError; @@ -154,6 +155,17 @@ impl BphxEvaporator { self } + /// Sets the channel frictional pressure-drop correlation (independent of HTC). + pub fn with_dp_correlation(mut self, correlation: BphxDpCorrelation) -> Self { + self.inner = self.inner.with_dp_correlation(correlation); + self + } + + /// Returns the configured pressure-drop correlation. + pub fn dp_correlation(&self) -> BphxDpCorrelation { + self.inner.dp_correlation() + } + /// Returns the component name. pub fn name(&self) -> &str { self.inner.name() diff --git a/crates/components/src/heat_exchanger/bphx_exchanger.rs b/crates/components/src/heat_exchanger/bphx_exchanger.rs index dcc4318..0effc09 100644 --- a/crates/components/src/heat_exchanger/bphx_exchanger.rs +++ b/crates/components/src/heat_exchanger/bphx_exchanger.rs @@ -33,6 +33,7 @@ use super::bphx_correlation::{ BphxCorrelation, CorrelationEvaluation, CorrelationParams, CorrelationResult, CorrelationSelector, ValidityStatus, }; +use super::bphx_dp::{evaluate_channel_pressure_drop, BphxDpCorrelation, BphxDpEval}; use super::bphx_geometry::{BphxGeometry, BphxType}; use super::correlation_registry::{ CorrelationSelectionError, ExchangerGeometryType, FlowRegime, SelectionOutcome, @@ -55,6 +56,7 @@ pub struct BphxExchanger { inner: HeatExchanger, geometry: BphxGeometry, correlation_selector: CorrelationSelector, + dp_correlation: BphxDpCorrelation, refrigerant_id: String, secondary_fluid_id: String, fluid_backend: Option>, @@ -70,6 +72,7 @@ impl std::fmt::Debug for BphxExchanger { .field("ua", &self.ua()) .field("geometry", &self.geometry) .field("correlation", &self.correlation_selector.correlation) + .field("dp_correlation", &self.dp_correlation) .field("refrigerant_id", &self.refrigerant_id) .field("secondary_fluid_id", &self.secondary_fluid_id) .field("has_fluid_backend", &self.fluid_backend.is_some()) @@ -108,6 +111,7 @@ impl BphxExchanger { inner: HeatExchanger::new(model, "BphxExchanger"), geometry, correlation_selector: CorrelationSelector::default(), + dp_correlation: BphxDpCorrelation::default(), refrigerant_id: String::new(), secondary_fluid_id: String::new(), fluid_backend: None, @@ -125,6 +129,7 @@ impl BphxExchanger { inner: HeatExchanger::new(model, "BphxExchanger"), geometry, correlation_selector: CorrelationSelector::default(), + dp_correlation: BphxDpCorrelation::default(), refrigerant_id: String::new(), secondary_fluid_id: String::new(), fluid_backend: None, @@ -151,6 +156,17 @@ impl BphxExchanger { self } + /// Sets the channel frictional pressure-drop correlation (independent of HTC). + pub fn with_dp_correlation(mut self, correlation: BphxDpCorrelation) -> Self { + self.dp_correlation = correlation; + self + } + + /// Returns the configured pressure-drop correlation. + pub fn dp_correlation(&self) -> BphxDpCorrelation { + self.dp_correlation + } + /// Sets the refrigerant fluid identifier. pub fn with_refrigerant(mut self, fluid: impl Into) -> Self { self.set_refrigerant_id(fluid); @@ -338,12 +354,10 @@ impl BphxExchanger { self.last_validity_warning.set(false); } - /// Computes the pressure drop using a simplified correlation. + /// Computes the channel frictional pressure drop for the configured correlation. /// - /// ΔP = f_dp × (2 × f × L × G²) / (ρ × d_h) - /// - /// where f is the friction factor, L is the plate length, G is mass flux. - /// The result is scaled by `calib().z_dp`. + /// ΔP = z_dp × (2 × f × L × G²) / (ρ × d_h) with Fanning-like `f` + /// (SimplifiedChannel default; Martin1996 when selected). /// /// # Arguments /// @@ -352,46 +366,158 @@ impl BphxExchanger { /// /// # Returns /// - /// Pressure drop in Pa, scaled by f_dp calibration factor. + /// Pressure drop in Pa, scaled by `calib().z_dp`. Returns 0 on invalid inputs + /// or (for Martin) invalid chevron — callers that need errors should use + /// [`evaluate_channel_pressure_drop`]. pub fn compute_pressure_drop(&self, mass_flux: f64, rho: f64) -> f64 { - if rho < 1e-10 || self.geometry.dh < 1e-10 { - return 0.0; + evaluate_channel_pressure_drop( + self.dp_correlation, + &self.geometry, + mass_flux, + rho, + self.calib().z_dp, + ) + .map(|e| e.delta_p_pa) + .unwrap_or(0.0) + } + + /// Live `z_dp`: free actuator from state when indexed, else calib value. + fn live_z_dp(&self, state: &StateSlice) -> f64 { + self.inner + .calib_indices_ref() + .z_dp + .and_then(|idx| state.get(idx).copied()) + .unwrap_or_else(|| self.calib().z_dp) + } + + fn side_channel_dp( + &self, + mass_flow: f64, + rho: f64, + z_dp: f64, + ) -> Result { + let g = self.geometry.mass_flux(mass_flow); + evaluate_channel_pressure_drop(self.dp_correlation, &self.geometry, g, rho, z_dp) + } + + /// Overwrites the two 4-port pressure rows with `P_out − P_in + ΔP = 0`. + fn overwrite_pressure_closures_with_dp( + &self, + state: &StateSlice, + residuals: &mut ResidualVector, + ) -> Result<(), ComponentError> { + let Some(edges) = self.inner.four_port_edges() else { + return Ok(()); + }; + let row0 = self.inner.model_n_equations(); + if residuals.len() < row0 + 2 { + return Err(ComponentError::InvalidResidualDimensions { + expected: row0 + 2, + actual: residuals.len(), + }); } - let re = mass_flux.abs() * self.geometry.dh / 0.0002; - // C¹ laminar→turbulent blend over [2300, 4000] (same rationale as the - // pipe friction factor and the Gnielinski correlation): a hard switch at - // Re = 2300 put a kink in the pressure-drop residual that the analytic - // Newton Jacobian could not see, hurting convergence near the transition. - // Reynolds uses |mass_flux| so transient reverse flow during iteration - // yields a physical friction factor instead of the flat `64` the old - // signed/`max(1.0)` form produced. The lower clamp only guards the exact - // div-by-zero at stagnation; at 1e-9 it is far below any reachable Re, so - // it introduces no visible kink (the previous `max(1.0)` kinked at Re = 1). - const RE_LAMINAR: f64 = 2300.0; - const RE_TURBULENT: f64 = 4000.0; - let f_laminar = 64.0 / re.max(1e-9); - let f = if re <= RE_LAMINAR { - f_laminar - } else { - let f_turbulent = 0.079 * re.powf(-0.25); - if re >= RE_TURBULENT { - f_turbulent - } else { - entropyk_core::smoothing::cubic_blend( - f_laminar, - f_turbulent, - re, - RE_LAMINAR, - RE_TURBULENT, - ) - } + let z_dp = self.live_z_dp(state); + let (m_h, p_h_in, h_h_in) = edges.hot_in; + let (_, p_h_out, _) = edges.hot_out; + let (m_c, p_c_in, h_c_in) = edges.cold_in; + let (_, p_c_out, _) = edges.cold_out; + + let max_idx = [m_h, p_h_in, h_h_in, p_h_out, m_c, p_c_in, h_c_in, p_c_out] + .into_iter() + .max() + .unwrap_or(0); + if max_idx >= state.len() { + return Err(ComponentError::InvalidStateDimensions { + expected: max_idx + 1, + actual: state.len(), + }); + } + + let rho_hot = self.inner.side_density( + "hot", + self.inner.hot_fluid_id_str(), + state[p_h_in], + state[h_h_in], + )?; + let rho_cold = self.inner.side_density( + "cold", + self.inner.cold_fluid_id_str(), + state[p_c_in], + state[h_c_in], + )?; + + let dp_hot = self.side_channel_dp(state[m_h], rho_hot, z_dp)?; + let dp_cold = self.side_channel_dp(state[m_c], rho_cold, z_dp)?; + + residuals[row0] = state[p_h_out] - state[p_h_in] + dp_hot.delta_p_pa; + residuals[row0 + 1] = state[p_c_out] - state[p_c_in] + dp_cold.delta_p_pa; + Ok(()) + } + + /// Adds analytic ∂r/∂ṁ and ∂r/∂z_dp for BPHX pressure rows (P entries already set). + fn append_pressure_dp_jacobian( + &self, + state: &StateSlice, + jacobian: &mut JacobianBuilder, + ) -> Result<(), ComponentError> { + let Some(edges) = self.inner.four_port_edges() else { + return Ok(()); }; + let row0 = self.inner.model_n_equations(); + let z_dp = self.live_z_dp(state); + let (m_h, p_h_in, h_h_in) = edges.hot_in; + let (_, p_h_out, _) = edges.hot_out; + let (m_c, p_c_in, h_c_in) = edges.cold_in; + let (_, p_c_out, _) = edges.cold_out; - let dp_base = - 2.0 * f * self.geometry.plate_length * mass_flux.powi(2) / (rho * self.geometry.dh); + let max_idx = [m_h, p_h_in, h_h_in, p_h_out, m_c, p_c_in, h_c_in, p_c_out] + .into_iter() + .max() + .unwrap_or(0); + if max_idx >= state.len() { + return Err(ComponentError::InvalidStateDimensions { + expected: max_idx + 1, + actual: state.len(), + }); + } - dp_base * self.calib().z_dp + let rho_hot = self.inner.side_density( + "hot", + self.inner.hot_fluid_id_str(), + state[p_h_in], + state[h_h_in], + )?; + let rho_cold = self.inner.side_density( + "cold", + self.inner.cold_fluid_id_str(), + state[p_c_in], + state[h_c_in], + )?; + + let dp_hot = self.side_channel_dp(state[m_h], rho_hot, z_dp)?; + let dp_cold = self.side_channel_dp(state[m_c], rho_cold, z_dp)?; + + let a_flow = self.geometry.channel_flow_area() * self.geometry.n_channels() as f64; + if a_flow <= 1e-30 { + return Err(ComponentError::InvalidState( + "BPHX channel flow area too small for pressure-drop Jacobian".into(), + )); + } + let d_g_d_m = 1.0 / a_flow; + + // r = P_out − P_in + ΔP(G(ṁ), z_dp) ⇒ ∂r/∂ṁ = ∂ΔP/∂G · ∂G/∂ṁ + jacobian.add_entry(row0, m_h, dp_hot.d_delta_p_d_g * d_g_d_m); + jacobian.add_entry(row0 + 1, m_c, dp_cold.d_delta_p_d_g * d_g_d_m); + + if let Some(z_idx) = self.inner.calib_indices_ref().z_dp { + // ΔP = z_dp · ΔP_corr ⇒ ∂r/∂z_dp = ΔP_corr = ΔP / z_dp + // Keep the column even near z_dp→0 so a free actuator stays visible to Newton. + let z_safe = z_dp.max(1e-30); + jacobian.add_entry(row0, z_idx, dp_hot.delta_p_pa / z_safe); + jacobian.add_entry(row0 + 1, z_idx, dp_cold.delta_p_pa / z_safe); + } + Ok(()) } /// Updates UA based on computed HTC. @@ -413,7 +539,9 @@ impl Component for BphxExchanger { state: &StateSlice, residuals: &mut ResidualVector, ) -> Result<(), ComponentError> { - self.inner.compute_residuals(state, residuals) + self.inner.compute_residuals(state, residuals)?; + // Replace isobaric P_out − P_in with channel ΔP on both sides. + self.overwrite_pressure_closures_with_dp(state, residuals) } fn jacobian_entries( @@ -421,7 +549,9 @@ impl Component for BphxExchanger { state: &StateSlice, jacobian: &mut JacobianBuilder, ) -> Result<(), ComponentError> { - self.inner.jacobian_entries(state, jacobian) + self.inner.jacobian_entries(state, jacobian)?; + // Inner already wrote ∂r/∂P_out=+1, ∂r/∂P_in=−1; add ṁ and z_dp terms. + self.append_pressure_dp_jacobian(state, jacobian) } fn get_ports(&self) -> &[ConnectedPort] { @@ -844,4 +974,131 @@ mod tests { "ΔP slope kink near Re = 1 ({s_below:.3e} vs {s_above:.3e})" ); } + + #[test] + fn test_bphx_both_sides_pressure_residuals_include_dp() { + use entropyk_fluids::TestBackend; + use std::sync::Arc; + + // State layout per port triple (ṁ, P, h): + // hot_in 0..3, hot_out 3..6, cold_in 6..9, cold_out 9..12 + let mut hx = BphxExchanger::new(test_geometry()) + .with_fluid_backend(Arc::new(TestBackend::new())); + hx.set_hot_fluid("Water"); + hx.set_cold_fluid("Water"); + hx.set_port_context(&[ + Some((0, 1, 2)), + Some((3, 4, 5)), + Some((6, 7, 8)), + Some((9, 10, 11)), + ]); + + let mut state = vec![0.0; 12]; + // Liquid water-ish P-h (TestBackend): ~2 bar, ~80 kJ/kg + state[0] = 0.25; // m_hot + state[1] = 200_000.0; + state[2] = 80_000.0; + state[3] = 0.25; + state[4] = 200_000.0; // isobaric guess — residual should show +ΔP + state[5] = 70_000.0; + state[6] = 0.20; // m_cold + state[7] = 200_000.0; + state[8] = 50_000.0; + state[9] = 0.20; + state[10] = 200_000.0; + state[11] = 60_000.0; + + let n = hx.n_equations(); + assert!(n >= 4, "4-port BPHX should include two pressure rows"); + let mut residuals = vec![0.0; n]; + hx.compute_residuals(&state, &mut residuals).unwrap(); + + let row_p = hx.inner.model_n_equations(); + assert!( + residuals[row_p].abs() > 1.0, + "hot pressure residual should include ΔP, got {}", + residuals[row_p] + ); + assert!( + residuals[row_p + 1].abs() > 1.0, + "cold pressure residual should include ΔP, got {}", + residuals[row_p + 1] + ); + + // z_dp scale: half calib → half residual contribution at same P_out=P_in + let mut hx_half = BphxExchanger::new(test_geometry()) + .with_fluid_backend(Arc::new(TestBackend::new())); + hx_half.set_hot_fluid("Water"); + hx_half.set_cold_fluid("Water"); + hx_half.set_port_context(&[ + Some((0, 1, 2)), + Some((3, 4, 5)), + Some((6, 7, 8)), + Some((9, 10, 11)), + ]); + let mut calib = Calib::default(); + calib.z_dp = 0.5; + hx_half.set_calib(calib); + let mut r_half = vec![0.0; n]; + hx_half.compute_residuals(&state, &mut r_half).unwrap(); + assert!( + (r_half[row_p] - 0.5 * residuals[row_p]).abs() < 1e-6 * residuals[row_p].abs().max(1.0), + "z_dp should scale hot ΔP residual" + ); + assert!( + (r_half[row_p + 1] - 0.5 * residuals[row_p + 1]).abs() + < 1e-6 * residuals[row_p + 1].abs().max(1.0), + "z_dp should scale cold ΔP residual" + ); + } + + #[test] + fn test_bphx_martin_dp_in_residuals_differs() { + use entropyk_fluids::TestBackend; + use std::sync::Arc; + + let ports = [ + Some((0, 1, 2)), + Some((3, 4, 5)), + Some((6, 7, 8)), + Some((9, 10, 11)), + ]; + let mut state = vec![0.0; 12]; + state[0] = 0.30; + state[1] = 200_000.0; + state[2] = 80_000.0; + state[3] = 0.30; + state[4] = 200_000.0; + state[5] = 70_000.0; + state[6] = 0.30; + state[7] = 200_000.0; + state[8] = 50_000.0; + state[9] = 0.30; + state[10] = 200_000.0; + state[11] = 60_000.0; + + let mut simp = BphxExchanger::new(test_geometry()) + .with_fluid_backend(Arc::new(TestBackend::new())); + simp.set_hot_fluid("Water"); + simp.set_cold_fluid("Water"); + simp.set_port_context(&ports); + + let mut martin = BphxExchanger::new(test_geometry()) + .with_dp_correlation(BphxDpCorrelation::Martin1996) + .with_fluid_backend(Arc::new(TestBackend::new())); + martin.set_hot_fluid("Water"); + martin.set_cold_fluid("Water"); + martin.set_port_context(&ports); + + let n = simp.n_equations(); + let mut r_s = vec![0.0; n]; + let mut r_m = vec![0.0; n]; + simp.compute_residuals(&state, &mut r_s).unwrap(); + martin.compute_residuals(&state, &mut r_m).unwrap(); + let row = simp.inner.model_n_equations(); + assert!( + (r_s[row] - r_m[row]).abs() > 1.0, + "Martin residual should differ from Simplified" + ); + } } diff --git a/crates/components/src/heat_exchanger/correlation_registry.rs b/crates/components/src/heat_exchanger/correlation_registry.rs index 01753bb..3e39479 100644 --- a/crates/components/src/heat_exchanger/correlation_registry.rs +++ b/crates/components/src/heat_exchanger/correlation_registry.rs @@ -38,6 +38,14 @@ pub enum CorrelationId { Friedel1979, /// Müller-Steinhagen-Heck (1986) two-phase pressure-drop correlation. MullerSteinhagenHeck1986, + /// BPHX SimplifiedChannel single-phase channel friction (default ΔP). + SimplifiedChannel, + /// Martin (1996/1999) chevron-plate single-phase friction. + Martin1996, + /// Amalfi–Thome (2016) boiling plate ΔP — metadata only (not implemented). + AmalfiThome2016, + /// Tao & Infante Ferreira condensation-in-BPHE DP survey — metadata only. + TaoInfanteFerreiraSurvey, } impl CorrelationId { @@ -58,6 +66,10 @@ impl CorrelationId { Self::Mostinski1963 => "mostinski-1963", Self::Friedel1979 => "friedel-1979", Self::MullerSteinhagenHeck1986 => "muller-steinhagen-heck-1986", + Self::SimplifiedChannel => "simplified-channel", + Self::Martin1996 => "martin-1996", + Self::AmalfiThome2016 => "amalfi-thome-2016", + Self::TaoInfanteFerreiraSurvey => "tao-infante-ferreira-survey", } } } @@ -337,6 +349,28 @@ const POOL_BOILING_GEOMETRIES: &[ExchangerGeometryType] = &[ ExchangerGeometryType::FinnedTube, ]; const DP_TUBE_GEOMETRIES: &[ExchangerGeometryType] = &[ExchangerGeometryType::SmoothTube]; +const BPHX_DP_BOUNDS: &[BoundedQuantity] = &[ + bound( + BoundedQuantityKind::Reynolds, + 1.0, + 100_000.0, + SINGLE_PHASE_REGIMES, + ), + bound( + BoundedQuantityKind::MassFlux, + 1.0, + 2000.0, + SINGLE_PHASE_REGIMES, + ), +]; +const AMALFI_BOUNDS: &[BoundedQuantity] = &[ + bound(BoundedQuantityKind::MassFlux, 5.0, 600.0, EVAPORATION), + bound(BoundedQuantityKind::Quality, 0.0, 1.0, EVAPORATION), +]; +const TAO_BOUNDS: &[BoundedQuantity] = &[ + bound(BoundedQuantityKind::MassFlux, 5.0, 600.0, CONDENSATION), + bound(BoundedQuantityKind::Quality, 0.0, 1.0, CONDENSATION), +]; const fn bound( kind: BoundedQuantityKind, @@ -481,6 +515,42 @@ pub fn correlation_metadata(id: CorrelationId) -> CorrelationMetadata { RefrigerantApplicability::EvidenceIncomplete, "Muller-Steinhagen, H., Heck, K. (1986). Chem. Eng. Process. 20, 297-308", ), + CorrelationId::SimplifiedChannel => ( + "BPHX SimplifiedChannel (default)", + CorrelationPurpose::PressureDrop, + PLATES, + SINGLE_PHASE_REGIMES, + BPHX_DP_BOUNDS, + RefrigerantApplicability::FluidAgnostic, + "Entropyk channel Blasius/laminar C¹ blend; ΔP = z_dp·2 f L G²/(ρ dh)", + ), + CorrelationId::Martin1996 => ( + "Martin (1996/1999) chevron plate friction", + CorrelationPurpose::PressureDrop, + PLATES, + SINGLE_PHASE_REGIMES, + BPHX_DP_BOUNDS, + RefrigerantApplicability::FluidAgnostic, + "Martin, H. (1996). Chem. Eng. Process. 35:301–310; Martin 1999 / Shah appendix form", + ), + CorrelationId::AmalfiThome2016 => ( + "Amalfi–Thome (2016) boiling plate ΔP [not implemented]", + CorrelationPurpose::PressureDrop, + PLATES, + EVAPORATION, + AMALFI_BOUNDS, + RefrigerantApplicability::EvidenceIncomplete, + "Amalfi, Vakili-Farahani, Thome (2016). Int. J. Refrigeration — Parts 1–2; deferred", + ), + CorrelationId::TaoInfanteFerreiraSurvey => ( + "Tao & Infante Ferreira condensation BPHE DP survey [not implemented]", + CorrelationPurpose::PressureDrop, + PLATES, + CONDENSATION, + TAO_BOUNDS, + RefrigerantApplicability::EvidenceIncomplete, + "Tao & Infante Ferreira — condensation-in-BPHE pressure-drop survey; deferred", + ), }; CorrelationMetadata { id, @@ -516,6 +586,10 @@ pub fn registered_correlations() -> Vec { CorrelationId::Mostinski1963, CorrelationId::Friedel1979, CorrelationId::MullerSteinhagenHeck1986, + CorrelationId::SimplifiedChannel, + CorrelationId::Martin1996, + CorrelationId::AmalfiThome2016, + CorrelationId::TaoInfanteFerreiraSurvey, ] .into_iter() .map(correlation_metadata) @@ -1369,16 +1443,23 @@ mod tests { #[test] fn registry_does_not_claim_unverified_fluid_agnostic_coverage() { - assert!(registered_correlations().iter().all(|metadata| matches!( - metadata.domain.refrigerants, - RefrigerantApplicability::EvidenceIncomplete - ))); + // Geometry-based single-phase BPHX channel DP may declare FluidAgnostic; + // HTC / two-phase entries must stay EvidenceIncomplete until verified. + assert!(registered_correlations().iter().all(|metadata| { + matches!( + metadata.domain.refrigerants, + RefrigerantApplicability::EvidenceIncomplete + ) || matches!( + metadata.id, + CorrelationId::SimplifiedChannel | CorrelationId::Martin1996 + ) + })); } #[test] fn central_registry_contains_heat_transfer_and_pressure_drop_entries() { let registered = registered_correlations(); - assert_eq!(registered.len(), 14); + assert_eq!(registered.len(), 18); assert!(registered .iter() .any(|entry| entry.id == CorrelationId::Friedel1979 @@ -1387,6 +1468,18 @@ mod tests { entry.id == CorrelationId::MullerSteinhagenHeck1986 && entry.domain.purpose == CorrelationPurpose::PressureDrop })); + assert!(registered.iter().any(|entry| { + entry.id == CorrelationId::SimplifiedChannel + && entry.domain.purpose == CorrelationPurpose::PressureDrop + })); + assert!(registered.iter().any(|entry| { + entry.id == CorrelationId::Martin1996 + && entry.domain.purpose == CorrelationPurpose::PressureDrop + })); + assert!(registered.iter().any(|entry| { + entry.id == CorrelationId::AmalfiThome2016 + && entry.domain.purpose == CorrelationPurpose::PressureDrop + })); assert!(registered .iter() .any(|entry| entry.id == CorrelationId::Cooper1984)); diff --git a/crates/components/src/heat_exchanger/exchanger.rs b/crates/components/src/heat_exchanger/exchanger.rs index 9572113..85b15ca 100644 --- a/crates/components/src/heat_exchanger/exchanger.rs +++ b/crates/components/src/heat_exchanger/exchanger.rs @@ -163,6 +163,15 @@ impl HxSideConditions { } } +/// Live four-port edge indices: each triple is `(ṁ_idx, P_idx, h_idx)`. +#[derive(Debug, Clone, Copy)] +pub(crate) struct FourPortEdgeIndices { + pub hot_in: (usize, usize, usize), + pub hot_out: (usize, usize, usize), + pub cold_in: (usize, usize, usize), + pub cold_out: (usize, usize, usize), +} + /// Generic heat exchanger component with 4 ports. /// /// Uses the Strategy Pattern for heat transfer calculations via the @@ -516,7 +525,7 @@ impl HeatExchanger { } /// `true` when all 4 edges are wired (Modelica-style 4-port mode). - fn edges_ready(&self) -> bool { + pub(crate) fn edges_ready(&self) -> bool { self.hot_in_idx.is_some() && self.hot_out_idx.is_some() && self.cold_in_idx.is_some() @@ -525,6 +534,60 @@ impl HeatExchanger { && !self.cold_fluid_id_str.is_empty() } + /// Live four-port edge index triples `(ṁ, P, h)` when wired. + pub(crate) fn four_port_edges(&self) -> Option { + if !self.edges_ready() { + return None; + } + Some(FourPortEdgeIndices { + hot_in: self.hot_in_idx?, + hot_out: self.hot_out_idx?, + cold_in: self.cold_in_idx?, + cold_out: self.cold_out_idx?, + }) + } + + /// Thermal-model equation count (excludes the two pressure-closure rows). + pub(crate) fn model_n_equations(&self) -> usize { + self.model.n_equations() + } + + /// Calibration state indices (e.g. free `z_dp` actuator). + pub(crate) fn calib_indices_ref(&self) -> &entropyk_core::CalibIndices { + &self.calib_indices + } + + /// Hot-side fluid id used for live property queries. + pub(crate) fn hot_fluid_id_str(&self) -> &str { + &self.hot_fluid_id_str + } + + /// Cold-side fluid id used for live property queries. + pub(crate) fn cold_fluid_id_str(&self) -> &str { + &self.cold_fluid_id_str + } + + /// Inlet density [kg/m³] from the attached fluid backend at `(P, h)`. + pub(crate) fn side_density( + &self, + side: &str, + fluid_id: &str, + p_pa: f64, + h_jkg: f64, + ) -> Result { + self.query_live_property(side, fluid_id, Property::Density, p_pa, h_jkg) + .and_then(|rho| { + if rho.is_finite() && rho > 1e-10 { + Ok(rho) + } else { + Err(ComponentError::DomainViolation(DomainViolation { + component: Some(self.name.clone()), + detail: format!("{} {}-side density is invalid: {}", self.name, side, rho), + })) + } + }) + } + fn live_state_required_error(&self) -> ComponentError { ComponentError::InvalidState(format!( "{} requires live four-port edge state (hot_inlet, hot_outlet, cold_inlet, cold_outlet); inlet-only boundary conditions cannot define outlet states", diff --git a/crates/components/src/heat_exchanger/mod.rs b/crates/components/src/heat_exchanger/mod.rs index 46d1f32..ec7239a 100644 --- a/crates/components/src/heat_exchanger/mod.rs +++ b/crates/components/src/heat_exchanger/mod.rs @@ -50,6 +50,7 @@ pub mod air_cooled_condenser; pub mod bphx_condenser; pub mod bphx_correlation; +pub mod bphx_dp; pub mod bphx_evaporator; pub mod bphx_exchanger; pub mod bphx_geometry; @@ -83,6 +84,7 @@ pub use bphx_correlation::{ BphxCorrelation, CorrelationEvaluation, CorrelationParams, CorrelationResult, CorrelationSelector, ValidityStatus, }; +pub use bphx_dp::{evaluate_channel_pressure_drop, BphxDpCorrelation, BphxDpEval}; pub use bphx_evaporator::BphxEvaporator; pub use bphx_exchanger::BphxExchanger; pub use bphx_geometry::{BphxGeometry, BphxGeometryBuilder, BphxGeometryError, BphxType}; diff --git a/crates/components/src/lib.rs b/crates/components/src/lib.rs index 3512c95..b134c19 100644 --- a/crates/components/src/lib.rs +++ b/crates/components/src/lib.rs @@ -132,11 +132,11 @@ pub use free_cooling_exchanger::{ }; pub use heat_exchanger::model::FluidState; pub use heat_exchanger::{ - AirCooledCondenser, CoilGeometry, Condenser, CondenserCoil, CondenserRating, Economizer, - EpsNtuModel, Evaporator, EvaporatorCoil, EvaporatorRating, ExchangerType, FanCoilUnit, - FinCoilCondenser, FinType, FloodedCondenser, FloodedEvaporator, FloodedPoolBoilingConfig, - FlowConfiguration, GasCooler, HeatExchanger, HeatExchangerBuilder, HeatTransferModel, - HxSideConditions, LmtdModel, MchxCondenserCoil, ShellAndTubeHx, UaMode, + AirCooledCondenser, BphxDpCorrelation, CoilGeometry, Condenser, CondenserCoil, CondenserRating, + Economizer, EpsNtuModel, Evaporator, EvaporatorCoil, EvaporatorRating, ExchangerType, + FanCoilUnit, FinCoilCondenser, FinType, FloodedCondenser, FloodedEvaporator, + FloodedPoolBoilingConfig, FlowConfiguration, GasCooler, HeatExchanger, HeatExchangerBuilder, + HeatTransferModel, HxSideConditions, LmtdModel, MchxCondenserCoil, ShellAndTubeHx, UaMode, }; pub use heat_source::HeatSource; pub use isenthalpic_expansion_valve::IsenthalpicExpansionValve; diff --git a/crates/entropyk/src/lib.rs b/crates/entropyk/src/lib.rs index 5ff4cd2..8b6dbe8 100644 --- a/crates/entropyk/src/lib.rs +++ b/crates/entropyk/src/lib.rs @@ -174,8 +174,9 @@ pub use builder::{SystemBuilder, SystemBuilderError}; mod result; pub use result::{ - extract_simulation_result, ComponentResult, ConvergenceSummary, EdgeResult, EnergyResult, - PortState, SimulationOutcome, SimulationResult, SystemSummary, + extract_simulation_result, extract_solved_variables, ComponentResult, ConvergenceSummary, + EdgeResult, EnergyResult, PortState, SimulationOutcome, SimulationResult, SolvedVariable, + SystemSummary, }; // ============================================================================= diff --git a/crates/entropyk/src/result.rs b/crates/entropyk/src/result.rs index 81f46b8..4ee11cf 100644 --- a/crates/entropyk/src/result.rs +++ b/crates/entropyk/src/result.rs @@ -173,6 +173,36 @@ pub struct SystemSummary { pub cop_heating: Option, } +// ───────────────────────────────────────────────────────────────────────────── +// Solved unknown (free actuator / calibration factor) +// ───────────────────────────────────────────────────────────────────────────── + +/// A named solver-computed unknown with its converged value and physical bounds. +/// +/// These correspond to bounded control variables that the Newton solve treats +/// as unknowns — free actuators (e.g. expansion-valve `opening`, condenser +/// `fan_speed`), hard inverse-control links, and saturated-controller +/// actuators (calibration factors like `z_ua`, `z_dp`, `z_flow`). They ride +/// inside the flat `raw_state_vector` returned by the solver; this struct +/// attaches a stable id, owning component, short label, and `[min, max]` +/// bounds so the UI can surface them without parsing the raw vector. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(rename_all = "camelCase")] +pub struct SolvedVariable { + /// Stable id from [`BoundedVariableId`] (e.g. `"exv__opening"`). + pub id: String, + /// Owning component name, when known (`None` for globals). + pub component: Option, + /// Short human label, e.g. `"opening"`, `"z_ua"`, `"z_dp"`. + pub variable: String, + /// Converged value of the unknown. + pub value: f64, + /// Lower bound (inclusive). + pub min: f64, + /// Upper bound (inclusive). + pub max: f64, +} + // ───────────────────────────────────────────────────────────────────────────── // Top-level SimulationResult // ───────────────────────────────────────────────────────────────────────────── @@ -195,6 +225,12 @@ pub struct SimulationResult { pub edges: Vec, /// Aggregated system performance summary. pub summary: SystemSummary, + /// Named solver-computed unknowns (free actuators + calibration factors). + /// + /// Empty for systems with no bounded control variables. Skipped on + /// serialization when empty so older consumers keep working. + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub solved_variables: Vec, } impl SimulationResult { @@ -422,6 +458,9 @@ pub fn extract_simulation_result(system: &System, converged: &ConvergedState) -> status, }; + // --- Solved unknowns (free actuators + calibration factors) --- + let solved_variables = extract_solved_variables(system, state); + SimulationResult { status, convergence, @@ -429,6 +468,93 @@ pub fn extract_simulation_result(system: &System, converged: &ConvergedState) -> components, edges, summary, + solved_variables, + } +} + +/// Extracts named solver unknowns (bounded variables / free actuators) from a +/// solved state vector. +/// +/// Every bounded variable registered on `system` whose value lives in the +/// Newton state vector — as a physical free actuator, a hard inverse-control +/// link, or a saturated-controller actuator — is materialized into a +/// [`SolvedVariable`] carrying its converged value and `[min, max]` bounds. +/// +/// Bounded variables with no resolvable state slot (purely registered, never +/// wired into the solve) are silently skipped; a single `tracing::debug!` +/// line summarizes how many were skipped so coverage gaps are diagnosable +/// without panic risk. Out-of-range indices (state shorter than expected, +/// e.g. a truncated warm-start vector) are likewise skipped. +/// +/// # Arguments +/// +/// * `system` - The solved system (must be finalized). +/// * `state` - Converged Newton state vector slice (length = +/// `system.full_state_vector_len()` when complete). +pub fn extract_solved_variables(system: &System, state: &[f64]) -> Vec { + let mut out = Vec::new(); + let mut skipped = 0usize; + + for bv in system.bounded_variables() { + let id = bv.id(); + // Free actuators are laid out in their own block at the tail of the + // state vector; everything else (hard-control links and saturated- + // controller actuators) is resolvable via control_variable_state_index. + let idx = system + .free_actuators() + .position(|fid| fid == id) + .map(|i| system.free_actuator_index(i)) + .or_else(|| system.control_variable_state_index(id)); + + match idx { + Some(i) if i < state.len() => { + let id_str = id.as_str(); + let component = bv.component_id().map(|s| s.to_string()); + let variable = derive_solved_variable_label(id_str, bv.component_id()); + out.push(SolvedVariable { + id: id_str.to_string(), + component, + variable, + value: state[i], + min: bv.min(), + max: bv.max(), + }); + } + _ => { + skipped += 1; + } + } + } + + if skipped > 0 { + tracing::debug!( + skipped, + total = system.bounded_variable_count(), + "Bounded variables skipped when extracting solved variables (no resolvable state slot)" + ); + } + + out +} + +/// Derives the short user-facing variable label from a bounded-variable id. +/// +/// Bounded-variable ids follow the convention `"{component}__{factor}"`, e.g. +/// `"exv__opening"`, `"cond__z_ua"`. When the component is known we strip the +/// `"{component}__"` prefix; otherwise we fall back to the segment after the +/// last `__`. The solver-internal `"actuator"` suffix (used for saturated +/// controllers driving an `opening`/`injection` factor) is remapped to the +/// user-facing `"opening"` label. +fn derive_solved_variable_label(id: &str, component_id: Option<&str>) -> String { + let stripped = match component_id { + Some(c) => id + .strip_prefix(&format!("{}__", c)) + .unwrap_or(id), + None => id.rfind("__").map(|i| &id[i + 2..]).unwrap_or(id), + }; + match stripped { + "actuator" => "opening".to_string(), + other => other.to_string(), } } @@ -572,6 +698,7 @@ mod tests { components: vec![], edges: vec![], summary: SystemSummary::default(), + solved_variables: vec![], }; let json = result.to_json().unwrap(); assert!(json.contains("\"status\": \"converged\"")); @@ -595,4 +722,65 @@ mod tests { assert_eq!(result.edges, de.edges); assert_eq!(result.summary, de.summary); } + + #[test] + fn test_solved_variable_serialization() { + let sv = SolvedVariable { + id: "exv__opening".to_string(), + component: Some("exv".to_string()), + variable: "opening".to_string(), + value: 0.62, + min: 0.02, + max: 1.0, + }; + let json = serde_json::to_string(&sv).unwrap(); + let de: SolvedVariable = serde_json::from_str(&json).unwrap(); + assert_eq!(sv, de); + assert!(json.contains("\"component\":\"exv\"")); + assert!(json.contains("\"variable\":\"opening\"")); + } + + #[test] + fn test_solved_variables_skipped_when_empty() { + // Empty solved_variables should be omitted from serialized JSON + // (skip_serializing_if = "Vec::is_empty") for backward compatibility. + let result = SimulationResult { + status: SimulationOutcome::Converged, + convergence: ConvergenceSummary { + iterations: 1, + final_residual: 0.0, + converged: true, + status: SimulationOutcome::Converged, + }, + metadata: entropyk_solver::SimulationMetadata::new("h".to_string()), + components: vec![], + edges: vec![], + summary: SystemSummary::default(), + solved_variables: vec![], + }; + let json = result.to_json().unwrap(); + assert!(!json.contains("solvedVariables")); + } + + #[test] + fn test_derive_solved_variable_label() { + // Component known → strip "{component}__" prefix. + assert_eq!( + derive_solved_variable_label("exv__opening", Some("exv")), + "opening" + ); + assert_eq!( + derive_solved_variable_label("cond__z_ua", Some("cond")), + "z_ua" + ); + // Solver-internal "actuator" suffix remapped to user-facing "opening". + assert_eq!( + derive_solved_variable_label("exv__actuator", Some("exv")), + "opening" + ); + // Component unknown → fall back to last "__"-separated segment. + assert_eq!(derive_solved_variable_label("glob__z_dp", None), "z_dp"); + // No separator and no component → id itself. + assert_eq!(derive_solved_variable_label("global_var", None), "global_var"); + } }