From 93bd1f4d38cee3bbe31d780716595f4803bfcaac Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sat, 4 Apr 2026 15:46:27 -0400 Subject: [PATCH 1/4] Some test fixes and adding initial simple version of matrix data struct --- src/lib.rs | 1 + src/matrix.rs | 135 ++++++++++++++++++++++++++++++++ src/mixtures/gas_mixture.rs | 20 ++--- src/properties/transport_fit.rs | 28 +++---- 4 files changed, 161 insertions(+), 23 deletions(-) create mode 100644 src/matrix.rs diff --git a/src/lib.rs b/src/lib.rs index eee60c5..95cdfcf 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,4 +1,5 @@ pub mod consts; +pub mod matrix; pub mod mixtures; pub mod properties; diff --git a/src/matrix.rs b/src/matrix.rs new file mode 100644 index 0000000..0fbe4c9 --- /dev/null +++ b/src/matrix.rs @@ -0,0 +1,135 @@ +use std::fmt; +use std::{ + fmt::Display, + ops::{Add, Div, Mul, Sub}, +}; + +/// Numeric trait as an alias to make the generics a little bit cleaner +pub trait Numeric: + Add + + Sub + + Mul + + Div + + Sized + + Clone + + Display +{ +} + +/// Blanket implementation for numeric types +impl Numeric for T where + T: Add + + Sub + + Mul + + Div + + Sized + + Clone + + Display +{ +} + +// Basic Error handling for Matrix operations +#[derive(Debug)] +pub enum MatrixError { + IndexError(usize, usize, usize, usize), +} + +impl fmt::Display for MatrixError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + MatrixError::IndexError(i, j, m, n) => write!( + f, + "Error accessing index [{i},{j}] for matrix with dimensions [{}, {}]", + m, n + ), + } + } +} + +fn make_index_error(i: usize, j: usize, m: &Matrix) -> MatrixError { + MatrixError::IndexError(i, j, m.m, m.n) +} + +pub struct Matrix { + data: Vec, + m: usize, + n: usize, +} + +impl Matrix { + pub fn new(m: usize, n: usize, init_val: T) -> Self { + let data = vec![init_val; n * m]; + Matrix { data, m, n } + } + + fn index(&self, i: usize, j: usize) -> usize { + i * self.n + j + } + pub fn get(&self, i: usize, j: usize) -> Result<&T, MatrixError> { + self.data + .get(self.index(i, j)) + .ok_or(make_index_error(i, j, &self)) + } + + pub fn set(&mut self, i: usize, j: usize, x: T) { + let index = self.index(i, j); + self.data[index] = x; + } +} + +impl Display for Matrix { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + let msg = (0..self.m) + .flat_map(|i| { + let mut parts = vec!["|".to_string()]; + parts.extend((0..self.n).map(|j| -> String { + if let Ok(val) = self.get(i, j) { + format!("{}", val) + } else { + "err".to_string() + } + })); + parts.push("|\n".to_string()); + parts + }) + .collect::>() + .join(" "); + write!(f, "\n{msg}") + } +} + +#[cfg(test)] +mod test { + use crate::{assert_delta, matrix::Matrix}; + fn gen_test_matrix() -> Matrix { + let mut m = Matrix::new(3, 4, 0.0); + m.set(0, 0, 1.0); + m.set(1, 1, 1.0); + m.set(2, 2, 1.0); + m.set(0, 3, 2.0); + m.set(1, 3, 2.0); + m.set(2, 3, 2.0); + m + } + + #[test] + fn test_matrix_basics() { + let m = gen_test_matrix(); + + // Validate that the matrix type returns the right values + assert_delta!(m.get(0, 0).unwrap(), 1.0, 1e-12); + assert_delta!(m.get(1, 1).unwrap(), 1.0, 1e-12); + assert_delta!(m.get(2, 2).unwrap(), 1.0, 1e-12); + assert_delta!(m.get(0, 3).unwrap(), 2.0, 1e-12); + assert_delta!(m.get(1, 3).unwrap(), 2.0, 1e-12); + assert_delta!(m.get(2, 3).unwrap(), 2.0, 1e-12); + + // Validate that the memory is laid out as expected + assert_delta!(m.data[0], 1.0, 1e-12); + assert_delta!(m.data[3], 2.0, 1e-12); + assert_delta!(m.data[5], 1.0, 1e-12); + assert_delta!(m.data[7], 2.0, 1e-12); + assert_delta!(m.data[10], 1.0, 1e-12); + assert_delta!(m.data[11], 2.0, 1e-12); + } +} diff --git a/src/mixtures/gas_mixture.rs b/src/mixtures/gas_mixture.rs index cc09169..651bef5 100644 --- a/src/mixtures/gas_mixture.rs +++ b/src/mixtures/gas_mixture.rs @@ -8,7 +8,7 @@ use crate::{ pub struct GasMixture { pub(crate) ns: Vec, - pub(crate) nsum: f64 + pub(crate) nsum: f64, pub(crate) species: Vec, pub(crate) transport_data: Vec, } @@ -27,14 +27,16 @@ impl GasMixture { 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() + (n/self.nsum).ln() + p.h_over_rt(temp) - p.s_over_r(temp) + + (pressure / P_REF).ln() + + (n / self.nsum).ln() } Phase::Condensed => todo!(), } }) .collect() } - + // Calculate the normalized entropy (S/R) for each mixture component // // Equations 2.17 from reference paper @@ -48,19 +50,19 @@ impl GasMixture { 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() + p.s_over_r(temp) - (n / self.nsum).ln() - (pressure / P_REF).ln() } Phase::Condensed => todo!(), } }) .collect() } - + // Calculate the normalized mixture enthalpy (H/RT) // 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{ - self.ns + pub fn mixture_h_over_rt(&self, temp: f64) -> Vec { + self.ns .iter() .zip(self.species.iter()) .map(|(n, s)| -> f64 { @@ -69,11 +71,11 @@ impl GasMixture { let p = s .polynomial_at(temp) .expect("Gas doesn't have a polynomial"); - n*p.h_over_rt(temp) + n * p.h_over_rt(temp) } Phase::Condensed => todo!(), } }) - .collect() + .collect() } } diff --git a/src/properties/transport_fit.rs b/src/properties/transport_fit.rs index 6145a2d..2e5229d 100644 --- a/src/properties/transport_fit.rs +++ b/src/properties/transport_fit.rs @@ -120,75 +120,75 @@ mod test { }; assert_delta!( - data.conductivity_at(0.5), + data.conductivity_at(0.5).unwrap(), data.conductivities[0].compute(0.5), 1e-12 ); assert_delta!( - data.conductivity_at(1.0), + data.conductivity_at(1.0).unwrap(), data.conductivities[0].compute(1.0), 1e-12 ); assert_delta!( - data.conductivity_at(50.0), + data.conductivity_at(50.0).unwrap(), data.conductivities[0].compute(50.0), 1e-12 ); assert_delta!( - data.conductivity_at(100.0), + data.conductivity_at(100.0).unwrap(), data.conductivities[0].compute(100.0), 1e-12 ); assert_delta!( - data.conductivity_at(100.0 + 1e-12), + data.conductivity_at(100.0 + 1e-12).unwrap(), data.conductivities[1].compute(100.0 + 1e-12), 1e-12 ); assert_delta!( - data.conductivity_at(200.0), + data.conductivity_at(200.0).unwrap(), data.conductivities[1].compute(200.0), 1e-12 ); assert_delta!( - data.conductivity_at(500.0), + data.conductivity_at(500.0).unwrap(), data.conductivities[1].compute(500.0), 1e-12 ); assert_delta!( - data.viscosity_at(0.5), + data.viscosity_at(0.5).unwrap(), data.viscosities[0].compute(0.5), 1e-12 ); assert_delta!( - data.viscosity_at(1.0), + data.viscosity_at(1.0).unwrap(), data.viscosities[0].compute(1.0), 1e-12 ); assert_delta!( - data.viscosity_at(50.0), + data.viscosity_at(50.0).unwrap(), data.viscosities[0].compute(50.0), 1e-12 ); assert_delta!( - data.viscosity_at(100.0), + data.viscosity_at(100.0).unwrap(), data.viscosities[0].compute(100.0), 1e-12 ); assert_delta!( - data.viscosity_at(100.0 + 1e-12), + data.viscosity_at(100.0 + 1e-12).unwrap(), data.viscosities[1].compute(100.0 + 1e-12), 1e-12 ); assert_delta!( - data.viscosity_at(200.0), + data.viscosity_at(200.0).unwrap(), data.viscosities[1].compute(200.0), 1e-12 ); assert_delta!( - data.viscosity_at(500.0), + data.viscosity_at(500.0).unwrap(), data.viscosities[1].compute(500.0), 1e-12 ); From df0e93735bce8eb94ed6ab4ca055b46c6ecd3568 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sat, 4 Apr 2026 21:50:06 -0400 Subject: [PATCH 2/4] Add basic matrix math --- src/matrix.rs | 102 ++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 99 insertions(+), 3 deletions(-) diff --git a/src/matrix.rs b/src/matrix.rs index 0fbe4c9..b196e3b 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -11,8 +11,9 @@ pub trait Numeric: + Mul + Div + Sized - + Clone + + Copy + Display + + Default { } @@ -25,6 +26,8 @@ impl Numeric for T where + Sized + Clone + Display + + Default + + Copy { } @@ -32,6 +35,9 @@ impl Numeric for T where #[derive(Debug)] pub enum MatrixError { IndexError(usize, usize, usize, usize), + AddError(usize, usize, usize, usize), + MultiplicationError(usize, usize, usize, usize), + FromVecError(usize, usize, usize), } impl fmt::Display for MatrixError { @@ -42,6 +48,19 @@ impl fmt::Display for MatrixError { "Error accessing index [{i},{j}] for matrix with dimensions [{}, {}]", m, n ), + + MatrixError::MultiplicationError(i, j, m, n) => write!( + f, + "Matrices with dimensions [{i},{j}] and [{m},{n}] cannot be multiplied", + ), + MatrixError::AddError(i, j, m, n) => write!( + f, + "Matrices with dimensions [{i},{j}] and [{m},{n}] cannot be added", + ), + MatrixError::FromVecError(i, j, len) => write!( + f, + "Matrices with dimensions [{i},{j}] cannot be created from vec with len={len}", + ), } } } @@ -62,13 +81,58 @@ impl Matrix { Matrix { data, m, n } } + pub fn from_vec(m: usize, n: usize, data: Vec) -> Result { + if m * n != data.len() { + return Err(MatrixError::FromVecError(m, n, data.len())); + } + Ok(Self { m, n, data }) + } + fn index(&self, i: usize, j: usize) -> usize { i * self.n + j } pub fn get(&self, i: usize, j: usize) -> Result<&T, MatrixError> { self.data .get(self.index(i, j)) - .ok_or(make_index_error(i, j, &self)) + .ok_or(make_index_error(i, j, self)) + } + + pub fn add(&self, other: &Matrix) -> Result, MatrixError> { + // Compatibility check + if self.m != other.m || self.n != other.n { + return Err(MatrixError::AddError(self.m, self.n, other.m, other.n)); + } + + let data = self + .data + .iter() + .zip(other.data.iter()) + .map(|(a, b)| *a + *b) + .collect(); + Self::from_vec(self.m, self.n, data) + } + + pub fn mul(&self, other: &Matrix) -> Result, MatrixError> { + let mut c = Matrix::new(self.m, other.n, T::default()); + // Compatibility check + if self.n != other.m { + return Err(MatrixError::MultiplicationError( + self.m, self.n, other.m, other.n, + )); + } + + for i in 0..self.m { + for k in 0..self.n { + for j in 0..other.n { + c.set( + i, + j, + *c.get(i, j)? + (*self.get(i, k)?) * (*other.get(k, j)?), + ); + } + } + } + Ok(c) } pub fn set(&mut self, i: usize, j: usize, x: T) { @@ -94,7 +158,7 @@ impl Display for Matrix { }) .collect::>() .join(" "); - write!(f, "\n{msg}") + write!(f, "\n {msg}") } } @@ -131,5 +195,37 @@ mod test { assert_delta!(m.data[7], 2.0, 1e-12); assert_delta!(m.data[10], 1.0, 1e-12); assert_delta!(m.data[11], 2.0, 1e-12); + + // Test from_vec + let data = vec![1, 2, 3, 4, 5, 6]; + let m = Matrix::from_vec(3, 2, data).unwrap(); + assert_eq!(*m.get(0, 0).unwrap(), 1); + assert_eq!(*m.get(0, 1).unwrap(), 2); + assert_eq!(*m.get(1, 0).unwrap(), 3); + assert_eq!(*m.get(1, 1).unwrap(), 4); + assert_eq!(*m.get(2, 0).unwrap(), 5); + assert_eq!(*m.get(2, 1).unwrap(), 6); + } + + #[test] + fn test_matrix_math() { + let a = Matrix::from_vec(3, 2, vec![1, 2, 3, 4, 5, 6]).unwrap(); + let b = Matrix::from_vec(3, 2, vec![0, 1, 0, 1, 0, 1]).unwrap(); + + let c = a.add(&b).unwrap(); + assert_eq!(*c.get(0, 0).unwrap(), 1); + assert_eq!(*c.get(0, 1).unwrap(), 3); + assert_eq!(*c.get(1, 0).unwrap(), 3); + assert_eq!(*c.get(1, 1).unwrap(), 5); + assert_eq!(*c.get(2, 0).unwrap(), 5); + assert_eq!(*c.get(2, 1).unwrap(), 7); + + assert!(a.mul(&b).is_err()); + + let a = Matrix::from_vec(3, 2, vec![1, 2, 3, 4, 5, 6]).unwrap(); + let b = Matrix::from_vec(2, 3, vec![0, 1, 0, 1, 0, 1]).unwrap(); + let c = a.mul(&b).unwrap(); + + assert_eq!(c.data, vec![2, 1, 2, 4, 3, 4, 6, 5, 6]); } } From dccd2885a2aec262b244ae32999a3107454afd4e Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sat, 4 Apr 2026 22:25:31 -0400 Subject: [PATCH 3/4] Add transpose function --- src/matrix.rs | 37 ++++++++++++++++++++++++++++++++++--- 1 file changed, 34 insertions(+), 3 deletions(-) diff --git a/src/matrix.rs b/src/matrix.rs index b196e3b..8d39ff5 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -92,11 +92,19 @@ impl Matrix { i * self.n + j } pub fn get(&self, i: usize, j: usize) -> Result<&T, MatrixError> { + if i >= self.m || j >= self.n { + return Err(make_index_error(i, j, self)); + } self.data .get(self.index(i, j)) .ok_or(make_index_error(i, j, self)) } + pub fn set(&mut self, i: usize, j: usize, x: T) { + let index = self.index(i, j); + self.data[index] = x; + } + pub fn add(&self, other: &Matrix) -> Result, MatrixError> { // Compatibility check if self.m != other.m || self.n != other.n { @@ -135,9 +143,15 @@ impl Matrix { Ok(c) } - pub fn set(&mut self, i: usize, j: usize, x: T) { - let index = self.index(i, j); - self.data[index] = x; + pub fn transpose(&self) -> Result { + let mut c = Self::new(self.n, self.m, T::default()); + + for i in 0..self.m { + for j in 0..self.n { + c.set(j, i, *self.get(i, j)?); + } + } + Ok(c) } } @@ -227,5 +241,22 @@ mod test { let c = a.mul(&b).unwrap(); assert_eq!(c.data, vec![2, 1, 2, 4, 3, 4, 6, 5, 6]); + + let a = Matrix::from_vec(3, 4, vec![1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]).unwrap(); + let atranspose = a.transpose().unwrap(); + + assert_eq!(*atranspose.get(0, 0).unwrap(), 1); + assert_eq!(*atranspose.get(0, 1).unwrap(), 5); + assert_eq!(*atranspose.get(0, 2).unwrap(), 9); + assert!(atranspose.get(0, 3).is_err()); + assert_eq!(*atranspose.get(1, 0).unwrap(), 2); + assert_eq!(*atranspose.get(1, 1).unwrap(), 6); + assert_eq!(*atranspose.get(1, 2).unwrap(), 10); + assert_eq!(*atranspose.get(2, 0).unwrap(), 3); + assert_eq!(*atranspose.get(2, 1).unwrap(), 7); + assert_eq!(*atranspose.get(2, 2).unwrap(), 11); + assert_eq!(*atranspose.get(3, 0).unwrap(), 4); + assert_eq!(*atranspose.get(3, 1).unwrap(), 8); + assert_eq!(*atranspose.get(3, 2).unwrap(), 12); } } From 6fa468ca6cb7a4c493bdaf717df7dccc2e17bd30 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Tue, 7 Apr 2026 22:16:38 -0400 Subject: [PATCH 4/4] Add base error struct as well as swap_rows function in matrix --- src/error.rs | 48 ++++++++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 2 ++ src/matrix.rs | 42 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 92 insertions(+) create mode 100644 src/error.rs diff --git a/src/error.rs b/src/error.rs new file mode 100644 index 0000000..b8b6439 --- /dev/null +++ b/src/error.rs @@ -0,0 +1,48 @@ +use std::fmt; + +use crate::{matrix::MatrixError, properties::PropertiesError, solvers::SolverError}; + +#[derive(Debug)] +pub enum CEAError { + Matrix(MatrixError), + Solver(SolverError), + Properties(PropertiesError), +} + +impl fmt::Display for CEAError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + CEAError::Matrix(e) => write!(f, "matrix error: {}", e), + CEAError::Solver(e) => write!(f, "solver error: {}", e), + CEAError::Properties(e) => write!(f, "properties error: {}", e), + } + } +} + +impl std::error::Error for CEAError { + fn source(&self) -> Option<&(dyn std::error::Error + 'static)> { + match self { + CEAError::Matrix(e) => Some(e), + CEAError::Solver(e) => Some(e), + CEAError::Properties(e) => Some(e), + } + } +} + +impl From for CEAError { + fn from(e: MatrixError) -> Self { + CEAError::Matrix(e) + } +} + +impl From for CEAError { + fn from(e: SolverError) -> Self { + CEAError::Solver(e) + } +} + +impl From for CEAError { + fn from(e: PropertiesError) -> Self { + CEAError::Properties(e) + } +} diff --git a/src/lib.rs b/src/lib.rs index 95cdfcf..3c96114 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,7 +1,9 @@ pub mod consts; +pub mod error; pub mod matrix; pub mod mixtures; pub mod properties; +pub mod solvers; #[macro_export] macro_rules! assert_delta { diff --git a/src/matrix.rs b/src/matrix.rs index 8d39ff5..5763240 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -65,10 +65,13 @@ impl fmt::Display for MatrixError { } } +impl std::error::Error for MatrixError {} + fn make_index_error(i: usize, j: usize, m: &Matrix) -> MatrixError { MatrixError::IndexError(i, j, m.m, m.n) } +#[derive(Clone)] pub struct Matrix { data: Vec, m: usize, @@ -153,6 +156,21 @@ impl Matrix { } Ok(c) } + + pub fn shape(&self) -> (usize, usize) { + (self.m, self.n) + } + + pub fn swap_rows(&mut self, i: usize, j: usize) { + if i == j { + return; + } + for col in 0..self.n { + let first_index = self.index(i, col); + let second_index = self.index(j, col); + self.data.swap(first_index, second_index); + } + } } impl Display for Matrix { @@ -259,4 +277,28 @@ mod test { assert_eq!(*atranspose.get(3, 1).unwrap(), 8); assert_eq!(*atranspose.get(3, 2).unwrap(), 12); } + + #[test] + fn test_swap_rows() { + let mut m = Matrix::from_vec(3, 3, vec![1, 2, 3, 4, 5, 6, 7, 8, 9]).unwrap(); + + // Swap rows 0 and 2 + m.swap_rows(0, 2); + assert_eq!(*m.get(0, 0).unwrap(), 7); + assert_eq!(*m.get(0, 1).unwrap(), 8); + assert_eq!(*m.get(0, 2).unwrap(), 9); + assert_eq!(*m.get(2, 0).unwrap(), 1); + assert_eq!(*m.get(2, 1).unwrap(), 2); + assert_eq!(*m.get(2, 2).unwrap(), 3); + // Row 1 unchanged + assert_eq!(*m.get(1, 0).unwrap(), 4); + assert_eq!(*m.get(1, 1).unwrap(), 5); + assert_eq!(*m.get(1, 2).unwrap(), 6); + + // Swapping a row with itself is a no-op + m.swap_rows(1, 1); + assert_eq!(*m.get(1, 0).unwrap(), 4); + assert_eq!(*m.get(1, 1).unwrap(), 5); + assert_eq!(*m.get(1, 2).unwrap(), 6); + } }