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};
15use num_dual::*;
16use quantity::*;
17use std::fmt;
18use std::ops::Sub;
19
20mod cache;
21mod composition;
22mod properties;
23mod residual_properties;
24mod statevec;
25pub(crate) use cache::Cache;
26pub use composition::Composition;
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    /// mole fractions
74    pub molefracs: OVector<D, N>,
75    /// partial number densities in Angstrom^-3
76    pub partial_density: OVector<D, N>,
77}
78
79impl<N: Dim, D: DualNum<f64> + Copy> StateHD<D, N>
80where
81    DefaultAllocator: Allocator<N>,
82{
83    /// Create a new `StateHD` for given temperature, volume and composition.
84    pub fn new(temperature: D, volume: D, moles: &OVector<D, N>) -> Self {
85        Self::new_density(temperature, &(moles / volume))
86    }
87
88    /// Create a new `StateHD` for given temperature and partial densities
89    pub fn new_density(temperature: D, partial_density: &OVector<D, N>) -> Self {
90        let molefracs = partial_density / partial_density.sum();
91
92        Self {
93            temperature,
94            molefracs,
95            partial_density: partial_density.clone(),
96        }
97    }
98
99    // Since the molefracs can not be reproduced from moles if the density is zero,
100    // this constructor exists specifically for these cases.
101    pub(crate) fn new_virial(temperature: D, density: D, molefracs: &OVector<D, N>) -> Self {
102        let partial_density = molefracs * density;
103        Self {
104            temperature,
105            molefracs: molefracs.clone(),
106            partial_density,
107        }
108    }
109}
110
111/// Thermodynamic state of the system.
112///
113/// The state is always specified by the variables of the Helmholtz energy: volume $V$,
114/// temperature $T$ and mole numbers $N_i$. Additional to these variables, the state saves
115/// properties like the density, that can be calculated directly from the basic variables.
116/// The state also contains a reference to the equation of state used to create the state.
117/// Therefore, it can be used directly to calculate all state properties.
118///
119/// Calculated partial derivatives are cached in the state. Therefore, the second evaluation
120/// of a property like the pressure, does not require a recalculation of the equation of state.
121/// This can be used in situations where both lower and higher order derivatives are required, as
122/// in a calculation of a derivative all lower derivatives have to be calculated internally as well.
123/// Since they are cached it is more efficient to calculate the highest derivatives first.
124/// For example during the calculation of the isochoric heat capacity $c_v$, the entropy and the
125/// Helmholtz energy are calculated as well.
126///
127/// `State` objects are meant to be immutable. If individual fields like `volume` are changed, the
128/// calculations are wrong as the internal fields of the state are not updated.
129///
130/// ## Contents
131///
132/// + [State properties](#state-properties)
133/// + [Mass specific state properties](#mass-specific-state-properties)
134/// + [Transport properties](#transport-properties)
135/// + [Critical points](#critical-points)
136/// + [State constructors](#state-constructors)
137/// + [Stability analysis](#stability-analysis)
138/// + [Flash calculations](#flash-calculations)
139#[derive(Debug, Clone)]
140pub struct State<E, N: Dim = Dyn, D: DualNum<f64> + Copy = f64>
141where
142    DefaultAllocator: Allocator<N>,
143{
144    /// Equation of state
145    pub eos: E,
146    /// Temperature $T$
147    pub temperature: Temperature<D>,
148    /// Molar volume $v=\frac{V}{N}$
149    pub molar_volume: MolarVolume<D>,
150    /// Total number of moles $N=\sum_iN_i$
151    pub total_moles: Option<Moles<D>>,
152    /// Total density $\rho=\frac{N}{V}=\sum_i\rho_i$
153    pub density: Density<D>,
154    /// Mole fractions $x_i=\frac{N_i}{N}=\frac{\rho_i}{\rho}$
155    pub molefracs: OVector<D, N>,
156    /// Cache
157    cache: Cache<D, N>,
158}
159
160impl<E, N: Dim, D: DualNum<f64> + Copy> State<E, N, D>
161where
162    DefaultAllocator: Allocator<N>,
163{
164    /// Set the total amount of substance to the given value.
165    ///
166    /// This method does not introduce inconsistencies, because the
167    /// total moles are the only field that stores information about
168    /// the size of the state.
169    pub fn set_total_moles(mut self, total_moles: Moles<D>) -> State<E, N, D> {
170        self.total_moles = Some(total_moles);
171        self
172    }
173
174    /// Partial densities $\rho_i=\frac{N_i}{V}$
175    pub fn partial_density(&self) -> Density<OVector<D, N>> {
176        Dimensionless::new(&self.molefracs) * self.density
177    }
178
179    /// Mole numbers $N_i$
180    pub fn moles(&self) -> FeosResult<Moles<OVector<D, N>>> {
181        Ok(Dimensionless::new(&self.molefracs) * self.total_moles()?)
182    }
183
184    /// Total moles $N=\sum_iN_i$
185    pub fn total_moles(&self) -> FeosResult<Moles<D>> {
186        self.total_moles.ok_or(FeosError::IntensiveState)
187    }
188
189    /// Volume $V$
190    pub fn volume(&self) -> FeosResult<Volume<D>> {
191        Ok(self.molar_volume * self.total_moles()?)
192    }
193}
194
195impl<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy> fmt::Display for State<E, N, D>
196where
197    DefaultAllocator: Allocator<N>,
198{
199    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
200        if self.eos.components() == 1 {
201            write!(
202                f,
203                "T = {:.5}, ρ = {:.5}",
204                self.temperature.re(),
205                self.density.re()
206            )
207        } else {
208            write!(
209                f,
210                "T = {:.5}, ρ = {:.5}, x = {:.5?}",
211                self.temperature.re(),
212                self.density.re(),
213                self.molefracs.map(|x| x.re()).as_slice()
214            )
215        }
216    }
217}
218
219impl<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy> State<E, N, D>
220where
221    DefaultAllocator: Allocator<N>,
222{
223    /// Return a new `State` given a temperature, an array of mole numbers and a volume.
224    ///
225    /// This function will perform a validation of the given properties, i.e. test for signs
226    /// and if values are finite. It will **not** validate physics, i.e. if the resulting
227    /// densities are below the maximum packing fraction.
228    pub fn new_nvt<X: Composition<D, N>>(
229        eos: &E,
230        temperature: Temperature<D>,
231        volume: Volume<D>,
232        composition: X,
233    ) -> FeosResult<Self> {
234        let (molefracs, total_moles) = composition.into_molefracs(eos)?;
235        let Some(total_moles) = total_moles else {
236            return Err(FeosError::UndeterminedState(
237                "Missing total mole number in the specification!".into(),
238            ));
239        };
240
241        let density = total_moles / volume;
242        Self::new(eos, temperature, density, (molefracs, total_moles))
243    }
244
245    /// Return a new `State` given a temperature and the partial density of all components.
246    ///
247    /// This function will perform a validation of the given properties, i.e. test for signs
248    /// and if values are finite. It will **not** validate physics, i.e. if the resulting
249    /// densities are below the maximum packing fraction.
250    pub fn new_density(
251        eos: &E,
252        temperature: Temperature<D>,
253        partial_density: Density<OVector<D, N>>,
254    ) -> FeosResult<Self> {
255        let density = partial_density.sum();
256        let molefracs = partial_density.convert_into(density);
257        Self::new(eos, temperature, density, molefracs)
258    }
259
260    /// Return a new `State` for a pure component given a temperature and a density.
261    ///
262    /// This function will perform a validation of the given properties, i.e. test for signs
263    /// and if values are finite. It will **not** validate physics, i.e. if the resulting
264    /// densities are below the maximum packing fraction.
265    pub fn new_pure(eos: &E, temperature: Temperature<D>, density: Density<D>) -> FeosResult<Self>
266    where
267        (): Composition<D, N>,
268    {
269        Self::new(eos, temperature, density, ())
270    }
271
272    /// Return a new `State` given a temperature, a density and the composition.
273    ///
274    /// This function will perform a validation of the given properties, i.e. test for signs
275    /// and if values are finite. It will **not** validate physics, i.e. if the resulting
276    /// densities are below the maximum packing fraction.
277    pub fn new<X: Composition<D, N>>(
278        eos: &E,
279        temperature: Temperature<D>,
280        density: Density<D>,
281        composition: X,
282    ) -> FeosResult<Self> {
283        let (molefracs, total_moles) = composition.into_molefracs(eos)?;
284        Self::_new(eos, temperature, density, molefracs, total_moles)
285    }
286
287    fn _new(
288        eos: &E,
289        temperature: Temperature<D>,
290        density: Density<D>,
291        molefracs: OVector<D, N>,
292        total_moles: Option<Moles<D>>,
293    ) -> FeosResult<Self> {
294        let molar_volume = density.inv();
295        validate(temperature, density, &molefracs)?;
296        Ok(State {
297            eos: eos.clone(),
298            temperature,
299            molar_volume,
300            density,
301            molefracs,
302            total_moles,
303            cache: Cache::new(),
304        })
305    }
306
307    /// Return a new `State` for the combination of inputs.
308    ///
309    /// # Errors
310    ///
311    /// When the state cannot be created using the combination of inputs (is over- or underdetermined).
312    pub fn build<X: Composition<D, N>>(
313        eos: &E,
314        temperature: Temperature<D>,
315        volume: Option<Volume<D>>,
316        density: Option<Density<D>>,
317        composition: X,
318        pressure: Option<Pressure<D>>,
319        density_initialization: Option<DensityInitialization>,
320    ) -> FeosResult<Self> {
321        Self::_build(
322            eos,
323            temperature,
324            volume,
325            density,
326            composition,
327            pressure,
328            density_initialization,
329        )?
330        .ok_or_else(|| FeosError::UndeterminedState(String::from("Missing input parameters.")))
331    }
332
333    fn _build<X: Composition<D, N>>(
334        eos: &E,
335        temperature: Temperature<D>,
336        volume: Option<Volume<D>>,
337        density: Option<Density<D>>,
338        composition: X,
339        pressure: Option<Pressure<D>>,
340        density_initialization: Option<DensityInitialization>,
341    ) -> FeosResult<Option<Self>> {
342        // unwrap composition
343        let (x, n) = composition.into_molefracs(eos)?;
344
345        let t = temperature;
346        let di = density_initialization;
347        // find the appropriate state constructor
348        match (volume, density, n, pressure) {
349            (None, None, None, None) => Ok(None),
350            (None, None, Some(_), None) => Ok(None),
351            (Some(_), None, None, None) => Ok(None),
352            (None, None, _, Some(p)) => State::new_npt(eos, t, p, (x, n), di).map(Some),
353            (None, Some(d), _, None) => State::new(eos, t, d, (x, n)).map(Some),
354            (Some(v), None, None, Some(p)) => State::new_tpvx(eos, t, p, v, x, di).map(Some),
355            (Some(v), None, Some(n), None) => State::new_nvt(eos, t, v, (x, n)).map(Some),
356            (Some(v), Some(d), None, None) => State::new_nvt(eos, t, v, (x, d * v)).map(Some),
357            (Some(_), Some(_), Some(_), _) => Err(FeosError::UndeterminedState(String::from(
358                "Density is overdetermined.",
359            ))),
360            (_, _, _, Some(_)) => Err(FeosError::UndeterminedState(String::from(
361                "Pressure is overdetermined.",
362            ))),
363        }
364    }
365
366    /// Return a new `State` using a density iteration. [DensityInitialization] is used to
367    /// influence the calculation with respect to the possible solutions.
368    pub fn new_npt<X: Composition<D, N>>(
369        eos: &E,
370        temperature: Temperature<D>,
371        pressure: Pressure<D>,
372        composition: X,
373        density_initialization: Option<DensityInitialization>,
374    ) -> FeosResult<Self> {
375        let (molefracs, total_moles) = composition.into_molefracs(eos)?;
376        density_iteration(
377            eos,
378            temperature,
379            pressure,
380            &molefracs,
381            density_initialization,
382        )
383        .and_then(|density| Self::_new(eos, temperature, density, molefracs, total_moles))
384    }
385
386    /// Return a new `State` for given pressure $p$, volume $V$, temperature $T$ and composition $x_i$.
387    pub fn new_tpvx(
388        eos: &E,
389        temperature: Temperature<D>,
390        pressure: Pressure<D>,
391        volume: Volume<D>,
392        molefracs: OVector<D, N>,
393        density_initialization: Option<DensityInitialization>,
394    ) -> FeosResult<Self> {
395        let density = density_iteration(
396            eos,
397            temperature,
398            pressure,
399            &molefracs,
400            density_initialization,
401        )?;
402        Self::new_nvt(eos, temperature, volume, (molefracs, density * volume))
403    }
404}
405
406impl<E: Total<N, D>, N: Gradients, D: DualNum<f64> + Copy> State<E, N, D>
407where
408    DefaultAllocator: Allocator<N>,
409{
410    /// Return a new `State` for the combination of inputs.
411    ///
412    /// # Errors
413    ///
414    /// When the state cannot be created using the combination of inputs (is over- or underdetermined).
415    #[expect(clippy::too_many_arguments)]
416    pub fn build_full<X: Composition<D, N> + Clone>(
417        eos: &E,
418        temperature: Option<Temperature<D>>,
419        volume: Option<Volume<D>>,
420        density: Option<Density<D>>,
421        composition: X,
422        pressure: Option<Pressure<D>>,
423        molar_enthalpy: Option<MolarEnergy<D>>,
424        molar_entropy: Option<MolarEntropy<D>>,
425        molar_internal_energy: Option<MolarEnergy<D>>,
426        density_initialization: Option<DensityInitialization>,
427        initial_temperature: Option<Temperature<D>>,
428    ) -> FeosResult<Self> {
429        let state = if let Some(temperature) = temperature {
430            Self::_build(
431                eos,
432                temperature,
433                volume,
434                density,
435                composition.clone(),
436                pressure,
437                density_initialization,
438            )?
439        } else {
440            None
441        };
442
443        let ti = initial_temperature;
444        match state {
445            Some(state) => Ok(state),
446            None => {
447                match (
448                    temperature,
449                    pressure,
450                    volume,
451                    molar_enthalpy,
452                    molar_entropy,
453                    molar_internal_energy,
454                ) {
455                    (Some(t), None, None, Some(h), None, None) => {
456                        State::new_nth(eos, t, h, composition, density_initialization)
457                    }
458                    (Some(t), None, None, None, Some(s), None) => {
459                        State::new_nts(eos, t, s, composition, density_initialization)
460                    }
461                    (None, Some(p), None, Some(h), None, None) => {
462                        State::new_nph(eos, p, h, composition, density_initialization, ti)
463                    }
464                    (None, Some(p), None, None, Some(s), None) => {
465                        State::new_nps(eos, p, s, composition, density_initialization, ti)
466                    }
467                    (None, None, Some(v), None, None, Some(u)) => {
468                        State::new_nvu(eos, v, u, composition, ti)
469                    }
470                    _ => Err(FeosError::UndeterminedState(String::from(
471                        "Missing input parameters.",
472                    ))),
473                }
474            }
475        }
476    }
477
478    /// Return a new `State` for given pressure $p$ and molar enthalpy $h$.
479    pub fn new_nph<X: Composition<D, N> + Clone>(
480        eos: &E,
481        pressure: Pressure<D>,
482        molar_enthalpy: MolarEnergy<D>,
483        composition: X,
484        density_initialization: Option<DensityInitialization>,
485        initial_temperature: Option<Temperature<D>>,
486    ) -> FeosResult<Self> {
487        let t0 = initial_temperature.unwrap_or(Temperature::from_reduced(D::from(298.15)));
488        let mut density = density_initialization;
489        let f = |x0| {
490            let s = State::new_npt(eos, x0, pressure, composition.clone(), density)?;
491            let dfx = s.molar_isobaric_heat_capacity(Contributions::Total);
492            let fx = s.molar_enthalpy(Contributions::Total) - molar_enthalpy;
493            density = Some(DensityInitialization::InitialDensity(s.density.re()));
494            Ok((fx, dfx, s))
495        };
496        newton(t0, f, Temperature::from_reduced(1.0e-8))
497    }
498
499    /// Return a new `State` for given temperature $T$ and molar enthalpy $h$.
500    pub fn new_nth<X: Composition<D, N> + Clone>(
501        eos: &E,
502        temperature: Temperature<D>,
503        molar_enthalpy: MolarEnergy<D>,
504        composition: X,
505        density_initialization: Option<DensityInitialization>,
506    ) -> FeosResult<Self> {
507        let (x, _) = composition.clone().into_molefracs(eos)?;
508        let rho0 = match density_initialization {
509            Some(DensityInitialization::InitialDensity(r)) => {
510                Density::from_reduced(D::from(r.into_reduced()))
511            }
512            Some(DensityInitialization::Liquid) => eos.max_density(&x)?,
513            Some(DensityInitialization::Vapor) => eos.max_density(&x)? * 1.0e-5,
514            None => eos.max_density(&x)? * 0.01,
515        };
516        let f = |rho| {
517            let s = State::new(eos, temperature, rho, composition.clone())?;
518            let dfx = -s.molar_volume
519                * s.molar_volume
520                * (s.molar_volume * s.dp_dv(Contributions::Total)
521                    + temperature * s.dp_dt(Contributions::Total));
522            let fx = s.molar_enthalpy(Contributions::Total) - molar_enthalpy;
523            Ok((fx, dfx, s))
524        };
525        newton(rho0, f, Density::from_reduced(1.0e-12))
526    }
527
528    /// Return a new `State` for given temperature $T$ and molar entropy $s$.
529    pub fn new_nts<X: Composition<D, N> + Clone>(
530        eos: &E,
531        temperature: Temperature<D>,
532        molar_entropy: MolarEntropy<D>,
533        composition: X,
534        density_initialization: Option<DensityInitialization>,
535    ) -> FeosResult<Self> {
536        let (x, _) = composition.clone().into_molefracs(eos)?;
537        let rho0 = match density_initialization {
538            Some(DensityInitialization::InitialDensity(r)) => {
539                Density::from_reduced(D::from(r.into_reduced()))
540            }
541            Some(DensityInitialization::Liquid) => eos.max_density(&x)?,
542            Some(DensityInitialization::Vapor) => eos.max_density(&x)? * 1.0e-5,
543            None => eos.max_density(&x)? * 0.01,
544        };
545        let f = |rho| {
546            let s = State::new(eos, temperature, rho, composition.clone())?;
547            let dfx = -s.molar_volume * s.molar_volume * s.dp_dt(Contributions::Total);
548            let fx = s.molar_entropy(Contributions::Total) - molar_entropy;
549            Ok((fx, dfx, s))
550        };
551        newton(rho0, f, Density::from_reduced(1.0e-12))
552    }
553
554    /// Return a new `State` for given pressure $p$ and molar entropy $s$.
555    pub fn new_nps<X: Composition<D, N> + Clone>(
556        eos: &E,
557        pressure: Pressure<D>,
558        molar_entropy: MolarEntropy<D>,
559        composition: X,
560        density_initialization: Option<DensityInitialization>,
561        initial_temperature: Option<Temperature<D>>,
562    ) -> FeosResult<Self> {
563        let t0 = initial_temperature.unwrap_or(Temperature::from_reduced(D::from(298.15)));
564        let mut density = density_initialization;
565        let f = |x0| {
566            let s = State::new_npt(eos, x0, pressure, composition.clone(), density)?;
567            let dfx = s.molar_isobaric_heat_capacity(Contributions::Total) / s.temperature;
568            let fx = s.molar_entropy(Contributions::Total) - molar_entropy;
569            density = Some(DensityInitialization::InitialDensity(s.density.re()));
570            Ok((fx, dfx, s))
571        };
572        newton(t0, f, Temperature::from_reduced(1.0e-8))
573    }
574
575    /// Return a new `State` for given volume $V$ and molar internal energy $u$.
576    pub fn new_nvu<X: Composition<D, N> + Clone>(
577        eos: &E,
578        volume: Volume<D>,
579        molar_internal_energy: MolarEnergy<D>,
580        composition: X,
581        initial_temperature: Option<Temperature<D>>,
582    ) -> FeosResult<Self> {
583        let t0 = initial_temperature.unwrap_or(Temperature::from_reduced(D::from(298.15)));
584        let f = |x0| {
585            let s = State::new_nvt(eos, x0, volume, composition.clone())?;
586            let fx = s.molar_internal_energy(Contributions::Total) - molar_internal_energy;
587            let dfx = s.molar_isochoric_heat_capacity(Contributions::Total);
588            Ok((fx, dfx, s))
589        };
590        newton(t0, f, Temperature::from_reduced(1.0e-8))
591    }
592}
593
594fn is_close<U: Copy>(
595    x: Quantity<f64, U>,
596    y: Quantity<f64, U>,
597    atol: Quantity<f64, U>,
598    rtol: f64,
599) -> bool {
600    (x - y).abs() <= atol + rtol * y.abs()
601}
602
603fn newton<E: Residual<N, D>, N: Dim, D: DualNum<f64> + Copy, F, X: Copy, Y>(
604    mut x0: Quantity<D, X>,
605    mut f: F,
606    atol: Quantity<f64, X>,
607) -> FeosResult<State<E, N, D>>
608where
609    DefaultAllocator: Allocator<N>,
610    Y: Sub<X> + Sub<<Y as Sub<X>>::Output, Output = X>,
611    F: FnMut(
612        Quantity<D, X>,
613    ) -> FeosResult<(
614        Quantity<D, Y>,
615        Quantity<D, <Y as Sub<X>>::Output>,
616        State<E, N, D>,
617    )>,
618{
619    let rtol = 1e-10;
620    let maxiter = 50;
621
622    for _ in 0..maxiter {
623        let (fx, dfx, mut state) = f(x0)?;
624        let x = x0 - fx / dfx;
625        if is_close(x.re(), x0.re(), atol, rtol) {
626            // Ensure that at least NDERIV iterations are performed (for implicit AD)
627            for _ in 0..D::NDERIV {
628                let (fx, dfx, s) = f(x0)?;
629                x0 -= fx / dfx;
630                state = s;
631            }
632            return Ok(state);
633        }
634        x0 = x;
635    }
636    Err(FeosError::NotConverged("newton".to_owned()))
637}
638
639/// Validate the given temperature, mole numbers and volume.
640///
641/// Properties are valid if
642/// * they are finite
643/// * they have a positive sign
644///
645/// There is no validation of the physical state, e.g.
646/// if resulting densities are below maximum packing fraction.
647fn validate<N: Dim, D: DualNum<f64>>(
648    temperature: Temperature<D>,
649    density: Density<D>,
650    molefracs: &OVector<D, N>,
651) -> FeosResult<()>
652where
653    DefaultAllocator: Allocator<N>,
654{
655    let t = temperature.re().to_reduced();
656    let rho = density.re().to_reduced();
657    if !t.is_finite() || t.is_sign_negative() {
658        return Err(FeosError::InvalidState(
659            String::from("validate"),
660            String::from("temperature"),
661            t,
662        ));
663    }
664    if !rho.is_finite() || rho.is_sign_negative() {
665        return Err(FeosError::InvalidState(
666            String::from("validate"),
667            String::from("density"),
668            rho,
669        ));
670    }
671    for n in molefracs.iter() {
672        if !n.re().is_finite() || n.re().is_sign_negative() {
673            return Err(FeosError::InvalidState(
674                String::from("validate"),
675                String::from("molefracs"),
676                n.re(),
677            ));
678        }
679    }
680    Ok(())
681}
682
683mod critical_point;
684
685#[cfg(test)]
686mod tests {
687    use super::*;
688    use nalgebra::dvector;
689
690    #[test]
691    fn test_validate() {
692        let temperature = 298.15 * KELVIN;
693        let density = 3000.0 * MOL / METER.powi::<3>();
694        let molefracs = dvector![0.03, 0.02, 0.05];
695        assert!(validate(temperature, density, &molefracs).is_ok());
696    }
697
698    #[test]
699    fn test_negative_temperature() {
700        let temperature = -298.15 * KELVIN;
701        let density = 3000.0 * MOL / METER.powi::<3>();
702        let molefracs = dvector![0.03, 0.02, 0.05];
703        assert!(validate(temperature, density, &molefracs).is_err());
704    }
705
706    #[test]
707    fn test_nan_temperature() {
708        let temperature = f64::NAN * KELVIN;
709        let density = 3000.0 * MOL / METER.powi::<3>();
710        let molefracs = dvector![0.03, 0.02, 0.05];
711        assert!(validate(temperature, density, &molefracs).is_err());
712    }
713
714    #[test]
715    fn test_negative_mole_number() {
716        let temperature = 298.15 * KELVIN;
717        let density = 3000.0 * MOL / METER.powi::<3>();
718        let molefracs = dvector![-0.03, 0.02, 0.05];
719        assert!(validate(temperature, density, &molefracs).is_err());
720    }
721
722    #[test]
723    fn test_nan_mole_number() {
724        let temperature = 298.15 * KELVIN;
725        let density = 3000.0 * MOL / METER.powi::<3>();
726        let molefracs = dvector![f64::NAN, 0.02, 0.05];
727        assert!(validate(temperature, density, &molefracs).is_err());
728    }
729
730    #[test]
731    fn test_negative_density() {
732        let temperature = 298.15 * KELVIN;
733        let density = -3000.0 * MOL / METER.powi::<3>();
734        let molefracs = dvector![0.01, 0.02, 0.05];
735        assert!(validate(temperature, density, &molefracs).is_err());
736    }
737}
feos_core/state/mod.rs

feos_core/state/
mod.rs
