Another gasmixture refactor to make the code less fortrany

src/mixtures/gas_mixture.rs

@@ -6,11 +6,18 @@ 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) species: Vec<SpeciesThermoData>,
+    pub(crate) gasses: Vec<MixtureComponent>,
-    pub(crate) coeffs: Matrix<f64>,
+    pub(crate) condensed: Vec<MixtureComponent>,
     pub(crate) elements: HashMap<String, usize>,
     pub(crate) binitial: Vec<f64>,
 }
@@ -40,23 +47,30 @@ 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,
-            species: species.to_vec(),
+            gasses,
-            coeffs,
+            condensed,
             elements,
             binitial,
         }
@@ -64,18 +78,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.ns
+        self.gasses
             .iter()
-            .zip(self.species.iter())
+            .chain(self.condensed.iter())
-            .map(|(n, s)| -> f64 {
+            .map(|c| -> f64 {
-                match s.phase {
+                match c.s.phase {
                     Phase::Gas => {
-                        let p = s
+                        let p =
-                            .polynomial_at(temp)
+                            c.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()
-                            + (n / self.nsum).ln()
+                            + (c.n / self.nsum).ln()
                     }
                     Phase::Condensed => todo!(),
                 }
@@ -87,16 +101,16 @@ impl GasMixture {
     //
     // Equations 2.17 from reference paper
     pub fn gas_entropies_over_rt(&self, temp: f64, pressure: f64) -> Vec<f64> {
-        self.ns
+        self.gasses
             .iter()
-            .zip(self.species.iter())
+            .chain(self.condensed.iter())
-            .map(|(n, s)| -> f64 {
+            .map(|c| -> f64 {
-                match s.phase {
+                match c.s.phase {
                     Phase::Gas => {
-                        let p = s
+                        let p =
-                            .polynomial_at(temp)
+                            c.s.polynomial_at(temp)
                                 .expect("Gas doesn't have a polynomial");
-                        p.s_over_r(temp) - (n / self.nsum).ln() - (pressure / P_REF).ln()
+                        p.s_over_r(temp) - (c.n / self.nsum).ln() - (pressure / P_REF).ln()
                     }
                     Phase::Condensed => todo!(),
                 }
@@ -108,16 +122,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.ns
+        self.gasses
             .iter()
-            .zip(self.species.iter())
+            .chain(self.condensed.iter())
-            .map(|(n, s)| -> f64 {
+            .map(|c| -> f64 {
-                match s.phase {
+                match c.s.phase {
                     Phase::Gas => {
-                        let p = s
+                        let p =
-                            .polynomial_at(temp)
+                            c.s.polynomial_at(temp)
                                 .expect("Gas doesn't have a polynomial");
-                        n * p.h_over_rt(temp)
+                        c.n * p.h_over_rt(temp)
                     }
                     Phase::Condensed => todo!(),
                 }
@@ -163,15 +177,17 @@ mod test {
             0.01665842352342649,
             0.02498763528513974,
         ];
-        assert_vec_delta!(expected_ns, gas.ns, 1e-12);
+        let ns: Vec<f64> = gas.gasses.iter().map(|c| c.n).collect();
+
+        assert_vec_delta!(expected_ns, 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.gasses[0].a[0], 2.0, 1e-12);
-        assert_delta!(gas.coeffs.get(0, 1).unwrap(), 0.0, 1e-12);
+        assert_delta!(gas.gasses[1].a[0], 0.0, 1e-12);
-        assert_delta!(gas.coeffs.get(0, 2).unwrap(), 2.0, 1e-12);
+        assert_delta!(gas.gasses[2].a[0], 2.0, 1e-12);
-        assert_delta!(gas.coeffs.get(1, 0).unwrap(), 0.0, 1e-12);
+        assert_delta!(gas.gasses[0].a[1], 0.0, 1e-12);
-        assert_delta!(gas.coeffs.get(1, 1).unwrap(), 2.0, 1e-12);
+        assert_delta!(gas.gasses[1].a[1], 2.0, 1e-12);
-        assert_delta!(gas.coeffs.get(1, 2).unwrap(), 1.0, 1e-12);
+        assert_delta!(gas.gasses[2].a[1], 1.0, 1e-12);
 
         let expected_b = [0.06663369409370597, 0.05830448233199272];
         assert_vec_delta!(gas.binitial, expected_b, 1e-12);