|
| 1 | +use feos_core::{Components, IdealGas, Residual, StateHD}; |
| 2 | +use nalgebra::{Const, SVector, U1}; |
| 3 | +use ndarray::{arr1, Array1, ScalarOperand}; |
| 4 | +use num_dual::{Derivative, DualNum, DualVec}; |
| 5 | +use std::sync::Arc; |
| 6 | + |
| 7 | +mod phase_equilibria; |
| 8 | +mod residual; |
| 9 | +mod state; |
| 10 | +mod total; |
| 11 | +pub use phase_equilibria::PhaseEquilibriumAD; |
| 12 | +pub use residual::{ParametersAD, ResidualHelmholtzEnergy}; |
| 13 | +pub use state::{Eigen, StateAD}; |
| 14 | +pub use total::{EquationOfStateAD, IdealGasAD, TotalHelmholtzEnergy}; |
| 15 | + |
| 16 | +/// Used internally to implement the [Residual] and [IdealGas] traits from FeOs. |
| 17 | +pub struct FeOsWrapper<E, const N: usize>(E); |
| 18 | + |
| 19 | +impl<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 | + |
| 29 | +impl<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(¶meters, temperature, &density) * volume |
| 45 | + / temperature; |
| 46 | + vec![(R::RESIDUAL.into(), a)] |
| 47 | + } |
| 48 | +} |
| 49 | + |
| 50 | +impl<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(¶meters, 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. |
| 62 | +pub 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 | + |
| 67 | +impl<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. |
| 81 | +pub 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 | + |
| 89 | +impl<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 | +} |
0 commit comments