Merge pull request #2 from aselimov/thermo-calcs

Thermo calcs

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,
+        &name,
-        polynomials,
         elements,
         phase,
+        polynomials,
         molecular_weight,
         h_formation,
-    })
+    ))
 }
 
 fn parse_polynomials_block<'a>(
@@ -295,9 +295,17 @@ mod test {
             -1.742171366e+01,
         ];
 
-        assert_vec_delta!(species.polynomials[0].a, real_coeff_1, 1e-9);
+        assert_vec_delta!(species.polynomial_at(650.0).unwrap().a, real_coeff_1, 1e-9);
-        assert_delta!(species.polynomials[0].temp_range.0, 300.000, 1e-3);
+        assert_delta!(
-        assert_delta!(species.polynomials[0].temp_range.1, 1000.000, 1e-3);
+            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.polynomials[1].a, real_coeff_2, 1e-9);
+        assert_vec_delta!(species.polynomial_at(3500.0).unwrap().a, real_coeff_2, 1e-9);
-        assert_delta!(species.polynomials[1].temp_range.0, 1000.000, 1e-3);
+        assert_delta!(
-        assert_delta!(species.polynomials[1].temp_range.1, 6000.000, 1e-3);
+            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]
@@ -363,10 +379,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).unwrap().a,
             [
                 7.750600970e+04,
                 -1.440779717e+03,
@@ -380,11 +396,19 @@ END REACTANTS
             ],
             1e-9
         );
-        assert_delta!(alcl3.polynomials[0].temp_range.0, 300.0, 1e-3);
+        assert_delta!(
-        assert_delta!(alcl3.polynomials[0].temp_range.1, 1000.0, 1e-3);
+            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.polynomials[1].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.polynomials[1].temp_range.0, 1000.0, 1e-3);
+        assert_delta!(
-        assert_delta!(alcl3.polynomials[1].temp_range.1, 6000.0, 1e-3);
+            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];
@@ -416,10 +448,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).unwrap().a,
             [
                 1.009950160e+04,
                 -1.968275610e+02,
@@ -433,11 +465,11 @@ END REACTANTS
             ],
             1e-9
         );
-        assert_delta!(air.polynomials[0].temp_range.0, 300.0, 1e-3);
+        assert_delta!(air.polynomial_at(650.0).unwrap().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).unwrap().temp_range.1, 1000.0, 1e-3);
 
         assert_vec_delta!(
-            air.polynomials[1].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.polynomials[1].temp_range.0, 1000.0, 1e-3);
+        assert_delta!(
-        assert_delta!(air.polynomials[1].temp_range.1, 6000.0, 1e-3);
+            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];
@@ -467,6 +507,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);
     }
 }

src/properties/thermo_fit.rs

@@ -7,7 +7,7 @@ pub struct SpeciesThermoData {
     pub name: String,
     pub elements: Vec<SpeciesElement>,
     pub phase: Phase,
-    pub polynomials: Vec<ThermoPolynomial>,
+    polynomials: Vec<ThermoPolynomial>,
     pub molecular_weight: f64,
     pub h_formation: f64,
 }
@@ -20,3 +20,193 @@ pub struct SpeciesElement {
     pub element: String,
     pub count: f64,
 }
+
+impl SpeciesThermoData {
+    pub fn new(
+        name: &str,
+        elements: Vec<SpeciesElement>,
+        phase: Phase,
+        polynomials: Vec<ThermoPolynomial>,
+        molecular_weight: f64,
+        h_formation: f64,
+    ) -> Self {
+        Self {
+            name: name.to_string(),
+            elements,
+            phase,
+            polynomials,
+            molecular_weight,
+            h_formation,
+        }
+    }
+
+    pub fn num_polynomials(&self) -> usize {
+        self.polynomials.len()
+    }
+
+    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);
+        Some(&self.polynomials[i_polynomial])
+    }
+}
+
+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 * 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 / 2.0 - 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 super::ThermoPolynomial;
+    use crate::{
+        assert_delta,
+        properties::thermo_fit::{Phase, SpeciesThermoData},
+    };
+
+    #[test]
+    fn test_cp_over_r() {
+        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(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],
+            temp_range: (0.0, 0.0),
+        }
+        .cp_over_r(2.0);
+        assert_delta!(result, 35.0, 1e-10);
+    }
+
+    #[test]
+    fn test_h_over_rt() {
+        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(2.0);
+        assert_delta!(result, 63.44314718055995, 1e-10);
+
+        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(100.0);
+        assert_delta!(result, 203.0196, 1e-4);
+    }
+
+    #[test]
+    fn test_s_over_r() {
+        let result = ThermoPolynomial {
+            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(100.0);
+        assert_delta!(result, 200691797.0572474, 1e-7);
+
+        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.5, 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).unwrap(),
+            &data.polynomials[0]
+        ));
+        assert!(std::ptr::eq(
+            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).unwrap(),
+            &data.polynomials[1]
+        ));
+        assert!(std::ptr::eq(
+            data.polynomial_at(500.0 + 1e-12).unwrap(),
+            &data.polynomials[1]
+        ));
+    }
+}

src/properties/transport_db.rs

@@ -50,11 +50,11 @@ fn parse_species_transport_block<'a>(
             (ViscosityOrConductivity::Conductivity, fit) => conductivities.push(fit),
         }
     }
-    Ok(SpeciesTransportData {
+    Ok(SpeciesTransportData::new(
-        name,
+        &name,
         viscosities,
         conductivities,
-    })
+    ))
 }
 
 fn parse_species_header_line(line: &str) -> Result<(String, usize, usize), PropertiesError> {

src/properties/transport_fit.rs

@@ -1,7 +1,45 @@
 pub struct SpeciesTransportData {
     pub name: String,
-    pub viscosities: Vec<TransportFit>,
+    pub(crate) viscosities: Vec<TransportFit>,
-    pub conductivities: Vec<TransportFit>,
+    pub(crate) conductivities: Vec<TransportFit>,
+}
+
+impl SpeciesTransportData {
+    pub fn new(
+        name: &str,
+        viscosities: Vec<TransportFit>,
+        conductivities: Vec<TransportFit>,
+    ) -> Self {
+        SpeciesTransportData {
+            name: name.to_string(),
+            viscosities,
+            conductivities,
+        }
+    }
+
+    pub fn viscosity_at(&self, temp: f64) -> Option<f64> {
+        if self.viscosities.is_empty() {
+            return None;
+        }
+        let i_viscosity = self
+            .viscosities
+            .iter()
+            .rposition(|viscosity| temp > viscosity.temp_range.0)
+            .unwrap_or(0);
+        Some(self.viscosities[i_viscosity].compute(temp))
+    }
+
+    pub fn conductivity_at(&self, temp: f64) -> Option<f64> {
+        if self.conductivities.is_empty() {
+            return None;
+        }
+        let i_conductivity = self
+            .conductivities
+            .iter()
+            .rposition(|conductivity| temp > conductivity.temp_range.0)
+            .unwrap_or(0);
+        Some(self.conductivities[i_conductivity].compute(temp))
+    }
 }
 
 pub struct TransportFit {
@@ -11,3 +49,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(4.0), 60.73794361119891, 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
+        );
+    }
+}