Fix some calculation issues

This commit is contained in:
Alex Selimov 2026-03-29 22:26:59 -04:00
parent cdb74a4463
commit 6f4bc45927
3 changed files with 16 additions and 29 deletions

View file

@ -99,7 +99,7 @@ fn parse_species<'a>(
}
Ok(SpeciesThermoData::new(
name,
&name,
elements,
phase,
polynomials,

View file

@ -23,7 +23,7 @@ pub struct SpeciesElement {
impl SpeciesThermoData {
pub fn new(
name: String,
name: &str,
elements: Vec<SpeciesElement>,
phase: Phase,
polynomials: Vec<ThermoPolynomial>,
@ -31,7 +31,7 @@ impl SpeciesThermoData {
h_formation: f64,
) -> Self {
Self {
name,
name: name.to_string(),
elements,
phase,
polynomials,
@ -77,7 +77,7 @@ impl ThermoPolynomial {
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 * inv_temp.ln()
+ self.a[1] * inv_temp * temp.ln()
+ self.a[2]
+ self.a[3] * temp / 2.0
+ self.a[4] * temp * temp / 3.0
@ -89,7 +89,7 @@ impl ThermoPolynomial {
/// 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 - self.a[1] * inv_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
@ -103,14 +103,12 @@ impl ThermoPolynomial {
mod test {
use super::ThermoPolynomial;
use crate::{
assert_delta, assert_vec_delta,
assert_delta,
properties::thermo_fit::{Phase, SpeciesThermoData},
};
#[test]
fn test_cp_over_r() {
// At T=1: a[0]/T^2 + a[1]/T + a[2] + a[3]*T + a[4]*T^2 + a[5]*T^3 + a[6]*T^4
// = 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28
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),
@ -118,7 +116,6 @@ mod test {
.cp_over_r(1.0);
assert_delta!(result, 28.0, 1e-10);
// At T=2: 4/4 + 2/2 + 1 + 2*2 + 1*4 + 1*8 + 1*16 = 35
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),
@ -129,46 +126,36 @@ mod test {
#[test]
fn test_h_over_rt() {
// At T=1: ln(1/T)=0, so the a[1] term vanishes
// -a[0] + 0 + a[2] + a[3]/2 + a[4]/3 + a[5]/4 + a[6]/5 + a[7]
// = -1 + 0 + 3 + 12/2 + 9/3 + 8/4 + 5/5 + 8 = -1+3+6+3+2+1+8 = 22
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(1.0);
assert_delta!(result, 22.0, 1e-10);
.h_over_rt(2.0);
assert_delta!(result, 63.44314718055995, 1e-10);
// At T=2: a[1]=0 so log term vanishes; a[4..6]=0 to avoid fractions
// -4/4 + 0 + 3 + 4*2/2 + 0 + 0 + 0 + 2/2 = -1+3+4+1 = 7
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(2.0);
assert_delta!(result, 7.0, 1e-10);
.h_over_rt(100.0);
assert_delta!(result, 203.0196, 1e-4);
}
#[test]
fn test_s_over_r() {
// At T=1: ln(T)=0, so a[2] term vanishes
// -a[0] - a[1] + 0 + a[3] + a[4]/2 + a[5]/3 + a[6]/4 + a[8]
// = -1-2+0+4+6/2+12/3+8/4+5 = -1-2+4+3+4+2+5 = 15
let result = ThermoPolynomial {
a: vec![1.0, 2.0, 3.0, 4.0, 6.0, 12.0, 8.0, 0.0, 5.0],
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(1.0);
assert_delta!(result, 15.0, 1e-10);
.s_over_r(100.0);
assert_delta!(result, 200691797.0572474, 1e-7);
// At T=2: a[2]=0 so log term vanishes; a[5..6]=0 to avoid fractions
// -4/4 - 2/2 + 0 + 2*2 + 2*4/2 + 0 + 0 + 3 = -1-1+4+4+3 = 9
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.0, 1e-10);
assert_delta!(result, 9.5, 1e-10);
}
#[test]

View file

@ -1,7 +1,7 @@
pub struct SpeciesTransportData {
pub name: String,
viscosities: Vec<TransportFit>,
conductivities: Vec<TransportFit>,
pub(crate) viscosities: Vec<TransportFit>,
pub(crate) conductivities: Vec<TransportFit>,
}
impl SpeciesTransportData {