From 8e467417b29d55711d4d70fd821d2626d2c24107 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sun, 29 Mar 2026 13:14:12 -0400 Subject: [PATCH 1/8] Add thermo polynomial calculations --- src/properties/thermo_fit.rs | 90 ++++++++++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) diff --git a/src/properties/thermo_fit.rs b/src/properties/thermo_fit.rs index 324762b..47e08c2 100644 --- a/src/properties/thermo_fit.rs +++ b/src/properties/thermo_fit.rs @@ -20,3 +20,93 @@ pub struct SpeciesElement { pub element: String, pub count: f64, } + +impl ThermoPolynomial { + /// Calculate using eq 4.9 from reference paper + /// NOTE: This is normalized and unitless + pub fn cp_over_r(self, temp: f64) -> f64 { + let inv_temp = 1.0 / temp; + self.a[0] * inv_temp * inv_temp + + self.a[1] * inv_temp + + self.a[2] + + self.a[3] * temp + + self.a[4] * temp * temp + + self.a[5] * temp * temp * temp + + self.a[6] * temp * temp * temp * temp + } + /// Calculate using eq 4.10 from reference paper + /// NOTE: This is normalized and unitless + pub fn h_over_rt(self, temp: f64) -> f64 { + let inv_temp = 1.0 / temp; + -self.a[0] * inv_temp * inv_temp + + self.a[1] * inv_temp * inv_temp.ln() + + self.a[2] + + self.a[3] * temp / 2.0 + + self.a[4] * temp * temp / 3.0 + + self.a[5] * temp * temp * temp / 4.0 + + self.a[6] * temp * temp * temp * temp / 5.0 + + self.a[7] * inv_temp + } + /// Calculate using eq 4.11 from reference paper + /// NOTE: This is normalized and unitless + pub fn s_over_r(self, temp: f64) -> f64 { + let inv_temp = 1.0 / temp; + -self.a[0] * inv_temp * inv_temp - self.a[1] * inv_temp + + self.a[2] * temp.ln() + + self.a[3] * temp + + self.a[4] * temp * temp / 2.0 + + self.a[5] * temp * temp * temp / 3.0 + + self.a[6] * temp * temp * temp * temp / 4.0 + + self.a[8] + } +} + +#[cfg(test)] +mod test { + use crate::assert_delta; + use super::ThermoPolynomial; + + fn poly(a: Vec) -> ThermoPolynomial { + ThermoPolynomial { a, temp_range: (0.0, 0.0) } + } + + #[test] + fn test_cp_over_r() { + // At T=1: a[0]/T^2 + a[1]/T + a[2] + a[3]*T + a[4]*T^2 + a[5]*T^3 + a[6]*T^4 + // = 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28 + let result = poly(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]).cp_over_r(1.0); + assert_delta!(result, 28.0, 1e-10); + + // At T=2: 4/4 + 2/2 + 1 + 2*2 + 1*4 + 1*8 + 1*16 = 35 + let result = poly(vec![4.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0]).cp_over_r(2.0); + assert_delta!(result, 35.0, 1e-10); + } + + #[test] + fn test_h_over_rt() { + // At T=1: ln(1/T)=0, so the a[1] term vanishes + // -a[0] + 0 + a[2] + a[3]/2 + a[4]/3 + a[5]/4 + a[6]/5 + a[7] + // = -1 + 0 + 3 + 12/2 + 9/3 + 8/4 + 5/5 + 8 = -1+3+6+3+2+1+8 = 22 + let result = poly(vec![1.0, 2.0, 3.0, 12.0, 9.0, 8.0, 5.0, 8.0]).h_over_rt(1.0); + assert_delta!(result, 22.0, 1e-10); + + // At T=2: a[1]=0 so log term vanishes; a[4..6]=0 to avoid fractions + // -4/4 + 0 + 3 + 4*2/2 + 0 + 0 + 0 + 2/2 = -1+3+4+1 = 7 + let result = poly(vec![4.0, 0.0, 3.0, 4.0, 0.0, 0.0, 0.0, 2.0]).h_over_rt(2.0); + assert_delta!(result, 7.0, 1e-10); + } + + #[test] + fn test_s_over_r() { + // At T=1: ln(T)=0, so a[2] term vanishes + // -a[0] - a[1] + 0 + a[3] + a[4]/2 + a[5]/3 + a[6]/4 + a[8] + // = -1-2+0+4+6/2+12/3+8/4+5 = -1-2+4+3+4+2+5 = 15 + let result = poly(vec![1.0, 2.0, 3.0, 4.0, 6.0, 12.0, 8.0, 0.0, 5.0]).s_over_r(1.0); + assert_delta!(result, 15.0, 1e-10); + + // At T=2: a[2]=0 so log term vanishes; a[5..6]=0 to avoid fractions + // -4/4 - 2/2 + 0 + 2*2 + 2*4/2 + 0 + 0 + 3 = -1-1+4+4+3 = 9 + let result = poly(vec![4.0, 2.0, 0.0, 2.0, 2.0, 0.0, 0.0, 0.0, 3.0]).s_over_r(2.0); + assert_delta!(result, 9.0, 1e-10); + } +} From 13e7f7be165d1fbb48329dc075415f7b3dac84ed Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sun, 29 Mar 2026 18:03:33 -0400 Subject: [PATCH 2/8] Add wrapper calculations to select correct polynomial --- src/properties/thermo_fit.rs | 57 ++++++++++++++++++++++++++++++++---- 1 file changed, 52 insertions(+), 5 deletions(-) diff --git a/src/properties/thermo_fit.rs b/src/properties/thermo_fit.rs index 47e08c2..a5eba68 100644 --- a/src/properties/thermo_fit.rs +++ b/src/properties/thermo_fit.rs @@ -21,10 +21,54 @@ pub struct SpeciesElement { pub count: f64, } +impl SpeciesThermoData { + pub fn cp_over_r(&self, temp: f64) -> f64 { + //TODO: Not the most efficient. Can refactor to pre-compute tables + //and do 1-d linear interpolation if needed + // + //TODO: I Think condensed species need to be treated differently. Verify how that works in + //the paper. + let i_polynomial = self + .polynomials + .iter() + .rposition(|polynomial| temp > polynomial.temp_range.0) + .unwrap_or(0); + self.polynomials[i_polynomial].cp_over_r(temp) + } + + pub fn h_over_rt(&self, temp: f64) -> f64 { + //TODO: Not the most efficient. Can refactor to pre-compute tables + //and do 1-d linear interpolation if needed + // + //TODO: I Think condensed species need to be treated differently. Verify how that works in + //the paper. + let i_polynomial = self + .polynomials + .iter() + .rposition(|polynomial| temp > polynomial.temp_range.0) + .unwrap_or(0); + self.polynomials[i_polynomial].h_over_rt(temp) + } + + pub fn s_over_r(&self, temp: f64) -> f64 { + //TODO: Not the most efficient. Can refactor to pre-compute tables + //and do 1-d linear interpolation if needed + // + //TODO: I Think condensed species need to be treated differently. Verify how that works in + //the paper. + let i_polynomial = self + .polynomials + .iter() + .rposition(|polynomial| temp > polynomial.temp_range.0) + .unwrap_or(0); + self.polynomials[i_polynomial].s_over_r(temp) + } +} + impl ThermoPolynomial { /// Calculate using eq 4.9 from reference paper /// NOTE: This is normalized and unitless - pub fn cp_over_r(self, temp: f64) -> f64 { + pub fn cp_over_r(&self, temp: f64) -> f64 { let inv_temp = 1.0 / temp; self.a[0] * inv_temp * inv_temp + self.a[1] * inv_temp @@ -36,7 +80,7 @@ impl ThermoPolynomial { } /// Calculate using eq 4.10 from reference paper /// NOTE: This is normalized and unitless - pub fn h_over_rt(self, temp: f64) -> f64 { + pub fn h_over_rt(&self, temp: f64) -> f64 { let inv_temp = 1.0 / temp; -self.a[0] * inv_temp * inv_temp + self.a[1] * inv_temp * inv_temp.ln() @@ -49,7 +93,7 @@ impl ThermoPolynomial { } /// Calculate using eq 4.11 from reference paper /// NOTE: This is normalized and unitless - pub fn s_over_r(self, temp: f64) -> f64 { + pub fn s_over_r(&self, temp: f64) -> f64 { let inv_temp = 1.0 / temp; -self.a[0] * inv_temp * inv_temp - self.a[1] * inv_temp + self.a[2] * temp.ln() @@ -63,11 +107,14 @@ impl ThermoPolynomial { #[cfg(test)] mod test { - use crate::assert_delta; use super::ThermoPolynomial; + use crate::assert_delta; fn poly(a: Vec) -> ThermoPolynomial { - ThermoPolynomial { a, temp_range: (0.0, 0.0) } + ThermoPolynomial { + a, + temp_range: (0.0, 0.0), + } } #[test] From 6d9175aa484c41d386866bdbece97afb5241a372 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sun, 29 Mar 2026 18:19:07 -0400 Subject: [PATCH 3/8] Reorganize the bad forwarder pattern and add tests --- src/properties/thermo_fit.rs | 122 ++++++++++++++++++++++------------- 1 file changed, 78 insertions(+), 44 deletions(-) diff --git a/src/properties/thermo_fit.rs b/src/properties/thermo_fit.rs index a5eba68..7ed1abb 100644 --- a/src/properties/thermo_fit.rs +++ b/src/properties/thermo_fit.rs @@ -22,7 +22,7 @@ pub struct SpeciesElement { } impl SpeciesThermoData { - pub fn cp_over_r(&self, temp: f64) -> f64 { + pub fn polynomial_at(&self, temp: f64) -> &ThermoPolynomial { //TODO: Not the most efficient. Can refactor to pre-compute tables //and do 1-d linear interpolation if needed // @@ -33,35 +33,7 @@ impl SpeciesThermoData { .iter() .rposition(|polynomial| temp > polynomial.temp_range.0) .unwrap_or(0); - self.polynomials[i_polynomial].cp_over_r(temp) - } - - pub fn h_over_rt(&self, temp: f64) -> f64 { - //TODO: Not the most efficient. Can refactor to pre-compute tables - //and do 1-d linear interpolation if needed - // - //TODO: I Think condensed species need to be treated differently. Verify how that works in - //the paper. - let i_polynomial = self - .polynomials - .iter() - .rposition(|polynomial| temp > polynomial.temp_range.0) - .unwrap_or(0); - self.polynomials[i_polynomial].h_over_rt(temp) - } - - pub fn s_over_r(&self, temp: f64) -> f64 { - //TODO: Not the most efficient. Can refactor to pre-compute tables - //and do 1-d linear interpolation if needed - // - //TODO: I Think condensed species need to be treated differently. Verify how that works in - //the paper. - let i_polynomial = self - .polynomials - .iter() - .rposition(|polynomial| temp > polynomial.temp_range.0) - .unwrap_or(0); - self.polynomials[i_polynomial].s_over_r(temp) + &self.polynomials[i_polynomial] } } @@ -108,24 +80,28 @@ impl ThermoPolynomial { #[cfg(test)] mod test { use super::ThermoPolynomial; - use crate::assert_delta; - - fn poly(a: Vec) -> ThermoPolynomial { - ThermoPolynomial { - a, - temp_range: (0.0, 0.0), - } - } + use crate::{ + assert_delta, assert_vec_delta, + properties::thermo_fit::{Phase, SpeciesThermoData}, + }; #[test] fn test_cp_over_r() { // At T=1: a[0]/T^2 + a[1]/T + a[2] + a[3]*T + a[4]*T^2 + a[5]*T^3 + a[6]*T^4 // = 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28 - let result = poly(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]).cp_over_r(1.0); + let result = ThermoPolynomial { + a: vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0], + temp_range: (0.0, 0.0), + } + .cp_over_r(1.0); assert_delta!(result, 28.0, 1e-10); // At T=2: 4/4 + 2/2 + 1 + 2*2 + 1*4 + 1*8 + 1*16 = 35 - let result = poly(vec![4.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0]).cp_over_r(2.0); + let result = ThermoPolynomial { + a: vec![4.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0], + temp_range: (0.0, 0.0), + } + .cp_over_r(2.0); assert_delta!(result, 35.0, 1e-10); } @@ -134,12 +110,20 @@ mod test { // At T=1: ln(1/T)=0, so the a[1] term vanishes // -a[0] + 0 + a[2] + a[3]/2 + a[4]/3 + a[5]/4 + a[6]/5 + a[7] // = -1 + 0 + 3 + 12/2 + 9/3 + 8/4 + 5/5 + 8 = -1+3+6+3+2+1+8 = 22 - let result = poly(vec![1.0, 2.0, 3.0, 12.0, 9.0, 8.0, 5.0, 8.0]).h_over_rt(1.0); + let result = ThermoPolynomial { + a: vec![1.0, 2.0, 3.0, 12.0, 9.0, 8.0, 5.0, 8.0], + temp_range: (0.0, 0.0), + } + .h_over_rt(1.0); assert_delta!(result, 22.0, 1e-10); // At T=2: a[1]=0 so log term vanishes; a[4..6]=0 to avoid fractions // -4/4 + 0 + 3 + 4*2/2 + 0 + 0 + 0 + 2/2 = -1+3+4+1 = 7 - let result = poly(vec![4.0, 0.0, 3.0, 4.0, 0.0, 0.0, 0.0, 2.0]).h_over_rt(2.0); + let result = ThermoPolynomial { + a: vec![4.0, 0.0, 3.0, 4.0, 0.0, 0.0, 0.0, 2.0], + temp_range: (0.0, 0.0), + } + .h_over_rt(2.0); assert_delta!(result, 7.0, 1e-10); } @@ -148,12 +132,62 @@ mod test { // At T=1: ln(T)=0, so a[2] term vanishes // -a[0] - a[1] + 0 + a[3] + a[4]/2 + a[5]/3 + a[6]/4 + a[8] // = -1-2+0+4+6/2+12/3+8/4+5 = -1-2+4+3+4+2+5 = 15 - let result = poly(vec![1.0, 2.0, 3.0, 4.0, 6.0, 12.0, 8.0, 0.0, 5.0]).s_over_r(1.0); + let result = ThermoPolynomial { + a: vec![1.0, 2.0, 3.0, 4.0, 6.0, 12.0, 8.0, 0.0, 5.0], + temp_range: (0.0, 0.0), + } + .s_over_r(1.0); assert_delta!(result, 15.0, 1e-10); // At T=2: a[2]=0 so log term vanishes; a[5..6]=0 to avoid fractions // -4/4 - 2/2 + 0 + 2*2 + 2*4/2 + 0 + 0 + 3 = -1-1+4+4+3 = 9 - let result = poly(vec![4.0, 2.0, 0.0, 2.0, 2.0, 0.0, 0.0, 0.0, 3.0]).s_over_r(2.0); + let result = ThermoPolynomial { + a: vec![4.0, 2.0, 0.0, 2.0, 2.0, 0.0, 0.0, 0.0, 3.0], + temp_range: (0.0, 0.0), + } + .s_over_r(2.0); assert_delta!(result, 9.0, 1e-10); } + + #[test] + fn test_polynomial_at() { + let polynomials = vec![ + ThermoPolynomial { + a: vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0], + temp_range: (1.0, 100.0), + }, + ThermoPolynomial { + a: vec![4.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0], + temp_range: (100.0, 300.0), + }, + ]; + + let data = SpeciesThermoData { + name: "".to_string(), + elements: vec![], + phase: Phase::Gas, + polynomials, + molecular_weight: 1.0, + h_formation: 1.0, + }; + + assert!(std::ptr::eq(data.polynomial_at(0.5), &data.polynomials[0])); + assert!(std::ptr::eq(data.polynomial_at(50.0), &data.polynomials[0])); + assert!(std::ptr::eq( + data.polynomial_at(100.0), + &data.polynomials[0] + )); + assert!(std::ptr::eq( + data.polynomial_at(100.0 + 1e-12), + &data.polynomials[1] + )); + assert!(std::ptr::eq( + data.polynomial_at(150.0 + 1e-12), + &data.polynomials[1] + )); + assert!(std::ptr::eq( + data.polynomial_at(500.0 + 1e-12), + &data.polynomials[1] + )); + } } From aa4e8707758a235fa0c4b2b101b0a3b30a2e64a7 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sun, 29 Mar 2026 19:17:33 -0400 Subject: [PATCH 4/8] Add transport property calculations --- src/properties/transport_fit.rs | 165 ++++++++++++++++++++++++++++++++ 1 file changed, 165 insertions(+) diff --git a/src/properties/transport_fit.rs b/src/properties/transport_fit.rs index 8489d3c..0f5ea20 100644 --- a/src/properties/transport_fit.rs +++ b/src/properties/transport_fit.rs @@ -4,6 +4,26 @@ pub struct SpeciesTransportData { pub conductivities: Vec, } +impl SpeciesTransportData { + pub fn viscosity_at(&self, temp: f64) -> f64 { + let i_viscosity = self + .viscosities + .iter() + .rposition(|viscosity| temp > viscosity.temp_range.0) + .unwrap_or(0); + self.viscosities[i_viscosity].compute(temp) + } + + pub fn conductivity_at(&self, temp: f64) -> f64 { + let i_conductivity = self + .conductivities + .iter() + .rposition(|conductivity| temp > conductivity.temp_range.0) + .unwrap_or(0); + self.conductivities[i_conductivity].compute(temp) + } +} + pub struct TransportFit { pub temp_range: (f64, f64), pub a: f64, @@ -11,3 +31,148 @@ pub struct TransportFit { pub c: f64, pub d: f64, } + +impl TransportFit { + pub fn compute(&self, temp: f64) -> f64 { + let inv_temp = 1.0 / temp; + self.a * temp.ln() + self.b * inv_temp + self.c * inv_temp * inv_temp + self.d + } +} + +#[cfg(test)] +mod test { + use crate::{ + assert_delta, + properties::transport_fit::{SpeciesTransportData, TransportFit}, + }; + + #[test] + fn test_transport_fit_compute() { + let fit = TransportFit { + temp_range: (1.0, 2.0), + a: 10.0, + b: 20.0, + c: 30.0, + d: 40.0, + }; + + assert_delta!(fit.compute(1.0), 90.0, 1e-12); + } + + #[test] + fn test_calc_transport_properties() { + let viscosities = vec![ + TransportFit { + temp_range: (1.0, 100.0), + a: 10.0, + b: 20.0, + c: 30.0, + d: 40.0, + }, + TransportFit { + temp_range: (100.0, 200.0), + a: 1.0, + b: 2.0, + c: 3.0, + d: 4.0, + }, + ]; + + let conductivities = vec![ + TransportFit { + temp_range: (1.0, 100.0), + a: 10.0, + b: 20.0, + c: 30.0, + d: 40.0, + }, + TransportFit { + temp_range: (100.0, 200.0), + a: 1.0, + b: 2.0, + c: 3.0, + d: 4.0, + }, + ]; + + let data = SpeciesTransportData { + viscosities, + conductivities, + name: "".to_string(), + }; + + assert_delta!( + data.conductivity_at(0.5), + data.conductivities[0].compute(0.5), + 1e-12 + ); + assert_delta!( + data.conductivity_at(1.0), + data.conductivities[0].compute(1.0), + 1e-12 + ); + assert_delta!( + data.conductivity_at(50.0), + data.conductivities[0].compute(50.0), + 1e-12 + ); + assert_delta!( + data.conductivity_at(100.0), + data.conductivities[0].compute(100.0), + 1e-12 + ); + assert_delta!( + data.conductivity_at(100.0 + 1e-12), + data.conductivities[1].compute(100.0 + 1e-12), + 1e-12 + ); + + assert_delta!( + data.conductivity_at(200.0), + data.conductivities[1].compute(200.0), + 1e-12 + ); + assert_delta!( + data.conductivity_at(500.0), + data.conductivities[1].compute(500.0), + 1e-12 + ); + + assert_delta!( + data.viscosity_at(0.5), + data.viscosities[0].compute(0.5), + 1e-12 + ); + assert_delta!( + data.viscosity_at(1.0), + data.viscosities[0].compute(1.0), + 1e-12 + ); + assert_delta!( + data.viscosity_at(50.0), + data.viscosities[0].compute(50.0), + 1e-12 + ); + assert_delta!( + data.viscosity_at(100.0), + data.viscosities[0].compute(100.0), + 1e-12 + ); + assert_delta!( + data.viscosity_at(100.0 + 1e-12), + data.viscosities[1].compute(100.0 + 1e-12), + 1e-12 + ); + + assert_delta!( + data.viscosity_at(200.0), + data.viscosities[1].compute(200.0), + 1e-12 + ); + assert_delta!( + data.viscosity_at(500.0), + data.viscosities[1].compute(500.0), + 1e-12 + ); + } +} From cdb74a4463a8fc4aa214c1f07e56291bf4212a09 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sun, 29 Mar 2026 20:50:39 -0400 Subject: [PATCH 5/8] Make field private and refactor --- src/properties/thermo_db.rs | 48 ++++++++++++++++----------------- src/properties/thermo_fit.rs | 24 ++++++++++++++++- src/properties/transport_db.rs | 6 ++--- src/properties/transport_fit.rs | 16 +++++++++-- 4 files changed, 64 insertions(+), 30 deletions(-) diff --git a/src/properties/thermo_db.rs b/src/properties/thermo_db.rs index 96cc957..4468874 100644 --- a/src/properties/thermo_db.rs +++ b/src/properties/thermo_db.rs @@ -98,14 +98,14 @@ fn parse_species<'a>( lines.next().ok_or(PropertiesError::InvalidFile)?; } - Ok(SpeciesThermoData { + Ok(SpeciesThermoData::new( name, - polynomials, elements, phase, + polynomials, molecular_weight, h_formation, - }) + )) } fn parse_polynomials_block<'a>( @@ -295,9 +295,9 @@ mod test { -1.742171366e+01, ]; - assert_vec_delta!(species.polynomials[0].a, real_coeff_1, 1e-9); - assert_delta!(species.polynomials[0].temp_range.0, 300.000, 1e-3); - assert_delta!(species.polynomials[0].temp_range.1, 1000.000, 1e-3); + assert_vec_delta!(species.polynomial_at(650.0).a, real_coeff_1, 1e-9); + assert_delta!(species.polynomial_at(650.0).temp_range.0, 300.000, 1e-3); + assert_delta!(species.polynomial_at(650.0).temp_range.1, 1000.000, 1e-3); let real_coeff_2 = [ -3.523782900e+05, @@ -311,9 +311,9 @@ mod test { -2.695610360e+00, ]; - assert_vec_delta!(species.polynomials[1].a, real_coeff_2, 1e-9); - assert_delta!(species.polynomials[1].temp_range.0, 1000.000, 1e-3); - assert_delta!(species.polynomials[1].temp_range.1, 6000.000, 1e-3); + assert_vec_delta!(species.polynomial_at(3500.0).a, real_coeff_2, 1e-9); + assert_delta!(species.polynomial_at(3500.0).temp_range.0, 1000.000, 1e-3); + assert_delta!(species.polynomial_at(3500.0).temp_range.1, 6000.000, 1e-3); } #[test] @@ -363,10 +363,10 @@ END REACTANTS assert!(matches!(alcl3.phase, Phase::Gas)); assert_delta!(alcl3.molecular_weight, 133.3405380, 1e-7); assert_delta!(alcl3.h_formation, -584678.863, 1e-3); - assert_eq!(alcl3.polynomials.len(), 2); + assert_eq!(alcl3.num_polynomials(), 2); assert_vec_delta!( - alcl3.polynomials[0].a, + alcl3.polynomial_at(650.0).a, [ 7.750600970e+04, -1.440779717e+03, @@ -380,11 +380,11 @@ END REACTANTS ], 1e-9 ); - assert_delta!(alcl3.polynomials[0].temp_range.0, 300.0, 1e-3); - assert_delta!(alcl3.polynomials[0].temp_range.1, 1000.0, 1e-3); + assert_delta!(alcl3.polynomial_at(650.0).temp_range.0, 300.0, 1e-3); + assert_delta!(alcl3.polynomial_at(650.0).temp_range.1, 1000.0, 1e-3); assert_vec_delta!( - alcl3.polynomials[1].a, + alcl3.polynomial_at(3500.0).a, [ -1.378630916e+05, -5.579207290e+01, @@ -398,8 +398,8 @@ END REACTANTS ], 1e-9 ); - assert_delta!(alcl3.polynomials[1].temp_range.0, 1000.0, 1e-3); - assert_delta!(alcl3.polynomials[1].temp_range.1, 6000.0, 1e-3); + assert_delta!(alcl3.polynomial_at(3500.0).temp_range.0, 1000.0, 1e-3); + assert_delta!(alcl3.polynomial_at(3500.0).temp_range.1, 6000.0, 1e-3); // --- Air (reactant 0) --- let air = &thermo_db.reactants[0]; @@ -416,10 +416,10 @@ END REACTANTS assert!(matches!(air.phase, Phase::Gas)); assert_delta!(air.molecular_weight, 28.9651159, 1e-7); assert_delta!(air.h_formation, -125.530, 1e-3); - assert_eq!(air.polynomials.len(), 2); + assert_eq!(air.num_polynomials(), 2); assert_vec_delta!( - air.polynomials[0].a, + air.polynomial_at(650.0).a, [ 1.009950160e+04, -1.968275610e+02, @@ -433,11 +433,11 @@ END REACTANTS ], 1e-9 ); - assert_delta!(air.polynomials[0].temp_range.0, 300.0, 1e-3); - assert_delta!(air.polynomials[0].temp_range.1, 1000.0, 1e-3); + assert_delta!(air.polynomial_at(650.0).temp_range.0, 300.0, 1e-3); + assert_delta!(air.polynomial_at(650.0).temp_range.1, 1000.0, 1e-3); assert_vec_delta!( - air.polynomials[1].a, + air.polynomial_at(3500.0).a, [ 2.415214430e+05, -1.257874600e+03, @@ -451,8 +451,8 @@ END REACTANTS ], 1e-9 ); - assert_delta!(air.polynomials[1].temp_range.0, 1000.0, 1e-3); - assert_delta!(air.polynomials[1].temp_range.1, 6000.0, 1e-3); + assert_delta!(air.polynomial_at(3500.0).temp_range.0, 1000.0, 1e-3); + assert_delta!(air.polynomial_at(3500.0).temp_range.1, 6000.0, 1e-3); // --- n-Butanol (reactant 1) --- let butanol = &thermo_db.reactants[1]; @@ -467,6 +467,6 @@ END REACTANTS assert!(matches!(butanol.phase, Phase::Condensed)); assert_delta!(butanol.molecular_weight, 74.1216000, 1e-7); assert_delta!(butanol.h_formation, -278510.000, 1e-3); - assert_eq!(butanol.polynomials.len(), 0); + assert_eq!(butanol.num_polynomials(), 0); } } diff --git a/src/properties/thermo_fit.rs b/src/properties/thermo_fit.rs index 7ed1abb..8a4c3fd 100644 --- a/src/properties/thermo_fit.rs +++ b/src/properties/thermo_fit.rs @@ -7,7 +7,7 @@ pub struct SpeciesThermoData { pub name: String, pub elements: Vec, pub phase: Phase, - pub polynomials: Vec, + polynomials: Vec, pub molecular_weight: f64, pub h_formation: f64, } @@ -22,6 +22,28 @@ pub struct SpeciesElement { } impl SpeciesThermoData { + pub fn new( + name: String, + elements: Vec, + phase: Phase, + polynomials: Vec, + molecular_weight: f64, + h_formation: f64, + ) -> Self { + Self { + name, + elements, + phase, + polynomials, + molecular_weight, + h_formation, + } + } + + pub fn num_polynomials(&self) -> usize { + self.polynomials.len() + } + pub fn polynomial_at(&self, temp: f64) -> &ThermoPolynomial { //TODO: Not the most efficient. Can refactor to pre-compute tables //and do 1-d linear interpolation if needed diff --git a/src/properties/transport_db.rs b/src/properties/transport_db.rs index 2ebb5f5..32bcc00 100644 --- a/src/properties/transport_db.rs +++ b/src/properties/transport_db.rs @@ -50,11 +50,11 @@ fn parse_species_transport_block<'a>( (ViscosityOrConductivity::Conductivity, fit) => conductivities.push(fit), } } - Ok(SpeciesTransportData { - name, + Ok(SpeciesTransportData::new( + &name, viscosities, conductivities, - }) + )) } fn parse_species_header_line(line: &str) -> Result<(String, usize, usize), PropertiesError> { diff --git a/src/properties/transport_fit.rs b/src/properties/transport_fit.rs index 0f5ea20..60a7213 100644 --- a/src/properties/transport_fit.rs +++ b/src/properties/transport_fit.rs @@ -1,10 +1,22 @@ pub struct SpeciesTransportData { pub name: String, - pub viscosities: Vec, - pub conductivities: Vec, + viscosities: Vec, + conductivities: Vec, } impl SpeciesTransportData { + pub fn new( + name: &str, + viscosities: Vec, + conductivities: Vec, + ) -> Self { + SpeciesTransportData { + name: name.to_string(), + viscosities, + conductivities, + } + } + pub fn viscosity_at(&self, temp: f64) -> f64 { let i_viscosity = self .viscosities From 6f4bc45927d65212d38c5ac3e81cda83c4c3a763 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sun, 29 Mar 2026 22:26:59 -0400 Subject: [PATCH 6/8] Fix some calculation issues --- src/properties/thermo_db.rs | 2 +- src/properties/thermo_fit.rs | 39 +++++++++++---------------------- src/properties/transport_fit.rs | 4 ++-- 3 files changed, 16 insertions(+), 29 deletions(-) diff --git a/src/properties/thermo_db.rs b/src/properties/thermo_db.rs index 4468874..5eecb2d 100644 --- a/src/properties/thermo_db.rs +++ b/src/properties/thermo_db.rs @@ -99,7 +99,7 @@ fn parse_species<'a>( } Ok(SpeciesThermoData::new( - name, + &name, elements, phase, polynomials, diff --git a/src/properties/thermo_fit.rs b/src/properties/thermo_fit.rs index 8a4c3fd..2d5bdf0 100644 --- a/src/properties/thermo_fit.rs +++ b/src/properties/thermo_fit.rs @@ -23,7 +23,7 @@ pub struct SpeciesElement { impl SpeciesThermoData { pub fn new( - name: String, + name: &str, elements: Vec, phase: Phase, polynomials: Vec, @@ -31,7 +31,7 @@ impl SpeciesThermoData { h_formation: f64, ) -> Self { Self { - name, + name: name.to_string(), elements, phase, polynomials, @@ -77,7 +77,7 @@ impl ThermoPolynomial { pub fn h_over_rt(&self, temp: f64) -> f64 { let inv_temp = 1.0 / temp; -self.a[0] * inv_temp * inv_temp - + self.a[1] * inv_temp * inv_temp.ln() + + self.a[1] * inv_temp * temp.ln() + self.a[2] + self.a[3] * temp / 2.0 + self.a[4] * temp * temp / 3.0 @@ -89,7 +89,7 @@ impl ThermoPolynomial { /// NOTE: This is normalized and unitless pub fn s_over_r(&self, temp: f64) -> f64 { let inv_temp = 1.0 / temp; - -self.a[0] * inv_temp * inv_temp - self.a[1] * inv_temp + -self.a[0] * inv_temp * inv_temp / 2.0 - self.a[1] * inv_temp + self.a[2] * temp.ln() + self.a[3] * temp + self.a[4] * temp * temp / 2.0 @@ -103,14 +103,12 @@ impl ThermoPolynomial { mod test { use super::ThermoPolynomial; use crate::{ - assert_delta, assert_vec_delta, + assert_delta, properties::thermo_fit::{Phase, SpeciesThermoData}, }; #[test] fn test_cp_over_r() { - // At T=1: a[0]/T^2 + a[1]/T + a[2] + a[3]*T + a[4]*T^2 + a[5]*T^3 + a[6]*T^4 - // = 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28 let result = ThermoPolynomial { a: vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0], temp_range: (0.0, 0.0), @@ -118,7 +116,6 @@ mod test { .cp_over_r(1.0); assert_delta!(result, 28.0, 1e-10); - // At T=2: 4/4 + 2/2 + 1 + 2*2 + 1*4 + 1*8 + 1*16 = 35 let result = ThermoPolynomial { a: vec![4.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0], temp_range: (0.0, 0.0), @@ -129,46 +126,36 @@ mod test { #[test] fn test_h_over_rt() { - // At T=1: ln(1/T)=0, so the a[1] term vanishes - // -a[0] + 0 + a[2] + a[3]/2 + a[4]/3 + a[5]/4 + a[6]/5 + a[7] - // = -1 + 0 + 3 + 12/2 + 9/3 + 8/4 + 5/5 + 8 = -1+3+6+3+2+1+8 = 22 let result = ThermoPolynomial { a: vec![1.0, 2.0, 3.0, 12.0, 9.0, 8.0, 5.0, 8.0], temp_range: (0.0, 0.0), } - .h_over_rt(1.0); - assert_delta!(result, 22.0, 1e-10); + .h_over_rt(2.0); + assert_delta!(result, 63.44314718055995, 1e-10); - // At T=2: a[1]=0 so log term vanishes; a[4..6]=0 to avoid fractions - // -4/4 + 0 + 3 + 4*2/2 + 0 + 0 + 0 + 2/2 = -1+3+4+1 = 7 let result = ThermoPolynomial { a: vec![4.0, 0.0, 3.0, 4.0, 0.0, 0.0, 0.0, 2.0], temp_range: (0.0, 0.0), } - .h_over_rt(2.0); - assert_delta!(result, 7.0, 1e-10); + .h_over_rt(100.0); + assert_delta!(result, 203.0196, 1e-4); } #[test] fn test_s_over_r() { - // At T=1: ln(T)=0, so a[2] term vanishes - // -a[0] - a[1] + 0 + a[3] + a[4]/2 + a[5]/3 + a[6]/4 + a[8] - // = -1-2+0+4+6/2+12/3+8/4+5 = -1-2+4+3+4+2+5 = 15 let result = ThermoPolynomial { - a: vec![1.0, 2.0, 3.0, 4.0, 6.0, 12.0, 8.0, 0.0, 5.0], + a: vec![2.0, 3.0, 4.0, 1.0, 5.0, 2.0, 8.0, 1.0, 12.0], temp_range: (0.0, 0.0), } - .s_over_r(1.0); - assert_delta!(result, 15.0, 1e-10); + .s_over_r(100.0); + assert_delta!(result, 200691797.0572474, 1e-7); - // At T=2: a[2]=0 so log term vanishes; a[5..6]=0 to avoid fractions - // -4/4 - 2/2 + 0 + 2*2 + 2*4/2 + 0 + 0 + 3 = -1-1+4+4+3 = 9 let result = ThermoPolynomial { a: vec![4.0, 2.0, 0.0, 2.0, 2.0, 0.0, 0.0, 0.0, 3.0], temp_range: (0.0, 0.0), } .s_over_r(2.0); - assert_delta!(result, 9.0, 1e-10); + assert_delta!(result, 9.5, 1e-10); } #[test] diff --git a/src/properties/transport_fit.rs b/src/properties/transport_fit.rs index 60a7213..d89887d 100644 --- a/src/properties/transport_fit.rs +++ b/src/properties/transport_fit.rs @@ -1,7 +1,7 @@ pub struct SpeciesTransportData { pub name: String, - viscosities: Vec, - conductivities: Vec, + pub(crate) viscosities: Vec, + pub(crate) conductivities: Vec, } impl SpeciesTransportData { From 063ec157956ff01f25825d91e4c160a8906ba932 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sun, 29 Mar 2026 22:31:09 -0400 Subject: [PATCH 7/8] Prevent polynomial_at from panicking --- src/properties/thermo_db.rs | 76 +++++++++++++++++++++++++++--------- src/properties/thermo_fit.rs | 30 +++++++++----- 2 files changed, 78 insertions(+), 28 deletions(-) diff --git a/src/properties/thermo_db.rs b/src/properties/thermo_db.rs index 5eecb2d..7db0c05 100644 --- a/src/properties/thermo_db.rs +++ b/src/properties/thermo_db.rs @@ -295,9 +295,17 @@ mod test { -1.742171366e+01, ]; - assert_vec_delta!(species.polynomial_at(650.0).a, real_coeff_1, 1e-9); - assert_delta!(species.polynomial_at(650.0).temp_range.0, 300.000, 1e-3); - assert_delta!(species.polynomial_at(650.0).temp_range.1, 1000.000, 1e-3); + assert_vec_delta!(species.polynomial_at(650.0).unwrap().a, real_coeff_1, 1e-9); + assert_delta!( + species.polynomial_at(650.0).unwrap().temp_range.0, + 300.000, + 1e-3 + ); + assert_delta!( + species.polynomial_at(650.0).unwrap().temp_range.1, + 1000.000, + 1e-3 + ); let real_coeff_2 = [ -3.523782900e+05, @@ -311,9 +319,17 @@ mod test { -2.695610360e+00, ]; - assert_vec_delta!(species.polynomial_at(3500.0).a, real_coeff_2, 1e-9); - assert_delta!(species.polynomial_at(3500.0).temp_range.0, 1000.000, 1e-3); - assert_delta!(species.polynomial_at(3500.0).temp_range.1, 6000.000, 1e-3); + assert_vec_delta!(species.polynomial_at(3500.0).unwrap().a, real_coeff_2, 1e-9); + assert_delta!( + species.polynomial_at(3500.0).unwrap().temp_range.0, + 1000.000, + 1e-3 + ); + assert_delta!( + species.polynomial_at(3500.0).unwrap().temp_range.1, + 6000.000, + 1e-3 + ); } #[test] @@ -366,7 +382,7 @@ END REACTANTS assert_eq!(alcl3.num_polynomials(), 2); assert_vec_delta!( - alcl3.polynomial_at(650.0).a, + alcl3.polynomial_at(650.0).unwrap().a, [ 7.750600970e+04, -1.440779717e+03, @@ -380,11 +396,19 @@ END REACTANTS ], 1e-9 ); - assert_delta!(alcl3.polynomial_at(650.0).temp_range.0, 300.0, 1e-3); - assert_delta!(alcl3.polynomial_at(650.0).temp_range.1, 1000.0, 1e-3); + assert_delta!( + alcl3.polynomial_at(650.0).unwrap().temp_range.0, + 300.0, + 1e-3 + ); + assert_delta!( + alcl3.polynomial_at(650.0).unwrap().temp_range.1, + 1000.0, + 1e-3 + ); assert_vec_delta!( - alcl3.polynomial_at(3500.0).a, + alcl3.polynomial_at(3500.0).unwrap().a, [ -1.378630916e+05, -5.579207290e+01, @@ -398,8 +422,16 @@ END REACTANTS ], 1e-9 ); - assert_delta!(alcl3.polynomial_at(3500.0).temp_range.0, 1000.0, 1e-3); - assert_delta!(alcl3.polynomial_at(3500.0).temp_range.1, 6000.0, 1e-3); + assert_delta!( + alcl3.polynomial_at(3500.0).unwrap().temp_range.0, + 1000.0, + 1e-3 + ); + assert_delta!( + alcl3.polynomial_at(3500.0).unwrap().temp_range.1, + 6000.0, + 1e-3 + ); // --- Air (reactant 0) --- let air = &thermo_db.reactants[0]; @@ -419,7 +451,7 @@ END REACTANTS assert_eq!(air.num_polynomials(), 2); assert_vec_delta!( - air.polynomial_at(650.0).a, + air.polynomial_at(650.0).unwrap().a, [ 1.009950160e+04, -1.968275610e+02, @@ -433,11 +465,11 @@ END REACTANTS ], 1e-9 ); - assert_delta!(air.polynomial_at(650.0).temp_range.0, 300.0, 1e-3); - assert_delta!(air.polynomial_at(650.0).temp_range.1, 1000.0, 1e-3); + assert_delta!(air.polynomial_at(650.0).unwrap().temp_range.0, 300.0, 1e-3); + assert_delta!(air.polynomial_at(650.0).unwrap().temp_range.1, 1000.0, 1e-3); assert_vec_delta!( - air.polynomial_at(3500.0).a, + air.polynomial_at(3500.0).unwrap().a, [ 2.415214430e+05, -1.257874600e+03, @@ -451,8 +483,16 @@ END REACTANTS ], 1e-9 ); - assert_delta!(air.polynomial_at(3500.0).temp_range.0, 1000.0, 1e-3); - assert_delta!(air.polynomial_at(3500.0).temp_range.1, 6000.0, 1e-3); + assert_delta!( + air.polynomial_at(3500.0).unwrap().temp_range.0, + 1000.0, + 1e-3 + ); + assert_delta!( + air.polynomial_at(3500.0).unwrap().temp_range.1, + 6000.0, + 1e-3 + ); // --- n-Butanol (reactant 1) --- let butanol = &thermo_db.reactants[1]; diff --git a/src/properties/thermo_fit.rs b/src/properties/thermo_fit.rs index 2d5bdf0..893bbcf 100644 --- a/src/properties/thermo_fit.rs +++ b/src/properties/thermo_fit.rs @@ -44,18 +44,22 @@ impl SpeciesThermoData { self.polynomials.len() } - pub fn polynomial_at(&self, temp: f64) -> &ThermoPolynomial { + pub fn polynomial_at(&self, temp: f64) -> Option<&ThermoPolynomial> { //TODO: Not the most efficient. Can refactor to pre-compute tables //and do 1-d linear interpolation if needed // //TODO: I Think condensed species need to be treated differently. Verify how that works in //the paper. + if self.polynomials.is_empty() { + return None; + } + let i_polynomial = self .polynomials .iter() .rposition(|polynomial| temp > polynomial.temp_range.0) .unwrap_or(0); - &self.polynomials[i_polynomial] + Some(&self.polynomials[i_polynomial]) } } @@ -113,8 +117,8 @@ mod test { a: vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0], temp_range: (0.0, 0.0), } - .cp_over_r(1.0); - assert_delta!(result, 28.0, 1e-10); + .cp_over_r(100.0); + assert_delta!(result, 706050403.0201, 1e-4); let result = ThermoPolynomial { a: vec![4.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0], @@ -180,22 +184,28 @@ mod test { h_formation: 1.0, }; - assert!(std::ptr::eq(data.polynomial_at(0.5), &data.polynomials[0])); - assert!(std::ptr::eq(data.polynomial_at(50.0), &data.polynomials[0])); assert!(std::ptr::eq( - data.polynomial_at(100.0), + data.polynomial_at(0.5).unwrap(), &data.polynomials[0] )); assert!(std::ptr::eq( - data.polynomial_at(100.0 + 1e-12), + data.polynomial_at(50.0).unwrap(), + &data.polynomials[0] + )); + assert!(std::ptr::eq( + data.polynomial_at(100.0).unwrap(), + &data.polynomials[0] + )); + assert!(std::ptr::eq( + data.polynomial_at(100.0 + 1e-12).unwrap(), &data.polynomials[1] )); assert!(std::ptr::eq( - data.polynomial_at(150.0 + 1e-12), + data.polynomial_at(150.0 + 1e-12).unwrap(), &data.polynomials[1] )); assert!(std::ptr::eq( - data.polynomial_at(500.0 + 1e-12), + data.polynomial_at(500.0 + 1e-12).unwrap(), &data.polynomials[1] )); } From 828b87dcc9361b18f3410f88bb71b5289d72d517 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sun, 29 Mar 2026 22:36:21 -0400 Subject: [PATCH 8/8] Additional code review comments --- src/properties/transport_fit.rs | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/properties/transport_fit.rs b/src/properties/transport_fit.rs index d89887d..6145a2d 100644 --- a/src/properties/transport_fit.rs +++ b/src/properties/transport_fit.rs @@ -17,22 +17,28 @@ impl SpeciesTransportData { } } - pub fn viscosity_at(&self, temp: f64) -> f64 { + pub fn viscosity_at(&self, temp: f64) -> Option { + if self.viscosities.is_empty() { + return None; + } let i_viscosity = self .viscosities .iter() .rposition(|viscosity| temp > viscosity.temp_range.0) .unwrap_or(0); - self.viscosities[i_viscosity].compute(temp) + Some(self.viscosities[i_viscosity].compute(temp)) } - pub fn conductivity_at(&self, temp: f64) -> f64 { + pub fn conductivity_at(&self, temp: f64) -> Option { + if self.conductivities.is_empty() { + return None; + } let i_conductivity = self .conductivities .iter() .rposition(|conductivity| temp > conductivity.temp_range.0) .unwrap_or(0); - self.conductivities[i_conductivity].compute(temp) + Some(self.conductivities[i_conductivity].compute(temp)) } } @@ -68,7 +74,7 @@ mod test { d: 40.0, }; - assert_delta!(fit.compute(1.0), 90.0, 1e-12); + assert_delta!(fit.compute(4.0), 60.73794361119891, 1e-12); } #[test]