src/mixtures/gas_mixture.rs

@@ -6,18 +6,11 @@ use crate::{
     properties::thermo_fit::{Phase, SpeciesThermoData},
 };
 
-pub struct MixtureComponent {
-    //kg-moles component/kg_mixture
-    n: f64,
-    // Coefficients
-    a: Vec<f64>,
-    s: SpeciesThermoData,
-}
-
 pub struct GasMixture {
+    pub(crate) ns: Vec<f64>,
     pub(crate) nsum: f64,
-    pub(crate) gasses: Vec<MixtureComponent>,
-    pub(crate) condensed: Vec<MixtureComponent>,
+    pub(crate) species: Vec<SpeciesThermoData>,
+    pub(crate) coeffs: Matrix<f64>,
     pub(crate) elements: HashMap<String, usize>,
     pub(crate) binitial: Vec<f64>,
 }
@@ -47,30 +40,23 @@ impl GasMixture {
             }
         }
 
+        // 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);
+            });
+        }
+
         let binitial = get_b_current(&elements, species, &ns);
 
-        // Now separate SpeciesThermoData in gas and condensed MixtureComponents
-        let (gasses, condensed) = ns
-            .iter()
-            .zip(species.iter())
-            .map(|(n, s)| {
-                let mut a = vec![0.0; elements.len()];
-                s.elements.iter().for_each(|e| {
-                    let i = elements.get(&e.element).unwrap();
-                    a[*i] += e.count;
-                });
-                MixtureComponent {
-                    n: *n,
-                    a,
-                    s: s.clone(),
-                }
-            })
-            .partition(|c| matches!(c.s.phase, Phase::Gas));
-
         GasMixture {
+            ns,
             nsum,
-            gasses,
-            condensed,
+            species: species.to_vec(),
+            coeffs,
             elements,
             binitial,
         }
@@ -78,18 +64,18 @@ impl GasMixture {
     // 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> {
-        self.gasses
+        self.ns
             .iter()
-            .chain(self.condensed.iter())
-            .map(|c| -> f64 {
-                match c.s.phase {
+            .zip(self.species.iter())
+            .map(|(n, s)| -> f64 {
+                match s.phase {
                     Phase::Gas => {
-                        let p =
-                            c.s.polynomial_at(temp)
-                                .expect("Gas doesn't have a polynomial");
+                        let p = s
+                            .polynomial_at(temp)
+                            .expect("Gas doesn't have a polynomial");
                         p.h_over_rt(temp) - p.s_over_r(temp)
                             + (pressure / P_REF).ln()
-                            + (c.n / self.nsum).ln()
+                            + (n / self.nsum).ln()
                     }
                     Phase::Condensed => todo!(),
                 }
@@ -101,16 +87,16 @@ impl GasMixture {
     //
     // Equations 2.17 from reference paper
     pub fn gas_entropies_over_rt(&self, temp: f64, pressure: f64) -> Vec<f64> {
-        self.gasses
+        self.ns
             .iter()
-            .chain(self.condensed.iter())
-            .map(|c| -> f64 {
-                match c.s.phase {
+            .zip(self.species.iter())
+            .map(|(n, s)| -> f64 {
+                match s.phase {
                     Phase::Gas => {
-                        let p =
-                            c.s.polynomial_at(temp)
-                                .expect("Gas doesn't have a polynomial");
-                        p.s_over_r(temp) - (c.n / self.nsum).ln() - (pressure / P_REF).ln()
+                        let p = s
+                            .polynomial_at(temp)
+                            .expect("Gas doesn't have a polynomial");
+                        p.s_over_r(temp) - (n / self.nsum).ln() - (pressure / P_REF).ln()
                     }
                     Phase::Condensed => todo!(),
                 }
@@ -122,16 +108,16 @@ impl GasMixture {
     // Note that the enthalpy doesn't have a dependence on the pressure.
     // Equation 2.14 from the paper
     pub fn mixture_h_over_rt(&self, temp: f64) -> Vec<f64> {
-        self.gasses
+        self.ns
             .iter()
-            .chain(self.condensed.iter())
-            .map(|c| -> f64 {
-                match c.s.phase {
+            .zip(self.species.iter())
+            .map(|(n, s)| -> f64 {
+                match s.phase {
                     Phase::Gas => {
-                        let p =
-                            c.s.polynomial_at(temp)
-                                .expect("Gas doesn't have a polynomial");
-                        c.n * p.h_over_rt(temp)
+                        let p = s
+                            .polynomial_at(temp)
+                            .expect("Gas doesn't have a polynomial");
+                        n * p.h_over_rt(temp)
                     }
                     Phase::Condensed => todo!(),
                 }
@@ -177,17 +163,15 @@ mod test {
             0.01665842352342649,
             0.02498763528513974,
         ];
-        let ns: Vec<f64> = gas.gasses.iter().map(|c| c.n).collect();
-
-        assert_vec_delta!(expected_ns, ns, 1e-12);
+        assert_vec_delta!(expected_ns, gas.ns, 1e-12);
         assert_delta!(gas.nsum, 0.04997527057027948, 1e-12);
 
-        assert_delta!(gas.gasses[0].a[0], 2.0, 1e-12);
-        assert_delta!(gas.gasses[1].a[0], 0.0, 1e-12);
-        assert_delta!(gas.gasses[2].a[0], 2.0, 1e-12);
-        assert_delta!(gas.gasses[0].a[1], 0.0, 1e-12);
-        assert_delta!(gas.gasses[1].a[1], 2.0, 1e-12);
-        assert_delta!(gas.gasses[2].a[1], 1.0, 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);
 
         let expected_b = [0.06663369409370597, 0.05830448233199272];
         assert_vec_delta!(gas.binitial, expected_b, 1e-12);