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

feos_core/state/
mod.rs
