/// Only compute the ideal gas contribution

IdealGas,

/// Only compute the difference between the total and the ideal gas contribution

Residual,

/// Compute ideal gas and residual contributions

Total,

/// Calculate a vapor phase by initializing using the ideal gas.

Vapor,

/// Calculate a liquid phase by using the `max_density`.

Liquid,

/// Use the given density as initial value.

InitialDensity(D),

pub fn into_reduced(self) -> DensityInitialization<f64> {

match self {

Self::Vapor => DensityInitialization::Vapor,

Self::Liquid => DensityInitialization::Liquid,

Self::InitialDensity(d) => DensityInitialization::InitialDensity(d.into_reduced()),

}

}

DefaultAllocator: Allocator<N>,

/// temperature in Kelvin

pub temperature: D,

/// mole fractions

pub molefracs: OVector<D, N>,

/// partial number densities in Angstrom^-3

pub partial_density: OVector<D, N>,

DefaultAllocator: Allocator<N>,

/// Create a new `StateHD` for given temperature, volume and composition.

pub fn new(temperature: D, volume: D, moles: &OVector<D, N>) -> Self {

Self::new_density(temperature, &(moles / volume))

}

/// Create a new `StateHD` for given temperature and partial densities

pub fn new_density(temperature: D, partial_density: &OVector<D, N>) -> Self {

let molefracs = partial_density / partial_density.sum();

Self {

temperature,

molefracs,

partial_density: partial_density.clone(),

}

}

// Since the molefracs can not be reproduced from moles if the density is zero,

// this constructor exists specifically for these cases.

pub(crate) fn new_virial(temperature: D, density: D, molefracs: &OVector<D, N>) -> Self {

let partial_density = molefracs * density;

Self {

temperature,

molefracs: molefracs.clone(),

partial_density,

}

}

DefaultAllocator: Allocator<N>,

/// Equation of state

pub eos: E,

/// Temperature $T$

pub temperature: Temperature<D>,

/// Molar volume $v=\frac{V}{N}$

pub molar_volume: MolarVolume<D>,

/// Total number of moles $N=\sum_iN_i$

pub total_moles: Option<Moles<D>>,

/// Total density $\rho=\frac{N}{V}=\sum_i\rho_i$

pub density: Density<D>,

/// Mole fractions $x_i=\frac{N_i}{N}=\frac{\rho_i}{\rho}$

pub molefracs: OVector<D, N>,

/// Cache

cache: Cache<D, N>,

DefaultAllocator: Allocator<N>,

/// Set the total amount of substance to the given value.

///

/// This method does not introduce inconsistencies, because the

/// total moles are the only field that stores information about

/// the size of the state.

pub fn set_total_moles(mut self, total_moles: Moles<D>) -> State<E, N, D> {

self.total_moles = Some(total_moles);

self

}

/// Partial densities $\rho_i=\frac{N_i}{V}$

pub fn partial_density(&self) -> Density<OVector<D, N>> {

Dimensionless::new(&self.molefracs) * self.density

}

/// Mole numbers $N_i$

pub fn moles(&self) -> FeosResult<Moles<OVector<D, N>>> {

Ok(Dimensionless::new(&self.molefracs) * self.total_moles()?)

}

/// Total moles $N=\sum_iN_i$

pub fn total_moles(&self) -> FeosResult<Moles<D>> {

self.total_moles.ok_or(FeosError::IntensiveState)

}

/// Volume $V$

pub fn volume(&self) -> FeosResult<Volume<D>> {

Ok(self.molar_volume * self.total_moles()?)

}

DefaultAllocator: Allocator<N>,

fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {

if self.eos.components() == 1 {

write!(

f,

"T = {:.5}, ρ = {:.5}",

self.temperature.re(),

self.density.re()

)

} else {

write!(

f,

"T = {:.5}, ρ = {:.5}, x = {:.5?}",

self.temperature.re(),

self.density.re(),

self.molefracs.map(|x| x.re()).as_slice()

)

}

}

DefaultAllocator: Allocator<N>,

/// Return a new `State` given a temperature, an array of mole numbers and a volume.

///

/// This function will perform a validation of the given properties, i.e. test for signs

/// and if values are finite. It will **not** validate physics, i.e. if the resulting

/// densities are below the maximum packing fraction.

pub fn new_nvt<X: Composition<D, N>>(

eos: &E,

temperature: Temperature<D>,

volume: Volume<D>,

composition: X,

) -> FeosResult<Self> {

let (molefracs, total_moles) = composition.into_molefracs(eos)?;

let Some(total_moles) = total_moles else {

return Err(FeosError::UndeterminedState(

"Missing total mole number in the specification!".into(),

));

};

let density = total_moles / volume;

Self::new(eos, temperature, density, (molefracs, total_moles))

}

/// Return a new `State` given a temperature and the partial density of all components.

///

/// This function will perform a validation of the given properties, i.e. test for signs

/// and if values are finite. It will **not** validate physics, i.e. if the resulting

/// densities are below the maximum packing fraction.

pub fn new_density(

eos: &E,

temperature: Temperature<D>,

partial_density: Density<OVector<D, N>>,

) -> FeosResult<Self> {

let density = partial_density.sum();

let molefracs = partial_density.convert_into(density);

Self::new(eos, temperature, density, molefracs)

}

/// Return a new `State` for a pure component given a temperature and a density.

///

/// This function will perform a validation of the given properties, i.e. test for signs

/// and if values are finite. It will **not** validate physics, i.e. if the resulting

/// densities are below the maximum packing fraction.

pub fn new_pure(eos: &E, temperature: Temperature<D>, density: Density<D>) -> FeosResult<Self>

where

(): Composition<D, N>,

{

Self::new(eos, temperature, density, ())

}

/// Return a new `State` given a temperature, a density and the composition.

///

/// This function will perform a validation of the given properties, i.e. test for signs

/// and if values are finite. It will **not** validate physics, i.e. if the resulting

/// densities are below the maximum packing fraction.

pub fn new<X: Composition<D, N>>(

eos: &E,

temperature: Temperature<D>,

density: Density<D>,

composition: X,

) -> FeosResult<Self> {

let (molefracs, total_moles) = composition.into_molefracs(eos)?;

Self::_new(eos, temperature, density, molefracs, total_moles)

}

fn _new(

eos: &E,

temperature: Temperature<D>,

density: Density<D>,

molefracs: OVector<D, N>,

total_moles: Option<Moles<D>>,

) -> FeosResult<Self> {

let molar_volume = density.inv();

validate(temperature, density, &molefracs)?;

Ok(State {

eos: eos.clone(),

temperature,

molar_volume,

density,

molefracs,

total_moles,

cache: Cache::new(),

})

}

/// Return a new `State` for the combination of inputs.

///

/// # Errors

///

/// When the state cannot be created using the combination of inputs (is over- or underdetermined).

pub fn build<X: Composition<D, N>>(

eos: &E,

temperature: Temperature<D>,

volume: Option<Volume<D>>,

density: Option<Density<D>>,

composition: X,

pressure: Option<Pressure<D>>,

density_initialization: Option<DensityInitialization>,

) -> FeosResult<Self> {

Self::_build(

eos,

temperature,

volume,

density,

composition,

pressure,

density_initialization,

)?

.ok_or_else(|| FeosError::UndeterminedState(String::from("Missing input parameters.")))

}

fn _build<X: Composition<D, N>>(

eos: &E,

temperature: Temperature<D>,

volume: Option<Volume<D>>,

density: Option<Density<D>>,

composition: X,

pressure: Option<Pressure<D>>,

density_initialization: Option<DensityInitialization>,

) -> FeosResult<Option<Self>> {

// unwrap composition

let (x, n) = composition.into_molefracs(eos)?;

let t = temperature;

let di = density_initialization;

// find the appropriate state constructor

match (volume, density, n, pressure) {

(None, None, None, None) => Ok(None),

(None, None, Some(_), None) => Ok(None),

(Some(_), None, None, None) => Ok(None),

(None, None, _, Some(p)) => State::new_npt(eos, t, p, (x, n), di).map(Some),

(None, Some(d), _, None) => State::new(eos, t, d, (x, n)).map(Some),

(Some(v), None, None, Some(p)) => State::new_tpvx(eos, t, p, v, x, di).map(Some),

(Some(v), None, Some(n), None) => State::new_nvt(eos, t, v, (x, n)).map(Some),

(Some(v), Some(d), None, None) => State::new_nvt(eos, t, v, (x, d * v)).map(Some),

(Some(_), Some(_), Some(_), _) => Err(FeosError::UndeterminedState(String::from(

"Density is overdetermined.",

))),

(_, _, _, Some(_)) => Err(FeosError::UndeterminedState(String::from(

"Pressure is overdetermined.",

))),

}

}

/// Return a new `State` using a density iteration. [DensityInitialization] is used to

/// influence the calculation with respect to the possible solutions.

pub fn new_npt<X: Composition<D, N>>(

eos: &E,

temperature: Temperature<D>,

pressure: Pressure<D>,

composition: X,

density_initialization: Option<DensityInitialization>,

) -> FeosResult<Self> {

let (molefracs, total_moles) = composition.into_molefracs(eos)?;

density_iteration(

eos,

temperature,

pressure,

&molefracs,

density_initialization,

)

.and_then(|density| Self::_new(eos, temperature, density, molefracs, total_moles))

}

/// Return a new `State` for given pressure $p$, volume $V$, temperature $T$ and composition $x_i$.

pub fn new_tpvx(

eos: &E,

temperature: Temperature<D>,

pressure: Pressure<D>,

volume: Volume<D>,

molefracs: OVector<D, N>,

density_initialization: Option<DensityInitialization>,

) -> FeosResult<Self> {

let density = density_iteration(

eos,

temperature,

pressure,

&molefracs,

density_initialization,

)?;

Self::new_nvt(eos, temperature, volume, (molefracs, density * volume))

}

DefaultAllocator: Allocator<N>,

/// Return a new `State` for the combination of inputs.

///

/// # Errors

///

/// When the state cannot be created using the combination of inputs (is over- or underdetermined).

#[expect(clippy::too_many_arguments)]

pub fn build_full<X: Composition<D, N> + Clone>(

eos: &E,

temperature: Option<Temperature<D>>,

volume: Option<Volume<D>>,

density: Option<Density<D>>,

composition: X,

pressure: Option<Pressure<D>>,

molar_enthalpy: Option<MolarEnergy<D>>,

molar_entropy: Option<MolarEntropy<D>>,

molar_internal_energy: Option<MolarEnergy<D>>,

density_initialization: Option<DensityInitialization>,

initial_temperature: Option<Temperature<D>>,

) -> FeosResult<Self> {

let state = if let Some(temperature) = temperature {

Self::_build(

eos,

temperature,

volume,

density,

composition.clone(),

pressure,

density_initialization,

)?

} else {

None

};

let ti = initial_temperature;

match state {

Some(state) => Ok(state),

None => {

match (

temperature,

pressure,

volume,

molar_enthalpy,

molar_entropy,

molar_internal_energy,

) {

(Some(t), None, None, Some(h), None, None) => {

State::new_nth(eos, t, h, composition, density_initialization)

}

(Some(t), None, None, None, Some(s), None) => {

State::new_nts(eos, t, s, composition, density_initialization)

}

(None, Some(p), None, Some(h), None, None) => {

State::new_nph(eos, p, h, composition, density_initialization, ti)

}

(None, Some(p), None, None, Some(s), None) => {

State::new_nps(eos, p, s, composition, density_initialization, ti)

}

(None, None, Some(v), None, None, Some(u)) => {

State::new_nvu(eos, v, u, composition, ti)

}

_ => Err(FeosError::UndeterminedState(String::from(

"Missing input parameters.",

))),

}

}

}

}

/// Return a new `State` for given pressure $p$ and molar enthalpy $h$.

pub fn new_nph<X: Composition<D, N> + Clone>(

eos: &E,

pressure: Pressure<D>,

molar_enthalpy: MolarEnergy<D>,

composition: X,

density_initialization: Option<DensityInitialization>,

initial_temperature: Option<Temperature<D>>,

) -> FeosResult<Self> {

let t0 = initial_temperature.unwrap_or(Temperature::from_reduced(D::from(298.15)));

let mut density = density_initialization;

let f = |x0| {

let s = State::new_npt(eos, x0, pressure, composition.clone(), density)?;

let dfx = s.molar_isobaric_heat_capacity(Contributions::Total);

let fx = s.molar_enthalpy(Contributions::Total) - molar_enthalpy;

density = Some(DensityInitialization::InitialDensity(s.density.re()));

Ok((fx, dfx, s))

};

newton(t0, f, Temperature::from_reduced(1.0e-8))

}

/// Return a new `State` for given temperature $T$ and molar enthalpy $h$.

pub fn new_nth<X: Composition<D, N> + Clone>(

eos: &E,

temperature: Temperature<D>,

molar_enthalpy: MolarEnergy<D>,

composition: X,

density_initialization: Option<DensityInitialization>,

) -> FeosResult<Self> {

let (x, _) = composition.clone().into_molefracs(eos)?;

let rho0 = match density_initialization {

Some(DensityInitialization::InitialDensity(r)) => {

Density::from_reduced(D::from(r.into_reduced()))

}

Some(DensityInitialization::Liquid) => eos.max_density(&x)?,

Some(DensityInitialization::Vapor) => eos.max_density(&x)? * 1.0e-5,

None => eos.max_density(&x)? * 0.01,

};

let f = |rho| {

let s = State::new(eos, temperature, rho, composition.clone())?;

let dfx = -s.molar_volume

* s.molar_volume

* (s.molar_volume * s.dp_dv(Contributions::Total)

+ temperature * s.dp_dt(Contributions::Total));

let fx = s.molar_enthalpy(Contributions::Total) - molar_enthalpy;

Ok((fx, dfx, s))

};

newton(rho0, f, Density::from_reduced(1.0e-12))

}

/// Return a new `State` for given temperature $T$ and molar entropy $s$.

pub fn new_nts<X: Composition<D, N> + Clone>(

eos: &E,

temperature: Temperature<D>,

molar_entropy: MolarEntropy<D>,

composition: X,

density_initialization: Option<DensityInitialization>,

) -> FeosResult<Self> {

let (x, _) = composition.clone().into_molefracs(eos)?;

let rho0 = match density_initialization {

Some(DensityInitialization::InitialDensity(r)) => {

Density::from_reduced(D::from(r.into_reduced()))

}

Some(DensityInitialization::Liquid) => eos.max_density(&x)?,

Some(DensityInitialization::Vapor) => eos.max_density(&x)? * 1.0e-5,

None => eos.max_density(&x)? * 0.01,

};

let f = |rho| {

let s = State::new(eos, temperature, rho, composition.clone())?;

let dfx = -s.molar_volume * s.molar_volume * s.dp_dt(Contributions::Total);

let fx = s.molar_entropy(Contributions::Total) - molar_entropy;

Ok((fx, dfx, s))

};

newton(rho0, f, Density::from_reduced(1.0e-12))

}

/// Return a new `State` for given pressure $p$ and molar entropy $s$.

pub fn new_nps<X: Composition<D, N> + Clone>(

eos: &E,

pressure: Pressure<D>,

molar_entropy: MolarEntropy<D>,

composition: X,

density_initialization: Option<DensityInitialization>,

initial_temperature: Option<Temperature<D>>,

) -> FeosResult<Self> {

let t0 = initial_temperature.unwrap_or(Temperature::from_reduced(D::from(298.15)));

let mut density = density_initialization;

let f = |x0| {

let s = State::new_npt(eos, x0, pressure, composition.clone(), density)?;

let dfx = s.molar_isobaric_heat_capacity(Contributions::Total) / s.temperature;

let fx = s.molar_entropy(Contributions::Total) - molar_entropy;

density = Some(DensityInitialization::InitialDensity(s.density.re()));

Ok((fx, dfx, s))

};

newton(t0, f, Temperature::from_reduced(1.0e-8))

}

/// Return a new `State` for given volume $V$ and molar internal energy $u$.

pub fn new_nvu<X: Composition<D, N> + Clone>(

eos: &E,

volume: Volume<D>,

molar_internal_energy: MolarEnergy<D>,

composition: X,

initial_temperature: Option<Temperature<D>>,

) -> FeosResult<Self> {

let t0 = initial_temperature.unwrap_or(Temperature::from_reduced(D::from(298.15)));

let f = |x0| {

let s = State::new_nvt(eos, x0, volume, composition.clone())?;

let fx = s.molar_internal_energy(Contributions::Total) - molar_internal_energy;

let dfx = s.molar_isochoric_heat_capacity(Contributions::Total);

Ok((fx, dfx, s))

};

newton(t0, f, Temperature::from_reduced(1.0e-8))

}

x: Quantity<f64, U>,

y: Quantity<f64, U>,

atol: Quantity<f64, U>,

rtol: f64,

) -> bool {

(x - y).abs() <= atol + rtol * y.abs()

mut x0: Quantity<D, X>,

mut f: F,

atol: Quantity<f64, X>,

) -> FeosResult<State<E, N, D>>

DefaultAllocator: Allocator<N>,

Y: Sub<X> + Sub<<Y as Sub<X>>::Output, Output = X>,

F: FnMut(

Quantity<D, X>,

) -> FeosResult<(

Quantity<D, Y>,

Quantity<D, <Y as Sub<X>>::Output>,

State<E, N, D>,

)>,

let rtol = 1e-10;

let maxiter = 50;

for _ in 0..maxiter {

let (fx, dfx, mut state) = f(x0)?;

let x = x0 - fx / dfx;

if is_close(x.re(), x0.re(), atol, rtol) {

// Ensure that at least NDERIV iterations are performed (for implicit AD)

for _ in 0..D::NDERIV {

let (fx, dfx, s) = f(x0)?;

x0 -= fx / dfx;

state = s;

}

return Ok(state);

}

x0 = x;

}

Err(FeosError::NotConverged("newton".to_owned()))

temperature: Temperature<D>,

density: Density<D>,

molefracs: &OVector<D, N>,

) -> FeosResult<()>

DefaultAllocator: Allocator<N>,

let t = temperature.re().to_reduced();

let rho = density.re().to_reduced();

if !t.is_finite() || t.is_sign_negative() {

return Err(FeosError::InvalidState(

String::from("validate"),

String::from("temperature"),

t,

));

}

if !rho.is_finite() || rho.is_sign_negative() {

return Err(FeosError::InvalidState(

String::from("validate"),

String::from("density"),

rho,

));

}

for n in molefracs.iter() {

if !n.re().is_finite() || n.re().is_sign_negative() {

return Err(FeosError::InvalidState(

String::from("validate"),

String::from("molefracs"),

n.re(),

));

}

}

Ok(())

use super::*;

use nalgebra::dvector;

#[test]

fn test_validate() {

let temperature = 298.15 * KELVIN;

let density = 3000.0 * MOL / METER.powi::<3>();

let molefracs = dvector![0.03, 0.02, 0.05];

assert!(validate(temperature, density, &molefracs).is_ok());

}

#[test]

fn test_negative_temperature() {

let temperature = -298.15 * KELVIN;

let density = 3000.0 * MOL / METER.powi::<3>();

let molefracs = dvector![0.03, 0.02, 0.05];

assert!(validate(temperature, density, &molefracs).is_err());

}

#[test]

fn test_nan_temperature() {

let temperature = f64::NAN * KELVIN;

let density = 3000.0 * MOL / METER.powi::<3>();

let molefracs = dvector![0.03, 0.02, 0.05];

assert!(validate(temperature, density, &molefracs).is_err());

}

#[test]

fn test_negative_mole_number() {

let temperature = 298.15 * KELVIN;

let density = 3000.0 * MOL / METER.powi::<3>();

let molefracs = dvector![-0.03, 0.02, 0.05];

assert!(validate(temperature, density, &molefracs).is_err());

}

#[test]

fn test_nan_mole_number() {

let temperature = 298.15 * KELVIN;

let density = 3000.0 * MOL / METER.powi::<3>();

let molefracs = dvector![f64::NAN, 0.02, 0.05];

assert!(validate(temperature, density, &molefracs).is_err());

}

#[test]

fn test_negative_density() {

let temperature = 298.15 * KELVIN;

let density = -3000.0 * MOL / METER.powi::<3>();

let molefracs = dvector![0.01, 0.02, 0.05];

assert!(validate(temperature, density, &molefracs).is_err());

}