feos_core/state/
mod.rs

1//! Description of a thermodynamic state.
2//!
3//! A thermodynamic state in SAFT is defined by
4//! * a temperature
5//! * an array of mole numbers
6//! * the volume
7//!
8//! Internally, all properties are computed using such states as input.
9use crate::density_iteration::density_iteration;
10use crate::equation_of_state::Residual;
11use crate::errors::{FeosError, FeosResult};
12use crate::{ReferenceSystem, Total};
13use nalgebra::allocator::Allocator;
14use nalgebra::{DefaultAllocator, Dim, Dyn, OVector, U1};
15use num_dual::*;
16use quantity::*;
17use std::fmt;
18use std::ops::Sub;
19
20mod builder;
21mod cache;
22mod properties;
23mod residual_properties;
24mod statevec;
25pub use builder::StateBuilder;
26pub(crate) use cache::Cache;
27pub use statevec::StateVec;
28
29/// Possible contributions that can be computed.
30#[derive(Clone, Copy, PartialEq)]
31pub enum Contributions {
32    /// Only compute the ideal gas contribution
33    IdealGas,
34    /// Only compute the difference between the total and the ideal gas contribution
35    Residual,
36    /// Compute ideal gas and residual contributions
37    Total,
38}
39
40/// Initial values in a density iteration.
41#[derive(Clone, Copy)]
42pub enum DensityInitialization<D = Density> {
43    /// Calculate a vapor phase by initializing using the ideal gas.
44    Vapor,
45    /// Calculate a liquid phase by using the `max_density`.
46    Liquid,
47    /// Use the given density as initial value.
48    InitialDensity(D),
49}
50
51impl DensityInitialization {
52    pub fn into_reduced(self) -> DensityInitialization<f64> {
53        match self {
54            Self::Vapor => DensityInitialization::Vapor,
55            Self::Liquid => DensityInitialization::Liquid,
56            Self::InitialDensity(d) => DensityInitialization::InitialDensity(d.into_reduced()),
57        }
58    }
59}
60
61/// Thermodynamic state of the system in reduced variables
62/// including their derivatives.
63///
64/// Properties are stored as generalized (hyper) dual numbers which allows
65/// for automatic differentiation.
66#[derive(Clone, Debug)]
67pub struct StateHD<D: DualNum<f64> + Copy, N: Dim = Dyn>
68where
69    DefaultAllocator: Allocator<N>,
70{
71    /// temperature in Kelvin
72    pub temperature: D,
73    // /// volume in Angstrom^3
74    // pub molar_volume: D,
75    /// mole fractions
76    pub molefracs: OVector<D, N>,
77    /// partial number densities in Angstrom^-3
78    pub partial_density: OVector<D, N>,
79}
80
81impl<N: Dim, D: DualNum<f64> + Copy> StateHD<D, N>
82where
83    DefaultAllocator: Allocator<N>,
84{
85    /// Create a new `StateHD` for given temperature, molar volume and composition.
86    pub fn new(temperature: D, molar_volume: D, molefracs: &OVector<D, N>) -> Self {
87        let partial_density = molefracs / molar_volume;
88
89        Self {
90            temperature,
91            molefracs: molefracs.clone(),
92            partial_density,
93        }
94    }
95
96    /// Create a new `StateHD` for given temperature and partial densities
97    pub fn new_density(temperature: D, partial_density: &OVector<D, N>) -> Self {
98        let molefracs = partial_density / partial_density.sum();
99
100        Self {
101            temperature,
102            molefracs,
103            partial_density: partial_density.clone(),
104        }
105    }
106
107    // Since the molefracs can not be reproduced from moles if the density is zero,
108    // this constructor exists specifically for these cases.
109    pub(crate) fn new_virial(temperature: D, density: D, molefracs: &OVector<D, N>) -> Self {
110        let partial_density = molefracs * density;
111        Self {
112            temperature,
113            molefracs: molefracs.clone(),
114            partial_density,
115        }
116    }
117}
118
119/// Thermodynamic state of the system.
120///
121/// The state is always specified by the variables of the Helmholtz energy: volume $V$,
122/// temperature $T$ and mole numbers $N_i$. Additional to these variables, the state saves
123/// properties like the density, that can be calculated directly from the basic variables.
124/// The state also contains a reference to the equation of state used to create the state.
125/// Therefore, it can be used directly to calculate all state properties.
126///
127/// Calculated partial derivatives are cached in the state. Therefore, the second evaluation
128/// of a property like the pressure, does not require a recalculation of the equation of state.
129/// This can be used in situations where both lower and higher order derivatives are required, as
130/// in a calculation of a derivative all lower derivatives have to be calculated internally as well.
131/// Since they are cached it is more efficient to calculate the highest derivatives first.
132/// For example during the calculation of the isochoric heat capacity $c_v$, the entropy and the
133/// Helmholtz energy are calculated as well.
134///
135/// `State` objects are meant to be immutable. If individual fields like `volume` are changed, the
136/// calculations are wrong as the internal fields of the state are not updated.
137///
138/// ## Contents
139///
140/// + [State properties](#state-properties)
141/// + [Mass specific state properties](#mass-specific-state-properties)
142/// + [Transport properties](#transport-properties)
143/// + [Critical points](#critical-points)
144/// + [State constructors](#state-constructors)
145/// + [Stability analysis](#stability-analysis)
146/// + [Flash calculations](#flash-calculations)
147#[derive(Debug)]
148pub struct State<E, N: Dim = Dyn, D: DualNum<f64> + Copy = f64>
149where
150    DefaultAllocator: Allocator<N>,
151{
152    /// Equation of state
153    pub eos: E,
154    /// Temperature $T$
155    pub temperature: Temperature<D>,
156    /// Volume $V$
157    pub volume: Volume<D>,
158    /// Mole numbers $N_i$
159    pub moles: Moles<OVector<D, N>>,
160    /// Total number of moles $N=\sum_iN_i$
161    pub total_moles: Moles<D>,
162    /// Partial densities $\rho_i=\frac{N_i}{V}$
163    pub partial_density: Density<OVector<D, N>>,
164    /// Total density $\rho=\frac{N}{V}=\sum_i\rho_i$
165    pub density: Density<D>,
166    /// Mole fractions $x_i=\frac{N_i}{N}=\frac{\rho_i}{\rho}$
167    pub molefracs: OVector<D, N>,
168    /// Cache
169    cache: Cache<D, N>,
170}
171
172impl<E: Clone, N: Dim, D: DualNum<f64> + Copy> Clone for State<E, N, D>
173where
174    DefaultAllocator: Allocator<N>,
175{
176    fn clone(&self) -> Self {
177        Self {
178            eos: self.eos.clone(),
179            temperature: self.temperature,
180            volume: self.volume,
181            moles: self.moles.clone(),
182            total_moles: self.total_moles,
183            partial_density: self.partial_density.clone(),
184            density: self.density,
185            molefracs: self.molefracs.clone(),
186            cache: self.cache.clone(),
187        }
188    }
189}
190
191impl<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy> fmt::Display for State<E, N, D>
192where
193    DefaultAllocator: Allocator<N>,
194{
195    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
196        if self.eos.components() == 1 {
197            write!(
198                f,
199                "T = {:.5}, ρ = {:.5}",
200                self.temperature.re(),
201                self.density.re()
202            )
203        } else {
204            write!(
205                f,
206                "T = {:.5}, ρ = {:.5}, x = {:.5?}",
207                self.temperature.re(),
208                self.density.re(),
209                self.molefracs.map(|x| x.re()).as_slice()
210            )
211        }
212    }
213}
214
215impl<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy> State<E, N, D>
216where
217    DefaultAllocator: Allocator<N>,
218{
219    /// Return a new `State` given a temperature, an array of mole numbers and a volume.
220    ///
221    /// This function will perform a validation of the given properties, i.e. test for signs
222    /// and if values are finite. It will **not** validate physics, i.e. if the resulting
223    /// densities are below the maximum packing fraction.
224    pub fn new_nvt(
225        eos: &E,
226        temperature: Temperature<D>,
227        volume: Volume<D>,
228        moles: &Moles<OVector<D, N>>,
229    ) -> FeosResult<Self> {
230        let total_moles = moles.sum();
231        let molefracs = (moles / total_moles).into_value();
232        let density = total_moles / volume;
233        validate(temperature, density, &molefracs)?;
234
235        Ok(Self::new_unchecked(
236            eos,
237            temperature,
238            density,
239            total_moles,
240            &molefracs,
241        ))
242    }
243
244    /// Return a new `State` for which the total amount of substance is unspecified.
245    ///
246    /// Internally the total number of moles will be set to 1 mol.
247    ///
248    /// This function will perform a validation of the given properties, i.e. test for signs
249    /// and if values are finite. It will **not** validate physics, i.e. if the resulting
250    /// densities are below the maximum packing fraction.
251    pub fn new_intensive(
252        eos: &E,
253        temperature: Temperature<D>,
254        density: Density<D>,
255        molefracs: &OVector<D, N>,
256    ) -> FeosResult<Self> {
257        validate(temperature, density, molefracs)?;
258        let total_moles = Moles::new(D::one());
259        Ok(Self::new_unchecked(
260            eos,
261            temperature,
262            density,
263            total_moles,
264            molefracs,
265        ))
266    }
267
268    fn new_unchecked(
269        eos: &E,
270        temperature: Temperature<D>,
271        density: Density<D>,
272        total_moles: Moles<D>,
273        molefracs: &OVector<D, N>,
274    ) -> Self {
275        let volume = total_moles / density;
276        let moles = Dimensionless::new(molefracs.clone()) * total_moles;
277        let partial_density = moles.clone() / volume;
278
279        State {
280            eos: eos.clone(),
281            temperature,
282            volume,
283            moles,
284            total_moles,
285            partial_density,
286            density,
287            molefracs: molefracs.clone(),
288            cache: Cache::new(),
289        }
290    }
291
292    /// Return a new `State` for a pure component given a temperature and a density. The moles
293    /// are set to the reference value for each component.
294    ///
295    /// This function will perform a validation of the given properties, i.e. test for signs
296    /// and if values are finite. It will **not** validate physics, i.e. if the resulting
297    /// densities are below the maximum packing fraction.
298    pub fn new_pure(eos: &E, temperature: Temperature<D>, density: Density<D>) -> FeosResult<Self> {
299        let molefracs = OVector::from_element_generic(N::from_usize(1), U1, D::one());
300        Self::new_intensive(eos, temperature, density, &molefracs)
301    }
302
303    /// Return a new `State` for the combination of inputs.
304    ///
305    /// The function attempts to create a new state using the given input values. If the state
306    /// is overdetermined, it will choose a method based on the following hierarchy.
307    /// 1. Create a state non-iteratively from the set of $T$, $V$, $\rho$, $\rho_i$, $N$, $N_i$ and $x_i$.
308    /// 2. Use a density iteration for a given pressure.
309    ///
310    /// The [StateBuilder] provides a convenient way of calling this function without the need to provide
311    /// all the optional input values.
312    ///
313    /// # Errors
314    ///
315    /// When the state cannot be created using the combination of inputs.
316    #[expect(clippy::too_many_arguments)]
317    pub fn new(
318        eos: &E,
319        temperature: Option<Temperature<D>>,
320        volume: Option<Volume<D>>,
321        density: Option<Density<D>>,
322        partial_density: Option<&Density<OVector<D, N>>>,
323        total_moles: Option<Moles<D>>,
324        moles: Option<&Moles<OVector<D, N>>>,
325        molefracs: Option<&OVector<D, N>>,
326        pressure: Option<Pressure<D>>,
327        density_initialization: Option<DensityInitialization>,
328    ) -> FeosResult<Self> {
329        Self::_new(
330            eos,
331            temperature,
332            volume,
333            density,
334            partial_density,
335            total_moles,
336            moles,
337            molefracs,
338            pressure,
339            density_initialization,
340        )?
341        .map_err(|_| FeosError::UndeterminedState(String::from("Missing input parameters.")))
342    }
343
344    #[expect(clippy::too_many_arguments)]
345    #[expect(clippy::type_complexity)]
346    fn _new(
347        eos: &E,
348        temperature: Option<Temperature<D>>,
349        volume: Option<Volume<D>>,
350        density: Option<Density<D>>,
351        partial_density: Option<&Density<OVector<D, N>>>,
352        total_moles: Option<Moles<D>>,
353        moles: Option<&Moles<OVector<D, N>>>,
354        molefracs: Option<&OVector<D, N>>,
355        pressure: Option<Pressure<D>>,
356        density_initialization: Option<DensityInitialization>,
357    ) -> FeosResult<Result<Self, Option<Moles<OVector<D, N>>>>> {
358        // check for density
359        if density.and(partial_density).is_some() {
360            return Err(FeosError::UndeterminedState(String::from(
361                "Both density and partial density given.",
362            )));
363        }
364        let rho = density.or_else(|| partial_density.map(|pd| pd.sum()));
365
366        // check for total moles
367        if moles.and(total_moles).is_some() {
368            return Err(FeosError::UndeterminedState(String::from(
369                "Both moles and total moles given.",
370            )));
371        }
372        let mut n = total_moles.or_else(|| moles.map(|m| m.sum()));
373
374        // check if total moles can be inferred from volume
375        if rho.and(n).and(volume).is_some() {
376            return Err(FeosError::UndeterminedState(String::from(
377                "Density is overdetermined.",
378            )));
379        }
380        n = n.or_else(|| rho.and_then(|d| volume.map(|v| v * d)));
381
382        // check for composition
383        if partial_density.and(moles).is_some() {
384            return Err(FeosError::UndeterminedState(String::from(
385                "Composition is overdetermined.",
386            )));
387        }
388        let x = partial_density
389            .map(|pd| pd / pd.sum())
390            .or_else(|| moles.map(|ms| ms / ms.sum()))
391            .map(Quantity::into_value);
392        let x_u = match (x, molefracs, eos.components()) {
393            (Some(_), Some(_), _) => {
394                return Err(FeosError::UndeterminedState(String::from(
395                    "Composition is overdetermined.",
396                )));
397            }
398            (Some(x), None, _) => x,
399            (None, Some(x), _) => x.clone(),
400            (None, None, 1) => OVector::from_element_generic(N::from_usize(1), U1, D::from(1.0)),
401            _ => {
402                return Err(FeosError::UndeterminedState(String::from(
403                    "Missing composition.",
404                )));
405            }
406        };
407        let x_u = &x_u / x_u.sum();
408
409        // If no extensive property is given, moles is set to the reference value.
410        if let (None, None) = (volume, n) {
411            n = Some(Moles::from_reduced(D::from(1.0)))
412        }
413        let n_i = n.map(|n| Dimensionless::new(&x_u) * n);
414        let v = volume.or_else(|| rho.and_then(|d| n.map(|n| n / d)));
415
416        // check if new state can be created using default constructor
417        if let (Some(v), Some(t), Some(n_i)) = (v, temperature, &n_i) {
418            return Ok(Ok(State::new_nvt(eos, t, v, n_i)?));
419        }
420
421        // Check if new state can be created using density iteration
422        if let (Some(p), Some(t), Some(n_i)) = (pressure, temperature, &n_i) {
423            return Ok(Ok(State::new_npt(eos, t, p, n_i, density_initialization)?));
424        }
425        if let (Some(p), Some(t), Some(v)) = (pressure, temperature, v) {
426            return Ok(Ok(State::new_npvx(
427                eos,
428                t,
429                p,
430                v,
431                &x_u,
432                density_initialization,
433            )?));
434        }
435        Ok(Err(n_i.to_owned()))
436    }
437
438    /// Return a new `State` using a density iteration. [DensityInitialization] is used to
439    /// influence the calculation with respect to the possible solutions.
440    pub fn new_npt(
441        eos: &E,
442        temperature: Temperature<D>,
443        pressure: Pressure<D>,
444        moles: &Moles<OVector<D, N>>,
445        density_initialization: Option<DensityInitialization>,
446    ) -> FeosResult<Self> {
447        let total_moles = moles.sum();
448        let molefracs = (moles / total_moles).into_value();
449        let density = Self::new_xpt(
450            eos,
451            temperature,
452            pressure,
453            &molefracs,
454            density_initialization,
455        )?
456        .density;
457        Ok(Self::new_unchecked(
458            eos,
459            temperature,
460            density,
461            total_moles,
462            &molefracs,
463        ))
464    }
465
466    /// Return a new `State` using a density iteration. [DensityInitialization] is used to
467    /// influence the calculation with respect to the possible solutions.
468    pub fn new_xpt(
469        eos: &E,
470        temperature: Temperature<D>,
471        pressure: Pressure<D>,
472        molefracs: &OVector<D, N>,
473        density_initialization: Option<DensityInitialization>,
474    ) -> FeosResult<Self> {
475        density_iteration(
476            eos,
477            temperature,
478            pressure,
479            molefracs,
480            density_initialization,
481        )
482        .and_then(|density| Self::new_intensive(eos, temperature, density, molefracs))
483    }
484
485    /// Return a new `State` for given pressure $p$, volume $V$, temperature $T$ and composition $x_i$.
486    pub fn new_npvx(
487        eos: &E,
488        temperature: Temperature<D>,
489        pressure: Pressure<D>,
490        volume: Volume<D>,
491        molefracs: &OVector<D, N>,
492        density_initialization: Option<DensityInitialization>,
493    ) -> FeosResult<Self> {
494        let density = Self::new_xpt(
495            eos,
496            temperature,
497            pressure,
498            molefracs,
499            density_initialization,
500        )?
501        .density;
502        let total_moles = density * volume;
503        Ok(Self::new_unchecked(
504            eos,
505            temperature,
506            density,
507            total_moles,
508            molefracs,
509        ))
510    }
511}
512
513impl<E: Total<N, D>, N: Gradients, D: DualNum<f64> + Copy> State<E, N, D>
514where
515    DefaultAllocator: Allocator<N>,
516{
517    /// Return a new `State` for the combination of inputs.
518    ///
519    /// The function attempts to create a new state using the given input values. If the state
520    /// is overdetermined, it will choose a method based on the following hierarchy.
521    /// 1. Create a state non-iteratively from the set of $T$, $V$, $\rho$, $\rho_i$, $N$, $N_i$ and $x_i$.
522    /// 2. Use a density iteration for a given pressure.
523    /// 3. Determine the state using a Newton iteration from (in this order): $(p, h)$, $(p, s)$, $(T, h)$, $(T, s)$, $(V, u)$
524    ///
525    /// The [StateBuilder] provides a convenient way of calling this function without the need to provide
526    /// all the optional input values.
527    ///
528    /// # Errors
529    ///
530    /// When the state cannot be created using the combination of inputs.
531    #[expect(clippy::too_many_arguments)]
532    pub fn new_full(
533        eos: &E,
534        temperature: Option<Temperature<D>>,
535        volume: Option<Volume<D>>,
536        density: Option<Density<D>>,
537        partial_density: Option<&Density<OVector<D, N>>>,
538        total_moles: Option<Moles<D>>,
539        moles: Option<&Moles<OVector<D, N>>>,
540        molefracs: Option<&OVector<D, N>>,
541        pressure: Option<Pressure<D>>,
542        molar_enthalpy: Option<MolarEnergy<D>>,
543        molar_entropy: Option<MolarEntropy<D>>,
544        molar_internal_energy: Option<MolarEnergy<D>>,
545        density_initialization: Option<DensityInitialization>,
546        initial_temperature: Option<Temperature<D>>,
547    ) -> FeosResult<Self> {
548        let state = Self::_new(
549            eos,
550            temperature,
551            volume,
552            density,
553            partial_density,
554            total_moles,
555            moles,
556            molefracs,
557            pressure,
558            density_initialization,
559        )?;
560
561        let ti = initial_temperature;
562        match state {
563            Ok(state) => Ok(state),
564            Err(n_i) => {
565                // Check if new state can be created using molar_enthalpy and temperature
566                if let (Some(p), Some(h), Some(n_i)) = (pressure, molar_enthalpy, &n_i) {
567                    return State::new_nph(eos, p, h, n_i, density_initialization, ti);
568                }
569                if let (Some(p), Some(s), Some(n_i)) = (pressure, molar_entropy, &n_i) {
570                    return State::new_nps(eos, p, s, n_i, density_initialization, ti);
571                }
572                if let (Some(t), Some(h), Some(n_i)) = (temperature, molar_enthalpy, &n_i) {
573                    return State::new_nth(eos, t, h, n_i, density_initialization);
574                }
575                if let (Some(t), Some(s), Some(n_i)) = (temperature, molar_entropy, &n_i) {
576                    return State::new_nts(eos, t, s, n_i, density_initialization);
577                }
578                if let (Some(u), Some(v), Some(n_i)) = (molar_internal_energy, volume, &n_i) {
579                    return State::new_nvu(eos, v, u, n_i, ti);
580                }
581                Err(FeosError::UndeterminedState(String::from(
582                    "Missing input parameters.",
583                )))
584            }
585        }
586    }
587
588    /// Return a new `State` for given pressure $p$ and molar enthalpy $h$.
589    pub fn new_nph(
590        eos: &E,
591        pressure: Pressure<D>,
592        molar_enthalpy: MolarEnergy<D>,
593        moles: &Moles<OVector<D, N>>,
594        density_initialization: Option<DensityInitialization>,
595        initial_temperature: Option<Temperature<D>>,
596    ) -> FeosResult<Self> {
597        let t0 = initial_temperature.unwrap_or(Temperature::from_reduced(D::from(298.15)));
598        let mut density = density_initialization;
599        let f = |x0| {
600            let s = State::new_npt(eos, x0, pressure, moles, density)?;
601            let dfx = s.molar_isobaric_heat_capacity(Contributions::Total);
602            let fx = s.molar_enthalpy(Contributions::Total) - molar_enthalpy;
603            density = Some(DensityInitialization::InitialDensity(s.density.re()));
604            Ok((fx, dfx, s))
605        };
606        newton(t0, f, Temperature::from_reduced(1.0e-8))
607    }
608
609    /// Return a new `State` for given temperature $T$ and molar enthalpy $h$.
610    pub fn new_nth(
611        eos: &E,
612        temperature: Temperature<D>,
613        molar_enthalpy: MolarEnergy<D>,
614        moles: &Moles<OVector<D, N>>,
615        density_initialization: Option<DensityInitialization>,
616    ) -> FeosResult<Self> {
617        let x = moles.convert_to(moles.sum());
618        let rho0 = match density_initialization {
619            Some(DensityInitialization::InitialDensity(r)) => {
620                Density::from_reduced(D::from(r.into_reduced()))
621            }
622            Some(DensityInitialization::Liquid) => eos.max_density(&Some(x))?,
623            Some(DensityInitialization::Vapor) => eos.max_density(&Some(x))? * 1.0e-5,
624            None => eos.max_density(&Some(x))? * 0.01,
625        };
626        let n_inv = moles.sum().inv();
627        let f = |x0| {
628            let s = State::new_nvt(eos, temperature, moles.sum() / x0, moles)?;
629            let dfx = -s.volume / s.density
630                * n_inv
631                * (s.volume * s.dp_dv(Contributions::Total)
632                    + temperature * s.dp_dt(Contributions::Total));
633            let fx = s.molar_enthalpy(Contributions::Total) - molar_enthalpy;
634            Ok((fx, dfx, s))
635        };
636        newton(rho0, f, Density::from_reduced(1.0e-12))
637    }
638
639    /// Return a new `State` for given temperature $T$ and molar entropy $s$.
640    pub fn new_nts(
641        eos: &E,
642        temperature: Temperature<D>,
643        molar_entropy: MolarEntropy<D>,
644        moles: &Moles<OVector<D, N>>,
645        density_initialization: Option<DensityInitialization>,
646    ) -> FeosResult<Self> {
647        let x = moles.convert_to(moles.sum());
648        let rho0 = match density_initialization {
649            Some(DensityInitialization::InitialDensity(r)) => {
650                Density::from_reduced(D::from(r.into_reduced()))
651            }
652            Some(DensityInitialization::Liquid) => eos.max_density(&Some(x))?,
653            Some(DensityInitialization::Vapor) => eos.max_density(&Some(x))? * 1.0e-5,
654            None => eos.max_density(&Some(x))? * 0.01,
655        };
656        let n_inv = moles.sum().inv();
657        let f = |x0| {
658            let s = State::new_nvt(eos, temperature, moles.sum() / x0, moles)?;
659            let dfx = -n_inv * s.volume / s.density * s.dp_dt(Contributions::Total);
660            let fx = s.molar_entropy(Contributions::Total) - molar_entropy;
661            Ok((fx, dfx, s))
662        };
663        newton(rho0, f, Density::from_reduced(1.0e-12))
664    }
665
666    /// Return a new `State` for given pressure $p$ and molar entropy $s$.
667    pub fn new_nps(
668        eos: &E,
669        pressure: Pressure<D>,
670        molar_entropy: MolarEntropy<D>,
671        moles: &Moles<OVector<D, N>>,
672        density_initialization: Option<DensityInitialization>,
673        initial_temperature: Option<Temperature<D>>,
674    ) -> FeosResult<Self> {
675        let t0 = initial_temperature.unwrap_or(Temperature::from_reduced(D::from(298.15)));
676        let mut density = density_initialization;
677        let f = |x0| {
678            let s = State::new_npt(eos, x0, pressure, moles, density)?;
679            let dfx = s.molar_isobaric_heat_capacity(Contributions::Total) / s.temperature;
680            let fx = s.molar_entropy(Contributions::Total) - molar_entropy;
681            density = Some(DensityInitialization::InitialDensity(s.density.re()));
682            Ok((fx, dfx, s))
683        };
684        newton(t0, f, Temperature::from_reduced(1.0e-8))
685    }
686
687    /// Return a new `State` for given volume $V$ and molar internal energy $u$.
688    pub fn new_nvu(
689        eos: &E,
690        volume: Volume<D>,
691        molar_internal_energy: MolarEnergy<D>,
692        moles: &Moles<OVector<D, N>>,
693        initial_temperature: Option<Temperature<D>>,
694    ) -> FeosResult<Self> {
695        let t0 = initial_temperature.unwrap_or(Temperature::from_reduced(D::from(298.15)));
696        let f = |x0| {
697            let s = State::new_nvt(eos, x0, volume, moles)?;
698            let fx = s.molar_internal_energy(Contributions::Total) - molar_internal_energy;
699            let dfx = s.molar_isochoric_heat_capacity(Contributions::Total);
700            Ok((fx, dfx, s))
701        };
702        newton(t0, f, Temperature::from_reduced(1.0e-8))
703    }
704}
705
706fn is_close<U: Copy>(
707    x: Quantity<f64, U>,
708    y: Quantity<f64, U>,
709    atol: Quantity<f64, U>,
710    rtol: f64,
711) -> bool {
712    (x - y).abs() <= atol + rtol * y.abs()
713}
714
715fn newton<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy, F, X: Copy, Y>(
716    mut x0: Quantity<D, X>,
717    mut f: F,
718    atol: Quantity<f64, X>,
719) -> FeosResult<State<E, N, D>>
720where
721    DefaultAllocator: Allocator<N>,
722    Y: Sub<X> + Sub<<Y as Sub<X>>::Output, Output = X>,
723    F: FnMut(
724        Quantity<D, X>,
725    ) -> FeosResult<(
726        Quantity<D, Y>,
727        Quantity<D, <Y as Sub<X>>::Output>,
728        State<E, N, D>,
729    )>,
730{
731    let rtol = 1e-10;
732    let maxiter = 50;
733
734    for _ in 0..maxiter {
735        let (fx, dfx, mut state) = f(x0)?;
736        let x = x0 - fx / dfx;
737        if is_close(x.re(), x0.re(), atol, rtol) {
738            // Ensure that at least NDERIV iterations are performed (for implicit AD)
739            for _ in 0..D::NDERIV {
740                let (fx, dfx, s) = f(x0)?;
741                x0 -= fx / dfx;
742                state = s;
743            }
744            return Ok(state);
745        }
746        x0 = x;
747    }
748    Err(FeosError::NotConverged("newton".to_owned()))
749}
750
751/// Validate the given temperature, mole numbers and volume.
752///
753/// Properties are valid if
754/// * they are finite
755/// * they have a positive sign
756///
757/// There is no validation of the physical state, e.g.
758/// if resulting densities are below maximum packing fraction.
759fn validate<N: Dim, D: DualNum<f64>>(
760    temperature: Temperature<D>,
761    density: Density<D>,
762    molefracs: &OVector<D, N>,
763) -> FeosResult<()>
764where
765    DefaultAllocator: Allocator<N>,
766{
767    let t = temperature.re().to_reduced();
768    let rho = density.re().to_reduced();
769    if !t.is_finite() || t.is_sign_negative() {
770        return Err(FeosError::InvalidState(
771            String::from("validate"),
772            String::from("temperature"),
773            t,
774        ));
775    }
776    if !rho.is_finite() || rho.is_sign_negative() {
777        return Err(FeosError::InvalidState(
778            String::from("validate"),
779            String::from("density"),
780            rho,
781        ));
782    }
783    for n in molefracs.iter() {
784        if !n.re().is_finite() || n.re().is_sign_negative() {
785            return Err(FeosError::InvalidState(
786                String::from("validate"),
787                String::from("molefracs"),
788                n.re(),
789            ));
790        }
791    }
792    Ok(())
793}
794
795mod critical_point;
796
797#[cfg(test)]
798mod tests {
799    use super::*;
800    use nalgebra::dvector;
801    use typenum::P3;
802
803    #[test]
804    fn test_validate() {
805        let temperature = 298.15 * KELVIN;
806        let density = 3000.0 * MOL / METER.powi::<P3>();
807        let molefracs = dvector![0.03, 0.02, 0.05];
808        assert!(validate(temperature, density, &molefracs).is_ok());
809    }
810
811    #[test]
812    fn test_negative_temperature() {
813        let temperature = -298.15 * KELVIN;
814        let density = 3000.0 * MOL / METER.powi::<P3>();
815        let molefracs = dvector![0.03, 0.02, 0.05];
816        assert!(validate(temperature, density, &molefracs).is_err());
817    }
818
819    #[test]
820    fn test_nan_temperature() {
821        let temperature = f64::NAN * KELVIN;
822        let density = 3000.0 * MOL / METER.powi::<P3>();
823        let molefracs = dvector![0.03, 0.02, 0.05];
824        assert!(validate(temperature, density, &molefracs).is_err());
825    }
826
827    #[test]
828    fn test_negative_mole_number() {
829        let temperature = 298.15 * KELVIN;
830        let density = 3000.0 * MOL / METER.powi::<P3>();
831        let molefracs = dvector![-0.03, 0.02, 0.05];
832        assert!(validate(temperature, density, &molefracs).is_err());
833    }
834
835    #[test]
836    fn test_nan_mole_number() {
837        let temperature = 298.15 * KELVIN;
838        let density = 3000.0 * MOL / METER.powi::<P3>();
839        let molefracs = dvector![f64::NAN, 0.02, 0.05];
840        assert!(validate(temperature, density, &molefracs).is_err());
841    }
842
843    #[test]
844    fn test_negative_density() {
845        let temperature = 298.15 * KELVIN;
846        let density = -3000.0 * MOL / METER.powi::<P3>();
847        let molefracs = dvector![0.01, 0.02, 0.05];
848        assert!(validate(temperature, density, &molefracs).is_err());
849    }
850}
feos_core/state/mod.rs

feos_core/state/
mod.rs
