feos_ad/core/
mod.rs

1use feos_core::{Components, IdealGas, Residual, StateHD};
2use nalgebra::{Const, SVector, U1};
3use ndarray::{arr1, Array1, ScalarOperand};
4use num_dual::{Derivative, DualNum, DualVec};
5use std::sync::Arc;
6
7mod phase_equilibria;
8mod residual;
9mod state;
10mod total;
11pub use phase_equilibria::PhaseEquilibriumAD;
12pub use residual::{ParametersAD, ResidualHelmholtzEnergy};
13pub use state::{Eigen, StateAD};
14pub use total::{EquationOfStateAD, IdealGasAD, TotalHelmholtzEnergy};
15
16/// Used internally to implement the [Residual] and [IdealGas] traits from FeOs.
17pub struct FeOsWrapper<E, const N: usize>(E);
18
19impl<R: ParametersAD, const N: usize> Components for FeOsWrapper<R, N> {
20    fn components(&self) -> usize {
21        N
22    }
23
24    fn subset(&self, _: &[usize]) -> Self {
25        panic!("Calculating properties of subsets of models is not possible with AD.")
26    }
27}
28
29impl<R: ResidualHelmholtzEnergy<N>, const N: usize> Residual for FeOsWrapper<R, N> {
30    fn compute_max_density(&self, moles: &Array1<f64>) -> f64 {
31        let total_moles = moles.sum();
32        let molefracs = SVector::from_fn(|i, _| moles[i] / total_moles);
33        self.0.compute_max_density(&molefracs)
34    }
35
36    fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy + ScalarOperand>(
37        &self,
38        state: &StateHD<D>,
39    ) -> Vec<(String, D)> {
40        let temperature = state.temperature;
41        let volume = state.volume;
42        let density = SVector::from_column_slice(state.partial_density.as_slice().unwrap());
43        let parameters = self.0.params();
44        let a = R::residual_helmholtz_energy_density(&parameters, temperature, &density) * volume
45            / temperature;
46        vec![(R::RESIDUAL.into(), a)]
47    }
48}
49
50impl<E: TotalHelmholtzEnergy<N>, const N: usize> IdealGas for FeOsWrapper<E, N> {
51    fn ln_lambda3<D: DualNum<f64> + Copy>(&self, temperature: D) -> Array1<D> {
52        let parameters = self.0.params();
53        arr1(&E::ln_lambda3(&parameters, temperature).data.0[0])
54    }
55
56    fn ideal_gas_model(&self) -> String {
57        E::IDEAL_GAS.into()
58    }
59}
60
61/// Struct that stores the reference to the equation of state together with the (possibly) dual parameters.
62pub struct HelmholtzEnergyWrapper<E: ParametersAD, D: DualNum<f64> + Copy, const N: usize> {
63    pub eos: Arc<FeOsWrapper<E, N>>,
64    pub parameters: E::Parameters<D>,
65}
66
67impl<E: ParametersAD, const N: usize> HelmholtzEnergyWrapper<E, f64, N> {
68    /// Manually set the parameters and their derivatives.
69    pub fn derivatives<D: DualNum<f64> + Copy>(
70        &self,
71        parameters: E::Parameters<D>,
72    ) -> HelmholtzEnergyWrapper<E, D, N> {
73        HelmholtzEnergyWrapper {
74            eos: self.eos.clone(),
75            parameters,
76        }
77    }
78}
79
80/// Models for which derivatives with respect to individual parameters can be calculated.
81pub trait NamedParameters: ParametersAD {
82    /// Return a mutable reference to the parameter named by `index` from the parameter set.
83    fn index_parameters_mut<'a, D: DualNum<f64> + Copy>(
84        parameters: &'a mut Self::Parameters<D>,
85        index: &str,
86    ) -> &'a mut D;
87}
88
89impl<E: NamedParameters, const N: usize> HelmholtzEnergyWrapper<E, f64, N> {
90    /// Initialize the parameters to calculate their derivatives.
91    pub fn named_derivatives<const P: usize>(
92        &self,
93        parameters: [&str; P],
94    ) -> HelmholtzEnergyWrapper<E, DualVec<f64, f64, Const<P>>, N> {
95        let mut params: E::Parameters<DualVec<f64, f64, Const<P>>> = self.eos.0.params();
96        for (i, p) in parameters.into_iter().enumerate() {
97            E::index_parameters_mut(&mut params, p).eps =
98                Derivative::derivative_generic(Const::<P>, U1, i)
99        }
100        self.derivatives(params)
101    }
102}
feos_ad/core/mod.rs

feos_ad/core/
mod.rs
