feos_core/ad/properties/
boiling_temperature.rs1use super::PropertyAD;
2use crate::ad::Gradient;
3use crate::{FeosResult, PhaseEquilibrium, ReferenceSystem, Residual};
4use nalgebra::allocator::Allocator;
5use nalgebra::{DefaultAllocator, U1};
6use num_dual::{DualNum, DualStruct, Gradients, first_derivative, partial2};
7use quantity::{_Temperature, KELVIN, PASCAL, Pressure, Temperature};
8
9pub struct BoilingTemperature(pub Pressure);
11
12impl<'a> From<&'a [f64]> for BoilingTemperature {
13 fn from(value: &'a [f64]) -> Self {
14 Self(value[0] * PASCAL)
15 }
16}
17
18impl<N: Gradients> PropertyAD<N> for BoilingTemperature
19where
20 DefaultAllocator: Allocator<N> + Allocator<U1, N> + Allocator<N, N>,
21{
22 type Unit = _Temperature;
23 const REFERENCE: Temperature = KELVIN;
24
25 fn evaluate<E: Residual<N, D>, D: DualNum<f64, Inner = f64> + Copy>(
26 &self,
27 eos: &E,
28 ) -> FeosResult<Temperature<D>> {
29 let p = Pressure::from_inner(&self.0);
30 PhaseEquilibrium::pure_p(eos, p, None, Default::default()).map(|(t, _)| t)
31 }
32
33 fn evaluate_gradient<E: Residual<N, Gradient<P>>, const P: usize>(
34 &self,
35 eos: &E,
36 ) -> FeosResult<quantity::Quantity<Gradient<P>, Self::Unit>> {
37 let eos_f64 = eos.re();
38 let (temperature, [vapor_density, liquid_density]) =
39 PhaseEquilibrium::pure_p(&eos_f64, self.0, None, Default::default())?;
40
41 let t = temperature.into_reduced();
42 let v1 = 1.0 / liquid_density.to_reduced();
43 let v2 = 1.0 / vapor_density.to_reduced();
44 let p = self.0.into_reduced();
45 let t = Gradient::from(t);
46 let t = t + {
47 let v1 = Gradient::from(v1);
48 let v2 = Gradient::from(v2);
49 let p = Gradient::from(p);
50 let x = E::pure_molefracs();
51
52 let residual_entropy = |v| {
53 let (a, s) = first_derivative(
54 partial2(
55 |t, &v, x| eos.lift().residual_helmholtz_energy(t, v, x),
56 &v,
57 &x,
58 ),
59 t,
60 );
61 (a, -s)
62 };
63 let (a1, s1) = residual_entropy(v1);
64 let (a2, s2) = residual_entropy(v2);
65
66 let ln_rho = (v1 / v2).ln();
67 (p * (v2 - v1) + (a2 - a1 + t * ln_rho)) / (s2 - s1 - ln_rho)
68 };
69 Ok(Temperature::from_reduced(t))
70 }
71}