Skip to content

Commit 55a40ef

Browse files
prehnerprehner
committed
Enable AD for pure-component VLEs with given pressure
1 parent 324d986 commit 55a40ef

File tree

8 files changed

+411
-221
lines changed

8 files changed

+411
-221
lines changed

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

Lines changed: 65 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ use crate::{FeosResult, PhaseEquilibrium, ReferenceSystem, Residual};
44
use nalgebra::{Const, SVector, U1, U2};
55
#[cfg(feature = "rayon")]
66
use ndarray::{Array1, Array2, ArrayView2, Zip};
7-
use num_dual::{Derivative, DualSVec, DualStruct};
7+
use num_dual::{Derivative, DualNum, DualSVec, DualStruct, first_derivative, partial2};
88
use quantity::{Density, Pressure, Temperature};
99
#[cfg(feature = "rayon")]
1010
use quantity::{KELVIN, KILO, METER, MOL, PASCAL};
@@ -66,6 +66,50 @@ pub trait PropertiesAD {
6666
Ok(Pressure::from_reduced(p))
6767
}
6868

69+
fn boiling_temperature<const P: usize>(
70+
&self,
71+
pressure: Pressure,
72+
) -> FeosResult<Temperature<Gradient<P>>>
73+
where
74+
Self: Residual<U1, Gradient<P>>,
75+
{
76+
let eos_f64 = self.re();
77+
let (temperature, [vapor_density, liquid_density]) =
78+
PhaseEquilibrium::pure_p(&eos_f64, pressure, None, Default::default())?;
79+
80+
// implicit differentiation is implemented here instead of just calling pure_t with dual
81+
// numbers, because for the first derivative, we can avoid calculating density derivatives.
82+
let t = temperature.into_reduced();
83+
let v1 = 1.0 / liquid_density.to_reduced();
84+
let v2 = 1.0 / vapor_density.to_reduced();
85+
let p = pressure.into_reduced();
86+
let t = Gradient::from(t);
87+
let t = t + {
88+
let v1 = Gradient::from(v1);
89+
let v2 = Gradient::from(v2);
90+
let p = Gradient::from(p);
91+
let x = Self::pure_molefracs();
92+
93+
let residual_entropy = |v| {
94+
let (a, s) = first_derivative(
95+
partial2(
96+
|t, &v, x| self.lift().residual_molar_helmholtz_energy(t, v, x),
97+
&v,
98+
&x,
99+
),
100+
t,
101+
);
102+
(a, -s)
103+
};
104+
let (a1, s1) = residual_entropy(v1);
105+
let (a2, s2) = residual_entropy(v2);
106+
107+
let ln_rho = (v1 / v2).ln();
108+
(p * (v2 - v1) + (a2 - a1 + t * ln_rho)) / (s2 - s1 - ln_rho)
109+
};
110+
Ok(Temperature::from_reduced(t))
111+
}
112+
69113
fn equilibrium_liquid_density<const P: usize>(
70114
&self,
71115
temperature: Temperature,
@@ -111,6 +155,26 @@ pub trait PropertiesAD {
111155
)
112156
}
113157

158+
#[cfg(feature = "rayon")]
159+
fn boiling_temperature_parallel<const P: usize>(
160+
parameter_names: [String; P],
161+
parameters: ArrayView2<f64>,
162+
input: ArrayView2<f64>,
163+
) -> (Array1<f64>, Array2<f64>, Array1<bool>)
164+
where
165+
Self: ParametersAD<1>,
166+
{
167+
parallelize::<_, Self, _, _>(
168+
parameter_names,
169+
parameters,
170+
input,
171+
|eos: &Self::Lifted<Gradient<P>>, inp| {
172+
eos.boiling_temperature(inp[0] * PASCAL)
173+
.map(|p| p.convert_into(KELVIN))
174+
},
175+
)
176+
}
177+
114178
#[cfg(feature = "rayon")]
115179
fn liquid_density_parallel<const P: usize>(
116180
parameter_names: [String; P],

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

Lines changed: 34 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,9 @@
11
use crate::{FeosError, FeosResult, ReferenceSystem, state::StateHD};
2+
use nalgebra::SVector;
23
use nalgebra::{DVector, DefaultAllocator, Dim, Dyn, OMatrix, OVector, U1, allocator::Allocator};
3-
use num_dual::{DualNum, Gradients, partial, partial2, second_derivative, third_derivative};
4+
use num_dual::{
5+
DualNum, Gradients, hessian, partial, partial2, second_derivative, third_derivative,
6+
};
47
use quantity::ad::first_derivative;
58
use quantity::*;
69
use std::ops::{Deref, Div};
@@ -313,12 +316,41 @@ where
313316
molar_volume,
314317
);
315318
(
316-
a * density,
319+
a,
317320
-da + temperature * density,
318321
molar_volume * molar_volume * d2a + temperature,
319322
)
320323
}
321324

325+
/// calculates a_res, p, s_res, dp_drho, dp_dt
326+
fn p_dpdrho_dpdt(
327+
&self,
328+
temperature: D,
329+
density: D,
330+
molefracs: &OVector<D, N>,
331+
) -> (D, D, D, D, D) {
332+
let molar_volume = density.recip();
333+
let (a, da, d2a) = hessian::<_, _, _, nalgebra::U2, _>(
334+
partial(
335+
|vt: SVector<_, 2>, x: &OVector<_, N>| {
336+
let [[v, t]] = vt.data.0;
337+
self.lift().residual_molar_helmholtz_energy(t, v, x)
338+
},
339+
molefracs,
340+
),
341+
&SVector::from([molar_volume, temperature]),
342+
);
343+
let [[da_dv, da_dt]] = da.data.0;
344+
let [[d2a_dv2, d2a_dvdt], _] = d2a.data.0;
345+
(
346+
a,
347+
-da_dv + temperature * density,
348+
-da_dt,
349+
molar_volume * molar_volume * d2a_dv2 + temperature,
350+
-d2a_dvdt + density,
351+
)
352+
}
353+
322354
/// calculates p, dp_drho, d2p_drho2
323355
fn p_dpdrho_d2pdrho2(
324356
&self,

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

Lines changed: 1 addition & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
use crate::equation_of_state::Residual;
2-
use crate::errors::{FeosError, FeosResult};
2+
use crate::errors::FeosResult;
33
use crate::state::{DensityInitialization, State};
44
use crate::{Contributions, ReferenceSystem};
55
use nalgebra::allocator::Allocator;
@@ -44,12 +44,6 @@ pub struct PhaseEquilibrium<E, const P: usize, N: Dim = Dyn, D: DualNum<f64> + C
4444
where
4545
DefaultAllocator: Allocator<N>;
4646

47-
// impl<E, const P: usize> Clone for PhaseEquilibrium<E, P> {
48-
// fn clone(&self) -> Self {
49-
// Self(self.0.clone())
50-
// }
51-
// }
52-
5347
impl<E: Residual, const P: usize> fmt::Display for PhaseEquilibrium<E, P> {
5448
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
5549
for (i, s) in self.0.iter().enumerate() {
@@ -224,14 +218,6 @@ impl<E: Residual<N>, N: Dim> PhaseEquilibrium<E, 2, N>
224218
where
225219
DefaultAllocator: Allocator<N>,
226220
{
227-
pub(super) fn check_trivial_solution(self) -> FeosResult<Self> {
228-
if Self::is_trivial_solution(self.vapor(), self.liquid()) {
229-
Err(FeosError::TrivialSolution)
230-
} else {
231-
Ok(self)
232-
}
233-
}
234-
235221
/// Check if the two states form a trivial solution
236222
pub fn is_trivial_solution(state1: &State<E, N>, state2: &State<E, N>) -> bool {
237223
let rho1 = state1.molefracs.clone() * state1.density.into_reduced();

0 commit comments

Comments
 (0)