Skip to content

Commit b37ef8e

Browse files
prehnerprehner
committed
Generic states to merge feos-core and feos-ad
1 parent 2b89f82 commit b37ef8e

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

43 files changed

+2259
-1587
lines changed

‎crates/feos-ad/Cargo.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@ approx = { workspace = true }
2323
quantity = { workspace = true, features = ["num-dual", "approx"] }
2424

2525
[features]
26-
default = []
26+
default = ["parameter_fit", "pcsaft"]
2727
pcsaft = ["feos/pcsaft"]
2828
gc_pcsaft = ["feos/gc_pcsaft"]
2929
parameter_fit = ["ndarray/rayon"]

‎crates/feos-ad/src/core/mod.rs‎

Lines changed: 51 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -1,78 +1,65 @@
1-
use feos_core::{Components, IdealGas, Residual, StateHD};
2-
use nalgebra::{Const, SVector, U1};
3-
use ndarray::{Array1, ScalarOperand, arr1};
4-
use num_dual::{Derivative, DualNum, DualVec};
5-
use std::sync::Arc;
1+
use std::ops::Deref;
62

7-
mod phase_equilibria;
3+
use nalgebra::{Const, U1};
4+
use num_dual::{Derivative, DualNum, DualSVec};
5+
6+
// mod phase_equilibria;
87
mod residual;
98
mod state;
109
mod total;
11-
pub use phase_equilibria::PhaseEquilibriumAD;
12-
pub use residual::{ParametersAD, ResidualHelmholtzEnergy};
13-
pub use state::{Eigen, StateAD};
10+
// pub use phase_equilibria::PhaseEquilibriumAD;
11+
pub use residual::ResidualHelmholtzEnergy;
12+
pub use state::{Eigen, PhaseEquilibriumAD, StateAD};
1413
pub use total::{EquationOfStateAD, IdealGasAD, TotalHelmholtzEnergy};
1514

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-
}
15+
/// A model that can be evaluated with derivatives of its parameters.
16+
pub trait ParametersAD: Sized + Deref<Target = Self::Parameters<f64>> {
17+
/// The type of the structure that stores the parameters internally.
18+
type Parameters<D: DualNum<f64> + Copy>: Clone;
2319

24-
fn subset(&self, _: &[usize]) -> Self {
25-
panic!("Calculating properties of subsets of models is not possible with AD.")
26-
}
27-
}
20+
// /// Return a reference to the parameters.
21+
// fn params(&self) -> &Self::Parameters<f64>;
2822

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-
}
23+
/// Lift the parameters to the given type of dual number.
24+
fn params_from_inner<D: DualNum<f64> + Copy, D2: DualNum<f64, Inner = D> + Copy>(
25+
parameters: &Self::Parameters<D>,
26+
) -> Self::Parameters<D2>;
3527

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-
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(&parameters, temperature).data.0[0])
28+
/// Wraps the model in the [HelmholtzEnergyWrapper] struct, so that it can be used
29+
/// as an argument to [StateAD](crate::StateAD) and [PhaseEquilibriumAD](crate::PhaseEquilibriumAD) constructors.
30+
fn wrap<'a, const N: usize>(&'a self) -> HelmholtzEnergyWrapper<'a, Self, f64, N> {
31+
HelmholtzEnergyWrapper {
32+
eos: self,
33+
parameters: self,
34+
}
5435
}
5536

56-
fn ideal_gas_model(&self) -> String {
57-
E::IDEAL_GAS.into()
37+
/// Manually set the parameters and their derivatives.
38+
fn derivatives<'a, D: DualNum<f64> + Copy, const N: usize>(
39+
&'a self,
40+
parameters: &'a Self::Parameters<D>,
41+
) -> HelmholtzEnergyWrapper<'a, Self, D, N> {
42+
HelmholtzEnergyWrapper {
43+
eos: self,
44+
parameters,
45+
}
5846
}
5947
}
6048

6149
/// 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>,
50+
#[derive(Copy)]
51+
pub struct HelmholtzEnergyWrapper<'a, E: ParametersAD, D: DualNum<f64> + Copy, const N: usize> {
52+
pub eos: &'a E,
53+
pub parameters: &'a E::Parameters<D>,
6554
}
6655

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,
56+
impl<'a, E: ParametersAD, D: DualNum<f64> + Copy, const N: usize> Clone
57+
for HelmholtzEnergyWrapper<'a, E, D, N>
58+
{
59+
fn clone(&self) -> Self {
60+
Self {
61+
eos: self.eos,
62+
parameters: self.parameters,
7663
}
7764
}
7865
}
@@ -84,19 +71,17 @@ pub trait NamedParameters: ParametersAD + for<'a> From<&'a [f64]> {
8471
parameters: &'a mut Self::Parameters<D>,
8572
index: &str,
8673
) -> &'a mut D;
87-
}
8874

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>(
75+
/// Return the parameters with the appropriate derivatives.
76+
fn named_derivatives<const P: usize>(
9277
&self,
9378
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();
79+
) -> Self::Parameters<DualSVec<f64, f64, P>> {
80+
let mut params: Self::Parameters<DualSVec<f64, f64, P>> = Self::params_from_inner(self);
9681
for (i, p) in parameters.into_iter().enumerate() {
97-
E::index_parameters_mut(&mut params, p).eps =
82+
Self::index_parameters_mut(&mut params, p).eps =
9883
Derivative::derivative_generic(Const::<P>, U1, i)
9984
}
100-
self.derivatives(params)
85+
params
10186
}
10287
}

‎crates/feos-ad/src/core/phase_equilibria.rs‎

Lines changed: 14 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,12 @@
1-
use super::{HelmholtzEnergyWrapper, ParametersAD, ResidualHelmholtzEnergy, StateAD};
1+
use super::{HelmholtzEnergyWrapper, ParametersAD, ResidualHelmholtzEnergy, StateAD2};
22
use feos_core::{Contributions, FeosResult, PhaseEquilibrium, ReferenceSystem};
33
use nalgebra::SVector;
44
use ndarray::{Array, arr1};
55
use num_dual::{DualNum, linalg::LU};
66
use quantity::{Dimensionless, Moles, Pressure, Temperature};
77

88
impl<'a, R: ResidualHelmholtzEnergy<N>, D: DualNum<f64> + Copy, const N: usize>
9-
StateAD<'a, R, D, N>
9+
StateAD2<'a, R, D, N>
1010
{
1111
/// Perform a Tp-flash calculation. Returns the [PhaseEquilibriumAD] and the vapor fraction.
1212
pub fn tp_flash(&self) -> FeosResult<(PhaseEquilibriumAD<'a, R, D, N>, Dimensionless<D>)> {
@@ -94,8 +94,8 @@ impl<'a, R: ResidualHelmholtzEnergy<N>, D: DualNum<f64> + Copy, const N: usize>
9494
let molefracs_v = rho_v * molar_volume_v;
9595
Ok((
9696
PhaseEquilibriumAD {
97-
liquid: StateAD::new(self.eos, t, molar_volume_l, molefracs_l),
98-
vapor: StateAD::new(self.eos, t, molar_volume_v, molefracs_v),
97+
liquid: StateAD2::new(self.eos, t, molar_volume_l, molefracs_l),
98+
vapor: StateAD2::new(self.eos, t, molar_volume_v, molefracs_v),
9999
},
100100
Dimensionless::from_reduced(v_v / molar_volume_v / self.molefracs.sum()),
101101
))
@@ -104,8 +104,8 @@ impl<'a, R: ResidualHelmholtzEnergy<N>, D: DualNum<f64> + Copy, const N: usize>
104104

105105
/// An equilibrium state consisting of a vapor and a liquid phase.
106106
pub struct PhaseEquilibriumAD<'a, E: ParametersAD, D: DualNum<f64> + Copy, const N: usize> {
107-
pub liquid: StateAD<'a, E, D, N>,
108-
pub vapor: StateAD<'a, E, D, N>,
107+
pub liquid: StateAD2<'a, E, D, N>,
108+
pub vapor: StateAD2<'a, E, D, N>,
109109
}
110110

111111
impl<'a, R: ResidualHelmholtzEnergy<1>, D: DualNum<f64> + Copy> PhaseEquilibriumAD<'a, R, D, 1> {
@@ -132,8 +132,8 @@ impl<'a, R: ResidualHelmholtzEnergy<1>, D: DualNum<f64> + Copy> PhaseEquilibrium
132132
}
133133
Ok((
134134
Self {
135-
liquid: StateAD::new(eos, t, density1.recip(), molefracs),
136-
vapor: StateAD::new(eos, t, density2.recip(), molefracs),
135+
liquid: StateAD2::new(eos, t, density1.recip(), molefracs),
136+
vapor: StateAD2::new(eos, t, density2.recip(), molefracs),
137137
},
138138
Pressure::from_reduced(p),
139139
))
@@ -207,7 +207,7 @@ impl<'a, R: ResidualHelmholtzEnergy<N>, D: DualNum<f64> + Copy, const N: usize>
207207
pressure: f64,
208208
density: f64,
209209
partial_density_other_phase: SVector<f64, N>,
210-
) -> FeosResult<(StateAD<'a, R, D, N>, StateAD<'a, R, D, N>, Pressure<D>)> {
210+
) -> FeosResult<(StateAD2<'a, R, D, N>, StateAD2<'a, R, D, N>, Pressure<D>)> {
211211
let mut rho = SVector::from_fn(|i, _| D::from(partial_density_other_phase[i]));
212212
let mut v = D::from(density.recip());
213213
let t = temperature.into_reduced();
@@ -253,8 +253,8 @@ impl<'a, R: ResidualHelmholtzEnergy<N>, D: DualNum<f64> + Copy, const N: usize>
253253
let v_o = rho.sum().recip();
254254
let molefracs_other_phase = rho * v_o;
255255
Ok((
256-
StateAD::new(eos, t, v, molefracs),
257-
StateAD::new(eos, t, v_o, molefracs_other_phase),
256+
StateAD2::new(eos, t, v, molefracs),
257+
StateAD2::new(eos, t, v_o, molefracs_other_phase),
258258
Pressure::from_reduced(p),
259259
))
260260
}
@@ -545,23 +545,23 @@ mod test_gc_pcsaft {
545545
let temperature = Temperature::from_reduced(Dual::from(300.0));
546546
let pressure = Pressure::from_reduced(Dual::from(0.005));
547547
let molefracs = SVector::from([Dual::from(0.5); 2]);
548-
let (vle_dual, phi_dual) = StateAD::new_tp(
548+
let (vle_dual, phi_dual) = StateAD2::new_tp(
549549
&eos_dual,
550550
temperature,
551551
pressure,
552552
molefracs,
553553
DensityInitialization::None,
554554
)?
555555
.tp_flash()?;
556-
let (vle, phi) = StateAD::new_tp(
556+
let (vle, phi) = StateAD2::new_tp(
557557
&eos,
558558
Temperature::from_reduced(300.0),
559559
Pressure::from_reduced(0.005),
560560
SVector::from([0.5, 0.5]),
561561
DensityInitialization::None,
562562
)?
563563
.tp_flash()?;
564-
let (vle_h, phi_h) = StateAD::new_tp(
564+
let (vle_h, phi_h) = StateAD2::new_tp(
565565
&eos_h,
566566
Temperature::from_reduced(300.0),
567567
Pressure::from_reduced(0.005),

‎crates/feos-ad/src/core/residual.rs‎

Lines changed: 7 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -1,42 +1,20 @@
1-
use super::{FeOsWrapper, HelmholtzEnergyWrapper};
1+
use crate::ParametersAD;
22
use nalgebra::{SMatrix, SVector};
33
use num_dual::{
4-
first_derivative, gradient, hessian, partial_hessian, second_derivative, Dual, Dual2, Dual2Vec,
5-
DualNum, DualVec, HyperDualVec,
4+
Dual, Dual2, Dual2Vec, DualNum, DualVec, HyperDualVec, first_derivative, gradient, hessian,
5+
partial_hessian, second_derivative,
66
};
7-
use std::sync::Arc;
8-
9-
/// A model that can be evaluated with derivatives of its parameters.
10-
pub trait ParametersAD: Send + Sync + Sized {
11-
/// The type of the structure that stores the parameters internally.
12-
type Parameters<D: DualNum<f64> + Copy>: Clone;
13-
14-
/// Return the parameters in the given data type.
15-
fn params<D: DualNum<f64> + Copy>(&self) -> Self::Parameters<D>;
16-
17-
/// Lift the parameters to the given type of dual number.
18-
fn params_from_inner<D: DualNum<f64> + Copy, D2: DualNum<f64, Inner = D> + Copy>(
19-
parameters: &Self::Parameters<D>,
20-
) -> Self::Parameters<D2>;
21-
22-
/// Wraps the model in the [HelmholtzEnergyWrapper] struct, so that it can be used
23-
/// as an argument to [StateAD](crate::StateAD) and [PhaseEquilibriumAD](crate::PhaseEquilibriumAD) constructors.
24-
fn wrap<const N: usize>(self) -> HelmholtzEnergyWrapper<Self, f64, N> {
25-
let parameters = self.params();
26-
HelmholtzEnergyWrapper {
27-
eos: Arc::new(FeOsWrapper(self)),
28-
parameters,
29-
}
30-
}
31-
}
327

338
/// Implementation of a residual Helmholtz energy model.
349
pub trait ResidualHelmholtzEnergy<const N: usize>: ParametersAD {
3510
/// The name of the model.
3611
const RESIDUAL: &str;
3712

3813
/// Return a density (in reduced units) that corresponds to a dense liquid phase.
39-
fn compute_max_density(&self, molefracs: &SVector<f64, N>) -> f64;
14+
fn compute_max_density<D: DualNum<f64> + Copy>(
15+
parameters: &Self::Parameters<D>,
16+
molefracs: &SVector<D, N>,
17+
) -> D;
4018

4119
fn residual_helmholtz_energy_density<D: DualNum<f64> + Copy>(
4220
parameters: &Self::Parameters<D>,

0 commit comments

Comments
 (0)