|
| 1 | +use num_dual::DualNum; |
| 2 | +use serde::Deserialize; |
| 3 | +use serde_json::Value; |
| 4 | +use std::collections::HashMap; |
| 5 | + |
| 6 | +#[derive(Clone, Deserialize)] |
| 7 | +pub struct IdealGasFunctionJson { |
| 8 | + #[serde(rename = "type")] |
| 9 | + pub ty: String, |
| 10 | + #[serde(flatten)] |
| 11 | + parameters: HashMap<String, Value>, |
| 12 | +} |
| 13 | + |
| 14 | +pub struct IdealGasFunctionIterator { |
| 15 | + inner: IdealGasFunctionJson, |
| 16 | + count: usize, |
| 17 | + index: usize, |
| 18 | +} |
| 19 | + |
| 20 | +impl Iterator for IdealGasFunctionIterator { |
| 21 | + type Item = IdealGasFunction; |
| 22 | + |
| 23 | + fn next(&mut self) -> Option<Self::Item> { |
| 24 | + if self.index == self.count { |
| 25 | + return None; |
| 26 | + } |
| 27 | + let mut parameters: HashMap<_, _> = self |
| 28 | + .inner |
| 29 | + .parameters |
| 30 | + .iter() |
| 31 | + .map(|(k, v)| { |
| 32 | + ( |
| 33 | + k.clone(), |
| 34 | + if self.inner.ty == "IdealGasHelmholtzCP0AlyLee" { |
| 35 | + // AlyLee parameters are stored as list instead of named parameters |
| 36 | + v.clone() |
| 37 | + } else { |
| 38 | + v.as_array().map_or(v, |v| &v[self.index]).clone() |
| 39 | + }, |
| 40 | + ) |
| 41 | + }) |
| 42 | + .collect(); |
| 43 | + if self.inner.ty == "IdealGasHelmholtzCP0AlyLee" { |
| 44 | + self.index += 5 |
| 45 | + } else { |
| 46 | + self.index += 1; |
| 47 | + } |
| 48 | + parameters.insert("type".into(), serde_json::to_value(&self.inner.ty).unwrap()); |
| 49 | + Some(serde_json::from_value(serde_json::to_value(parameters).unwrap()).unwrap()) |
| 50 | + } |
| 51 | +} |
| 52 | + |
| 53 | +impl IntoIterator for IdealGasFunctionJson { |
| 54 | + type Item = IdealGasFunction; |
| 55 | + type IntoIter = IdealGasFunctionIterator; |
| 56 | + |
| 57 | + fn into_iter(self) -> Self::IntoIter { |
| 58 | + let count = self |
| 59 | + .parameters |
| 60 | + .values() |
| 61 | + .map(|e| e.as_array().map_or(1, Vec::len)) |
| 62 | + .max() |
| 63 | + .unwrap(); |
| 64 | + IdealGasFunctionIterator { |
| 65 | + count, |
| 66 | + inner: self, |
| 67 | + index: 0, |
| 68 | + } |
| 69 | + } |
| 70 | +} |
| 71 | + |
| 72 | +#[derive(Clone, Copy, Deserialize)] |
| 73 | +#[serde(tag = "type")] |
| 74 | +#[expect(non_snake_case)] |
| 75 | +pub enum IdealGasFunction { |
| 76 | + IdealGasHelmholtzLead { a1: f64, a2: f64 }, |
| 77 | + IdealGasHelmholtzLogTau { a: f64 }, |
| 78 | + IdealGasHelmholtzPower { n: f64, t: f64 }, |
| 79 | + IdealGasHelmholtzPlanckEinstein { n: f64, t: f64 }, |
| 80 | + IdealGasHelmholtzPlanckEinsteinFunctionT { Tcrit: f64, n: f64, v: f64 }, |
| 81 | + IdealGasHelmholtzPlanckEinsteinGeneralized { c: f64, d: f64, n: f64, t: f64 }, |
| 82 | + IdealGasHelmholtzCP0AlyLee { T0: f64, Tc: f64, c: [f64; 5] }, |
| 83 | + IdealGasHelmholtzCP0Constant { T0: f64, Tc: f64, cp_over_R: f64 }, |
| 84 | + IdealGasHelmholtzCP0PolyT { T0: f64, Tc: f64, c: f64, t: f64 }, |
| 85 | + IdealGasHelmholtzEnthalpyEntropyOffset { a1: f64, a2: f64 }, |
| 86 | +} |
| 87 | + |
| 88 | +impl IdealGasFunction { |
| 89 | + pub fn evaluate<D: DualNum<f64> + Copy>(&self, delta: D, tau: D) -> D { |
| 90 | + match *self { |
| 91 | + IdealGasFunction::IdealGasHelmholtzLead { a1, a2 } => delta.ln() + a1 + tau * a2, |
| 92 | + IdealGasFunction::IdealGasHelmholtzLogTau { a } => tau.ln() * a, |
| 93 | + IdealGasFunction::IdealGasHelmholtzPower { n, t } => tau.powf(t) * n, |
| 94 | + IdealGasFunction::IdealGasHelmholtzPlanckEinstein { n, t } => { |
| 95 | + (-(-tau * t).exp()).ln_1p() * n |
| 96 | + } |
| 97 | + IdealGasFunction::IdealGasHelmholtzPlanckEinsteinFunctionT { Tcrit, n, v } => { |
| 98 | + (-(-tau * v / Tcrit).exp()).ln_1p() * n |
| 99 | + } |
| 100 | + IdealGasFunction::IdealGasHelmholtzPlanckEinsteinGeneralized { c, d, n, t } => { |
| 101 | + ((tau * t).exp() * d + c).ln() * n |
| 102 | + } |
| 103 | + IdealGasFunction::IdealGasHelmholtzCP0AlyLee { T0, Tc, c } => { |
| 104 | + let [a, b, c, d, e] = c; |
| 105 | + let mut res = D::zero(); |
| 106 | + if a > 0.0 { |
| 107 | + let tau0 = Tc / T0; |
| 108 | + res += (-tau / tau0 + 1.0 + (tau / tau0).ln()) * a; |
| 109 | + } |
| 110 | + if b > 0.0 { |
| 111 | + res += (-(-tau * 2.0 * c / Tc).exp()).ln_1p() * b; |
| 112 | + } |
| 113 | + if d > 0.0 { |
| 114 | + res -= ((-tau * 2.0 * e / Tc).exp()).ln_1p() * d; |
| 115 | + } |
| 116 | + res |
| 117 | + } |
| 118 | + IdealGasFunction::IdealGasHelmholtzCP0Constant { T0, Tc, cp_over_R } => { |
| 119 | + let tau0 = Tc / T0; |
| 120 | + (-tau / tau0 + 1.0 + (tau / tau0).ln()) * cp_over_R |
| 121 | + } |
| 122 | + IdealGasFunction::IdealGasHelmholtzCP0PolyT { T0, Tc, c, t } => { |
| 123 | + // unfortunately some models use floats and other models use 0 or -1 which needs to be treated separately... |
| 124 | + if t.abs() < 10.0 * f64::EPSILON { |
| 125 | + let tau0 = Tc / T0; |
| 126 | + (-tau / tau0 + 1.0 + (tau / tau0).ln()) * c |
| 127 | + } else if (t + 1.0).abs() < 10.0 * f64::EPSILON { |
| 128 | + let tau0 = Tc / T0; |
| 129 | + (-tau / Tc * (tau / tau0).ln() + (tau - tau0) / Tc) * c |
| 130 | + } else { |
| 131 | + (-tau.powf(-t) * Tc.powf(t) / (t * (t + 1.0)) |
| 132 | + - tau * T0.powf(t + 1.0) / (Tc * (t + 1.0)) |
| 133 | + + T0.powf(t) / t) |
| 134 | + * c |
| 135 | + } |
| 136 | + } |
| 137 | + IdealGasFunction::IdealGasHelmholtzEnthalpyEntropyOffset { a1, a2 } => tau * a2 + a1, |
| 138 | + } |
| 139 | + } |
| 140 | +} |
0 commit comments