Update gas mixture with correct n calculation

src/mixtures/gas_mixture.rs

@@ -1,19 +1,62 @@
+use std::collections::HashMap;
+
 use crate::{
     consts::P_REF,
-    properties::{
+    matrix::Matrix,
-        thermo_fit::{Phase, SpeciesThermoData},
+    properties::thermo_fit::{Phase, SpeciesThermoData},
-        transport_fit::SpeciesTransportData,
-    },
 };
 
 pub struct GasMixture {
     pub(crate) ns: Vec<f64>,
     pub(crate) nsum: f64,
     pub(crate) species: Vec<SpeciesThermoData>,
-    pub(crate) transport_data: Vec<SpeciesTransportData>,
+    pub(crate) coeffs: Matrix<f64>,
+    pub(crate) elements: HashMap<String, usize>,
 }
 
 impl GasMixture {
+    pub fn new(ns_kg_moles: &[f64], species: &[SpeciesThermoData]) -> Self {
+        // First calculate the total and per species kg-moles/kg mixture
+        let (nsum_kg_moles, mass_sum) = ns_kg_moles
+            .iter()
+            .zip(species.iter())
+            .fold((0.0, 0.0), |acc, (n, s)| {
+                (acc.0 + n, acc.1 + n * s.molecular_weight)
+            });
+
+        let nsum = nsum_kg_moles / mass_sum;
+        let ns = ns_kg_moles.iter().map(|n| n / mass_sum).collect();
+
+        // Now build out the element list
+        let mut elements = HashMap::new();
+        let mut ei = 0;
+        for s in species.iter() {
+            for element in s.elements.iter() {
+                if !elements.contains_key(&element.element) {
+                    elements.insert(element.element.clone(), ei);
+                    ei += 1
+                }
+            }
+        }
+
+        // Now build the coefficients
+        let mut coeffs = Matrix::new(ei, species.len(), 0.0);
+        for (j, s) in species.iter().enumerate() {
+            s.elements.iter().for_each(|e| {
+                // Safe to unwrap because elements should always have e
+                let i = *elements.get(&e.element).unwrap();
+                coeffs.set(i, j, e.count);
+            });
+        }
+
+        GasMixture {
+            ns,
+            nsum,
+            species: species.to_vec(),
+            coeffs,
+            elements,
+        }
+    }
     // Calculate the normalized chemical potential (μ/RT) for each component in the mixture.
     // Equations 2.11 from reference paper
     pub fn gas_chem_potentials_over_rt(&self, temp: f64, pressure: f64) -> Vec<f64> {
@@ -78,3 +121,35 @@ impl GasMixture {
             .collect()
     }
 }
+
+#[cfg(test)]
+mod test {
+    use crate::{
+        assert_delta, assert_vec_delta, mixtures::gas_mixture::GasMixture,
+        properties::test_helpers::h2_o2_h2o_thermo_data,
+    };
+
+    #[test]
+    pub fn test_gas_mixture_creation() {
+        let species = h2_o2_h2o_thermo_data();
+
+        let ns = [1.0, 2.0, 3.0];
+
+        let gas = GasMixture::new(&ns, &species);
+
+        let expected_ns = [
+            0.008329211761713246,
+            0.01665842352342649,
+            0.02498763528513974,
+        ];
+        assert_vec_delta!(expected_ns, gas.ns, 1e-12);
+        assert_delta!(gas.nsum, 0.04997527057027948, 1e-12);
+
+        assert_delta!(gas.coeffs.get(0, 0).unwrap(), 2.0, 1e-12);
+        assert_delta!(gas.coeffs.get(0, 1).unwrap(), 0.0, 1e-12);
+        assert_delta!(gas.coeffs.get(0, 2).unwrap(), 2.0, 1e-12);
+        assert_delta!(gas.coeffs.get(1, 0).unwrap(), 0.0, 1e-12);
+        assert_delta!(gas.coeffs.get(1, 1).unwrap(), 2.0, 1e-12);
+        assert_delta!(gas.coeffs.get(1, 2).unwrap(), 1.0, 1e-12);
+    }
+}

src/properties/mod.rs

@@ -5,5 +5,7 @@ pub mod thermo_fit;
 pub mod transport_db;
 pub mod transport_fit;
 mod utils;
+#[cfg(test)]
+pub mod test_helpers;
 
 pub use error::PropertiesError;

src/properties/test_helpers.rs

@@ -0,0 +1,45 @@
+use super::thermo_fit::{Phase, SpeciesElement, SpeciesThermoData, ThermoPolynomial};
+
+pub fn simple_thermo_polynomial(
+    cp_over_r: f64,
+    h_b: f64,
+    s_b: f64,
+    temp_range: (f64, f64),
+) -> ThermoPolynomial {
+    ThermoPolynomial {
+        a: vec![0.0, 0.0, cp_over_r, 0.0, 0.0, 0.0, 0.0, h_b, s_b],
+        temp_range,
+    }
+}
+
+pub fn h2_o2_h2o_thermo_data() -> Vec<SpeciesThermoData> {
+    vec![
+        SpeciesThermoData::new(
+            "H2",
+            vec![SpeciesElement { element: "H".to_string(), count: 2.0 }],
+            Phase::Gas,
+            vec![simple_thermo_polynomial(3.5, -1012.0, 29.11, (200.0, 6000.0))],
+            2.01594,
+            0.0,
+        ),
+        SpeciesThermoData::new(
+            "O2",
+            vec![SpeciesElement { element: "O".to_string(), count: 2.0 }],
+            Phase::Gas,
+            vec![simple_thermo_polynomial(3.5, -1005.0, 30.03, (200.0, 6000.0))],
+            31.9988,
+            0.0,
+        ),
+        SpeciesThermoData::new(
+            "H2O",
+            vec![
+                SpeciesElement { element: "H".to_string(), count: 2.0 },
+                SpeciesElement { element: "O".to_string(), count: 1.0 },
+            ],
+            Phase::Gas,
+            vec![simple_thermo_polynomial(4.0, -29192.0, 23.03, (200.0, 6000.0))],
+            18.01528,
+            -241826.0,
+        ),
+    ]
+}

src/properties/thermo_fit.rs

@@ -1,8 +1,10 @@
+#[derive(Clone)]
 pub enum Phase {
     Gas,
     Condensed,
 }
 
+#[derive(Clone)]
 pub struct SpeciesThermoData {
     pub name: String,
     pub elements: Vec<SpeciesElement>,
@@ -11,11 +13,14 @@ pub struct SpeciesThermoData {
     pub molecular_weight: f64,
     pub h_formation: f64,
 }
+
+#[derive(Clone)]
 pub struct ThermoPolynomial {
     pub a: Vec<f64>,
     pub temp_range: (f64, f64),
 }
 
+#[derive(Clone)]
 pub struct SpeciesElement {
     pub element: String,
     pub count: f64,