feos_ad/core/
state.rs

1use super::{
2    FeOsWrapper, HelmholtzEnergyWrapper, ParametersAD, ResidualHelmholtzEnergy,
3    TotalHelmholtzEnergy,
4};
5use feos_core::{DensityInitialization, EosResult, ReferenceSystem, State};
6use nalgebra::{Const, SMatrix, SVector};
7use ndarray::arr1;
8use num_dual::{hessian, jacobian, third_derivative, Dual2Vec, Dual3, DualNum, DualVec};
9use quantity::{MolarEnergy, MolarEntropy, MolarVolume, Moles, Pressure, Temperature};
10
11/// An (intensive) thermodynamic state representing a single phase.
12pub struct StateAD<'a, E: ParametersAD, D: DualNum<f64> + Copy, const N: usize> {
13    pub eos: &'a HelmholtzEnergyWrapper<E, D, N>,
14    pub temperature: Temperature<D>,
15    pub molar_volume: MolarVolume<D>,
16    pub molefracs: SVector<D, N>,
17    pub reduced_temperature: D,
18    pub reduced_molar_volume: D,
19}
20
21impl<'a, E: ParametersAD, D: DualNum<f64> + Copy, const N: usize> StateAD<'a, E, D, N> {
22    /// Crate a state from its thermodynamic variables (temperature, molar volume, composition)
23    pub fn new(
24        eos: &'a HelmholtzEnergyWrapper<E, D, N>,
25        temperature: D,
26        molar_volume: D,
27        molefracs: SVector<D, N>,
28    ) -> Self {
29        Self {
30            eos,
31            temperature: Temperature::from_reduced(temperature),
32            molar_volume: MolarVolume::from_reduced(molar_volume),
33            molefracs,
34            reduced_temperature: temperature,
35            reduced_molar_volume: molar_volume,
36        }
37    }
38
39    fn from_state<
40        F: Fn(
41            &E::Parameters<DualVec<D, f64, Const<2>>>,
42            DualVec<D, f64, Const<2>>,
43            DualVec<D, f64, Const<2>>,
44            &SVector<DualVec<D, f64, Const<2>>, N>,
45        ) -> SVector<DualVec<D, f64, Const<2>>, 2>,
46    >(
47        eos: &'a HelmholtzEnergyWrapper<E, D, N>,
48        state: State<FeOsWrapper<E, N>>,
49        molefracs: SVector<D, N>,
50        f: F,
51        rhs: SVector<D, 2>,
52    ) -> Self {
53        let mut vars = SVector::from([
54            D::from(state.temperature.to_reduced()),
55            D::from(state.density.to_reduced().recip()),
56        ]);
57        let x = molefracs.map(DualVec::from_re);
58        let params = E::params_from_inner(&eos.parameters);
59        for _ in 0..D::NDERIV {
60            let (mut f, jac) = jacobian(|vars| f(&params, vars[0], vars[1], &x), vars);
61            f -= rhs;
62            let det = (jac[(0, 0)] * jac[(1, 1)] - jac[(0, 1)] * jac[(1, 0)]).recip();
63            vars[0] -= det * (jac[(1, 1)] * f[0] - jac[(0, 1)] * f[1]);
64            vars[1] -= det * (jac[(0, 0)] * f[1] - jac[(1, 0)] * f[0]);
65        }
66        let [temperature, molar_volume] = vars.data.0[0];
67        Self::new(eos, temperature, molar_volume, molefracs)
68    }
69}
70
71impl<'a, E: ResidualHelmholtzEnergy<N>, D: DualNum<f64> + Copy, const N: usize>
72    StateAD<'a, E, D, N>
73{
74    /// Calculate a state from given temperature, pressure and composition.
75    pub fn new_tp(
76        eos: &'a HelmholtzEnergyWrapper<E, D, N>,
77        temperature: Temperature<D>,
78        pressure: Pressure<D>,
79        molefracs: SVector<D, N>,
80        density_initialization: DensityInitialization,
81    ) -> EosResult<Self> {
82        let t = temperature.re();
83        let p = pressure.re();
84        let moles = Moles::from_reduced(arr1(&molefracs.data.0[0].map(|x| x.re())));
85        let state = State::new_npt(&eos.eos, t, p, &moles, density_initialization)?;
86        let mut density = D::from(state.density.to_reduced());
87        let t = temperature.into_reduced();
88        for _ in 0..D::NDERIV {
89            let (_, p, dp_drho) = E::dp_drho(&eos.parameters, t, density.recip(), &molefracs);
90            density -= (p - pressure.into_reduced()) / dp_drho;
91        }
92        Ok(Self::new(eos, t, density.recip(), molefracs))
93    }
94
95    pub fn pressure(&self) -> Pressure<D> {
96        Pressure::from_reduced(E::pressure(
97            &self.eos.parameters,
98            self.reduced_temperature,
99            self.reduced_molar_volume,
100            &self.molefracs,
101        ))
102    }
103}
104
105impl<'a, E: TotalHelmholtzEnergy<N>, D: DualNum<f64> + Copy, const N: usize> StateAD<'a, E, D, N> {
106    /// Calculate a state from given pressure, molar entropy and composition.
107    pub fn new_ps(
108        eos: &'a HelmholtzEnergyWrapper<E, D, N>,
109        pressure: Pressure<D>,
110        molar_entropy: MolarEntropy<D>,
111        molefracs: SVector<D, N>,
112        density_initialization: DensityInitialization,
113        initial_temperature: Option<Temperature>,
114    ) -> EosResult<Self> {
115        let moles = Moles::from_reduced(arr1(&molefracs.data.0[0].map(|x| x.re())));
116        let state = State::new_nps(
117            &eos.eos,
118            pressure.re(),
119            molar_entropy.re(),
120            &moles,
121            density_initialization,
122            initial_temperature,
123        )?;
124        Ok(Self::from_state(
125            eos,
126            state,
127            molefracs,
128            E::pressure_entropy,
129            SVector::from([pressure.into_reduced(), molar_entropy.into_reduced()]),
130        ))
131    }
132
133    /// Calculate a state from given pressure, molar enthalpy and composition.
134    pub fn new_ph(
135        eos: &'a HelmholtzEnergyWrapper<E, D, N>,
136        pressure: Pressure<D>,
137        molar_enthalpy: MolarEnergy<D>,
138        molefracs: SVector<D, N>,
139        density_initialization: DensityInitialization,
140        initial_temperature: Option<Temperature>,
141    ) -> EosResult<Self> {
142        let moles = Moles::from_reduced(arr1(&molefracs.data.0[0].map(|x| x.re())));
143        let state = State::new_nph(
144            &eos.eos,
145            pressure.re(),
146            molar_enthalpy.re(),
147            &moles,
148            density_initialization,
149            initial_temperature,
150        )?;
151        Ok(Self::from_state(
152            eos,
153            state,
154            molefracs,
155            E::pressure_enthalpy,
156            SVector::from([pressure.into_reduced(), molar_enthalpy.into_reduced()]),
157        ))
158    }
159
160    pub fn molar_entropy(&self) -> MolarEntropy<D> {
161        MolarEntropy::from_reduced(E::molar_entropy(
162            &self.eos.parameters,
163            self.reduced_temperature,
164            self.reduced_molar_volume,
165            &self.molefracs,
166        ))
167    }
168
169    pub fn molar_enthalpy(&self) -> MolarEnergy<D> {
170        MolarEnergy::from_reduced(E::molar_enthalpy(
171            &self.eos.parameters,
172            self.reduced_temperature,
173            self.reduced_molar_volume,
174            &self.molefracs,
175        ))
176    }
177
178    pub fn molar_isochoric_heat_capacity(&self) -> MolarEntropy<D> {
179        MolarEntropy::from_reduced(E::molar_isochoric_heat_capacity(
180            &self.eos.parameters,
181            self.reduced_temperature,
182            self.reduced_molar_volume,
183            &self.molefracs,
184        ))
185    }
186
187    pub fn molar_isobaric_heat_capacity(&self) -> MolarEntropy<D> {
188        MolarEntropy::from_reduced(E::molar_isobaric_heat_capacity(
189            &self.eos.parameters,
190            self.reduced_temperature,
191            self.reduced_molar_volume,
192            &self.molefracs,
193        ))
194    }
195}
196
197impl<'a, E: ResidualHelmholtzEnergy<1>, D: DualNum<f64> + Copy> StateAD<'a, E, D, 1> {
198    /// Calculate the critical point of a pure component.
199    pub fn critical_point_pure(eos: &'a HelmholtzEnergyWrapper<E, D, 1>) -> EosResult<Self> {
200        Self::critical_point(eos, SVector::from([D::one()]))
201    }
202}
203
204impl<'a, E: ResidualHelmholtzEnergy<N>, D: DualNum<f64> + Copy, const N: usize>
205    StateAD<'a, E, D, N>
206{
207    /// Calculate the critical point of a mixture with given composition.
208    pub fn critical_point(
209        eos: &'a HelmholtzEnergyWrapper<E, D, N>,
210        molefracs: SVector<D, N>,
211    ) -> EosResult<Self>
212    where
213        Const<N>: Eigen<N>,
214    {
215        let moles = Moles::from_reduced(arr1(molefracs.map(|x| x.re()).as_slice()));
216        let state = State::critical_point(&eos.eos, Some(&moles), None, Default::default())?;
217        Ok(Self::from_state(
218            eos,
219            state,
220            molefracs,
221            Self::criticality_conditions,
222            SVector::from([D::from(0.0); 2]),
223        ))
224    }
225}
226
227impl<E: ResidualHelmholtzEnergy<N>, D: DualNum<f64> + Copy, const N: usize> StateAD<'_, E, D, N> {
228    fn criticality_conditions(
229        parameters: &E::Parameters<DualVec<D, f64, Const<2>>>,
230        temperature: DualVec<D, f64, Const<2>>,
231        molar_volume: DualVec<D, f64, Const<2>>,
232        molefracs: &SVector<DualVec<D, f64, Const<2>>, N>,
233    ) -> SVector<DualVec<D, f64, Const<2>>, 2>
234    where
235        Const<N>: Eigen<N>,
236    {
237        // calculate M
238        let sqrt_z = molefracs.map(|z| z.sqrt());
239        let z_mix = sqrt_z * sqrt_z.transpose();
240        let (_, _, m) = hessian(
241            |x| {
242                let params = E::params_from_inner(parameters);
243                let t = Dual2Vec::from_re(temperature);
244                let v = Dual2Vec::from_re(molar_volume);
245                E::residual_molar_helmholtz_energy(&params, t, v, &x)
246            },
247            *molefracs,
248        );
249        let m = m.component_mul(&z_mix) / temperature + SMatrix::identity();
250
251        // calculate smallest eigenvalue and corresponding eigenvector
252        let (l, u) = <Const<N> as Eigen<N>>::eigen(m);
253
254        let (_, _, _, c2) = third_derivative(
255            |s| {
256                let x = molefracs.map(Dual3::from_re);
257                let x = x + sqrt_z.component_mul(&u).map(Dual3::from_re) * s;
258                let params = E::params_from_inner(parameters);
259                let t = Dual3::from_re(temperature);
260                let v = Dual3::from_re(molar_volume);
261                let ig = x.component_mul(&x.map(|x| (x / v).ln() - 1.0)).sum();
262                E::residual_molar_helmholtz_energy(&params, t, v, &x) / t + ig
263            },
264            DualVec::from_re(D::zero()),
265        );
266
267        SVector::from([l, c2])
268    }
269}
270
271pub trait Eigen<const N: usize> {
272    fn eigen<D: DualNum<f64> + Copy>(matrix: SMatrix<D, N, N>) -> (D, SVector<D, N>);
273}
274
275impl Eigen<1> for Const<1> {
276    fn eigen<D: DualNum<f64> + Copy>(matrix: SMatrix<D, 1, 1>) -> (D, SVector<D, 1>) {
277        let [[l]] = matrix.data.0;
278        (l, SVector::from([D::one()]))
279    }
280}
281
282impl Eigen<2> for Const<2> {
283    fn eigen<D: DualNum<f64> + Copy>(matrix: SMatrix<D, 2, 2>) -> (D, SVector<D, 2>) {
284        let [[a, b], [_, c]] = matrix.data.0;
285        let l = (a + c - ((a - c).powi(2) + b * b * 4.0).sqrt()) * 0.5;
286        let u = SVector::from([D::one(), (l - a) / b]);
287        let u = u / (u[0] * u[0] + u[1] * u[1]).sqrt();
288        (l, u)
289    }
290}
291
292#[cfg(test)]
293mod test {
294    use super::*;
295    use crate::eos::ideal_gas::test::joback;
296    use crate::eos::pcsaft::test::pcsaft;
297    use crate::eos::PcSaftPure;
298    use crate::EquationOfStateAD;
299    use approx::assert_relative_eq;
300    use feos_core::{Contributions, EosResult, EquationOfState, PhaseEquilibrium};
301    use num_dual::{Dual, Dual64};
302    use quantity::{BAR, JOULE, KELVIN, MOL, PASCAL};
303    use std::sync::Arc;
304
305    #[test]
306    fn test_critical_point() -> EosResult<()> {
307        let (pcsaft, eos) = pcsaft()?;
308        let mut params = pcsaft.params();
309        let mut params_dual = pcsaft.params::<Dual64>();
310        params_dual[0] = params_dual[0].derivative();
311        let pcsaft = pcsaft.wrap().derivatives(params_dual);
312        let cp = State::critical_point(&eos, None, None, Default::default())?;
313        let state = StateAD::critical_point_pure(&pcsaft)?;
314        let t = state.temperature.re();
315        let rho = 1.0 / state.molar_volume.re();
316        println!("{:.5} {:.5}", t, rho);
317        println!("{:.5} {:.5}", cp.temperature, cp.density);
318        assert_relative_eq!(t, cp.temperature, max_relative = 1e-10);
319        assert_relative_eq!(rho, cp.density, max_relative = 1e-10);
320
321        let h = 1e-8;
322        params[0] += h;
323        let eos = PcSaftPure(params).wrap();
324        let cp_h = State::critical_point(&eos.eos, None, None, Default::default())?;
325        let dt = (cp_h.temperature - cp.temperature).to_reduced() / h;
326        let drho = (cp_h.density - cp.density).to_reduced() / h;
327
328        println!(
329            "{:.5e} {:.5e}",
330            state.reduced_temperature.eps,
331            state.reduced_molar_volume.recip().eps
332        );
333        println!("{:.5e} {:.5e}", dt, drho);
334        assert_relative_eq!(state.reduced_temperature.eps, dt, max_relative = 1e-6);
335        assert_relative_eq!(
336            state.reduced_molar_volume.recip().eps,
337            drho,
338            max_relative = 1e-6
339        );
340        Ok(())
341    }
342
343    #[test]
344    fn test_state_tp() -> EosResult<()> {
345        let (pcsaft, eos) = pcsaft()?;
346        let pcsaft = pcsaft.wrap();
347        let p = BAR;
348        let t = 300.0 * KELVIN;
349        let state_feos = State::new_npt(
350            &eos,
351            t,
352            p,
353            &(arr1(&[1.0]) * MOL),
354            DensityInitialization::Liquid,
355        )?;
356        let state = StateAD::new_tp(
357            &pcsaft,
358            t,
359            p,
360            SVector::from([1.0]),
361            DensityInitialization::Liquid,
362        )?;
363        let density = 1.0 / state.molar_volume;
364        println!("{:.5}", density);
365        println!("{:.5}", state_feos.density);
366        assert_relative_eq!(density, state_feos.density, max_relative = 1e-10);
367        Ok(())
368    }
369
370    #[test]
371    fn test_state_tp_derivative() -> EosResult<()> {
372        let (pcsaft, residual) = pcsaft()?;
373        let p = BAR;
374        let t = 300.0 * KELVIN;
375        let state_feos = State::new_npt(
376            &residual,
377            t,
378            p,
379            &(arr1(&[1.0]) * MOL),
380            DensityInitialization::Liquid,
381        )?;
382        let h = 1e2 * PASCAL;
383        let state_h = State::new_npt(
384            &residual,
385            t,
386            p + h,
387            &(arr1(&[1.0]) * MOL),
388            DensityInitialization::Liquid,
389        )?;
390        let params: [Dual64; 8] = pcsaft.params();
391        let eos_ad = pcsaft.wrap().derivatives(params);
392        let t = Temperature::from_reduced(Dual::from(t.to_reduced()));
393        let p = Pressure::from_reduced(Dual::from(p.to_reduced()).derivative());
394        let state = StateAD::new_tp(
395            &eos_ad,
396            t,
397            p,
398            SVector::from([Dual::from(1.0)]),
399            DensityInitialization::Liquid,
400        )?;
401        let density = state.molar_volume.into_reduced().recip();
402        println!("{:.5} {:.5}", density.re, density.eps);
403        let density_h = ((state_h.density - state_feos.density) / h).into_reduced();
404        println!("{:.5} {:.5}", state_feos.density.into_reduced(), density_h);
405        assert_relative_eq!(density.eps, density_h, max_relative = 1e-6);
406        Ok(())
407    }
408
409    #[test]
410    fn test_state_ps() -> EosResult<()> {
411        let (joback, ideal_gas) = joback()?;
412        let (pcsaft, residual) = pcsaft()?;
413        let eos_ad = EquationOfStateAD::new([joback], pcsaft).wrap();
414        let eos = Arc::new(EquationOfState::new(ideal_gas, residual));
415        let vle = PhaseEquilibrium::pure(&eos, 250.0 * KELVIN, None, Default::default())?;
416        let p = vle.liquid().pressure(Contributions::Total);
417        let s = vle.liquid().molar_entropy(Contributions::Total);
418        let t = vle.liquid().temperature;
419        let state = StateAD::new_ps(
420            &eos_ad,
421            p,
422            s,
423            SVector::from([1.0]),
424            DensityInitialization::Liquid,
425            Some(t),
426        )?;
427        let density = 1.0 / state.molar_volume;
428        println!("{:.5} {:.5}", state.temperature, density);
429        println!(
430            "{:.5} {:.5}",
431            vle.liquid().temperature,
432            vle.liquid().density,
433        );
434        assert_relative_eq!(
435            state.temperature,
436            vle.liquid().temperature,
437            max_relative = 1e-10
438        );
439        assert_relative_eq!(density, vle.liquid().density, max_relative = 1e-10);
440        Ok(())
441    }
442
443    #[test]
444    fn test_state_ps_derivative() -> EosResult<()> {
445        let (joback, ideal_gas) = joback()?;
446        let (pcsaft, residual) = pcsaft()?;
447        let eos = Arc::new(EquationOfState::new(ideal_gas, residual));
448        let vle = PhaseEquilibrium::pure(&eos, 250.0 * KELVIN, None, Default::default())?;
449        let h = 1e-3 * JOULE / KELVIN / MOL;
450        let state_h = State::new_nps(
451            &eos,
452            vle.liquid().pressure(Contributions::Total),
453            vle.liquid().molar_entropy(Contributions::Total) + h,
454            &vle.liquid().moles,
455            DensityInitialization::Liquid,
456            Some(vle.liquid().temperature),
457        )?;
458        let p = vle.liquid().pressure(Contributions::Total).to_reduced();
459        let s = vle
460            .liquid()
461            .molar_entropy(Contributions::Total)
462            .to_reduced();
463        let t = vle.liquid().temperature;
464        let eos_ad = EquationOfStateAD::new([joback], pcsaft)
465            .wrap()
466            .derivatives(([joback.params()], pcsaft.params()));
467        let p: Pressure<Dual64> = Pressure::from_reduced(Dual::from(p));
468        let s = MolarEntropy::from_reduced(Dual::from(s).derivative());
469        let state = StateAD::new_ps(
470            &eos_ad,
471            p,
472            s,
473            SVector::from([Dual::from(1.0)]),
474            DensityInitialization::Liquid,
475            Some(t),
476        )?;
477        println!(
478            "{:.5e} {:.5e}",
479            state.reduced_temperature.eps,
480            state.reduced_molar_volume.recip().eps
481        );
482        println!(
483            "{:.5e} {:.5e}",
484            ((state_h.temperature - vle.liquid().temperature) / h).to_reduced(),
485            ((state_h.density - vle.liquid().density) / h).to_reduced(),
486        );
487        assert_relative_eq!(
488            state.reduced_temperature.eps,
489            ((state_h.temperature - vle.liquid().temperature) / h).to_reduced(),
490            max_relative = 1e-6
491        );
492        assert_relative_eq!(
493            state.reduced_molar_volume.recip().eps,
494            ((state_h.density - vle.liquid().density) / h).to_reduced(),
495            max_relative = 1e-6
496        );
497        Ok(())
498    }
499
500    #[test]
501    fn test_state_ph() -> EosResult<()> {
502        let (joback, ideal_gas) = joback()?;
503        let (pcsaft, residual) = pcsaft()?;
504        let eos_ad = EquationOfStateAD::new([joback], pcsaft).wrap();
505        let eos = Arc::new(EquationOfState::new(ideal_gas, residual));
506        let vle = PhaseEquilibrium::pure(&eos, 250.0 * KELVIN, None, Default::default())?;
507        let p = vle.liquid().pressure(Contributions::Total);
508        let h = vle.liquid().molar_enthalpy(Contributions::Total);
509        let t = vle.liquid().temperature;
510        let state = StateAD::new_ph(
511            &eos_ad,
512            p,
513            h,
514            SVector::from([1.0]),
515            DensityInitialization::Liquid,
516            Some(t),
517        )?;
518        let density = 1.0 / state.molar_volume;
519        println!("{:.5} {:.5}", state.temperature, density);
520        println!(
521            "{:.5} {:.5}",
522            vle.liquid().temperature,
523            vle.liquid().density,
524        );
525        assert_relative_eq!(
526            state.temperature,
527            vle.liquid().temperature,
528            max_relative = 1e-10
529        );
530        assert_relative_eq!(density, vle.liquid().density, max_relative = 1e-10);
531        Ok(())
532    }
533
534    #[test]
535    fn test_state_ph_derivative() -> EosResult<()> {
536        let (joback, ideal_gas) = joback()?;
537        let (pcsaft, residual) = pcsaft()?;
538        let eos = Arc::new(EquationOfState::new(ideal_gas, residual));
539        let vle = PhaseEquilibrium::pure(&eos, 250.0 * KELVIN, None, Default::default())?;
540        let delta = 1e-1 * JOULE / MOL;
541        let state_h = State::new_nph(
542            &eos,
543            vle.liquid().pressure(Contributions::Total),
544            vle.liquid().molar_enthalpy(Contributions::Total) + delta,
545            &vle.liquid().moles,
546            DensityInitialization::Liquid,
547            Some(vle.liquid().temperature),
548        )?;
549        let p = vle.liquid().pressure(Contributions::Total).to_reduced();
550        let h = vle
551            .liquid()
552            .molar_enthalpy(Contributions::Total)
553            .to_reduced();
554        let t = vle.liquid().temperature;
555        let eos_ad = EquationOfStateAD::new([joback], pcsaft)
556            .wrap()
557            .derivatives(([joback.params()], pcsaft.params()));
558        let p: Pressure<Dual64> = Pressure::from_reduced(Dual::from(p));
559        let h = MolarEnergy::from_reduced(Dual::from(h).derivative());
560        let state = StateAD::new_ph(
561            &eos_ad,
562            p,
563            h,
564            SVector::from([Dual::from(1.0)]),
565            DensityInitialization::Liquid,
566            Some(t),
567        )?;
568        println!(
569            "{:.5e} {:.5e}",
570            state.reduced_temperature.eps,
571            state.reduced_molar_volume.recip().eps
572        );
573        println!(
574            "{:.5e} {:.5e}",
575            ((state_h.temperature - vle.liquid().temperature) / delta).to_reduced(),
576            ((state_h.density - vle.liquid().density) / delta).to_reduced(),
577        );
578        assert_relative_eq!(
579            state.reduced_temperature.eps,
580            ((state_h.temperature - vle.liquid().temperature) / delta).to_reduced(),
581            max_relative = 1e-6
582        );
583        assert_relative_eq!(
584            state.reduced_molar_volume.recip().eps,
585            ((state_h.density - vle.liquid().density) / delta).to_reduced(),
586            max_relative = 1e-6
587        );
588        Ok(())
589    }
590
591    #[test]
592    fn test_heat_capacities() -> EosResult<()> {
593        let (joback, ideal_gas) = joback()?;
594        let (pcsaft, residual) = pcsaft()?;
595        let eos = Arc::new(EquationOfState::new(ideal_gas, residual));
596        let eos_ad = EquationOfStateAD::new([joback], pcsaft)
597            .wrap()
598            .derivatives(([joback.params()], pcsaft.params()));
599
600        let temperature = 300.0 * KELVIN;
601        let pressure = 5.0 * BAR;
602
603        let state = State::new_npt(
604            &eos,
605            temperature,
606            pressure,
607            &(arr1(&[1.0]) * MOL),
608            DensityInitialization::None,
609        )?;
610        let state_ad = StateAD::new_tp(
611            &eos_ad,
612            temperature,
613            pressure,
614            SVector::from([1.0]),
615            DensityInitialization::None,
616        )?;
617
618        let c_v = state.molar_isochoric_heat_capacity(Contributions::Total);
619        let c_p = state.molar_isobaric_heat_capacity(Contributions::Total);
620        let c_v_ad = state_ad.molar_isochoric_heat_capacity();
621        let c_p_ad = state_ad.molar_isobaric_heat_capacity();
622
623        println!("{c_v} {c_p}");
624        println!("{c_v_ad} {c_p_ad}");
625
626        assert_relative_eq!(c_v, c_v_ad, max_relative = 1e-10);
627        assert_relative_eq!(c_p, c_p_ad, max_relative = 1e-10);
628
629        Ok(())
630    }
631}
feos_ad/core/state.rs

feos_ad/core/
state.rs
