diff --git a/crates/feos-core/src/ad/mod.rs b/crates/feos-core/src/ad/mod.rs index ce92d05d3..aae9fc198 100644 --- a/crates/feos-core/src/ad/mod.rs +++ b/crates/feos-core/src/ad/mod.rs @@ -4,7 +4,7 @@ use crate::{FeosResult, PhaseEquilibrium, ReferenceSystem, Residual}; use nalgebra::{Const, SVector, U1, U2}; #[cfg(feature = "rayon")] use ndarray::{Array1, Array2, ArrayView2, Zip}; -use num_dual::{Derivative, DualSVec, DualStruct}; +use num_dual::{Derivative, DualNum, DualSVec, DualStruct, first_derivative, partial2}; use quantity::{Density, Pressure, Temperature}; #[cfg(feature = "rayon")] use quantity::{KELVIN, KILO, METER, MOL, PASCAL}; @@ -57,8 +57,8 @@ pub trait PropertiesAD { let v2 = Gradient::from(v2); let x = Self::pure_molefracs(); - let a1 = self.residual_molar_helmholtz_energy(t, v1, &x); - let a2 = self.residual_molar_helmholtz_energy(t, v2, &x); + let a1 = self.residual_helmholtz_energy(t, v1, &x); + let a2 = self.residual_helmholtz_energy(t, v2, &x); (a1, a2) }; @@ -66,6 +66,50 @@ pub trait PropertiesAD { Ok(Pressure::from_reduced(p)) } + fn boiling_temperature( + &self, + pressure: Pressure, + ) -> FeosResult>> + where + Self: Residual>, + { + let eos_f64 = self.re(); + let (temperature, [vapor_density, liquid_density]) = + PhaseEquilibrium::pure_p(&eos_f64, pressure, None, Default::default())?; + + // implicit differentiation is implemented here instead of just calling pure_t with dual + // numbers, because for the first derivative, we can avoid calculating density derivatives. + let t = temperature.into_reduced(); + let v1 = 1.0 / liquid_density.to_reduced(); + let v2 = 1.0 / vapor_density.to_reduced(); + let p = pressure.into_reduced(); + let t = Gradient::from(t); + let t = t + { + let v1 = Gradient::from(v1); + let v2 = Gradient::from(v2); + let p = Gradient::from(p); + let x = Self::pure_molefracs(); + + let residual_entropy = |v| { + let (a, s) = first_derivative( + partial2( + |t, &v, x| self.lift().residual_helmholtz_energy(t, v, x), + &v, + &x, + ), + t, + ); + (a, -s) + }; + let (a1, s1) = residual_entropy(v1); + let (a2, s2) = residual_entropy(v2); + + let ln_rho = (v1 / v2).ln(); + (p * (v2 - v1) + (a2 - a1 + t * ln_rho)) / (s2 - s1 - ln_rho) + }; + Ok(Temperature::from_reduced(t)) + } + fn equilibrium_liquid_density( &self, temperature: Temperature, @@ -111,6 +155,26 @@ pub trait PropertiesAD { ) } + #[cfg(feature = "rayon")] + fn boiling_temperature_parallel( + parameter_names: [String; P], + parameters: ArrayView2, + input: ArrayView2, + ) -> (Array1, Array2, Array1) + where + Self: ParametersAD<1>, + { + parallelize::<_, Self, _, _>( + parameter_names, + parameters, + input, + |eos: &Self::Lifted>, inp| { + eos.boiling_temperature(inp[0] * PASCAL) + .map(|p| p.convert_into(KELVIN)) + }, + ) + } + #[cfg(feature = "rayon")] fn liquid_density_parallel( parameter_names: [String; P], @@ -184,16 +248,16 @@ pub trait PropertiesAD { let y = y.map(Gradient::from); let x = liquid_molefracs.map(Gradient::from); - let a_v = self.residual_molar_helmholtz_energy(t, v_v, &y); + let a_v = self.residual_helmholtz_energy(t, v_v, &y); let (p_l, mu_res_l, dp_l, dmu_l) = self.dmu_dv(t, v_l, &x); let vi_l = dmu_l / dp_l; let v_l = vi_l.dot(&y); let a_l = (mu_res_l - vi_l * p_l).dot(&y); (a_l, a_v, v_l, v_v) }; - let rho_l = vle.liquid().partial_density.to_reduced(); + let rho_l = vle.liquid().partial_density().to_reduced(); let rho_l = [rho_l[0], rho_l[1]]; - let rho_v = vle.vapor().partial_density.to_reduced(); + let rho_v = vle.vapor().partial_density().to_reduced(); let rho_v = [rho_v[0], rho_v[1]]; let p = -(a_v - a_l + t * (y[0] * (rho_v[0] / rho_l[0]).ln() + y[1] * (rho_v[1] / rho_l[1]).ln() - 1.0)) @@ -234,16 +298,16 @@ pub trait PropertiesAD { let x = x.map(Gradient::from); let y = vapor_molefracs.map(Gradient::from); - let a_l = self.residual_molar_helmholtz_energy(t, v_l, &x); + let a_l = self.residual_helmholtz_energy(t, v_l, &x); let (p_v, mu_res_v, dp_v, dmu_v) = self.dmu_dv(t, v_v, &y); let vi_v = dmu_v / dp_v; let v_v = vi_v.dot(&x); let a_v = (mu_res_v - vi_v * p_v).dot(&x); (a_l, a_v, v_l, v_v) }; - let rho_l = vle.liquid().partial_density.to_reduced(); + let rho_l = vle.liquid().partial_density().to_reduced(); let rho_l = [rho_l[0], rho_l[1]]; - let rho_v = vle.vapor().partial_density.to_reduced(); + let rho_v = vle.vapor().partial_density().to_reduced(); let rho_v = [rho_v[0], rho_v[1]]; let p = -(a_l - a_v + t * (x[0] * (rho_l[0] / rho_v[0]).ln() + x[1] * (rho_l[1] / rho_v[1]).ln() - 1.0)) diff --git a/crates/feos-core/src/cubic.rs b/crates/feos-core/src/cubic.rs index 562422b05..a12f615e5 100644 --- a/crates/feos-core/src/cubic.rs +++ b/crates/feos-core/src/cubic.rs @@ -221,7 +221,7 @@ mod tests { let parameters = PengRobinsonParameters::new_pure(propane)?; let pr = PengRobinson::new(parameters); let options = SolverOptions::new().verbosity(Verbosity::Iter); - let cp = State::critical_point(&&pr, None, None, None, options)?; + let cp = State::critical_point(&&pr, (), None, None, options)?; println!("{} {}", cp.temperature, cp.pressure(Contributions::Total)); assert_relative_eq!(cp.temperature, tc * KELVIN, max_relative = 1e-4); assert_relative_eq!( diff --git a/crates/feos-core/src/density_iteration.rs b/crates/feos-core/src/density_iteration.rs index a54cfd8b9..aaee68ec1 100644 --- a/crates/feos-core/src/density_iteration.rs +++ b/crates/feos-core/src/density_iteration.rs @@ -87,7 +87,7 @@ where let (a_res, da_res) = first_derivative( |molar_volume| { eos.lift() - .residual_molar_helmholtz_energy(t, molar_volume, &x) + .residual_helmholtz_energy(t, molar_volume, &x) }, molar_volume, ); diff --git a/crates/feos-core/src/equation_of_state/mod.rs b/crates/feos-core/src/equation_of_state/mod.rs index fa72f6997..dc849068f 100644 --- a/crates/feos-core/src/equation_of_state/mod.rs +++ b/crates/feos-core/src/equation_of_state/mod.rs @@ -1,9 +1,10 @@ -use crate::{ReferenceSystem, StateHD}; +use crate::ReferenceSystem; +use crate::state::StateHD; use nalgebra::{ Const, DVector, DefaultAllocator, Dim, Dyn, OVector, SVector, U1, allocator::Allocator, }; use num_dual::DualNum; -use quantity::{Energy, MolarEnergy, Moles, Temperature, Volume}; +use quantity::{Dimensionless, MolarEnergy, MolarVolume, Temperature}; use std::ops::Deref; mod residual; @@ -164,17 +165,17 @@ where fn ideal_gas_helmholtz_energy + Copy>( &self, temperature: Temperature, - volume: Volume, - moles: &Moles>, - ) -> Energy { + volume: MolarVolume, + moles: &OVector, + ) -> MolarEnergy { let total_moles = moles.sum(); let molefracs = moles / total_moles; - let molar_volume = volume / total_moles; + let molar_volume = volume.into_reduced() / total_moles; MolarEnergy::from_reduced(self.ideal_gas_molar_helmholtz_energy( temperature.into_reduced(), - molar_volume.into_reduced(), + molar_volume, &molefracs, - )) * total_moles + )) * Dimensionless::new(total_moles) } } diff --git a/crates/feos-core/src/equation_of_state/residual.rs b/crates/feos-core/src/equation_of_state/residual.rs index 326aff883..af5fc88f6 100644 --- a/crates/feos-core/src/equation_of_state/residual.rs +++ b/crates/feos-core/src/equation_of_state/residual.rs @@ -1,6 +1,10 @@ -use crate::{FeosError, FeosResult, ReferenceSystem, state::StateHD}; -use nalgebra::{DVector, DefaultAllocator, Dim, Dyn, OMatrix, OVector, U1, allocator::Allocator}; -use num_dual::{DualNum, Gradients, partial, partial2, second_derivative, third_derivative}; +use crate::state::StateHD; +use crate::{Composition, FeosResult, ReferenceSystem}; +use nalgebra::allocator::Allocator; +use nalgebra::{DVector, DefaultAllocator, Dim, Dyn, OMatrix, OVector, SVector, U1, U2}; +use num_dual::{ + DualNum, Gradients, hessian, partial, partial2, second_derivative, third_derivative, +}; use quantity::ad::first_derivative; use quantity::*; use std::ops::{Deref, Div}; @@ -140,74 +144,42 @@ where .fold(D::zero(), |acc, (_, a)| acc + a) } - /// Evaluate the molar Helmholtz energy of each individual contribution - /// and return them together with a string representation of the contribution. - fn molar_helmholtz_energy_contributions( + /// Evaluate the Helmholtz energy of each individual contribution and return them + /// together with a string representation of the contribution. + fn helmholtz_energy_contributions( &self, temperature: D, - molar_volume: D, - molefracs: &OVector, + volume: D, + moles: &OVector, ) -> Vec<(&'static str, D)> { - let state = StateHD::new(temperature, molar_volume, molefracs); + let state = StateHD::new(temperature, volume, moles); self.reduced_helmholtz_energy_density_contributions(&state) .into_iter() - .map(|(n, f)| (n, f * temperature * molar_volume)) + .map(|(n, f)| (n, f * temperature * volume)) .collect() } - /// Evaluate the residual molar Helmholtz energy $a^\mathrm{res}$. - fn residual_molar_helmholtz_energy( - &self, - temperature: D, - molar_volume: D, - molefracs: &OVector, - ) -> D { - let state = StateHD::new(temperature, molar_volume, molefracs); - self.reduced_residual_helmholtz_energy_density(&state) * temperature * molar_volume - } - /// Evaluate the residual Helmholtz energy $A^\mathrm{res}$. fn residual_helmholtz_energy(&self, temperature: D, volume: D, moles: &OVector) -> D { - let state = StateHD::new_density(temperature, &(moles / volume)); + let state = StateHD::new(temperature, volume, moles); self.reduced_residual_helmholtz_energy_density(&state) * temperature * volume } - /// Evaluate the residual Helmholtz energy $A^\mathrm{res}$. - fn residual_helmholtz_energy_unit( + /// Evaluate the residual molar Helmholtz energy $a^\mathrm{res}$. + /// + /// The molefracs are treated as independently variable in order to + /// calculate derivatives like the chemical potential. + fn residual_molar_helmholtz_energy( &self, temperature: Temperature, - volume: Volume, - moles: &Moles>, - ) -> Energy { - let temperature = temperature.into_reduced(); - let total_moles = moles.sum(); - let molar_volume = (volume / total_moles).into_reduced(); - let molefracs = moles / total_moles; - let state = StateHD::new(temperature, molar_volume, &molefracs); - Pressure::from_reduced(self.reduced_residual_helmholtz_energy_density(&state) * temperature) - * volume - } - - /// Check if the provided optional molar concentration is consistent with the - /// equation of state. - /// - /// In general, the number of elements in `molefracs` needs to match the number - /// of components of the equation of state. For a pure component, however, - /// no molefracs need to be provided. - fn validate_molefracs(&self, molefracs: &Option>) -> FeosResult> { - let l = molefracs.as_ref().map_or(1, |m| m.len()); - if self.components() == l { - match molefracs { - Some(m) => Ok(m.clone()), - None => Ok(OVector::from_element_generic( - N::from_usize(1), - U1, - D::one(), - )), - } - } else { - Err(FeosError::IncompatibleComponents(self.components(), l)) - } + molar_volume: MolarVolume, + molefracs: &OVector, + ) -> MolarEnergy { + MolarEnergy::from_reduced(self.residual_helmholtz_energy( + temperature.into_reduced(), + molar_volume.into_reduced(), + molefracs, + )) } /// Calculate the maximum density. @@ -216,18 +188,18 @@ where /// equilibria and other iterations. It is not explicitly meant to /// be a mathematical limit for the density (if those exist in the /// equation of state anyways). - fn max_density(&self, molefracs: &Option>) -> FeosResult> { - let x = self.validate_molefracs(molefracs)?; + fn max_density>(&self, composition: X) -> FeosResult> { + let (x, _) = composition.into_molefracs(self)?; Ok(Density::from_reduced(self.compute_max_density(&x))) } /// Calculate the second virial coefficient $B(T)$ - fn second_virial_coefficient( + fn second_virial_coefficient>( &self, temperature: Temperature, - molefracs: &Option>, - ) -> MolarVolume { - let x = self.validate_molefracs(molefracs).unwrap(); + composition: X, + ) -> FeosResult> { + let (x, _) = composition.into_molefracs(self)?; let (_, _, d2f) = second_derivative( partial2( |rho, &t, x| { @@ -241,16 +213,16 @@ where D::from(0.0), ); - Quantity::from_reduced(d2f * 0.5) + Ok(Quantity::from_reduced(d2f * 0.5)) } /// Calculate the third virial coefficient $C(T)$ - fn third_virial_coefficient( + fn third_virial_coefficient>( &self, temperature: Temperature, - molefracs: &Option>, - ) -> Quot, Density> { - let x = self.validate_molefracs(molefracs).unwrap(); + composition: X, + ) -> FeosResult, Density>> { + let (x, _) = composition.into_molefracs(self)?; let (_, _, _, d3f) = third_derivative( partial2( |rho, &t, x| { @@ -264,36 +236,43 @@ where D::from(0.0), ); - Quantity::from_reduced(d3f / 3.0) + Ok(Quantity::from_reduced(d3f / 3.0)) } /// Calculate the temperature derivative of the second virial coefficient $B'(T)$ - fn second_virial_coefficient_temperature_derivative( + fn second_virial_coefficient_temperature_derivative>( &self, temperature: Temperature, - molefracs: &Option>, - ) -> Quot, Temperature> { + composition: X, + ) -> FeosResult, Temperature>> { + let (molefracs, _) = composition.into_molefracs(self)?; let (_, db_dt) = first_derivative( partial( - |t, x| self.lift().second_virial_coefficient(t, x), - molefracs, + |t, x: &OVector<_, _>| self.lift().second_virial_coefficient(t, x).unwrap(), + &molefracs, ), temperature, ); - db_dt + Ok(db_dt) } /// Calculate the temperature derivative of the third virial coefficient $C'(T)$ - fn third_virial_coefficient_temperature_derivative( + #[expect(clippy::type_complexity)] + fn third_virial_coefficient_temperature_derivative>( &self, temperature: Temperature, - molefracs: &Option>, - ) -> Quot, Density>, Temperature> { + composition: X, + ) -> FeosResult, Density>, Temperature>> { + let (molefracs, _) = composition.into_molefracs(self)?; let (_, dc_dt) = first_derivative( - partial(|t, x| self.lift().third_virial_coefficient(t, x), molefracs), + partial( + // TODO: Fallible partial would be nice here... + |t, x: &OVector<_, _>| self.lift().third_virial_coefficient(t, x).unwrap(), + &molefracs, + ), temperature, ); - dc_dt + Ok(dc_dt) } // The following methods are used in phase equilibrium algorithms @@ -303,22 +282,48 @@ where let molar_volume = density.recip(); let (a, da, d2a) = second_derivative( partial2( - |molar_volume, &t, x| { - self.lift() - .residual_molar_helmholtz_energy(t, molar_volume, x) - }, + |molar_volume, &t, x| self.lift().residual_helmholtz_energy(t, molar_volume, x), &temperature, molefracs, ), molar_volume, ); ( - a * density, + a, -da + temperature * density, molar_volume * molar_volume * d2a + temperature, ) } + /// calculates a_res, p, s_res, dp_drho, dp_dt + fn p_dpdrho_dpdt( + &self, + temperature: D, + density: D, + molefracs: &OVector, + ) -> (D, D, D, D, D) { + let molar_volume = density.recip(); + let (a, da, d2a) = hessian::<_, _, _, U2, _>( + partial( + |vt: SVector<_, 2>, x: &OVector<_, N>| { + let [[v, t]] = vt.data.0; + self.lift().residual_helmholtz_energy(t, v, x) + }, + molefracs, + ), + &SVector::from([molar_volume, temperature]), + ); + let [[da_dv, da_dt]] = da.data.0; + let [[d2a_dv2, d2a_dvdt], _] = d2a.data.0; + ( + a, + -da_dv + temperature * density, + -da_dt, + molar_volume * molar_volume * d2a_dv2 + temperature, + -d2a_dvdt + density, + ) + } + /// calculates p, dp_drho, d2p_drho2 fn p_dpdrho_d2pdrho2( &self, @@ -329,10 +334,7 @@ where let molar_volume = density.recip(); let (_, da, d2a, d3a) = third_derivative( partial2( - |molar_volume, &t, x| { - self.lift() - .residual_molar_helmholtz_energy(t, molar_volume, x) - }, + |molar_volume, &t, x| self.lift().residual_helmholtz_energy(t, molar_volume, x), &temperature, molefracs, ), @@ -423,22 +425,22 @@ where fn viscosity_reference( &self, temperature: Temperature, - volume: Volume, - moles: &Moles>, + molar_volume: MolarVolume, + molefracs: &OVector, ) -> Viscosity; fn viscosity_correlation(&self, s_res: D, x: &OVector) -> D; fn diffusion_reference( &self, temperature: Temperature, - volume: Volume, - moles: &Moles>, + molar_volume: MolarVolume, + molefracs: &OVector, ) -> Diffusivity; fn diffusion_correlation(&self, s_res: D, x: &OVector) -> D; fn thermal_conductivity_reference( &self, temperature: Temperature, - volume: Volume, - moles: &Moles>, + molar_volume: MolarVolume, + molefracs: &OVector, ) -> ThermalConductivity; fn thermal_conductivity_correlation(&self, s_res: D, x: &OVector) -> D; } @@ -451,10 +453,11 @@ where fn viscosity_reference( &self, temperature: Temperature, - volume: Volume, - moles: &Moles>, + molar_volume: MolarVolume, + molefracs: &OVector, ) -> Viscosity { - self.deref().viscosity_reference(temperature, volume, moles) + self.deref() + .viscosity_reference(temperature, molar_volume, molefracs) } fn viscosity_correlation(&self, s_res: D, x: &OVector) -> D { self.deref().viscosity_correlation(s_res, x) @@ -462,10 +465,11 @@ where fn diffusion_reference( &self, temperature: Temperature, - volume: Volume, - moles: &Moles>, + molar_volume: MolarVolume, + molefracs: &OVector, ) -> Diffusivity { - self.deref().diffusion_reference(temperature, volume, moles) + self.deref() + .diffusion_reference(temperature, molar_volume, molefracs) } fn diffusion_correlation(&self, s_res: D, x: &OVector) -> D { self.deref().diffusion_correlation(s_res, x) @@ -473,11 +477,11 @@ where fn thermal_conductivity_reference( &self, temperature: Temperature, - volume: Volume, - moles: &Moles>, + molar_volume: MolarVolume, + molefracs: &OVector, ) -> ThermalConductivity { self.deref() - .thermal_conductivity_reference(temperature, volume, moles) + .thermal_conductivity_reference(temperature, molar_volume, molefracs) } fn thermal_conductivity_correlation(&self, s_res: D, x: &OVector) -> D { self.deref().thermal_conductivity_correlation(s_res, x) diff --git a/crates/feos-core/src/errors.rs b/crates/feos-core/src/errors.rs index 24834480f..f02946bc0 100644 --- a/crates/feos-core/src/errors.rs +++ b/crates/feos-core/src/errors.rs @@ -22,7 +22,7 @@ pub enum FeosError { IncompatibleComponents(usize, usize), #[error("Invalid state in {0}: {1} = {2}.")] InvalidState(String, String, f64), - #[error("Undetermined state: {0}.")] + #[error("Undetermined state: {0}")] UndeterminedState(String), #[error("System is supercritical.")] SuperCritical, diff --git a/crates/feos-core/src/lib.rs b/crates/feos-core/src/lib.rs index b51ae0653..a2552007d 100644 --- a/crates/feos-core/src/lib.rs +++ b/crates/feos-core/src/lib.rs @@ -41,7 +41,7 @@ pub use errors::{FeosError, FeosResult}; #[cfg(feature = "ndarray")] pub use phase_equilibria::{PhaseDiagram, PhaseDiagramHetero}; pub use phase_equilibria::{PhaseEquilibrium, TemperatureOrPressure}; -pub use state::{Contributions, DensityInitialization, State, StateBuilder, StateHD, StateVec}; +pub use state::{Composition, Contributions, DensityInitialization, State, StateHD, StateVec}; /// Level of detail in the iteration output. #[derive(Copy, Clone, PartialOrd, PartialEq, Eq, Default)] @@ -205,7 +205,7 @@ impl< mod tests { use crate::Contributions; use crate::FeosResult; - use crate::StateBuilder; + use crate::State; use crate::cubic::*; use crate::equation_of_state::{EquationOfState, IdealGas}; use crate::parameter::*; @@ -268,20 +268,12 @@ mod tests { let parameters = PengRobinsonParameters::new_pure(propane.clone())?; let residual = PengRobinson::new(parameters); - let sr = StateBuilder::new(&&residual) - .temperature(300.0 * KELVIN) - .pressure(1.0 * BAR) - .total_moles(2.0 * MOL) - .build()?; + let sr = State::new_npt(&&residual, 300.0 * KELVIN, 1.0 * BAR, 2.0 * MOL, None)?; let parameters = PengRobinsonParameters::new_pure(propane.clone())?; let residual = PengRobinson::new(parameters); let eos = EquationOfState::new(vec![NoIdealGas], residual); - let s = StateBuilder::new(&&eos) - .temperature(300.0 * KELVIN) - .pressure(1.0 * BAR) - .total_moles(2.0 * MOL) - .build()?; + let s = State::new_npt(&&eos, 300.0 * KELVIN, 1.0 * BAR, 2.0 * MOL, None)?; // pressure assert_relative_eq!( @@ -348,7 +340,7 @@ mod tests { ); assert_relative_eq!( s.gibbs_energy(Contributions::Residual) - - s.total_moles + - s.total_moles() * RGAS * s.temperature * s.compressibility(Contributions::Total).ln(), @@ -424,13 +416,13 @@ mod tests { max_relative = 1e-15 ); assert_relative_eq!( - s.dp_dni(Contributions::Total), - sr.dp_dni(Contributions::Total), + s.n_dp_dni(Contributions::Total), + sr.n_dp_dni(Contributions::Total), max_relative = 1e-15 ); assert_relative_eq!( - s.dp_dni(Contributions::Residual), - sr.dp_dni(Contributions::Residual), + s.n_dp_dni(Contributions::Residual), + sr.n_dp_dni(Contributions::Residual), max_relative = 1e-15 ); @@ -448,8 +440,8 @@ mod tests { max_relative = 1e-15 ); assert_relative_eq!( - s.dmu_dni(Contributions::Residual), - sr.dmu_dni(Contributions::Residual), + s.n_dmu_dni(Contributions::Residual), + sr.n_dmu_dni(Contributions::Residual), max_relative = 1e-15 ); assert_relative_eq!( @@ -462,7 +454,7 @@ mod tests { assert_relative_eq!(s.ln_phi(), sr.ln_phi(), max_relative = 1e-15); assert_relative_eq!(s.dln_phi_dt(), sr.dln_phi_dt(), max_relative = 1e-15); assert_relative_eq!(s.dln_phi_dp(), sr.dln_phi_dp(), max_relative = 1e-15); - assert_relative_eq!(s.dln_phi_dnj(), sr.dln_phi_dnj(), max_relative = 1e-15); + assert_relative_eq!(s.n_dln_phi_dnj(), sr.n_dln_phi_dnj(), max_relative = 1e-15); assert_relative_eq!( s.thermodynamic_factor(), sr.thermodynamic_factor(), diff --git a/crates/feos-core/src/phase_equilibria/bubble_dew.rs b/crates/feos-core/src/phase_equilibria/bubble_dew.rs index da8a12b60..126b314b8 100644 --- a/crates/feos-core/src/phase_equilibria/bubble_dew.rs +++ b/crates/feos-core/src/phase_equilibria/bubble_dew.rs @@ -11,7 +11,7 @@ use nalgebra::{DMatrix, DVector, DefaultAllocator, Dim, Dyn, OVector, U1}; use ndarray::Array1; use num_dual::linalg::LU; use num_dual::{DualNum, DualStruct, Gradients}; -use quantity::{Density, Dimensionless, Moles, Pressure, Quantity, RGAS, SIUnit, Temperature}; +use quantity::{Density, Dimensionless, Pressure, Quantity, RGAS, SIUnit, Temperature}; const MAX_ITER_INNER: usize = 5; const TOL_INNER: f64 = 1e-9; @@ -190,18 +190,28 @@ where bubble: bool, options: (SolverOptions, SolverOptions), ) -> FeosResult { - let (temperature, pressure, iterate_p) = - temperature_or_pressure.temperature_pressure(tp_init); - Self::bubble_dew_point_tp( - eos, - temperature, - pressure, - vapor_molefracs, - liquid_molefracs, - bubble, - iterate_p, - options, - ) + if eos.components() == 1 { + Self::pure_with( + eos, + temperature_or_pressure, + D::from(if bubble { 0.0 } else { 1.0 }), + None, + options.1, + ) + } else { + let (temperature, pressure, iterate_p) = + temperature_or_pressure.temperature_pressure(tp_init); + Self::bubble_dew_point_tp( + eos, + temperature, + pressure, + vapor_molefracs, + liquid_molefracs, + bubble, + iterate_p, + options, + ) + } } #[expect(clippy::too_many_arguments)] @@ -324,7 +334,7 @@ where ) }; } - let state1 = State::new_intensive( + let state1 = State::new( eos, Temperature::from_reduced(t), Density::from_reduced(molar_volume.recip()), @@ -332,18 +342,18 @@ where )?; let rho2_total = rho2.sum(); let x2 = rho2 / rho2_total; - let state2 = State::new_intensive( + let state2 = State::new( eos, Temperature::from_reduced(t), Density::from_reduced(rho2_total), - &x2, + x2, )?; - Ok(PhaseEquilibrium(if bubble { - [state2, state1] + Ok(if bubble { + PhaseEquilibrium([state2, state1], [D::from(0.0), D::from(1.0)]) } else { - [state1, state2] - })) + PhaseEquilibrium([state1, state2], [D::from(1.0), D::from(0.0)]) + }) } fn newton_step_t( @@ -569,8 +579,8 @@ where } else { let mut t = temperature.into_reduced(); let mut p = pressure.into_reduced(); - let mut molar_volume = state1.density.into_reduced().recip(); - let mut rho2 = state2.partial_density.to_reduced(); + let mut molar_volume = state1.molar_volume.into_reduced(); + let mut rho2 = state2.partial_density().to_reduced(); let err = if iterate_p { Self::newton_step_t( &state1.eos, @@ -594,8 +604,19 @@ where }; *temperature = Temperature::from_reduced(t); *pressure = Pressure::from_reduced(p); - state1.density = Density::from_reduced(molar_volume.recip()); - state2.partial_density = Density::from_reduced(rho2); + state1 = State::new( + &state1.eos, + *temperature, + Density::from_reduced(molar_volume.recip()), + molefracs_spec, + )?; + let density = rho2.sum(); + state2 = State::new( + &state2.eos, + *temperature, + Density::from_reduced(density), + rho2 / density, + )?; Ok(err) }?; @@ -618,7 +639,7 @@ where ); Ok(( state1.density.into_reduced().recip(), - state2.partial_density.to_reduced(), + state2.partial_density().to_reduced(), )) } else { // not converged, return error @@ -743,7 +764,7 @@ where liquid_molefracs: &OVector, ) -> FeosResult<(Pressure, OVector)> { let density = 0.75 * Density::from_reduced(eos.compute_max_density(liquid_molefracs)); - let liquid = State::new_intensive(eos, temperature, density, liquid_molefracs)?; + let liquid = State::new(eos, temperature, density, liquid_molefracs)?; let v_l = liquid.partial_molar_volume(); let p_l = liquid.pressure(Contributions::Total); let mu_l = liquid.residual_chemical_potential(); @@ -767,7 +788,7 @@ where let mut x = vapor_molefracs.clone(); for _ in 0..5 { let density = Density::from_reduced(0.75 * eos.compute_max_density(&x)); - let liquid = State::new_intensive(eos, temperature, density, &x)?; + let liquid = State::new(eos, temperature, density, x)?; let v_l = liquid.partial_molar_volume(); let p_l = liquid.pressure(Contributions::Total); let mu_l = liquid.residual_chemical_potential(); @@ -795,7 +816,7 @@ where temperature: Temperature, molefracs: &OVector, ) -> FeosResult { - let [sp_v, sp_l] = State::spinodal(eos, temperature, Some(molefracs), Default::default())?; + let [sp_v, sp_l] = State::spinodal(eos, temperature, molefracs, Default::default())?; let pv = sp_v.pressure(Contributions::Total); let pl = sp_l.pressure(Contributions::Total); Ok(0.5 * (Pressure::from_reduced(0.0).max(pl) + pv)) @@ -809,7 +830,7 @@ where vapor_molefracs: Option<&OVector>, ) -> FeosResult<[State; 2]> { let liquid_state = - State::new_xpt(eos, temperature, pressure, liquid_molefracs, Some(Liquid))?; + State::new_npt(eos, temperature, pressure, liquid_molefracs, Some(Liquid))?; let xv = match vapor_molefracs { Some(xv) => xv.clone(), None => liquid_state @@ -817,7 +838,7 @@ where .map(f64::exp) .component_mul(liquid_molefracs), }; - let vapor_state = State::new_xpt(eos, temperature, pressure, &xv, Some(Vapor))?; + let vapor_state = State::new_npt(eos, temperature, pressure, xv, Some(Vapor))?; Ok([liquid_state, vapor_state]) } @@ -828,13 +849,7 @@ where vapor_molefracs: &OVector, liquid_molefracs: Option<&OVector>, ) -> FeosResult<[State; 2]> { - let vapor_state = State::new_npt( - eos, - temperature, - pressure, - &Moles::from_reduced(vapor_molefracs.clone()), - Some(Vapor), - )?; + let vapor_state = State::new_npt(eos, temperature, pressure, vapor_molefracs, Some(Vapor))?; let xl = match liquid_molefracs { Some(xl) => xl.clone(), None => { @@ -842,13 +857,13 @@ where .ln_phi() .map(f64::exp) .component_mul(vapor_molefracs); - let liquid_state = State::new_xpt(eos, temperature, pressure, &xl, Some(Liquid))?; + let liquid_state = State::new_npt(eos, temperature, pressure, xl, Some(Liquid))?; (vapor_state.ln_phi() - liquid_state.ln_phi()) .map(f64::exp) .component_mul(vapor_molefracs) } }; - let liquid_state = State::new_xpt(eos, temperature, pressure, &xl, Some(Liquid))?; + let liquid_state = State::new_npt(eos, temperature, pressure, xl, Some(Liquid))?; Ok([vapor_state, liquid_state]) } @@ -857,20 +872,20 @@ where pressure: Pressure, state1: &mut State, state2: &mut State, - moles_state2: Option<&Moles>>, + molefracs_state2: Option<&OVector>, ) -> FeosResult<()> { *state1 = State::new_npt( &state1.eos, temperature, pressure, - &state1.moles, + &state1.molefracs, Some(InitialDensity(state1.density)), )?; *state2 = State::new_npt( &state2.eos, temperature, pressure, - moles_state2.unwrap_or(&state2.moles), + molefracs_state2.unwrap_or(&state2.molefracs), Some(InitialDensity(state2.density)), )?; Ok(()) @@ -899,11 +914,11 @@ where "", "" ); - *state2 = State::new_xpt( + *state2 = State::new_npt( &state2.eos, state2.temperature, state2.pressure(Contributions::Total), - &x2, + x2, Some(InitialDensity(state2.density)), )?; Ok(err_out) diff --git a/crates/feos-core/src/phase_equilibria/mod.rs b/crates/feos-core/src/phase_equilibria/mod.rs index 9683e49a1..bd3ef6439 100644 --- a/crates/feos-core/src/phase_equilibria/mod.rs +++ b/crates/feos-core/src/phase_equilibria/mod.rs @@ -1,11 +1,11 @@ use crate::equation_of_state::Residual; -use crate::errors::{FeosError, FeosResult}; +use crate::errors::FeosResult; use crate::state::{DensityInitialization, State}; use crate::{Contributions, ReferenceSystem}; use nalgebra::allocator::Allocator; use nalgebra::{DefaultAllocator, Dim, Dyn, OVector}; -use num_dual::{DualNum, DualStruct, Gradients}; -use quantity::{Energy, Moles, Pressure, RGAS, Temperature}; +use num_dual::{DualNum, Gradients}; +use quantity::{MolarEnergy, Pressure, RGAS, Temperature}; use std::fmt; use std::fmt::Write; @@ -40,16 +40,11 @@ pub use phase_diagram_pure::PhaseDiagram; #[derive(Debug, Clone)] pub struct PhaseEquilibrium + Copy = f64>( pub [State; P], + pub [D; P], ) where DefaultAllocator: Allocator; -// impl Clone for PhaseEquilibrium { -// fn clone(&self) -> Self { -// Self(self.0.clone()) -// } -// } - impl fmt::Display for PhaseEquilibrium { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { for (i, s) in self.0.iter().enumerate() { @@ -119,53 +114,74 @@ impl PhaseEquilibrium { } } -impl, N: Dim> PhaseEquilibrium +impl, N: Dim, D: DualNum + Copy> PhaseEquilibrium where DefaultAllocator: Allocator, { - pub(super) fn from_states(state1: State, state2: State) -> Self { - let (vapor, liquid) = if state1.density.re() < state2.density.re() { - (state1, state2) - } else { - (state2, state1) - }; - Self([vapor, liquid]) + pub(super) fn single_phase(state: &State) -> Self { + Self::two_phase(state.clone(), state.clone()) } - /// Creates a new PhaseEquilibrium that contains two states at the - /// specified temperature, pressure and molefracs. - /// - /// The constructor can be used in custom phase equilibrium solvers or, - /// e.g., to generate initial guesses for an actual VLE solver. - /// In general, the two states generated are NOT in an equilibrium. - pub fn new_xpt( - eos: &E, - temperature: Temperature, - pressure: Pressure, - vapor_molefracs: &OVector, - liquid_molefracs: &OVector, - ) -> FeosResult { - let liquid = State::new_xpt( - eos, - temperature, - pressure, - liquid_molefracs, - Some(DensityInitialization::Liquid), - )?; - let vapor = State::new_xpt( - eos, - temperature, - pressure, - vapor_molefracs, - Some(DensityInitialization::Vapor), - )?; - Ok(Self([vapor, liquid])) + pub(super) fn two_phase(vapor: State, liquid: State) -> Self { + Self::with_vapor_phase_fraction(vapor, liquid, D::from(1.0)) } - pub(super) fn vapor_phase_fraction(&self) -> f64 { - (self.vapor().total_moles / (self.vapor().total_moles + self.liquid().total_moles)) - .into_value() + pub(super) fn with_vapor_phase_fraction( + vapor: State, + liquid: State, + vapor_phase_fraction: D, + ) -> Self { + Self( + [vapor, liquid], + [vapor_phase_fraction, -vapor_phase_fraction + 1.0], + ) } + + // fn from_states(state1: State, state2: State, beta: f64) -> Self { + // let (vapor, liquid) = if state1.density.re() < state2.density.re() { + // (state1, state2) + // } else { + // (state2, state1) + // }; + // Self([vapor, liquid], [beta, 1.0 - beta]) + // } + + // /// Creates a new PhaseEquilibrium that contains two states at the + // /// specified temperature, pressure and molefracs. + // /// + // /// The constructor can be used in custom phase equilibrium solvers or, + // /// e.g., to generate initial guesses for an actual VLE solver. + // /// In general, the two states generated are NOT in an equilibrium. + // pub fn new_xpt( + // eos: &E, + // temperature: Temperature, + // pressure: Pressure, + // vapor_molefracs: &OVector, + // liquid_molefracs: &OVector, + // ) -> FeosResult { + // let liquid = State::new_xpt( + // eos, + // temperature, + // pressure, + // liquid_molefracs, + // Some(DensityInitialization::Liquid), + // )?; + // let vapor = State::new_xpt( + // eos, + // temperature, + // pressure, + // vapor_molefracs, + // Some(DensityInitialization::Vapor), + // )?; + // Ok(Self([vapor, liquid])) + // } + + // pub(super) fn vapor_phase_fraction(&self) -> Option { + // self.vapor() + // .total_moles + // .zip(self.liquid().total_moles) + // .map(|(v, l)| (v / (l + v)).into_value()) + // } } impl, N: Gradients, const P: usize> PhaseEquilibrium @@ -182,38 +198,42 @@ where &s.eos, temperature, pressure, - &s.moles, + &*s, Some(DensityInitialization::InitialDensity(s.density)), )?; } Ok(self) } - pub(super) fn update_moles( + pub(super) fn update_composition( &mut self, pressure: Pressure, - moles: [&Moles>; P], + molefracs: [&OVector; P], + phase_fractions: [f64; P], ) -> FeosResult<()> { for (i, s) in self.0.iter_mut().enumerate() { *s = State::new_npt( &s.eos, s.temperature, pressure, - moles[i], + molefracs[i], Some(DensityInitialization::InitialDensity(s.density)), )?; } + self.1 = phase_fractions; Ok(()) } - // Total Gibbs energy excluding the constant contribution RT sum_i N_i ln(\Lambda_i^3) - pub(super) fn total_gibbs_energy(&self) -> Energy { - self.0.iter().fold(Energy::from_reduced(0.0), |acc, s| { - let ln_rho_m1 = s.partial_density.to_reduced().map(|r| r.ln() - 1.0); - acc + s.residual_helmholtz_energy() - + s.pressure(Contributions::Total) * s.volume - + RGAS * s.temperature * s.total_moles * s.molefracs.dot(&ln_rho_m1) - }) + // Total molar Gibbs energy excluding the constant contribution RT sum_i x_i ln(\Lambda_i^3) + pub(super) fn total_molar_gibbs_energy(&self) -> MolarEnergy { + self.0 + .iter() + .fold(MolarEnergy::from_reduced(0.0), |acc, s| { + let ln_rho_m1 = s.partial_density().to_reduced().map(|r| r.ln() - 1.0); + acc + s.residual_molar_helmholtz_energy() + + s.pressure(Contributions::Total) * s.molar_volume + + RGAS * s.temperature * s.molefracs.dot(&ln_rho_m1) + }) } } @@ -224,14 +244,6 @@ impl, N: Dim> PhaseEquilibrium where DefaultAllocator: Allocator, { - pub(super) fn check_trivial_solution(self) -> FeosResult { - if Self::is_trivial_solution(self.vapor(), self.liquid()) { - Err(FeosError::TrivialSolution) - } else { - Ok(self) - } - } - /// Check if the two states form a trivial solution pub fn is_trivial_solution(state1: &State, state2: &State) -> bool { let rho1 = state1.molefracs.clone() * state1.density.into_reduced(); diff --git a/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs b/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs index 2ebe07f38..0cba6b3b0 100644 --- a/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs +++ b/crates/feos-core/src/phase_equilibria/phase_diagram_binary.rs @@ -1,7 +1,7 @@ use super::bubble_dew::TemperatureOrPressure; use super::{PhaseDiagram, PhaseEquilibrium}; use crate::errors::{FeosError, FeosResult}; -use crate::state::{Contributions, DensityInitialization::Vapor, State, StateBuilder}; +use crate::state::{Contributions, DensityInitialization::Vapor, State}; use crate::{ReferenceSystem, Residual, SolverOptions, Subset}; use nalgebra::{DVector, dvector, stack, vector}; use ndarray::{Array1, s}; @@ -63,7 +63,7 @@ impl PhaseDiagram { None, SolverOptions::default(), )?; - let cp_vle = PhaseEquilibrium::from_states(cp.clone(), cp.clone()); + let cp_vle = PhaseEquilibrium::single_phase(&cp); ([0.0, cp.molefracs[0]], (vle2, cp_vle), bubble) } [None, Some(vle1)] => { @@ -75,7 +75,7 @@ impl PhaseDiagram { None, SolverOptions::default(), )?; - let cp_vle = PhaseEquilibrium::from_states(cp.clone(), cp.clone()); + let cp_vle = PhaseEquilibrium::single_phase(&cp); ([1.0, cp.molefracs[0]], (vle1, cp_vle), bubble) } [Some(vle2), Some(vle1)] => ([0.0, 1.0], (vle2, vle1), true), @@ -391,21 +391,21 @@ impl PhaseEquilibrium { let p0 = (vle1.vapor().pressure(Contributions::Total) + vle2.vapor().pressure(Contributions::Total)) * 0.5; - let nv0 = (&vle1.vapor().moles + &vle2.vapor().moles) * 0.5; - let mut v = State::new_npt(eos, temperature, p0, &nv0, Some(Vapor))?; + let y0 = (&vle1.vapor().molefracs + &vle2.vapor().molefracs) * 0.5; + let mut v = State::new_npt(eos, temperature, p0, y0, Some(Vapor))?; for _ in 0..options.max_iter.unwrap_or(MAX_ITER_HETERO) { // calculate properties - let dmu_drho_l1 = (l1.dmu_dni(Contributions::Total) * l1.volume).to_reduced(); - let dmu_drho_l2 = (l2.dmu_dni(Contributions::Total) * l2.volume).to_reduced(); - let dmu_drho_v = (v.dmu_dni(Contributions::Total) * v.volume).to_reduced(); - let dp_drho_l1 = (l1.dp_dni(Contributions::Total) * l1.volume) + let dmu_drho_l1 = (l1.n_dmu_dni(Contributions::Total) * l1.molar_volume).to_reduced(); + let dmu_drho_l2 = (l2.n_dmu_dni(Contributions::Total) * l2.molar_volume).to_reduced(); + let dmu_drho_v = (v.n_dmu_dni(Contributions::Total) * v.molar_volume).to_reduced(); + let dp_drho_l1 = (l1.n_dp_dni(Contributions::Total) * l1.molar_volume) .to_reduced() .transpose(); - let dp_drho_l2 = (l2.dp_dni(Contributions::Total) * l2.volume) + let dp_drho_l2 = (l2.n_dp_dni(Contributions::Total) * l2.molar_volume) .to_reduced() .transpose(); - let dp_drho_v = (v.dp_dni(Contributions::Total) * v.volume) + let dp_drho_v = (v.n_dp_dni(Contributions::Total) * v.molar_volume) .to_reduced() .transpose(); let mu_l1_res = l1.residual_chemical_potential().to_reduced(); @@ -418,15 +418,15 @@ impl PhaseEquilibrium { // calculate residual let delta_l1v_mu_ig = (RGAS * v.temperature).to_reduced() * (l1 - .partial_density + .partial_density() .to_reduced() - .component_div(&v.partial_density.to_reduced())) + .component_div(&v.partial_density().to_reduced())) .map(f64::ln); let delta_l2v_mu_ig = (RGAS * v.temperature).to_reduced() * (l2 - .partial_density + .partial_density() .to_reduced() - .component_div(&v.partial_density.to_reduced())) + .component_div(&v.partial_density().to_reduced())) .map(f64::ln); let res = stack![ mu_l1_res - &mu_v_res + delta_l1v_mu_ig; @@ -437,7 +437,7 @@ impl PhaseEquilibrium { // check for convergence if res.norm() < options.tol.unwrap_or(TOL_HETERO) { - return Ok(Self([v, l1, l2])); + return Ok(Self([v, l1, l2], [1.0, 0.0, 0.0])); } // calculate Jacobian @@ -453,11 +453,11 @@ impl PhaseEquilibrium { // apply Newton step let rho_l1 = - &l1.partial_density - &Density::from_reduced(dx.rows_range(0..2).into_owned()); + &l1.partial_density() - &Density::from_reduced(dx.rows_range(0..2).into_owned()); let rho_l2 = - &l2.partial_density - &Density::from_reduced(dx.rows_range(2..4).into_owned()); + &l2.partial_density() - &Density::from_reduced(dx.rows_range(2..4).into_owned()); let rho_v = - &v.partial_density - &Density::from_reduced(dx.rows_range(4..6).into_owned()); + &v.partial_density() - &Density::from_reduced(dx.rows_range(4..6).into_owned()); // check for negative densities for i in 0..2 { @@ -472,18 +472,9 @@ impl PhaseEquilibrium { } // update states - l1 = StateBuilder::new(eos) - .temperature(temperature) - .partial_density(&rho_l1) - .build()?; - l2 = StateBuilder::new(eos) - .temperature(temperature) - .partial_density(&rho_l2) - .build()?; - v = StateBuilder::new(eos) - .temperature(temperature) - .partial_density(&rho_v) - .build()?; + l1 = State::new_density(eos, temperature, rho_l1)?; + l2 = State::new_density(eos, temperature, rho_l2)?; + v = State::new_density(eos, temperature, rho_v)?; } Err(FeosError::NotConverged(String::from( "PhaseEquilibrium::heteroazeotrope_t", @@ -513,24 +504,24 @@ impl PhaseEquilibrium { let mut l1 = vle1.liquid().clone(); let mut l2 = vle2.liquid().clone(); let t0 = (vle1.vapor().temperature + vle2.vapor().temperature) * 0.5; - let nv0 = (&vle1.vapor().moles + &vle2.vapor().moles) * 0.5; - let mut v = State::new_npt(eos, t0, pressure, &nv0, Some(Vapor))?; + let y0 = (&vle1.vapor().molefracs + &vle2.vapor().molefracs) * 0.5; + let mut v = State::new_npt(eos, t0, pressure, y0, Some(Vapor))?; for _ in 0..options.max_iter.unwrap_or(MAX_ITER_HETERO) { // calculate properties - let dmu_drho_l1 = (l1.dmu_dni(Contributions::Total) * l1.volume).to_reduced(); - let dmu_drho_l2 = (l2.dmu_dni(Contributions::Total) * l2.volume).to_reduced(); - let dmu_drho_v = (v.dmu_dni(Contributions::Total) * v.volume).to_reduced(); + let dmu_drho_l1 = (l1.n_dmu_dni(Contributions::Total) * l1.molar_volume).to_reduced(); + let dmu_drho_l2 = (l2.n_dmu_dni(Contributions::Total) * l2.molar_volume).to_reduced(); + let dmu_drho_v = (v.n_dmu_dni(Contributions::Total) * v.molar_volume).to_reduced(); let dmu_res_dt_l1 = (l1.dmu_res_dt()).to_reduced(); let dmu_res_dt_l2 = (l2.dmu_res_dt()).to_reduced(); let dmu_res_dt_v = (v.dmu_res_dt()).to_reduced(); - let dp_drho_l1 = (l1.dp_dni(Contributions::Total) * l1.volume) + let dp_drho_l1 = (l1.n_dp_dni(Contributions::Total) * l1.molar_volume) .to_reduced() .transpose(); - let dp_drho_l2 = (l2.dp_dni(Contributions::Total) * l2.volume) + let dp_drho_l2 = (l2.n_dp_dni(Contributions::Total) * l2.molar_volume) .to_reduced() .transpose(); - let dp_drho_v = (v.dp_dni(Contributions::Total) * v.volume) + let dp_drho_v = (v.n_dp_dni(Contributions::Total) * v.molar_volume) .to_reduced() .transpose(); let dp_dt_l1 = (l1.dp_dt(Contributions::Total)).to_reduced(); @@ -545,14 +536,14 @@ impl PhaseEquilibrium { // calculate residual let delta_l1v_dmu_ig_dt = l1 - .partial_density + .partial_density() .to_reduced() - .component_div(&v.partial_density.to_reduced()) + .component_div(&v.partial_density().to_reduced()) .map(f64::ln); let delta_l2v_dmu_ig_dt = l2 - .partial_density + .partial_density() .to_reduced() - .component_div(&v.partial_density.to_reduced()) + .component_div(&v.partial_density().to_reduced()) .map(f64::ln); let delta_l1v_mu_ig = (RGAS * v.temperature).to_reduced() * &delta_l1v_dmu_ig_dt; let delta_l2v_mu_ig = (RGAS * v.temperature).to_reduced() * &delta_l2v_dmu_ig_dt; @@ -566,7 +557,7 @@ impl PhaseEquilibrium { // check for convergence if res.norm() < options.tol.unwrap_or(TOL_HETERO) { - return Ok(Self([v, l1, l2])); + return Ok(Self([v, l1, l2], [1.0, 0.0, 0.0])); } let jacobian = stack![ @@ -582,10 +573,11 @@ impl PhaseEquilibrium { // apply Newton step let rho_l1 = - l1.partial_density - Density::from_reduced(dx.rows_range(0..2).into_owned()); + l1.partial_density() - Density::from_reduced(dx.rows_range(0..2).into_owned()); let rho_l2 = - l2.partial_density - Density::from_reduced(dx.rows_range(2..4).into_owned()); - let rho_v = v.partial_density - Density::from_reduced(dx.rows_range(4..6).into_owned()); + l2.partial_density() - Density::from_reduced(dx.rows_range(2..4).into_owned()); + let rho_v = + v.partial_density() - Density::from_reduced(dx.rows_range(4..6).into_owned()); let t = v.temperature - Temperature::from_reduced(dx[6]); // check for negative densities and temperatures @@ -602,18 +594,9 @@ impl PhaseEquilibrium { } // update states - l1 = StateBuilder::new(eos) - .temperature(t) - .partial_density(&rho_l1) - .build()?; - l2 = StateBuilder::new(eos) - .temperature(t) - .partial_density(&rho_l2) - .build()?; - v = StateBuilder::new(eos) - .temperature(t) - .partial_density(&rho_v) - .build()?; + l1 = State::new_density(eos, t, rho_l1)?; + l2 = State::new_density(eos, t, rho_l2)?; + v = State::new_density(eos, t, rho_v)?; } Err(FeosError::NotConverged(String::from( "PhaseEquilibrium::heteroazeotrope_p", diff --git a/crates/feos-core/src/phase_equilibria/phase_diagram_pure.rs b/crates/feos-core/src/phase_equilibria/phase_diagram_pure.rs index daeb12132..fc6b12a17 100644 --- a/crates/feos-core/src/phase_equilibria/phase_diagram_pure.rs +++ b/crates/feos-core/src/phase_equilibria/phase_diagram_pure.rs @@ -37,7 +37,7 @@ impl PhaseDiagram { let sc = State::critical_point( eos, - None, + (), critical_temperature, None, SolverOptions::default(), @@ -54,7 +54,7 @@ impl PhaseDiagram { states.push(vle.clone()); } } - states.push(PhaseEquilibrium::from_states(sc.clone(), sc)); + states.push(PhaseEquilibrium::single_phase(&sc)); Ok(PhaseDiagram::new(states)) } @@ -104,7 +104,7 @@ impl PhaseDiagram { { let sc = State::critical_point( eos, - None, + (), critical_temperature, None, SolverOptions::default(), diff --git a/crates/feos-core/src/phase_equilibria/phase_envelope.rs b/crates/feos-core/src/phase_equilibria/phase_envelope.rs index 89610e538..c2e1ce060 100644 --- a/crates/feos-core/src/phase_equilibria/phase_envelope.rs +++ b/crates/feos-core/src/phase_equilibria/phase_envelope.rs @@ -20,7 +20,7 @@ impl PhaseDiagram { let sc = State::critical_point( eos, - Some(molefracs), + molefracs, critical_temperature, None, SolverOptions::default(), @@ -51,7 +51,7 @@ impl PhaseDiagram { states.push(vle.clone()); } } - states.push(PhaseEquilibrium::from_states(sc.clone(), sc)); + states.push(PhaseEquilibrium::single_phase(&sc)); Ok(PhaseDiagram::new(states)) } @@ -69,7 +69,7 @@ impl PhaseDiagram { let sc = State::critical_point( eos, - Some(molefracs), + molefracs, critical_temperature, None, SolverOptions::default(), @@ -116,7 +116,7 @@ impl PhaseDiagram { } } - states.push(PhaseEquilibrium::from_states(sc.clone(), sc)); + states.push(PhaseEquilibrium::single_phase(&sc)); Ok(PhaseDiagram::new(states)) } @@ -134,7 +134,7 @@ impl PhaseDiagram { let sc = State::critical_point( eos, - Some(molefracs), + molefracs, critical_temperature, None, SolverOptions::default(), @@ -145,12 +145,12 @@ impl PhaseDiagram { let temperatures = Temperature::linspace(min_temperature, max_temperature, npoints - 1); for ti in &temperatures { - let spinodal = State::spinodal(eos, ti, Some(molefracs), options).ok(); - if let Some(spinodal) = spinodal { - states.push(PhaseEquilibrium(spinodal)); + let spinodal = State::spinodal(eos, ti, molefracs, options).ok(); + if let Some([sp_v, sp_l]) = spinodal { + states.push(PhaseEquilibrium::two_phase(sp_v, sp_l)); } } - states.push(PhaseEquilibrium::from_states(sc.clone(), sc)); + states.push(PhaseEquilibrium::single_phase(&sc)); Ok(PhaseDiagram::new(states)) } diff --git a/crates/feos-core/src/phase_equilibria/stability_analysis.rs b/crates/feos-core/src/phase_equilibria/stability_analysis.rs index 69988f8c2..950f05026 100644 --- a/crates/feos-core/src/phase_equilibria/stability_analysis.rs +++ b/crates/feos-core/src/phase_equilibria/stability_analysis.rs @@ -91,7 +91,7 @@ where &self.eos, self.temperature, self.pressure(Contributions::Total), - &Moles::from_reduced(x_trial), + Moles::from_reduced(x_trial), Some(phase), ) } @@ -122,7 +122,7 @@ where &trial.eos, trial.temperature, trial.pressure(Contributions::Total), - &Moles::from_reduced(y), + Moles::from_reduced(y), Some(DensityInitialization::InitialDensity(trial.density)), )?; if (i > 4 && error > scaled_tol) || (tpd > tpd_old + 1E-05 && i > 2) { @@ -166,9 +166,9 @@ where let (n, _) = di.shape_generic(); // calculate residual and ideal hesse matrix - let mut hesse = (self.dln_phi_dnj() * Moles::from_reduced(1.0)).into_value(); + let mut hesse = self.n_dln_phi_dnj() / self.total_moles().into_reduced(); let lnphi = self.ln_phi(); - let y = self.moles.to_reduced(); + let y = self.moles().into_reduced(); let ln_y = y.map(|y| if y > f64::EPSILON { y.ln() } else { 0.0 }); let sq_y = y.map(f64::sqrt); let gradient = (&ln_y + &lnphi - di).component_mul(&sq_y); @@ -228,7 +228,7 @@ where &self.eos, self.temperature, self.pressure(Contributions::Total), - &Moles::from_reduced(y), + Moles::from_reduced(y), Some(DensityInitialization::InitialDensity(self.density)), )?; } diff --git a/crates/feos-core/src/phase_equilibria/tp_flash.rs b/crates/feos-core/src/phase_equilibria/tp_flash.rs index cd9d9bea3..91ca49201 100644 --- a/crates/feos-core/src/phase_equilibria/tp_flash.rs +++ b/crates/feos-core/src/phase_equilibria/tp_flash.rs @@ -8,7 +8,7 @@ use nalgebra::{DefaultAllocator, Dim, Matrix3, OVector, SVector, U1, U2, vector} use num_dual::{ Dual, Dual2Vec, DualNum, DualStruct, Gradients, first_derivative, implicit_derivative_sp, }; -use quantity::{Dimensionless, MOL, MolarVolume, Moles, Pressure, Quantity, Temperature}; +use quantity::{MOL, MolarVolume, Moles, Pressure, Quantity, Temperature}; const MAX_ITER_TP: usize = 400; const TOL_TP: f64 = 1e-8; @@ -54,8 +54,7 @@ impl, D: DualNum + Copy> PhaseEquilibrium { ) -> FeosResult { let z = feed.get(0).convert_into(feed.get(0) + feed.get(1)); let total_moles = feed.sum(); - let moles = vector![z.re(), 1.0 - z.re()] * MOL; - let vle_re = State::new_npt(&eos.re(), temperature.re(), pressure.re(), &moles, None)? + let vle_re = State::new_npt(&eos.re(), temperature.re(), pressure.re(), z.re(), None)? .tp_flash(None, options, None)?; // implicit differentiation @@ -82,7 +81,7 @@ impl, D: DualNum + Copy> PhaseEquilibrium { let eos = eos.lift(); let molar_gibbs_energy = |x: Dual2Vec<_, _, _>, v| { let molefracs = vector![x, -x + 1.0]; - let a_res = eos.residual_molar_helmholtz_energy(t, v, &molefracs); + let a_res = eos.residual_helmholtz_energy(t, v, &molefracs); let a_ig = (x * (x / v).ln() - (x - 1.0) * ((-x + 1.0) / v).ln() - 1.0) * t; a_res + a_ig + v * p }; @@ -99,9 +98,12 @@ impl, D: DualNum + Copy> PhaseEquilibrium { let state = |x: D, v, phi| { let volume = MolarVolume::from_reduced(v * phi) * total_moles; let moles = Quantity::new(vector![x, -x + 1.0] * phi * total_moles.convert_into(MOL)); - State::new_nvt(eos, temperature, volume, &moles) + State::new_nvt(eos, temperature, volume, moles) }; - Ok(Self([state(y, v_v, beta)?, state(x, v_l, -beta + 1.0)?])) + Ok(Self( + [state(y, v_v, beta)?, state(x, v_l, -beta + 1.0)?], + [beta, -beta + 1.0], + )) } } @@ -184,12 +186,14 @@ where )?; // check convergence - let beta = new_vle_state.vapor_phase_fraction(); + // unwrap is safe here, because after the first successive substitution step the + // phase amounts in new_vle_state are known. let tpd = [ self.tangent_plane_distance(new_vle_state.vapor()), self.tangent_plane_distance(new_vle_state.liquid()), ]; - let dg = (1.0 - beta) * tpd[1] + beta * tpd[0]; + let b = new_vle_state.1; + let dg = b[0] * tpd[0] + b[1] * tpd[1]; // fix if only tpd[1] is positive if tpd[0] < 0.0 && dg >= 0.0 { @@ -287,7 +291,7 @@ where } // calculate total Gibbs energy before the extrapolation - let gibbs = self.total_gibbs_energy(); + let gibbs = self.total_molar_gibbs_energy(); // extrapolate K values let delta_vec = [ @@ -315,7 +319,7 @@ where // calculate new states let mut trial_vle_state = self.clone(); trial_vle_state.update_states(feed_state, &k)?; - if trial_vle_state.total_gibbs_energy() < gibbs { + if trial_vle_state.total_molar_gibbs_energy() < gibbs { *self = trial_vle_state; } } @@ -377,17 +381,21 @@ where fn update_states(&mut self, feed_state: &State, k: &OVector) -> FeosResult<()> { // calculate vapor phase fraction using Rachford-Rice algorithm - let mut beta = self.vapor_phase_fraction(); - beta = rachford_rice(&feed_state.molefracs, k, Some(beta))?; + let beta = self.1[0]; + let b = rachford_rice(&feed_state.molefracs, k, Some(beta))?; // update VLE - let v = feed_state.moles.clone().component_mul(&Dimensionless::new( - k.map(|k| beta * k / (1.0 - beta + beta * k)), - )); - let l = feed_state.moles.clone().component_mul(&Dimensionless::new( - k.map(|k| (1.0 - beta) / (1.0 - beta + beta * k)), - )); - self.update_moles(feed_state.pressure(Contributions::Total), [&v, &l])?; + let v = feed_state + .molefracs + .component_mul(&k.map(|k| b * k / (1.0 - b + b * k))); + let l = feed_state + .molefracs + .component_mul(&k.map(|k| (1.0 - b) / (1.0 - b + b * k))); + self.update_composition( + feed_state.pressure(Contributions::Total), + [&v, &l], + [b, 1.0 - b], + )?; Ok(()) } @@ -396,9 +404,10 @@ where let state1 = stable_states.pop(); let state2 = stable_states.pop(); if let Some(s1) = state1 { - let init1 = Self::from_states(s1.clone(), feed_state.clone()); + // TODO: sort + let init1 = Self::two_phase(s1.clone(), feed_state.clone()); if let Some(s2) = state2 { - Ok((Self::from_states(s1, s2), Some(init1))) + Ok((Self::two_phase(s1, s2), Some(init1))) } else { Ok((init1, None)) } diff --git a/crates/feos-core/src/phase_equilibria/vle_pure.rs b/crates/feos-core/src/phase_equilibria/vle_pure.rs index 5e58601af..6b2822437 100644 --- a/crates/feos-core/src/phase_equilibria/vle_pure.rs +++ b/crates/feos-core/src/phase_equilibria/vle_pure.rs @@ -5,40 +5,65 @@ use crate::errors::{FeosError, FeosResult}; use crate::state::{Contributions, DensityInitialization, State}; use crate::{ReferenceSystem, SolverOptions, TemperatureOrPressure, Verbosity}; use nalgebra::allocator::Allocator; -use nalgebra::{DVector, DefaultAllocator, Dim, dvector}; -use num_dual::{DualNum, DualStruct, Gradients}; -use quantity::{Density, Pressure, RGAS, Temperature}; +use nalgebra::{DVector, DefaultAllocator, Dim, SVector, U1, U2}; +use num_dual::{DualNum, DualStruct, Gradients, gradient, partial}; +use quantity::{Density, Pressure, Temperature}; const SCALE_T_NEW: f64 = 0.7; const MAX_ITER_PURE: usize = 50; const TOL_PURE: f64 = 1e-12; /// # Pure component phase equilibria -impl PhaseEquilibrium { +impl, N: Gradients, D: DualNum + Copy> PhaseEquilibrium +where + DefaultAllocator: Allocator + Allocator + Allocator, +{ /// Calculate a phase equilibrium for a pure component. - pub fn pure( + pub fn pure>( eos: &E, temperature_or_pressure: TP, initial_state: Option<&Self>, options: SolverOptions, ) -> FeosResult { - if let Some(t) = temperature_or_pressure.temperature() { + let (t, [rho_v, rho_l]) = if let Some(t) = temperature_or_pressure.temperature() { let (_, rho) = Self::pure_t(eos, t, initial_state, options)?; - Ok(Self(rho.map(|r| { - State::new_intensive(eos, t, r, &dvector![1.0]).unwrap() - }))) + (t, rho) } else if let Some(p) = temperature_or_pressure.pressure() { - Self::pure_p(eos, p, initial_state, options) + Self::pure_p(eos, p, initial_state, options)? } else { unreachable!() - } + }; + let x = E::pure_molefracs(); + Ok(Self::two_phase( + State::new(eos, t, rho_v, &x)?, + State::new(eos, t, rho_l, x)?, + )) + } + + /// Calculate a phase equilibrium for a pure component. + pub fn pure_with>( + eos: &E, + temperature_or_pressure: TP, + vapor_phase_fraction: D, + initial_state: Option<&Self>, + options: SolverOptions, + ) -> FeosResult { + let (t, [rho_v, rho_l]) = if let Some(t) = temperature_or_pressure.temperature() { + let (_, rho) = Self::pure_t(eos, t, initial_state, options)?; + (t, rho) + } else if let Some(p) = temperature_or_pressure.pressure() { + Self::pure_p(eos, p, initial_state, options)? + } else { + unreachable!() + }; + let x = E::pure_molefracs(); + Ok(Self::with_vapor_phase_fraction( + State::new(eos, t, rho_v, &x)?, + State::new(eos, t, rho_l, x)?, + vapor_phase_fraction, + )) } -} -impl, N: Gradients, D: DualNum + Copy> PhaseEquilibrium -where - DefaultAllocator: Allocator, -{ /// Calculate a phase equilibrium for a pure component /// and given temperature. pub fn pure_t( @@ -89,9 +114,9 @@ where for _ in 0..D::NDERIV { let v_l = liquid_density.recip(); let v_v = vapor_density.recip(); - let (f_l, p_l, dp_l) = eos.p_dpdrho(t, liquid_density, &x); - let (f_v, p_v, dp_v) = eos.p_dpdrho(t, vapor_density, &x); - pressure = -(f_l * v_l - f_v * v_v + t * (v_v / v_l).ln()) / (v_l - v_v); + let (a_l, p_l, dp_l) = eos.p_dpdrho(t, liquid_density, &x); + let (a_v, p_v, dp_v) = eos.p_dpdrho(t, vapor_density, &x); + pressure = -(a_l - a_v + t * (v_v / v_l).ln()) / (v_l - v_v); liquid_density += (pressure - p_l) / dp_l; vapor_density += (pressure - p_v) / dp_v; } @@ -133,15 +158,14 @@ where for i in 1..=max_iter { // calculate properties - let (f_l_res, p_l, p_rho_l) = eos.p_dpdrho(temperature, liquid_density, &x); - let (f_v_res, p_v, p_rho_v) = eos.p_dpdrho(temperature, vapor_density, &x); + let (a_l_res, p_l, p_rho_l) = eos.p_dpdrho(temperature, liquid_density, &x); + let (a_v_res, p_v, p_rho_v) = eos.p_dpdrho(temperature, vapor_density, &x); // Estimate the new pressure let v_v = vapor_density.recip(); let v_l = liquid_density.recip(); let delta_v = v_v - v_l; - let delta_a = - f_v_res * v_v - f_l_res * v_l + temperature * (vapor_density / liquid_density).ln(); + let delta_a = a_v_res - a_l_res + temperature * (vapor_density / liquid_density).ln(); let mut p_new = -delta_a / delta_v; // If the pressure becomes negative, assume the gas phase is ideal. The @@ -211,7 +235,7 @@ where let x = E::pure_molefracs(); let v = (0.75 * eos.compute_max_density(&x)).recip(); let t = temperature.into_reduced(); - let a_res = eos.residual_molar_helmholtz_energy(t, v, &x); + let a_res = eos.residual_helmholtz_energy(t, v, &x); let p = t / v * (a_res / t - 1.0).exp(); let rho_v = p / t; let rho_l = v.recip(); @@ -238,198 +262,248 @@ where Ok((p, [rho_v, rho_l])) } -impl PhaseEquilibrium { - fn new_pt(eos: &E, temperature: Temperature, pressure: Pressure) -> FeosResult { - let liquid = State::new_xpt( - eos, - temperature, - pressure, - &dvector![1.0], - Some(DensityInitialization::Liquid), - )?; - let vapor = State::new_xpt( - eos, - temperature, - pressure, - &dvector![1.0], - Some(DensityInitialization::Vapor), - )?; - Ok(Self([vapor, liquid])) - } - +impl, N: Gradients, D: DualNum + Copy> PhaseEquilibrium +where + DefaultAllocator: Allocator + Allocator + Allocator, +{ /// Calculate a phase equilibrium for a pure component /// and given pressure. - fn pure_p( + pub fn pure_p( eos: &E, - pressure: Pressure, + pressure: Pressure, initial_state: Option<&Self>, options: SolverOptions, - ) -> FeosResult { - let (max_iter, tol, verbosity) = options.unwrap_or(MAX_ITER_PURE, TOL_PURE); + ) -> FeosResult<(Temperature, [Density; 2])> { + let eos_f64 = eos.re(); + let p = pressure.into_reduced(); // Initialize the phase equilibrium - let mut vle = match initial_state { - Some(init) => init - .clone() - .update_pressure(init.vapor().temperature, pressure)?, - None => PhaseEquilibrium::init_pure_p(eos, pressure)?, + let vle = match initial_state { + Some(init) => ( + init.vapor().temperature.into_reduced().re(), + [ + init.vapor().density.into_reduced().re(), + init.liquid().density.into_reduced().re(), + ], + ), + None => init_pure_p(&eos_f64, pressure.re())?, }; + let (t, [rho_v, rho_l]) = iterate_pure_p(&eos_f64, p.re(), vle, options)?; + // Implicit differentiation + let mut temperature = D::from(t); + let mut vapor_density = D::from(rho_v); + let mut liquid_density = D::from(rho_l); + let x = E::pure_molefracs(); + for _ in 0..D::NDERIV { + let v_l = liquid_density.recip(); + let v_v = vapor_density.recip(); + let (a_l, p_l, s_l, p_rho_l, p_t_l) = + eos.p_dpdrho_dpdt(temperature, liquid_density, &x); + let (a_v, p_v, s_v, p_rho_v, p_t_v) = eos.p_dpdrho_dpdt(temperature, vapor_density, &x); + let ln_rho = (v_l / v_v).ln(); + let delta_t = + (p * (v_v - v_l) + (a_v - a_l + temperature * ln_rho)) / (s_v - s_l - ln_rho); + temperature += delta_t; + liquid_density += (p - p_l - p_t_l * delta_t) / p_rho_l; + vapor_density += (p - p_v - p_t_v * delta_t) / p_rho_v; + } + Ok(( + Temperature::from_reduced(temperature), + [ + Density::from_reduced(vapor_density), + Density::from_reduced(liquid_density), + ], + )) + } +} + +/// Calculate a phase equilibrium for a pure component +/// and given pressure. +fn iterate_pure_p, N: Dim>( + eos: &E, + pressure: f64, + (mut temperature, [mut vapor_density, mut liquid_density]): (f64, [f64; 2]), + options: SolverOptions, +) -> FeosResult<(f64, [f64; 2])> +where + DefaultAllocator: Allocator, +{ + let (max_iter, tol, verbosity) = options.unwrap_or(MAX_ITER_PURE, TOL_PURE); + let x = E::pure_molefracs(); + + log_iter!( + verbosity, + " iter | residual | temperature | liquid density | vapor density " + ); + log_iter!(verbosity, "{:-<89}", ""); + log_iter!( + verbosity, + " {:4} | | {:13.8} | {:12.8} | {:12.8}", + 0, + Temperature::from_reduced(temperature), + Density::from_reduced(liquid_density), + Density::from_reduced(vapor_density) + ); + for i in 1..=max_iter { + // calculate properties + let (a_l_res, p_l, s_l_res, p_rho_l, p_t_l) = + eos.p_dpdrho_dpdt(temperature, liquid_density, &x); + let (a_v_res, p_v, s_v_res, p_rho_v, p_t_v) = + eos.p_dpdrho_dpdt(temperature, vapor_density, &x); + + // calculate the molar volumes + let v_l = liquid_density.recip(); + let v_v = vapor_density.recip(); + + // estimate the temperature steps + let ln_rho = (v_l / v_v).ln(); + let delta_t = (pressure * (v_v - v_l) + (a_v_res - a_l_res + temperature * ln_rho)) + / (s_v_res - s_l_res - ln_rho); + temperature += delta_t; + + // calculate Newton steps for the densities and update state. + let rho_l = liquid_density + (pressure - p_l - p_t_l * delta_t) / p_rho_l; + let rho_v = vapor_density + (pressure - p_v - p_t_v * delta_t) / p_rho_v; + + if rho_l.is_sign_negative() || rho_v.is_sign_negative() || delta_t.abs() > 1.0 { + // if densities are negative or the temperature step is large use density iteration instead + liquid_density = _density_iteration( + eos, + temperature, + pressure, + &x, + DensityInitialization::InitialDensity(liquid_density), + )?; + vapor_density = _density_iteration( + eos, + temperature, + pressure, + &x, + DensityInitialization::InitialDensity(vapor_density), + )?; + } else { + liquid_density = rho_l; + vapor_density = rho_v; + } + + // check for trivial solution + if (vapor_density / liquid_density - 1.0).abs() < TRIVIAL_REL_DEVIATION { + return Err(FeosError::TrivialSolution); + } + + // check for convergence + let res = delta_t.abs(); log_iter!( verbosity, - " iter | residual | temperature | liquid density | vapor density " - ); - log_iter!(verbosity, "{:-<89}", ""); - log_iter!( - verbosity, - " {:4} | | {:13.8} | {:12.8} | {:12.8}", - 0, - vle.vapor().temperature, - vle.liquid().density, - vle.vapor().density + " {:4} | {:14.8e} | {:13.8} | {:12.8} | {:12.8}", + i, + res, + Temperature::from_reduced(temperature), + Density::from_reduced(liquid_density), + Density::from_reduced(vapor_density) ); - for i in 1..=max_iter { - // calculate the pressures and derivatives - let (p_l, p_rho_l) = vle.liquid().p_dpdrho(); - let (p_v, p_rho_v) = vle.vapor().p_dpdrho(); - let p_t_l = vle.liquid().dp_dt(Contributions::Total); - let p_t_v = vle.vapor().dp_dt(Contributions::Total); - - // calculate the residual molar entropies (already cached) - let s_l_res = vle.liquid().residual_molar_entropy(); - let s_v_res = vle.vapor().residual_molar_entropy(); - - // calculate the residual molar Helmholtz energies (already cached) - let a_l_res = vle.liquid().residual_molar_helmholtz_energy(); - let a_v_res = vle.vapor().residual_molar_helmholtz_energy(); - - // calculate the molar volumes - let v_l = 1.0 / vle.liquid().density; - let v_v = 1.0 / vle.vapor().density; - - // estimate the temperature steps - let kt = RGAS * vle.vapor().temperature; - let ln_rho = (v_l / v_v).into_value().ln(); - let delta_t = (pressure * (v_v - v_l) + (a_v_res - a_l_res + kt * ln_rho)) - / (s_v_res - s_l_res - RGAS * ln_rho); - let t_new = vle.vapor().temperature + delta_t; - - // calculate Newton steps for the densities and update state. - let rho_l = vle.liquid().density + (pressure - p_l - p_t_l * delta_t) / p_rho_l; - let rho_v = vle.vapor().density + (pressure - p_v - p_t_v * delta_t) / p_rho_v; - - if rho_l.is_sign_negative() - || rho_v.is_sign_negative() - || delta_t.abs() > Temperature::from_reduced(1.0) - { - // if densities are negative or the temperature step is large use density iteration instead - vle = vle - .update_pressure(t_new, pressure)? - .check_trivial_solution()?; - } else { - // update state - vle = Self([ - State::new_pure(eos, t_new, rho_v)?, - State::new_pure(eos, t_new, rho_l)?, - ]); - } - - // check for convergence - let res = delta_t.abs(); - log_iter!( + if res < temperature * tol { + log_result!( verbosity, - " {:4} | {:14.8e} | {:13.8} | {:12.8} | {:12.8}", - i, - res, - vle.vapor().temperature, - vle.liquid().density, - vle.vapor().density + "PhaseEquilibrium::pure_p: calculation converged in {} step(s)\n", + i ); - if res < vle.vapor().temperature * tol { - log_result!( - verbosity, - "PhaseEquilibrium::pure_p: calculation converged in {} step(s)\n", - i - ); - return Ok(vle); - } + return Ok((temperature, [vapor_density, liquid_density])); } - Err(FeosError::NotConverged("pure_p".to_owned())) } + Err(FeosError::NotConverged("pure_p".to_owned())) +} - /// Initialize a new VLE for a pure substance for a given pressure. - fn init_pure_p(eos: &E, pressure: Pressure) -> FeosResult { - let trial_temperatures = [ - Temperature::from_reduced(300.0), - Temperature::from_reduced(500.0), - Temperature::from_reduced(200.0), - ]; - let x = dvector![1.0]; - let mut vle = None; - let mut t0 = Temperature::from_reduced(1.0); - for t in trial_temperatures.iter() { - t0 = *t; - let _vle = PhaseEquilibrium::new_pt(eos, *t, pressure)?; - if !Self::is_trivial_solution(_vle.vapor(), _vle.liquid()) { - return Ok(_vle); +/// Initialize a new VLE for a pure substance for a given pressure. +fn init_pure_p, N: Gradients>( + eos: &E, + pressure: Pressure, +) -> FeosResult<(f64, [f64; 2])> +where + DefaultAllocator: Allocator + Allocator + Allocator, +{ + let trial_temperatures = [300.0, 500.0, 200.0]; + let p = pressure.into_reduced(); + let x = E::pure_molefracs(); + let mut vle = None; + for t in trial_temperatures { + let liquid_density = _density_iteration(eos, t, p, &x, DensityInitialization::Liquid)?; + let vapor_density = _density_iteration(eos, t, p, &x, DensityInitialization::Vapor)?; + let _vle = (t, [vapor_density, liquid_density]); + if (vapor_density / liquid_density - 1.0).abs() >= TRIVIAL_REL_DEVIATION { + return Ok(_vle); + } + vle = Some(_vle); + } + let Some((t0, [mut rho_v, mut rho_l])) = vle else { + unreachable!() + }; + let [mut t_v, mut t_l] = [t0, t0]; + + let cp = State::critical_point(eos, &x, None, None, SolverOptions::default())?; + let cp_density = cp.density.into_reduced(); + if pressure > cp.pressure(Contributions::Total) { + return Err(FeosError::SuperCritical); + }; + + if rho_v < cp_density { + // reduce temperature of liquid phase... + for _ in 0..8 { + t_l *= SCALE_T_NEW; + rho_l = _density_iteration(eos, t_l, p, &x, DensityInitialization::Liquid)?; + if rho_l > cp_density { + break; } - vle = Some(_vle); } - - let cp = State::critical_point(eos, None, None, None, SolverOptions::default())?; - if pressure > cp.pressure(Contributions::Total) { - return Err(FeosError::SuperCritical); - }; - if let Some(mut e) = vle { - if e.vapor().density < cp.density { - for _ in 0..8 { - t0 *= SCALE_T_NEW; - e.0[1] = - State::new_xpt(eos, t0, pressure, &x, Some(DensityInitialization::Liquid))?; - if e.liquid().density > cp.density { - break; - } - } - } else { - for _ in 0..8 { - t0 /= SCALE_T_NEW; - e.0[0] = - State::new_xpt(eos, t0, pressure, &x, Some(DensityInitialization::Vapor))?; - if e.vapor().density < cp.density { - break; - } - } + } else { + // ...or increase temperature of vapor phase + for _ in 0..8 { + t_v /= SCALE_T_NEW; + rho_v = _density_iteration(eos, t_v, p, &x, DensityInitialization::Vapor)?; + if rho_v < cp_density { + break; } + } + } - for _ in 0..20 { - let h = |s: &State<_>| s.residual_enthalpy() + s.total_moles * RGAS * s.temperature; - t0 = (h(e.vapor()) - h(e.liquid())) - / (e.vapor().residual_entropy() - - e.liquid().residual_entropy() - - RGAS - * e.vapor().total_moles - * ((e.vapor().density / e.liquid().density).into_value().ln())); - let trial_state = - State::new_xpt(eos, t0, pressure, &x, Some(DensityInitialization::Vapor))?; - if trial_state.density < cp.density { - e.0[0] = trial_state; - } - let trial_state = - State::new_xpt(eos, t0, pressure, &x, Some(DensityInitialization::Liquid))?; - if trial_state.density > cp.density { - e.0[1] = trial_state; - } - if e.liquid().temperature == e.vapor().temperature { - return Ok(e); - } - } - Err(FeosError::IterationFailed( - "new_init_p: could not find proper initial state".to_owned(), - )) - } else { - unreachable!() + // determine new temperatures and assign them to either the liquid or the vapor phase until + // both phases have the same temperature + for _ in 0..20 { + let h_s = |t, v| { + let (a_res, da_res) = gradient::<_, _, _, U2, _>( + partial( + |t_v: SVector<_, _>, x| { + let [[t, v]] = t_v.data.0; + eos.lift().residual_helmholtz_energy(t, v, x) + }, + &x, + ), + &SVector::from([t, v]), + ); + let [[da_res_dt, da_res_dv]] = da_res.data.0; + (a_res - t * da_res_dt - v * da_res_dv + t, -da_res_dt) + }; + let (h_l, s_l_res) = h_s(t_l, rho_l.recip()); + let (h_v, s_v_res) = h_s(t_v, rho_v.recip()); + let t = (h_v - h_l) / (s_v_res - s_l_res - (rho_v / rho_l).ln()); + let trial_density = _density_iteration(eos, t, p, &x, DensityInitialization::Vapor)?; + if trial_density < cp_density { + rho_v = trial_density; + t_v = t; + } + let trial_density = _density_iteration(eos, t, p, &x, DensityInitialization::Liquid)?; + if trial_density > cp_density { + rho_l = trial_density; + t_l = t; + } + if t_l == t_v { + return Ok((t_l, [rho_v, rho_l])); } } + Err(FeosError::IterationFailed( + "new_init_p: could not find proper initial state".to_owned(), + )) } impl PhaseEquilibrium { @@ -453,7 +527,7 @@ impl PhaseEquilibrium { .map(|i| { let pure_eos = eos.subset(&[i]); PhaseEquilibrium::pure_p(&pure_eos, pressure, None, SolverOptions::default()) - .map(|vle| vle.vapor().temperature) + .map(|(t, _)| t) .ok() }) .collect() @@ -479,19 +553,19 @@ impl PhaseEquilibrium { let mut molefracs_liquid = molefracs_vapor.clone(); molefracs_vapor[i] = 1.0; molefracs_liquid[i] = 1.0; - let vapor = State::new_intensive( + let vapor = State::new( eos, vle_pure.vapor().temperature, vle_pure.vapor().density, - &molefracs_vapor, + molefracs_vapor, )?; - let liquid = State::new_intensive( + let liquid = State::new( eos, vle_pure.liquid().temperature, vle_pure.liquid().density, - &molefracs_liquid, + molefracs_liquid, )?; - Ok(PhaseEquilibrium::from_states(vapor, liquid)) + Ok(PhaseEquilibrium::two_phase(vapor, liquid)) }) .ok() }) diff --git a/crates/feos-core/src/state/builder.rs b/crates/feos-core/src/state/builder.rs deleted file mode 100644 index c418a7c1f..000000000 --- a/crates/feos-core/src/state/builder.rs +++ /dev/null @@ -1,251 +0,0 @@ -use super::{DensityInitialization, State}; -use crate::Total; -use crate::equation_of_state::Residual; -use crate::errors::FeosResult; -use nalgebra::DVector; -use quantity::*; - -/// A simple tool to construct [State]s with arbitrary input parameters. -/// -/// # Examples -/// ``` -/// # use feos_core::{FeosResult, StateBuilder}; -/// # use feos_core::cubic::{PengRobinson, PengRobinsonParameters}; -/// # use quantity::*; -/// # use nalgebra::dvector; -/// # use approx::assert_relative_eq; -/// # fn main() -> FeosResult<()> { -/// // Create a state for given T,V,N -/// let eos = &PengRobinson::new(PengRobinsonParameters::new_simple(&[369.8], &[41.9 * 1e5], &[0.15], &[15.0])?); -/// let state = StateBuilder::new(&eos) -/// .temperature(300.0 * KELVIN) -/// .volume(12.5 * METER.powi::<3>()) -/// .moles(&(dvector![2.5] * MOL)) -/// .build()?; -/// assert_eq!(state.density, 0.2 * MOL / METER.powi::<3>()); -/// -/// // For a pure component, the composition does not need to be specified. -/// let eos = &PengRobinson::new(PengRobinsonParameters::new_simple(&[369.8], &[41.9 * 1e5], &[0.15], &[15.0])?); -/// let state = StateBuilder::new(&eos) -/// .temperature(300.0 * KELVIN) -/// .volume(12.5 * METER.powi::<3>()) -/// .total_moles(2.5 * MOL) -/// .build()?; -/// assert_eq!(state.density, 0.2 * MOL / METER.powi::<3>()); -/// -/// // The state can be constructed without providing any extensive property. -/// let eos = &PengRobinson::new( -/// PengRobinsonParameters::new_simple( -/// &[369.8, 305.4], -/// &[41.9 * 1e5, 48.2 * 1e5], -/// &[0.15, 0.10], -/// &[15.0, 30.0] -/// )? -/// ); -/// let state = StateBuilder::new(&eos) -/// .temperature(300.0 * KELVIN) -/// .partial_density(&(dvector![0.2, 0.6] * MOL / METER.powi::<3>())) -/// .build()?; -/// assert_relative_eq!(state.molefracs, dvector![0.25, 0.75]); -/// assert_relative_eq!(state.density, 0.8 * MOL / METER.powi::<3>()); -/// # Ok(()) -/// # } -/// ``` -#[derive(Clone)] -pub struct StateBuilder<'a, E, const IG: bool> { - eos: &'a E, - temperature: Option, - volume: Option, - density: Option, - partial_density: Option<&'a Density>>, - total_moles: Option, - moles: Option<&'a Moles>>, - molefracs: Option<&'a DVector>, - pressure: Option, - molar_enthalpy: Option, - molar_entropy: Option, - molar_internal_energy: Option, - density_initialization: Option, - initial_temperature: Option, -} - -impl<'a, E: Residual> StateBuilder<'a, E, false> { - /// Create a new `StateBuilder` for the given equation of state. - pub fn new(eos: &'a E) -> Self { - StateBuilder { - eos, - temperature: None, - volume: None, - density: None, - partial_density: None, - total_moles: None, - moles: None, - molefracs: None, - pressure: None, - molar_enthalpy: None, - molar_entropy: None, - molar_internal_energy: None, - density_initialization: None, - initial_temperature: None, - } - } -} - -impl<'a, E: Residual, const IG: bool> StateBuilder<'a, E, IG> { - /// Provide the temperature for the new state. - pub fn temperature(mut self, temperature: Temperature) -> Self { - self.temperature = Some(temperature); - self - } - - /// Provide the volume for the new state. - pub fn volume(mut self, volume: Volume) -> Self { - self.volume = Some(volume); - self - } - - /// Provide the density for the new state. - pub fn density(mut self, density: Density) -> Self { - self.density = Some(density); - self - } - - /// Provide partial densities for the new state. - pub fn partial_density(mut self, partial_density: &'a Density>) -> Self { - self.partial_density = Some(partial_density); - self - } - - /// Provide the total moles for the new state. - pub fn total_moles(mut self, total_moles: Moles) -> Self { - self.total_moles = Some(total_moles); - self - } - - /// Provide the moles for the new state. - pub fn moles(mut self, moles: &'a Moles>) -> Self { - self.moles = Some(moles); - self - } - - /// Provide the molefracs for the new state. - pub fn molefracs(mut self, molefracs: &'a DVector) -> Self { - self.molefracs = Some(molefracs); - self - } - - /// Provide the pressure for the new state. - pub fn pressure(mut self, pressure: Pressure) -> Self { - self.pressure = Some(pressure); - self - } - - /// Specify a vapor state. - pub fn vapor(mut self) -> Self { - self.density_initialization = Some(DensityInitialization::Vapor); - self - } - - /// Specify a liquid state. - pub fn liquid(mut self) -> Self { - self.density_initialization = Some(DensityInitialization::Liquid); - self - } - - /// Provide an initial density used in density iterations. - pub fn initial_density(mut self, initial_density: Density) -> Self { - self.density_initialization = Some(DensityInitialization::InitialDensity(initial_density)); - self - } -} - -impl<'a, E: Total, const IG: bool> StateBuilder<'a, E, IG> { - /// Provide the molar enthalpy for the new state. - pub fn molar_enthalpy(mut self, molar_enthalpy: MolarEnergy) -> StateBuilder<'a, E, true> { - self.molar_enthalpy = Some(molar_enthalpy); - self.convert() - } - - /// Provide the molar entropy for the new state. - pub fn molar_entropy(mut self, molar_entropy: MolarEntropy) -> StateBuilder<'a, E, true> { - self.molar_entropy = Some(molar_entropy); - self.convert() - } - - /// Provide the molar internal energy for the new state. - pub fn molar_internal_energy( - mut self, - molar_internal_energy: MolarEnergy, - ) -> StateBuilder<'a, E, true> { - self.molar_internal_energy = Some(molar_internal_energy); - self.convert() - } - - /// Provide an initial temperature used in the Newton solver. - pub fn initial_temperature( - mut self, - initial_temperature: Temperature, - ) -> StateBuilder<'a, E, true> { - self.initial_temperature = Some(initial_temperature); - self.convert() - } - - fn convert(self) -> StateBuilder<'a, E, true> { - StateBuilder { - eos: self.eos, - temperature: self.temperature, - volume: self.volume, - density: self.density, - partial_density: self.partial_density, - total_moles: self.total_moles, - moles: self.moles, - molefracs: self.molefracs, - pressure: self.pressure, - molar_enthalpy: self.molar_enthalpy, - molar_entropy: self.molar_entropy, - molar_internal_energy: self.molar_internal_energy, - density_initialization: self.density_initialization, - initial_temperature: self.initial_temperature, - } - } -} - -impl StateBuilder<'_, E, false> { - /// Try to build the state with the given inputs. - pub fn build(self) -> FeosResult> { - State::new( - self.eos, - self.temperature, - self.volume, - self.density, - self.partial_density, - self.total_moles, - self.moles, - self.molefracs, - self.pressure, - self.density_initialization, - ) - } -} - -impl StateBuilder<'_, E, true> { - /// Try to build the state with the given inputs. - pub fn build(self) -> FeosResult> { - State::new_full( - self.eos, - self.temperature, - self.volume, - self.density, - self.partial_density, - self.total_moles, - self.moles, - self.molefracs, - self.pressure, - self.molar_enthalpy, - self.molar_entropy, - self.molar_internal_energy, - self.density_initialization, - self.initial_temperature, - ) - } -} diff --git a/crates/feos-core/src/state/cache.rs b/crates/feos-core/src/state/cache.rs index b9a1326f7..cefc798d8 100644 --- a/crates/feos-core/src/state/cache.rs +++ b/crates/feos-core/src/state/cache.rs @@ -12,17 +12,17 @@ pub struct Cache where DefaultAllocator: Allocator, { - pub a: OnceLock>, - pub da_dt: OnceLock>, + pub a: OnceLock>, + pub da_dt: OnceLock>, pub da_dv: OnceLock>, pub da_dn: OnceLock>>, - pub d2a_dt2: OnceLock>>, - pub d2a_dv2: OnceLock>>, + pub d2a_dt2: OnceLock>>, + pub d2a_dv2: OnceLock>>, pub d2a_dtdv: OnceLock>>, pub d2a_dndt: OnceLock, Diff<_MolarEnergy, _Temperature>>>, - pub d2a_dndv: OnceLock, Diff<_MolarEnergy, _Volume>>>, - pub d3a_dt3: OnceLock, _Temperature>>>, - pub d3a_dv3: OnceLock, _Volume>>>, + pub d2a_dndv: OnceLock, Diff<_Energy, _Volume>>>, + pub d3a_dt3: OnceLock, _Temperature>>>, + pub d3a_dv3: OnceLock, _MolarVolume>>>, } impl Cache diff --git a/crates/feos-core/src/state/composition.rs b/crates/feos-core/src/state/composition.rs new file mode 100644 index 000000000..0ad47ac9e --- /dev/null +++ b/crates/feos-core/src/state/composition.rs @@ -0,0 +1,304 @@ +use super::State; +use crate::equation_of_state::Residual; +use crate::{FeosError, FeosResult}; +use nalgebra::allocator::Allocator; +use nalgebra::{DefaultAllocator, Dim, Dyn, OVector, U1, U2, dvector, vector}; +use num_dual::{DualNum, DualStruct}; +use quantity::{Density, Moles, Quantity, SIUnit}; + +pub trait Composition + Copy, N: Dim> +where + DefaultAllocator: Allocator, +{ + #[expect(clippy::type_complexity)] + fn into_molefracs>( + self, + eos: &E, + ) -> FeosResult<(OVector, Option>)>; + fn density(&self) -> Option> { + None + } +} + +pub trait FullComposition + Copy, N: Dim>: Composition +where + DefaultAllocator: Allocator, +{ + fn into_moles>(self, eos: &E) -> FeosResult<(OVector, Moles)>; +} + +// trivial implementations +impl + Copy, N: Dim> Composition for (OVector, Moles) +where + DefaultAllocator: Allocator, +{ + fn into_molefracs>( + self, + _: &E, + ) -> FeosResult<(OVector, Option>)> { + Ok((self.0, Some(self.1))) + } +} + +impl + Copy, N: Dim> FullComposition for (OVector, Moles) +where + DefaultAllocator: Allocator, +{ + fn into_moles>(self, _: &E) -> FeosResult<(OVector, Moles)> { + Ok((self.0, self.1)) + } +} + +impl + Copy, N: Dim> Composition for (OVector, Option>) +where + DefaultAllocator: Allocator, +{ + fn into_molefracs>( + self, + _: &E, + ) -> FeosResult<(OVector, Option>)> { + Ok((self.0, self.1)) + } +} + +// copy the composition from a given state +impl + Copy, N: Dim> Composition for &State +where + DefaultAllocator: Allocator, +{ + fn into_molefracs>( + self, + _: &E1, + ) -> FeosResult<(OVector, Option>)> { + Ok(((self.molefracs.clone()), self.total_moles)) + } +} + +// a pure component needs no specification +impl + Copy> Composition for () { + fn into_molefracs>( + self, + _: &E, + ) -> FeosResult<(OVector, Option>)> { + Ok(((vector![D::one()]), None)) + } +} +impl + Copy> Composition for () { + fn into_molefracs>( + self, + eos: &E, + ) -> FeosResult<(OVector, Option>)> { + if eos.components() == 1 { + Ok(((dvector![D::one()]), None)) + } else { + Err(FeosError::UndeterminedState( + "The composition needs to be specified for a system with more than one component." + .into(), + )) + } + } +} + +// a binary mixture can be specified by a scalar (x1) +impl + Copy> Composition for D { + fn into_molefracs>( + self, + _: &E, + ) -> FeosResult<(OVector, Option>)> { + Ok(((vector![self, -self + 1.0]), None)) + } +} + +// this cannot be implemented generically for D due to mising specialization +impl Composition for f64 { + fn into_molefracs( + self, + eos: &E, + ) -> FeosResult<(OVector, Option)> { + if eos.components() == 2 { + Ok(((dvector![self, 1.0 - self]), None)) + } else { + Err(FeosError::UndeterminedState(format!( + "A scalar ({}) can only be used to specify a binary mixture!", + self + ))) + } + } +} + +// a pure component can be specified by the total mole number +impl + Copy> Composition for Moles { + fn into_molefracs>( + self, + _: &E, + ) -> FeosResult<(OVector, Option>)> { + Ok(((vector![D::one()]), Some(self))) + } +} + +impl + Copy> FullComposition for Moles { + fn into_moles>(self, _: &E) -> FeosResult<(OVector, Moles)> { + Ok(((vector![D::one()]), self)) + } +} + +impl + Copy> Composition for Moles { + fn into_molefracs>( + self, + eos: &E, + ) -> FeosResult<(OVector, Option>)> { + if eos.components() == 1 { + Ok(((dvector![D::one()]), Some(self))) + } else { + Err(FeosError::UndeterminedState(format!( + "A single mole number ({}) can only be used to specify a pure component!", + self.re() + ))) + } + } +} + +impl + Copy> FullComposition for Moles { + fn into_moles>(self, eos: &E) -> FeosResult<(OVector, Moles)> { + if eos.components() == 1 { + Ok(((dvector![D::one()]), self)) + } else { + Err(FeosError::UndeterminedState(format!( + "A single mole number ({}) can only be used to specify a pure component!", + self.re() + ))) + } + } +} + +// the mixture can be specified by its molefractions +// +// for a dynamic number of components, it is also possible to specify only the +// N-1 first components +impl + Copy, N: Dim> Composition for OVector +where + DefaultAllocator: Allocator, +{ + fn into_molefracs>( + self, + eos: &E, + ) -> FeosResult<(OVector, Option>)> { + (&self).into_molefracs(eos) + } +} + +impl + Copy, N: Dim> Composition for &OVector +where + DefaultAllocator: Allocator, +{ + fn into_molefracs>( + self, + eos: &E, + ) -> FeosResult<(OVector, Option>)> { + let sum = self.sum(); + if eos.components() == self.len() { + Ok(((self.clone() / sum), None)) + } else if eos.components() == self.len() + 1 { + let mut x = OVector::zeros_generic(N::from_usize(eos.components()), U1); + for i in 0..self.len() { + x[i] = self[i]; + } + x[self.len()] = -sum + 1.0; + Ok(((x), None)) + } else { + Err(FeosError::UndeterminedState(format!( + "The length of the composition vector ({}) does not match the number of components ({})!", + self.len(), + eos.components() + ))) + } + } +} + +// the mixture can be specified by its moles +impl + Copy, N: Dim> Composition for Moles> +where + DefaultAllocator: Allocator, +{ + fn into_molefracs>( + self, + eos: &E, + ) -> FeosResult<(OVector, Option>)> { + (&self).into_molefracs(eos) + } +} + +impl + Copy, N: Dim> FullComposition for Moles> +where + DefaultAllocator: Allocator, +{ + fn into_moles>(self, eos: &E) -> FeosResult<(OVector, Moles)> { + (&self).into_moles(eos) + } +} + +impl + Copy, N: Dim> Composition for &Moles> +where + DefaultAllocator: Allocator, +{ + fn into_molefracs>( + self, + eos: &E, + ) -> FeosResult<(OVector, Option>)> { + if eos.components() == self.len() { + let total_moles = self.sum(); + Ok(((self.convert_to(total_moles)), Some(total_moles))) + } else { + Err(FeosError::UndeterminedState(format!( + "The length of the composition vector ({}) does not match the number of components ({})!", + self.len(), + eos.components() + ))) + } + } +} + +impl + Copy, N: Dim> FullComposition for &Moles> +where + DefaultAllocator: Allocator, +{ + fn into_moles>(self, eos: &E) -> FeosResult<(OVector, Moles)> { + if eos.components() == self.len() { + let total_moles = self.sum(); + Ok(((self.convert_to(total_moles)), total_moles)) + } else { + Err(FeosError::UndeterminedState(format!( + "The length of the composition vector ({}) does not match the number of components ({})!", + self.len(), + eos.components() + ))) + } + } +} + +// the mixture can be specified with the partial density +impl + Copy, N: Dim> Composition + for Quantity, SIUnit<0, -3, 0, 0, 0, 1, 0>> +where + DefaultAllocator: Allocator, +{ + fn into_molefracs>( + self, + eos: &E, + ) -> FeosResult<(OVector, Option>)> { + if eos.components() == self.len() { + let density = self.sum(); + Ok(((self.convert_to(density)), None)) + } else { + panic!( + "The length of the composition vector ({}) does not match the number of components ({})!", + self.len(), + eos.components() + ) + } + } + + fn density(&self) -> Option> { + Some(self.sum()) + } +} diff --git a/crates/feos-core/src/state/critical_point.rs b/crates/feos-core/src/state/critical_point.rs index e77b95299..4de9f9562 100644 --- a/crates/feos-core/src/state/critical_point.rs +++ b/crates/feos-core/src/state/critical_point.rs @@ -1,4 +1,4 @@ -use super::{DensityInitialization, State}; +use super::{Composition, DensityInitialization, State}; use crate::equation_of_state::Residual; use crate::errors::{FeosError, FeosResult}; use crate::{ReferenceSystem, SolverOptions, Subset, TemperatureOrPressure, Verbosity}; @@ -31,14 +31,14 @@ impl State { let pure_eos = eos.subset(&[i]); let cp = State::critical_point( &pure_eos, - None, + (), initial_temperature, initial_density, options, )?; let mut molefracs = DVector::zeros(eos.components()); molefracs[i] = 1.0; - State::new_intensive(eos, cp.temperature, cp.density, &molefracs) + State::new(eos, cp.temperature, cp.density, molefracs) }) .collect() } @@ -81,7 +81,7 @@ where ); let density = rho[0] + rho[1]; let molefracs = OVector::from_fn_generic(n, U1, |i, _| rho[i] / density); - Self::new_intensive(eos, t, Density::from_reduced(density), &molefracs) + Self::new(eos, t, Density::from_reduced(density), molefracs) } else if let Some(p) = temperature_or_pressure.pressure() { let x = critical_point_binary_p( &eos_re, @@ -103,22 +103,22 @@ where let density = trho[1] + trho[2]; let molefracs = OVector::from_fn_generic(n, U1, |i, _| trho[i + 1] / density); let t = Temperature::from_reduced(trho[0]); - Self::new_intensive(eos, t, Density::from_reduced(density), &molefracs) + Self::new(eos, t, Density::from_reduced(density), molefracs) } else { unreachable!() } } /// Calculate the critical point of a system for given moles. - pub fn critical_point( + pub fn critical_point>( eos: &E, - molefracs: Option<&OVector>, + composition: X, initial_temperature: Option, initial_density: Option, options: SolverOptions, ) -> FeosResult { let eos_re = eos.re(); - let molefracs = molefracs.map_or_else(E::pure_molefracs, |x| x.clone()); + let (molefracs, total_moles) = composition.into_molefracs(eos)?; let x = &molefracs.map(|x| x.re()); let rho_init = initial_density.map(|r| r.into_reduced()); let trial_temperatures = [300.0, 700.0, 500.0]; @@ -144,11 +144,11 @@ where t_rho[1], &molefracs, ); - Self::new_intensive( + Self::new( eos, Temperature::from_reduced(temperature), Density::from_reduced(density), - &molefracs, + (molefracs, total_moles), ) } } @@ -411,18 +411,18 @@ impl, N: Gradients> State where DefaultAllocator: Allocator + Allocator + Allocator, { - pub fn spinodal( + pub fn spinodal + Clone>( eos: &E, temperature: Temperature, - molefracs: Option<&OVector>, + composition: X, options: SolverOptions, ) -> FeosResult<[Self; 2]> { - let critical_point = Self::critical_point(eos, molefracs, None, None, options)?; - let molefracs = molefracs.map_or_else(E::pure_molefracs, |x| x.clone()); + let critical_point = Self::critical_point(eos, composition, None, None, options)?; + let molefracs = &critical_point.molefracs; let spinodal_vapor = Self::calculate_spinodal( eos, temperature, - &molefracs, + molefracs, DensityInitialization::Vapor, options, )?; @@ -430,7 +430,7 @@ where let spinodal_liquid = Self::calculate_spinodal( eos, temperature, - &molefracs, + molefracs, DensityInitialization::InitialDensity(rho), options, )?; @@ -500,12 +500,7 @@ where "Spinodal calculation converged in {} step(s)\n", i ); - return Self::new_intensive( - eos, - temperature, - Density::from_reduced(rho), - molefracs, - ); + return Self::new(eos, temperature, Density::from_reduced(rho), molefracs); } } Err(FeosError::SuperCritical) @@ -571,7 +566,7 @@ where // calculate pressure let a = partial2( - |v, &t, x| eos.lift().residual_molar_helmholtz_energy(t, v, x), + |v, &t, x| eos.lift().residual_helmholtz_energy(t, v, x), &temperature, &molefracs, ); diff --git a/crates/feos-core/src/state/mod.rs b/crates/feos-core/src/state/mod.rs index 5da7a43c5..ec99d4f19 100644 --- a/crates/feos-core/src/state/mod.rs +++ b/crates/feos-core/src/state/mod.rs @@ -11,19 +11,19 @@ use crate::equation_of_state::Residual; use crate::errors::{FeosError, FeosResult}; use crate::{ReferenceSystem, Total}; use nalgebra::allocator::Allocator; -use nalgebra::{DefaultAllocator, Dim, Dyn, OVector, U1}; +use nalgebra::{DefaultAllocator, Dim, Dyn, OVector}; use num_dual::*; use quantity::*; use std::fmt; use std::ops::Sub; -mod builder; mod cache; +mod composition; mod properties; mod residual_properties; mod statevec; -pub use builder::StateBuilder; pub(crate) use cache::Cache; +pub use composition::{Composition, FullComposition}; pub use statevec::StateVec; /// Possible contributions that can be computed. @@ -70,8 +70,6 @@ where { /// temperature in Kelvin pub temperature: D, - // /// volume in Angstrom^3 - // pub molar_volume: D, /// mole fractions pub molefracs: OVector, /// partial number densities in Angstrom^-3 @@ -82,15 +80,9 @@ impl + Copy> StateHD where DefaultAllocator: Allocator, { - /// Create a new `StateHD` for given temperature, molar volume and composition. - pub fn new(temperature: D, molar_volume: D, molefracs: &OVector) -> Self { - let partial_density = molefracs / molar_volume; - - Self { - temperature, - molefracs: molefracs.clone(), - partial_density, - } + /// Create a new `StateHD` for given temperature, volume and composition. + pub fn new(temperature: D, volume: D, moles: &OVector) -> Self { + Self::new_density(temperature, &(moles / volume)) } /// Create a new `StateHD` for given temperature and partial densities @@ -144,7 +136,7 @@ where /// + [State constructors](#state-constructors) /// + [Stability analysis](#stability-analysis) /// + [Flash calculations](#flash-calculations) -#[derive(Debug)] +#[derive(Debug, Clone)] pub struct State + Copy = f64> where DefaultAllocator: Allocator, @@ -153,14 +145,10 @@ where pub eos: E, /// Temperature $T$ pub temperature: Temperature, - /// Volume $V$ - pub volume: Volume, - /// Mole numbers $N_i$ - pub moles: Moles>, + /// Molar volume $v=\frac{V}{N}$ + pub molar_volume: MolarVolume, /// Total number of moles $N=\sum_iN_i$ - pub total_moles: Moles, - /// Partial densities $\rho_i=\frac{N_i}{V}$ - pub partial_density: Density>, + pub total_moles: Option>, /// Total density $\rho=\frac{N}{V}=\sum_i\rho_i$ pub density: Density, /// Mole fractions $x_i=\frac{N_i}{N}=\frac{\rho_i}{\rho}$ @@ -169,22 +157,38 @@ where cache: Cache, } -impl + Copy> Clone for State +impl + Copy> State where DefaultAllocator: Allocator, { - fn clone(&self) -> Self { - Self { - eos: self.eos.clone(), - temperature: self.temperature, - volume: self.volume, - moles: self.moles.clone(), - total_moles: self.total_moles, - partial_density: self.partial_density.clone(), - density: self.density, - molefracs: self.molefracs.clone(), - cache: self.cache.clone(), - } + /// 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) -> State { + self.total_moles = Some(total_moles); + self + } + + /// Partial densities $\rho_i=\frac{N_i}{V}$ + pub fn partial_density(&self) -> Density> { + Dimensionless::new(&self.molefracs) * self.density + } + + /// Mole numbers $N_i$ + pub fn moles(&self) -> Moles> { + Dimensionless::new(&self.molefracs) * self.total_moles() + } + + /// Total moles $N=\sum_iN_i$ + pub fn total_moles(&self) -> Moles { + self.total_moles.expect("Extensive properties can only be evaluated for states that are initialized with extensive properties!") + } + + /// Volume $V$ + pub fn volume(&self) -> Volume { + self.molar_volume * self.total_moles() } } @@ -221,83 +225,76 @@ where /// 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( + pub fn new_nvt>( eos: &E, temperature: Temperature, volume: Volume, - moles: &Moles>, + composition: X, ) -> FeosResult { - let total_moles = moles.sum(); - let molefracs = (moles / total_moles).into_value(); + let (molefracs, total_moles) = composition.into_moles(eos)?; let density = total_moles / volume; - validate(temperature, density, &molefracs)?; + Self::new(eos, temperature, density, (molefracs, total_moles)) + } - Ok(Self::new_unchecked( - eos, - temperature, - density, - total_moles, - &molefracs, - )) + /// 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, + partial_density: Density>, + ) -> FeosResult { + let density = partial_density.sum(); + Self::new(eos, temperature, density, partial_density) } - /// Return a new `State` for which the total amount of substance is unspecified. + /// Return a new `State` for a pure component given a temperature and a density. /// - /// Internally the total number of moles will be set to 1 mol. + /// 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, density: Density) -> FeosResult + where + (): Composition, + { + 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_intensive( + pub fn new>( eos: &E, temperature: Temperature, density: Density, - molefracs: &OVector, + composition: X, ) -> FeosResult { - validate(temperature, density, molefracs)?; - let total_moles = Moles::new(D::one()); - Ok(Self::new_unchecked( - eos, - temperature, - density, - total_moles, - molefracs, - )) + let (molefracs, total_moles) = composition.into_molefracs(eos)?; + Self::_new(eos, temperature, density, molefracs, total_moles) } - fn new_unchecked( + fn _new( eos: &E, temperature: Temperature, density: Density, - total_moles: Moles, - molefracs: &OVector, - ) -> Self { - let volume = total_moles / density; - let moles = Dimensionless::new(molefracs.clone()) * total_moles; - let partial_density = moles.clone() / volume; - - State { + molefracs: OVector, + total_moles: Option>, + ) -> FeosResult { + let molar_volume = density.inv(); + validate(temperature, density, &molefracs)?; + Ok(State { eos: eos.clone(), temperature, - volume, - moles, - total_moles, - partial_density, + molar_volume, density, - molefracs: molefracs.clone(), + molefracs, + total_moles, cache: Cache::new(), - } - } - - /// Return a new `State` for a pure component given a temperature and a density. The moles - /// are set to the reference value for each component. - /// - /// 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, density: Density) -> FeosResult { - let molefracs = OVector::from_element_generic(N::from_usize(1), U1, D::one()); - Self::new_intensive(eos, temperature, density, &molefracs) + }) } /// Return a new `State` for the combination of inputs. @@ -313,200 +310,105 @@ where /// # Errors /// /// When the state cannot be created using the combination of inputs. - #[expect(clippy::too_many_arguments)] - pub fn new( + pub fn build>( eos: &E, - temperature: Option>, + temperature: Temperature, volume: Option>, density: Option>, - partial_density: Option<&Density>>, - total_moles: Option>, - moles: Option<&Moles>>, - molefracs: Option<&OVector>, + composition: X, pressure: Option>, density_initialization: Option, ) -> FeosResult { - Self::_new( + Self::_build( eos, temperature, volume, density, - partial_density, - total_moles, - moles, - molefracs, + composition, pressure, density_initialization, )? - .map_err(|_| FeosError::UndeterminedState(String::from("Missing input parameters."))) + .ok_or_else(|| FeosError::UndeterminedState(String::from("Missing input parameters."))) } - #[expect(clippy::too_many_arguments)] - #[expect(clippy::type_complexity)] - fn _new( + fn _build>( eos: &E, - temperature: Option>, + temperature: Temperature, volume: Option>, density: Option>, - partial_density: Option<&Density>>, - total_moles: Option>, - moles: Option<&Moles>>, - molefracs: Option<&OVector>, + composition: X, pressure: Option>, density_initialization: Option, - ) -> FeosResult>>>> { - // check for density - if density.and(partial_density).is_some() { + ) -> FeosResult> { + // check if density is given twice + if density.and(composition.density()).is_some() { return Err(FeosError::UndeterminedState(String::from( "Both density and partial density given.", ))); } - let rho = density.or_else(|| partial_density.map(|pd| pd.sum())); - - // check for total moles - if moles.and(total_moles).is_some() { - return Err(FeosError::UndeterminedState(String::from( - "Both moles and total moles given.", - ))); - } - let mut n = total_moles.or_else(|| moles.map(|m| m.sum())); - - // check if total moles can be inferred from volume - if rho.and(n).and(volume).is_some() { - return Err(FeosError::UndeterminedState(String::from( + let density = density.or_else(|| composition.density()); + + // 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.", + ))), } - n = n.or_else(|| rho.and_then(|d| volume.map(|v| v * d))); - - // check for composition - if partial_density.and(moles).is_some() { - return Err(FeosError::UndeterminedState(String::from( - "Composition is overdetermined.", - ))); - } - let x = partial_density - .map(|pd| pd / pd.sum()) - .or_else(|| moles.map(|ms| ms / ms.sum())) - .map(Quantity::into_value); - let x_u = match (x, molefracs, eos.components()) { - (Some(_), Some(_), _) => { - return Err(FeosError::UndeterminedState(String::from( - "Composition is overdetermined.", - ))); - } - (Some(x), None, _) => x, - (None, Some(x), _) => x.clone(), - (None, None, 1) => OVector::from_element_generic(N::from_usize(1), U1, D::from(1.0)), - _ => { - return Err(FeosError::UndeterminedState(String::from( - "Missing composition.", - ))); - } - }; - let x_u = &x_u / x_u.sum(); - - // If no extensive property is given, moles is set to the reference value. - if let (None, None) = (volume, n) { - n = Some(Moles::from_reduced(D::from(1.0))) - } - let n_i = n.map(|n| Dimensionless::new(&x_u) * n); - let v = volume.or_else(|| rho.and_then(|d| n.map(|n| n / d))); - - // check if new state can be created using default constructor - if let (Some(v), Some(t), Some(n_i)) = (v, temperature, &n_i) { - return Ok(Ok(State::new_nvt(eos, t, v, n_i)?)); - } - - // Check if new state can be created using density iteration - if let (Some(p), Some(t), Some(n_i)) = (pressure, temperature, &n_i) { - return Ok(Ok(State::new_npt(eos, t, p, n_i, density_initialization)?)); - } - if let (Some(p), Some(t), Some(v)) = (pressure, temperature, v) { - return Ok(Ok(State::new_npvx( - eos, - t, - p, - v, - &x_u, - density_initialization, - )?)); - } - Ok(Err(n_i.to_owned())) } /// 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( + pub fn new_npt>( eos: &E, temperature: Temperature, pressure: Pressure, - moles: &Moles>, - density_initialization: Option, - ) -> FeosResult { - let total_moles = moles.sum(); - let molefracs = (moles / total_moles).into_value(); - let density = Self::new_xpt( - eos, - temperature, - pressure, - &molefracs, - density_initialization, - )? - .density; - Ok(Self::new_unchecked( - eos, - temperature, - density, - total_moles, - &molefracs, - )) - } - - /// Return a new `State` using a density iteration. [DensityInitialization] is used to - /// influence the calculation with respect to the possible solutions. - pub fn new_xpt( - eos: &E, - temperature: Temperature, - pressure: Pressure, - molefracs: &OVector, + composition: X, density_initialization: Option, ) -> FeosResult { + let (molefracs, total_moles) = composition.into_molefracs(eos)?; density_iteration( eos, temperature, pressure, - molefracs, + &molefracs, density_initialization, ) - .and_then(|density| Self::new_intensive(eos, temperature, density, molefracs)) + .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_npvx( + pub fn new_tpvx( eos: &E, temperature: Temperature, pressure: Pressure, volume: Volume, - molefracs: &OVector, + molefracs: OVector, density_initialization: Option, ) -> FeosResult { - let density = Self::new_xpt( + let density = density_iteration( eos, temperature, pressure, - molefracs, + &molefracs, density_initialization, - )? - .density; - let total_moles = density * volume; - Ok(Self::new_unchecked( - eos, - temperature, - density, - total_moles, - molefracs, - )) + )?; + Self::new_nvt(eos, temperature, volume, (molefracs, density * volume)) } } @@ -529,15 +431,12 @@ where /// /// When the state cannot be created using the combination of inputs. #[expect(clippy::too_many_arguments)] - pub fn new_full( + pub fn build_full + Clone>( eos: &E, temperature: Option>, volume: Option>, density: Option>, - partial_density: Option<&Density>>, - total_moles: Option>, - moles: Option<&Moles>>, - molefracs: Option<&OVector>, + composition: X, pressure: Option>, molar_enthalpy: Option>, molar_entropy: Option>, @@ -545,38 +444,42 @@ where density_initialization: Option, initial_temperature: Option>, ) -> FeosResult { - let state = Self::_new( - eos, - temperature, - volume, - density, - partial_density, - total_moles, - moles, - molefracs, - pressure, - density_initialization, - )?; + 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 { - Ok(state) => Ok(state), - Err(n_i) => { + Some(state) => Ok(state), + None => { // Check if new state can be created using molar_enthalpy and temperature - if let (Some(p), Some(h), Some(n_i)) = (pressure, molar_enthalpy, &n_i) { - return State::new_nph(eos, p, h, n_i, density_initialization, ti); + if let (Some(p), Some(h)) = (pressure, molar_enthalpy) { + return State::new_nph(eos, p, h, composition, density_initialization, ti); } - if let (Some(p), Some(s), Some(n_i)) = (pressure, molar_entropy, &n_i) { - return State::new_nps(eos, p, s, n_i, density_initialization, ti); + if let (Some(p), Some(s)) = (pressure, molar_entropy) { + return State::new_nps(eos, p, s, composition, density_initialization, ti); } - if let (Some(t), Some(h), Some(n_i)) = (temperature, molar_enthalpy, &n_i) { - return State::new_nth(eos, t, h, n_i, density_initialization); + if let (Some(t), Some(h)) = (temperature, molar_enthalpy) { + return State::new_nth(eos, t, h, composition, density_initialization); } - if let (Some(t), Some(s), Some(n_i)) = (temperature, molar_entropy, &n_i) { - return State::new_nts(eos, t, s, n_i, density_initialization); + if let (Some(t), Some(s)) = (temperature, molar_entropy) { + return State::new_nts(eos, t, s, composition, density_initialization); } - if let (Some(u), Some(v), Some(n_i)) = (molar_internal_energy, volume, &n_i) { - return State::new_nvu(eos, v, u, n_i, ti); + if let (Some(u), Some(v)) = (molar_internal_energy, volume) { + let (molefracs, total_moles) = composition.into_molefracs(eos)?; + if let Some(n) = total_moles { + return State::new_nvu(eos, v, u, (molefracs, n), ti); + } } Err(FeosError::UndeterminedState(String::from( "Missing input parameters.", @@ -586,18 +489,18 @@ where } /// Return a new `State` for given pressure $p$ and molar enthalpy $h$. - pub fn new_nph( + pub fn new_nph + Clone>( eos: &E, pressure: Pressure, molar_enthalpy: MolarEnergy, - moles: &Moles>, + composition: X, density_initialization: Option, initial_temperature: Option>, ) -> FeosResult { 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, moles, density)?; + 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())); @@ -607,28 +510,27 @@ where } /// Return a new `State` for given temperature $T$ and molar enthalpy $h$. - pub fn new_nth( + pub fn new_nth + Clone>( eos: &E, temperature: Temperature, molar_enthalpy: MolarEnergy, - moles: &Moles>, + composition: X, density_initialization: Option, ) -> FeosResult { - let x = moles.convert_to(moles.sum()); + 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(&Some(x))?, - Some(DensityInitialization::Vapor) => eos.max_density(&Some(x))? * 1.0e-5, - None => eos.max_density(&Some(x))? * 0.01, + 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 n_inv = moles.sum().inv(); - let f = |x0| { - let s = State::new_nvt(eos, temperature, moles.sum() / x0, moles)?; - let dfx = -s.volume / s.density - * n_inv - * (s.volume * s.dp_dv(Contributions::Total) + 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)) @@ -637,26 +539,25 @@ where } /// Return a new `State` for given temperature $T$ and molar entropy $s$. - pub fn new_nts( + pub fn new_nts + Clone>( eos: &E, temperature: Temperature, molar_entropy: MolarEntropy, - moles: &Moles>, + composition: X, density_initialization: Option, ) -> FeosResult { - let x = moles.convert_to(moles.sum()); + 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(&Some(x))?, - Some(DensityInitialization::Vapor) => eos.max_density(&Some(x))? * 1.0e-5, - None => eos.max_density(&Some(x))? * 0.01, + 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 n_inv = moles.sum().inv(); - let f = |x0| { - let s = State::new_nvt(eos, temperature, moles.sum() / x0, moles)?; - let dfx = -n_inv * s.volume / s.density * s.dp_dt(Contributions::Total); + 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)) }; @@ -664,18 +565,18 @@ where } /// Return a new `State` for given pressure $p$ and molar entropy $s$. - pub fn new_nps( + pub fn new_nps + Clone>( eos: &E, pressure: Pressure, molar_entropy: MolarEntropy, - moles: &Moles>, + composition: X, density_initialization: Option, initial_temperature: Option>, ) -> FeosResult { 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, moles, density)?; + 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())); @@ -685,16 +586,16 @@ where } /// Return a new `State` for given volume $V$ and molar internal energy $u$. - pub fn new_nvu( + pub fn new_nvu + Clone>( eos: &E, volume: Volume, molar_internal_energy: MolarEnergy, - moles: &Moles>, + composition: X, initial_temperature: Option>, ) -> FeosResult { let t0 = initial_temperature.unwrap_or(Temperature::from_reduced(D::from(298.15))); let f = |x0| { - let s = State::new_nvt(eos, x0, volume, moles)?; + 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)) diff --git a/crates/feos-core/src/state/properties.rs b/crates/feos-core/src/state/properties.rs index f45bfb751..c93c1e77d 100644 --- a/crates/feos-core/src/state/properties.rs +++ b/crates/feos-core/src/state/properties.rs @@ -20,11 +20,11 @@ where let ideal_gas = || { quantity::ad::gradient_copy( partial2( - |n, &t, &v| self.eos.ideal_gas_helmholtz_energy(t, v, &n), + |n: Dimensionless<_>, &t, &v| self.eos.ideal_gas_helmholtz_energy(t, v, &n), &self.temperature, - &self.volume, + &self.molar_volume, ), - &self.moles, + &Dimensionless::new(self.molefracs.clone()), ) .1 }; @@ -37,10 +37,15 @@ where let ideal_gas = || { quantity::ad::partial_hessian_copy( partial( - |(n, t), &v| self.eos.ideal_gas_helmholtz_energy(t, v, &n), - &self.volume, + |(n, t): (Dimensionless<_>, _), &v| { + self.eos.ideal_gas_helmholtz_energy(t, v, &n) + }, + &self.molar_volume, + ), + ( + &Dimensionless::new(self.molefracs.clone()), + self.temperature, ), - (&self.moles, self.temperature), ) .3 }; @@ -49,7 +54,7 @@ where /// Molar isochoric heat capacity: $c_v=\left(\frac{\partial u}{\partial T}\right)_{V,N_i}$ pub fn molar_isochoric_heat_capacity(&self, contributions: Contributions) -> MolarEntropy { - self.temperature * self.ds_dt(contributions) / self.total_moles + self.temperature * self.ds_dt(contributions) } /// Partial derivative of the molar isochoric heat capacity w.r.t. temperature: $\left(\frac{\partial c_V}{\partial T}\right)_{V,N_i}$ @@ -57,8 +62,7 @@ where &self, contributions: Contributions, ) -> as Div>>::Output { - (self.temperature * self.d2s_dt2(contributions) + self.ds_dt(contributions)) - / self.total_moles + self.temperature * self.d2s_dt2(contributions) + self.ds_dt(contributions) } /// Molar isobaric heat capacity: $c_p=\left(\frac{\partial h}{\partial T}\right)_{p,N_i}$ @@ -66,7 +70,7 @@ where match contributions { Contributions::Residual => self.residual_molar_isobaric_heat_capacity(), _ => { - self.temperature / self.total_moles + self.temperature * (self.ds_dt(contributions) - (self.dp_dt(contributions) * self.dp_dt(contributions)) / self.dp_dv(contributions)) @@ -76,13 +80,18 @@ where /// Entropy: $S=-\left(\frac{\partial A}{\partial T}\right)_{V,N_i}$ pub fn entropy(&self, contributions: Contributions) -> Entropy { - let residual = || self.residual_entropy(); + self.molar_entropy(contributions) * self.total_moles() + } + + /// Molar entropy: $s=\frac{S}{N}$ + pub fn molar_entropy(&self, contributions: Contributions) -> MolarEntropy { + let residual = || self.residual_molar_entropy(); let ideal_gas = || { -quantity::ad::first_derivative( partial2( |t, &v, n| self.eos.ideal_gas_helmholtz_energy(t, v, n), - &self.volume, - &self.moles, + &self.molar_volume, + &self.molefracs, ), self.temperature, ) @@ -91,29 +100,24 @@ where Self::contributions(ideal_gas, residual, contributions) } - /// Molar entropy: $s=\frac{S}{N}$ - pub fn molar_entropy(&self, contributions: Contributions) -> MolarEntropy { - self.entropy(contributions) / self.total_moles - } - /// Partial molar entropy: $s_i=\left(\frac{\partial S}{\partial N_i}\right)_{T,p,N_j}$ pub fn partial_molar_entropy(&self) -> MolarEntropy> { let c = Contributions::Total; - -(self.dmu_dt(c) + self.dp_dni(c) * (self.dp_dt(c) / self.dp_dv(c))) + -(self.dmu_dt(c) + self.n_dp_dni(c) * (self.dp_dt(c) / self.dp_dv(c))) } - /// Partial derivative of the entropy w.r.t. temperature: $\left(\frac{\partial S}{\partial T}\right)_{V,N_i}$ + /// Partial derivative of the molar entropy w.r.t. temperature: $\left(\frac{\partial s}{\partial T}\right)_{V,N_i}$ pub fn ds_dt( &self, contributions: Contributions, - ) -> as Div>>::Output { + ) -> as Div>>::Output { let residual = || self.ds_res_dt(); let ideal_gas = || { -quantity::ad::second_derivative( partial2( |t, &v, n| self.eos.ideal_gas_helmholtz_energy(t, v, n), - &self.volume, - &self.moles, + &self.molar_volume, + &self.molefracs, ), self.temperature, ) @@ -122,18 +126,18 @@ where Self::contributions(ideal_gas, residual, contributions) } - /// Second partial derivative of the entropy w.r.t. temperature: $\left(\frac{\partial^2 S}{\partial T^2}\right)_{V,N_i}$ + /// Second partial derivative of the molar entropy w.r.t. temperature: $\left(\frac{\partial^2 s}{\partial T^2}\right)_{V,N_i}$ pub fn d2s_dt2( &self, contributions: Contributions, - ) -> < as Div>>::Output as Div>>::Output { + ) -> < as Div>>::Output as Div>>::Output { let residual = || self.d2s_res_dt2(); let ideal_gas = || { -quantity::ad::third_derivative( partial2( |t, &v, n| self.eos.ideal_gas_helmholtz_energy(t, v, n), - &self.volume, - &self.moles, + &self.molar_volume, + &self.molefracs, ), self.temperature, ) @@ -144,14 +148,14 @@ where /// Enthalpy: $H=A+TS+pV$ pub fn enthalpy(&self, contributions: Contributions) -> Energy { - self.temperature * self.entropy(contributions) - + self.helmholtz_energy(contributions) - + self.pressure(contributions) * self.volume + self.molar_enthalpy(contributions) * self.total_moles() } /// Molar enthalpy: $h=\frac{H}{N}$ pub fn molar_enthalpy(&self, contributions: Contributions) -> MolarEnergy { - self.enthalpy(contributions) / self.total_moles + self.temperature * self.molar_entropy(contributions) + + self.molar_helmholtz_energy(contributions) + + self.pressure(contributions) * self.molar_volume } /// Partial molar enthalpy: $h_i=\left(\frac{\partial H}{\partial N_i}\right)_{T,p,N_j}$ @@ -163,13 +167,18 @@ where /// Helmholtz energy: $A$ pub fn helmholtz_energy(&self, contributions: Contributions) -> Energy { - let residual = || self.residual_helmholtz_energy(); + self.molar_helmholtz_energy(contributions) * self.total_moles() + } + + /// Molar Helmholtz energy: $a=\frac{A}{N}$ + pub fn molar_helmholtz_energy(&self, contributions: Contributions) -> MolarEnergy { + let residual = || self.residual_molar_helmholtz_energy(); let ideal_gas = || { quantity::ad::zeroth_derivative( partial2( |t, &v, n| self.eos.ideal_gas_helmholtz_energy(t, v, n), - &self.volume, - &self.moles, + &self.molar_volume, + &self.molefracs, ), self.temperature, ) @@ -177,43 +186,40 @@ where Self::contributions(ideal_gas, residual, contributions) } - /// Molar Helmholtz energy: $a=\frac{A}{N}$ - pub fn molar_helmholtz_energy(&self, contributions: Contributions) -> MolarEnergy { - self.helmholtz_energy(contributions) / self.total_moles - } - /// Internal energy: $U=A+TS$ pub fn internal_energy(&self, contributions: Contributions) -> Energy { - self.temperature * self.entropy(contributions) + self.helmholtz_energy(contributions) + self.molar_internal_energy(contributions) * self.total_moles() } /// Molar internal energy: $u=\frac{U}{N}$ pub fn molar_internal_energy(&self, contributions: Contributions) -> MolarEnergy { - self.internal_energy(contributions) / self.total_moles + self.temperature * self.molar_entropy(contributions) + + self.molar_helmholtz_energy(contributions) } /// Gibbs energy: $G=A+pV$ pub fn gibbs_energy(&self, contributions: Contributions) -> Energy { - self.pressure(contributions) * self.volume + self.helmholtz_energy(contributions) + self.molar_gibbs_energy(contributions) * self.total_moles() } /// Molar Gibbs energy: $g=\frac{G}{N}$ pub fn molar_gibbs_energy(&self, contributions: Contributions) -> MolarEnergy { - self.gibbs_energy(contributions) / self.total_moles + self.pressure(contributions) * self.molar_volume + + self.molar_helmholtz_energy(contributions) } /// Joule Thomson coefficient: $\mu_{JT}=\left(\frac{\partial T}{\partial p}\right)_{H,N_i}$ pub fn joule_thomson(&self) -> as Div>>::Output { let c = Contributions::Total; - -(self.volume + self.temperature * self.dp_dt(c) / self.dp_dv(c)) - / (self.total_moles * self.molar_isobaric_heat_capacity(c)) + -(self.molar_volume + self.temperature * self.dp_dt(c) / self.dp_dv(c)) + / self.molar_isobaric_heat_capacity(c) } /// Isentropic compressibility: $\kappa_s=-\frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{S,N_i}$ pub fn isentropic_compressibility(&self) -> InvP { let c = Contributions::Total; -self.molar_isochoric_heat_capacity(c) - / (self.molar_isobaric_heat_capacity(c) * self.dp_dv(c) * self.volume) + / (self.molar_isobaric_heat_capacity(c) * self.dp_dv(c) * self.molar_volume) } /// Isenthalpic compressibility: $\kappa_H=-\frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{H,N_i}$ @@ -224,7 +230,7 @@ where /// Thermal expansivity: $\alpha_p=-\frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_{p,N_i}$ pub fn thermal_expansivity(&self) -> InvT { let c = Contributions::Total; - -self.dp_dt(c) / self.dp_dv(c) / self.volume + -self.dp_dt(c) / (self.dp_dv(c) * self.molar_volume) } /// Grueneisen parameter: $\phi=V\left(\frac{\partial p}{\partial U}\right)_{V,n_i}=\frac{v}{c_v}\left(\frac{\partial p}{\partial T}\right)_{v,n_i}=\frac{\rho}{T}\left(\frac{\partial T}{\partial \rho}\right)_{s, n_i}$ @@ -251,11 +257,7 @@ where )); } if let Contributions::Residual | Contributions::Total = contributions { - res.extend( - self.eos - .lift() - .molar_helmholtz_energy_contributions(t, v, &x), - ); + res.extend(self.eos.lift().helmholtz_energy_contributions(t, v, &x)); } res.into_iter() .map(|(s, v)| (s, MolarEnergy::from_reduced(v.eps))) diff --git a/crates/feos-core/src/state/residual_properties.rs b/crates/feos-core/src/state/residual_properties.rs index 179ff6d19..883438c6d 100644 --- a/crates/feos-core/src/state/residual_properties.rs +++ b/crates/feos-core/src/state/residual_properties.rs @@ -7,11 +7,8 @@ use num_dual::{Dual, DualNum, Gradients, partial, partial2}; use quantity::*; use std::ops::{Add, Div, Neg, Sub}; -type DpDn = Quantity>::Output>; -type DeDn = Quantity>::Output>; type InvT = Quantity::Output>; type InvP = Quantity::Output>; -type InvM = Quantity::Output>; type POverT = Quantity>::Output>; /// # State properties @@ -38,25 +35,33 @@ where /// Residual Helmholtz energy $A^\text{res}$ pub fn residual_helmholtz_energy(&self) -> Energy { - *self.cache.a.get_or_init(|| { - self.eos - .residual_helmholtz_energy_unit(self.temperature, self.volume, &self.moles) - }) + self.residual_molar_helmholtz_energy() * self.total_moles() } /// Residual molar Helmholtz energy $a^\text{res}$ pub fn residual_molar_helmholtz_energy(&self) -> MolarEnergy { - self.residual_helmholtz_energy() / self.total_moles + *self.cache.a.get_or_init(|| { + self.eos.residual_molar_helmholtz_energy( + self.temperature, + self.molar_volume, + &self.molefracs, + ) + }) } /// Residual entropy $S^\text{res}=\left(\frac{\partial A^\text{res}}{\partial T}\right)_{V,N_i}$ pub fn residual_entropy(&self) -> Entropy { + self.residual_molar_entropy() * self.total_moles() + } + + /// Residual molar entropy $s^\text{res}=\left(\frac{\partial a^\text{res}}{\partial T}\right)_{V,N_i}$ + pub fn residual_molar_entropy(&self) -> MolarEntropy { -*self.cache.da_dt.get_or_init(|| { let (a, da_dt) = quantity::ad::first_derivative( partial2( - |t, &v, n| self.eos.lift().residual_helmholtz_energy_unit(t, v, n), - &self.volume, - &self.moles, + |t, &v, n| self.eos.lift().residual_molar_helmholtz_energy(t, v, n), + &self.molar_volume, + &self.molefracs, ), self.temperature, ); @@ -65,11 +70,6 @@ where }) } - /// Residual entropy $s^\text{res}=\left(\frac{\partial a^\text{res}}{\partial T}\right)_{V,N_i}$ - pub fn residual_molar_entropy(&self) -> MolarEntropy { - self.residual_entropy() / self.total_moles - } - /// Pressure: $p=-\left(\frac{\partial A}{\partial V}\right)_{T,N_i}$ pub fn pressure(&self, contributions: Contributions) -> Pressure { let ideal_gas = || self.density * RGAS * self.temperature; @@ -77,11 +77,11 @@ where -*self.cache.da_dv.get_or_init(|| { let (a, da_dv) = quantity::ad::first_derivative( partial2( - |v, &t, n| self.eos.lift().residual_helmholtz_energy_unit(t, v, n), + |v, &t, n| self.eos.lift().residual_molar_helmholtz_energy(t, v, n), &self.temperature, - &self.moles, + &self.molefracs, ), - self.volume, + self.molar_volume, ); let _ = self.cache.a.set(a); da_dv @@ -97,11 +97,13 @@ where .get_or_init(|| { let (a, mu) = quantity::ad::gradient_copy( partial2( - |n, &t, &v| self.eos.lift().residual_helmholtz_energy_unit(t, v, &n), + |n: Dimensionless<_>, &t, &v| { + self.eos.lift().residual_molar_helmholtz_energy(t, v, &n) + }, &self.temperature, - &self.volume, + &self.molar_volume, ), - &self.moles, + &Dimensionless::new(self.molefracs.clone()), ); let _ = self.cache.a.set(a); mu @@ -116,18 +118,21 @@ where // pressure derivatives - /// Partial derivative of pressure w.r.t. volume: $\left(\frac{\partial p}{\partial V}\right)_{T,N_i}$ - pub fn dp_dv(&self, contributions: Contributions) -> as Div>>::Output { - let ideal_gas = || -self.density * RGAS * self.temperature / self.volume; + /// Partial derivative of pressure w.r.t. molar volume: $\left(\frac{\partial p}{\partial v}\right)_{T,N_i}$ + pub fn dp_dv( + &self, + contributions: Contributions, + ) -> as Div>>::Output { + let ideal_gas = || -self.density * RGAS * self.temperature / self.molar_volume; let residual = || { -*self.cache.d2a_dv2.get_or_init(|| { let (a, da_dv, d2a_dv2) = quantity::ad::second_derivative( partial2( - |v, &t, n| self.eos.lift().residual_helmholtz_energy_unit(t, v, n), + |v, &t, n| self.eos.lift().residual_molar_helmholtz_energy(t, v, n), &self.temperature, - &self.moles, + &self.molefracs, ), - self.volume, + self.molar_volume, ); let _ = self.cache.a.set(a); let _ = self.cache.da_dv.set(da_dv); @@ -142,7 +147,7 @@ where &self, contributions: Contributions, ) -> as Div>>::Output { - -self.volume / self.density * self.dp_dv(contributions) + -self.molar_volume / self.density * self.dp_dv(contributions) } /// Partial derivative of pressure w.r.t. temperature: $\left(\frac{\partial p}{\partial T}\right)_{V,N_i}$ @@ -152,10 +157,10 @@ where -*self.cache.d2a_dtdv.get_or_init(|| { let (a, da_dt, da_dv, d2a_dtdv) = quantity::ad::second_partial_derivative( partial( - |(t, v), n| self.eos.lift().residual_helmholtz_energy_unit(t, v, n), - &self.moles, + |(t, v), n| self.eos.lift().residual_molar_helmholtz_energy(t, v, n), + &self.molefracs, ), - (self.temperature, self.volume), + (self.temperature, self.molar_volume), ); let _ = self.cache.a.set(a); let _ = self.cache.da_dt.set(da_dt); @@ -166,18 +171,23 @@ where Self::contributions(ideal_gas, residual, contributions) } - /// Partial derivative of pressure w.r.t. moles: $\left(\frac{\partial p}{\partial N_i}\right)_{T,V,N_j}$ - pub fn dp_dni(&self, contributions: Contributions) -> DpDn> { + /// Partial derivative of pressure w.r.t. moles: $N\left(\frac{\partial p}{\partial N_i}\right)_{T,V,N_j}$ + pub fn n_dp_dni(&self, contributions: Contributions) -> Pressure> { let residual = -self .cache .d2a_dndv .get_or_init(|| { let (a, da_dn, da_dv, dmu_dv) = quantity::ad::partial_hessian_copy( partial( - |(n, v), &t| self.eos.lift().residual_helmholtz_energy_unit(t, v, &n), + |(n, v): (Dimensionless<_>, _), &t| { + self.eos.lift().residual_molar_helmholtz_energy(t, v, &n) + }, &self.temperature, ), - (&self.moles, self.volume), + ( + &Dimensionless::new(self.molefracs.clone()), + self.molar_volume, + ), ); let _ = self.cache.a.set(a); let _ = self.cache.da_dn.set(da_dn); @@ -186,7 +196,7 @@ where }) .clone(); let (r, c) = residual.shape_generic(); - let ideal_gas = || self.temperature / self.volume * RGAS; + let ideal_gas = || self.temperature * self.density * RGAS; Quantity::from_fn_generic(r, c, |i, _| { Self::contributions(ideal_gas, || residual.get(i), contributions) }) @@ -196,18 +206,19 @@ where pub fn d2p_dv2( &self, contributions: Contributions, - ) -> < as Div>>::Output as Div>>::Output { - let ideal_gas = - || self.density * RGAS * self.temperature / (self.volume * self.volume) * 2.0; + ) -> < as Div>>::Output as Div>>::Output { + let ideal_gas = || { + self.density * RGAS * self.temperature / (self.molar_volume * self.molar_volume) * 2.0 + }; let residual = || { -*self.cache.d3a_dv3.get_or_init(|| { let (a, da_dv, d2a_dv2, d3a_dv3) = quantity::ad::third_derivative( partial2( - |v, &t, n| self.eos.lift().residual_helmholtz_energy_unit(t, v, n), + |v, &t, n| self.eos.lift().residual_molar_helmholtz_energy(t, v, n), &self.temperature, - &self.moles, + &self.molefracs, ), - self.volume, + self.molar_volume, ); let _ = self.cache.a.set(a); let _ = self.cache.da_dv.set(da_dv); @@ -223,68 +234,62 @@ where &self, contributions: Contributions, ) -> < as Div>>::Output as Div>>::Output { - self.volume / (self.density * self.density) - * (self.volume * self.d2p_dv2(contributions) + self.dp_dv(contributions) * 2.0) + self.molar_volume.powi::<3>() + * (self.molar_volume * self.d2p_dv2(contributions) + self.dp_dv(contributions) * 2.0) } /// Structure factor: $S(0)=k_BT\left(\frac{\partial\rho}{\partial p}\right)_{T,N_i}$ pub fn structure_factor(&self) -> D { - -(self.temperature * self.density * RGAS / (self.volume * self.dp_dv(Contributions::Total))) - .into_value() - } - - // This function is designed specifically for use in density iterations - pub(crate) fn p_dpdrho(&self) -> (Pressure, as Div>>::Output) { - let dp_dv = self.dp_dv(Contributions::Total); - ( - self.pressure(Contributions::Total), - (-self.volume * dp_dv / self.density), - ) + -(self.temperature * self.density * RGAS + / (self.molar_volume * self.dp_dv(Contributions::Total))) + .into_value() } /// Partial molar volume: $v_i=\left(\frac{\partial V}{\partial N_i}\right)_{T,p,N_j}$ pub fn partial_molar_volume(&self) -> MolarVolume> { - -self.dp_dni(Contributions::Total) / self.dp_dv(Contributions::Total) + -self.n_dp_dni(Contributions::Total) / self.dp_dv(Contributions::Total) } - /// Partial derivative of chemical potential w.r.t. moles: $\left(\frac{\partial\mu_i}{\partial N_j}\right)_{T,V,N_k}$ - pub fn dmu_dni(&self, contributions: Contributions) -> DeDn> + /// Partial derivative of chemical potential w.r.t. moles: $N\left(\frac{\partial\mu_i}{\partial N_j}\right)_{T,V,N_k}$ + pub fn n_dmu_dni(&self, contributions: Contributions) -> MolarEnergy> where DefaultAllocator: Allocator, { let (a, da_dn, d2a_dn2) = quantity::ad::hessian_copy( partial2( - |n, &t, &v| self.eos.lift().residual_helmholtz_energy_unit(t, v, &n), + |n: Dimensionless<_>, &t, &v| { + self.eos.lift().residual_molar_helmholtz_energy(t, v, &n) + }, &self.temperature, - &self.volume, + &self.molar_volume, ), - &self.moles, + &Dimensionless::new(self.molefracs.clone()), ); let _ = self.cache.a.set(a); let _ = self.cache.da_dn.set(da_dn); let residual = || d2a_dn2; let ideal_gas = || { Dimensionless::new(OMatrix::from_diagonal(&self.molefracs.map(|x| x.recip()))) - * (self.temperature * RGAS / self.total_moles) + * (self.temperature * RGAS) }; Self::contributions(ideal_gas, residual, contributions) } /// Isothermal compressibility: $\kappa_T=-\frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{T,N_i}$ pub fn isothermal_compressibility(&self) -> InvP { - (self.dp_dv(Contributions::Total) * self.volume).inv() + (self.dp_dv(Contributions::Total) * self.molar_volume).inv() } // entropy derivatives - /// Partial derivative of the residual entropy w.r.t. temperature: $\left(\frac{\partial S^\text{res}}{\partial T}\right)_{V,N_i}$ - pub fn ds_res_dt(&self) -> as Div>>::Output { + /// Partial derivative of the residual molar entropy w.r.t. temperature: $\left(\frac{\partial s^\text{res}}{\partial T}\right)_{V,N_i}$ + pub fn ds_res_dt(&self) -> as Div>>::Output { -*self.cache.d2a_dt2.get_or_init(|| { let (a, da_dt, d2a_dt2) = quantity::ad::second_derivative( partial2( - |t, &v, n| self.eos.lift().residual_helmholtz_energy_unit(t, v, n), - &self.volume, - &self.moles, + |t, &v, n| self.eos.lift().residual_molar_helmholtz_energy(t, v, n), + &self.molar_volume, + &self.molefracs, ), self.temperature, ); @@ -294,16 +299,16 @@ where }) } - /// Second partial derivative of the residual entropy w.r.t. temperature: $\left(\frac{\partial^2S^\text{res}}{\partial T^2}\right)_{V,N_i}$ + /// Second partial derivative of the residual molar entropy w.r.t. temperature: $\left(\frac{\partial^2s^\text{res}}{\partial T^2}\right)_{V,N_i}$ pub fn d2s_res_dt2( &self, - ) -> < as Div>>::Output as Div>>::Output { + ) -> < as Div>>::Output as Div>>::Output { -*self.cache.d3a_dt3.get_or_init(|| { let (a, da_dt, d2a_dt2, d3a_dt3) = quantity::ad::third_derivative( partial2( - |t, &v, n| self.eos.lift().residual_helmholtz_energy_unit(t, v, n), - &self.volume, - &self.moles, + |t, &v, n| self.eos.lift().residual_molar_helmholtz_energy(t, v, n), + &self.molar_volume, + &self.molefracs, ), self.temperature, ); @@ -321,10 +326,15 @@ where .get_or_init(|| { let (a, da_dn, da_dt, d2a_dndt) = quantity::ad::partial_hessian_copy( partial( - |(n, t), &v| self.eos.lift().residual_helmholtz_energy_unit(t, v, &n), - &self.volume, + |(n, t): (Dimensionless<_>, _), &v| { + self.eos.lift().residual_molar_helmholtz_energy(t, v, &n) + }, + &self.molar_volume, + ), + ( + &Dimensionless::new(self.molefracs.clone()), + self.temperature, ), - (&self.moles, self.temperature), ); let _ = self.cache.a.set(a); let _ = self.cache.da_dn.set(da_dn); @@ -359,17 +369,19 @@ where .add_scalar(-self.pressure(Contributions::Total).inv()) } - /// Partial derivative of the logarithm of the fugacity coefficient w.r.t. moles: $\left(\frac{\partial\ln\varphi_i}{\partial N_j}\right)_{T,p,N_k}$ - pub fn dln_phi_dnj(&self) -> InvM> + /// Partial derivative of the logarithm of the fugacity coefficient w.r.t. moles: $N\left(\frac{\partial\ln\varphi_i}{\partial N_j}\right)_{T,p,N_k}$ + pub fn n_dln_phi_dnj(&self) -> OMatrix where DefaultAllocator: Allocator, { - let dmu_dni = self.dmu_dni(Contributions::Residual); - let dp_dni = self.dp_dni(Contributions::Total); + let dmu_dni = self.n_dmu_dni(Contributions::Residual); + let dp_dni = self.n_dp_dni(Contributions::Total); let dp_dv = self.dp_dv(Contributions::Total); let (r, c) = dmu_dni.shape_generic(); let dp_dn_2 = Quantity::from_fn_generic(r, c, |i, j| dp_dni.get(i) * dp_dni.get(j)); - ((dmu_dni + dp_dn_2 / dp_dv) / (self.temperature * RGAS)).add_scalar(self.total_moles.inv()) + ((dmu_dni + dp_dn_2 / dp_dv) / (self.temperature * RGAS)) + .into_value() + .add_scalar(D::from(1.0)) } } @@ -380,11 +392,11 @@ impl State { (0..self.eos.components()) .map(|i| { let eos = self.eos.subset(&[i]); - let state = State::new_xpt( + let state = State::new_npt( &eos, self.temperature, pressure, - &dvector![1.0], + dvector![1.0], Some(crate::DensityInitialization::Liquid), )?; Ok(state.ln_phi()[0]) @@ -419,26 +431,21 @@ impl State { .unzip(); let solvent_molefracs = DVector::from_vec(solvent_molefracs); let solvent = eos.subset(&solvent_comps); - let vle = if solvent_comps.len() == 1 { - PhaseEquilibrium::pure(&solvent, temperature, None, Default::default()) - } else { - PhaseEquilibrium::bubble_point( - &solvent, - temperature, - &solvent_molefracs, - None, - None, - Default::default(), - ) - }?; - - // Calculate the liquid state including the Henry components - let liquid = State::new_nvt( - eos, + // let vle = if solvent_comps.len() == 1 { + // PhaseEquilibrium::pure(&solvent, temperature, None, Default::default()) + // } else { + let vle = PhaseEquilibrium::bubble_point( + &solvent, temperature, - vle.liquid().volume, - &(molefracs * vle.liquid().total_moles), + &solvent_molefracs, + None, + None, + Default::default(), )?; + // }?; + + // Calculate the liquid state including the Henry components + let liquid = State::new(eos, temperature, vle.liquid().density, molefracs.clone())?; // Calculate the vapor state including the Henry components let mut molefracs_vapor = molefracs.clone(); @@ -446,12 +453,7 @@ impl State { .into_iter() .zip(&vle.vapor().molefracs) .for_each(|(i, &y)| molefracs_vapor[i] = y); - let vapor = State::new_nvt( - eos, - temperature, - vle.vapor().volume, - &(molefracs_vapor * vle.vapor().total_moles), - )?; + let vapor = State::new(eos, temperature, vle.vapor().density, molefracs.clone())?; // Determine the Henry's law coefficients and return only those of the Henry components let p = vle.vapor().pressure(Contributions::Total).into_reduced(); @@ -474,11 +476,10 @@ impl State { impl State { /// Thermodynamic factor: $\Gamma_{ij}=\delta_{ij}+x_i\left(\frac{\partial\ln\varphi_i}{\partial x_j}\right)_{T,p,\Sigma}$ pub fn thermodynamic_factor(&self) -> DMatrix { - let dln_phi_dnj = (self.dln_phi_dnj() * Moles::from_reduced(1.0)).into_value(); - let moles = &self.molefracs * self.total_moles.into_reduced(); + let dln_phi_dnj = self.n_dln_phi_dnj(); let n = self.eos.components() - 1; DMatrix::from_fn(n, n, |i, j| { - moles[i] * (dln_phi_dnj[(i, j)] - dln_phi_dnj[(i, n)]) + if i == j { 1.0 } else { 0.0 } + dln_phi_dnj[(i, j)] - dln_phi_dnj[(i, n)] + if i == j { 1.0 } else { 0.0 } }) } } @@ -489,63 +490,62 @@ where { /// Residual molar isochoric heat capacity: $c_v^\text{res}=\left(\frac{\partial u^\text{res}}{\partial T}\right)_{V,N_i}$ pub fn residual_molar_isochoric_heat_capacity(&self) -> MolarEntropy { - self.ds_res_dt() * self.temperature / self.total_moles + self.ds_res_dt() * self.temperature } /// Partial derivative of the residual molar isochoric heat capacity w.r.t. temperature: $\left(\frac{\partial c_V^\text{res}}{\partial T}\right)_{V,N_i}$ pub fn dc_v_res_dt(&self) -> as Div>>::Output { - (self.temperature * self.d2s_res_dt2() + self.ds_res_dt()) / self.total_moles + self.temperature * self.d2s_res_dt2() + self.ds_res_dt() } /// Residual molar isobaric heat capacity: $c_p^\text{res}=\left(\frac{\partial h^\text{res}}{\partial T}\right)_{p,N_i}$ pub fn residual_molar_isobaric_heat_capacity(&self) -> MolarEntropy { let dp_dt = self.dp_dt(Contributions::Total); - self.temperature / self.total_moles - * (self.ds_res_dt() - dp_dt * dp_dt / self.dp_dv(Contributions::Total)) + self.temperature * (self.ds_res_dt() - dp_dt * dp_dt / self.dp_dv(Contributions::Total)) - RGAS } /// Residual enthalpy: $H^\text{res}(T,p,\mathbf{n})=A^\text{res}+TS^\text{res}+p^\text{res}V$ pub fn residual_enthalpy(&self) -> Energy { - self.temperature * self.residual_entropy() - + self.residual_helmholtz_energy() - + self.pressure(Contributions::Residual) * self.volume + self.residual_molar_enthalpy() * self.total_moles() } /// Residual molar enthalpy: $h^\text{res}(T,p,\mathbf{n})=a^\text{res}+Ts^\text{res}+p^\text{res}v$ pub fn residual_molar_enthalpy(&self) -> MolarEnergy { - self.residual_enthalpy() / self.total_moles + self.temperature * self.residual_molar_entropy() + + self.residual_molar_helmholtz_energy() + + self.pressure(Contributions::Residual) * self.molar_volume } /// Residual internal energy: $U^\text{res}(T,V,\mathbf{n})=A^\text{res}+TS^\text{res}$ pub fn residual_internal_energy(&self) -> Energy { - self.temperature * self.residual_entropy() + self.residual_helmholtz_energy() + self.residual_molar_internal_energy() * self.total_moles() } /// Residual molar internal energy: $u^\text{res}(T,V,\mathbf{n})=a^\text{res}+Ts^\text{res}$ pub fn residual_molar_internal_energy(&self) -> MolarEnergy { - self.residual_internal_energy() / self.total_moles + self.temperature * self.residual_molar_entropy() + self.residual_molar_helmholtz_energy() } /// Residual Gibbs energy: $G^\text{res}(T,p,\mathbf{n})=A^\text{res}+p^\text{res}V-NRT \ln Z$ pub fn residual_gibbs_energy(&self) -> Energy { - self.pressure(Contributions::Residual) * self.volume + self.residual_helmholtz_energy() - - self.total_moles - * RGAS - * self.temperature - * Dimensionless::new(self.compressibility(Contributions::Total).ln()) + self.residual_molar_gibbs_energy() * self.total_moles() } /// Residual Gibbs energy: $g^\text{res}(T,p,\mathbf{n})=a^\text{res}+p^\text{res}v-RT \ln Z$ pub fn residual_molar_gibbs_energy(&self) -> MolarEnergy { - self.residual_gibbs_energy() / self.total_moles + self.pressure(Contributions::Residual) * self.molar_volume + + self.residual_molar_helmholtz_energy() + - self.temperature + * RGAS + * Dimensionless::new(self.compressibility(Contributions::Total).ln()) } /// Molar Helmholtz energy $a^\text{res}$ evaluated for each residual contribution of the equation of state. pub fn residual_molar_helmholtz_energy_contributions( &self, ) -> Vec<(&'static str, MolarEnergy)> { - let residual_contributions = self.eos.molar_helmholtz_energy_contributions( + let residual_contributions = self.eos.helmholtz_energy_contributions( self.temperature.into_reduced(), self.density.into_reduced().recip(), &self.molefracs, @@ -566,10 +566,7 @@ where let v = Dual::from_re(self.temperature.into_reduced()); let mut x = self.molefracs.map(Dual::from_re); x[component].eps = D::one(); - let contributions = self - .eos - .lift() - .molar_helmholtz_energy_contributions(t, v, &x); + let contributions = self.eos.lift().helmholtz_energy_contributions(t, v, &x); let mut res = Vec::with_capacity(contributions.len()); for (s, v) in contributions { res.push((s, MolarEnergy::from_reduced(v.eps))); @@ -582,10 +579,7 @@ where let t = Dual::from_re(self.temperature.into_reduced()); let v = Dual::from_re(self.density.into_reduced().recip()).derivative(); let x = self.molefracs.map(Dual::from_re); - let contributions = self - .eos - .lift() - .molar_helmholtz_energy_contributions(t, v, &x); + let contributions = self.eos.lift().helmholtz_energy_contributions(t, v, &x); let mut res = Vec::with_capacity(contributions.len() + 1); res.push(("Ideal gas", self.density * RGAS * self.temperature)); for (s, v) in contributions { @@ -608,12 +602,15 @@ where /// Mass of each component: $m_i=n_iMW_i$ pub fn mass(&self) -> Mass> { - self.eos.molar_weight().component_mul(&self.moles) + self.eos + .molar_weight() + .component_mul(&Dimensionless::new(self.molefracs.clone())) + * self.total_moles() } /// Total mass: $m=\sum_im_i=nMW$ pub fn total_mass(&self) -> Mass { - self.total_moles * self.total_molar_weight() + self.total_molar_weight() * self.total_moles() } /// Mass density: $\rho^{(m)}=\frac{m}{V}$ @@ -623,7 +620,10 @@ where /// Mass fractions: $w_i=\frac{m_i}{m}$ pub fn massfracs(&self) -> OVector { - (self.mass() / self.total_mass()).into_value() + self.eos + .molar_weight() + .convert_into(self.total_molar_weight()) + .component_mul(&self.molefracs) } } @@ -639,7 +639,7 @@ where pub fn viscosity(&self) -> Viscosity { let s = self.residual_molar_entropy().into_reduced(); self.eos - .viscosity_reference(self.temperature, self.volume, &self.moles) + .viscosity_reference(self.temperature, self.molar_volume, &self.molefracs) * Dimensionless::new(self.eos.viscosity_correlation(s, &self.molefracs).exp()) } @@ -655,14 +655,14 @@ where /// Return the viscosity reference as used in entropy scaling. pub fn viscosity_reference(&self) -> Viscosity { self.eos - .viscosity_reference(self.temperature, self.volume, &self.moles) + .viscosity_reference(self.temperature, self.molar_volume, &self.molefracs) } /// Return the diffusion via entropy scaling. pub fn diffusion(&self) -> Diffusivity { let s = self.residual_molar_entropy().into_reduced(); self.eos - .diffusion_reference(self.temperature, self.volume, &self.moles) + .diffusion_reference(self.temperature, self.molar_volume, &self.molefracs) * Dimensionless::new(self.eos.diffusion_correlation(s, &self.molefracs).exp()) } @@ -678,19 +678,21 @@ where /// Return the diffusion reference as used in entropy scaling. pub fn diffusion_reference(&self) -> Diffusivity { self.eos - .diffusion_reference(self.temperature, self.volume, &self.moles) + .diffusion_reference(self.temperature, self.molar_volume, &self.molefracs) } /// Return the thermal conductivity via entropy scaling. pub fn thermal_conductivity(&self) -> ThermalConductivity { let s = self.residual_molar_entropy().into_reduced(); - self.eos - .thermal_conductivity_reference(self.temperature, self.volume, &self.moles) - * Dimensionless::new( - self.eos - .thermal_conductivity_correlation(s, &self.molefracs) - .exp(), - ) + self.eos.thermal_conductivity_reference( + self.temperature, + self.molar_volume, + &self.molefracs, + ) * Dimensionless::new( + self.eos + .thermal_conductivity_correlation(s, &self.molefracs) + .exp(), + ) } /// Return the logarithm of the reduced thermal conductivity. @@ -705,7 +707,10 @@ where /// Return the thermal conductivity reference as used in entropy scaling. pub fn thermal_conductivity_reference(&self) -> ThermalConductivity { - self.eos - .thermal_conductivity_reference(self.temperature, self.volume, &self.moles) + self.eos.thermal_conductivity_reference( + self.temperature, + self.molar_volume, + &self.molefracs, + ) } } diff --git a/crates/feos-core/src/state/statevec.rs b/crates/feos-core/src/state/statevec.rs index a62fe9413..63d51b3fd 100644 --- a/crates/feos-core/src/state/statevec.rs +++ b/crates/feos-core/src/state/statevec.rs @@ -63,7 +63,7 @@ impl StateVec<'_, E> { pub fn moles(&self) -> Moles> { Moles::from_shape_fn((self.0.len(), self.0[0].eos.components()), |(i, j)| { - self.0[i].moles.get(j) + self.0[i].moles().get(j) }) } diff --git a/crates/feos-derive/src/residual.rs b/crates/feos-derive/src/residual.rs index 12e4fe862..df3f5367f 100644 --- a/crates/feos-derive/src/residual.rs +++ b/crates/feos-derive/src/residual.rs @@ -1,4 +1,4 @@ -use super::{implement, OPT_IMPLS}; +use super::{OPT_IMPLS, implement}; use quote::quote; use syn::DeriveInput; @@ -201,19 +201,19 @@ fn impl_entropy_scaling( let name = &v.ident; if implement("entropy_scaling", v, &OPT_IMPLS)? { etar.push(quote! { - Self::#name(eos) => eos.viscosity_reference(temperature, volume, moles) + Self::#name(eos) => eos.viscosity_reference(temperature, molar_volume, molefracs) }); etac.push(quote! { Self::#name(eos) => eos.viscosity_correlation(s_res, x) }); dr.push(quote! { - Self::#name(eos) => eos.diffusion_reference(temperature, volume, moles) + Self::#name(eos) => eos.diffusion_reference(temperature, molar_volume, molefracs) }); dc.push(quote! { Self::#name(eos) => eos.diffusion_correlation(s_res, x) }); thcr.push(quote! { - Self::#name(eos) => eos.thermal_conductivity_reference(temperature, volume, moles) + Self::#name(eos) => eos.thermal_conductivity_reference(temperature, molar_volume, molefracs) }); thcc.push(quote! { Self::#name(eos) => eos.thermal_conductivity_correlation(s_res, x) @@ -245,8 +245,8 @@ fn impl_entropy_scaling( fn viscosity_reference( &self, temperature: Temperature, - volume: Volume, - moles: &Moles>, + molar_volume: MolarVolume, + molefracs: &DVector, ) -> Viscosity { match self { #(#etar,)* @@ -262,8 +262,8 @@ fn impl_entropy_scaling( fn diffusion_reference( &self, temperature: Temperature, - volume: Volume, - moles: &Moles>, + molar_volume: MolarVolume, + molefracs: &DVector, ) -> Diffusivity { match self { #(#dr,)* @@ -279,8 +279,8 @@ fn impl_entropy_scaling( fn thermal_conductivity_reference( &self, temperature: Temperature, - volume: Volume, - moles: &Moles>, + molar_volume: MolarVolume, + molefracs: &DVector, ) -> ThermalConductivity { match self { #(#thcr,)* diff --git a/crates/feos-dft/src/adsorption/mod.rs b/crates/feos-dft/src/adsorption/mod.rs index 9962f27a5..bf9ddabe0 100644 --- a/crates/feos-dft/src/adsorption/mod.rs +++ b/crates/feos-dft/src/adsorption/mod.rs @@ -1,11 +1,12 @@ //! Adsorption profiles and isotherms. use super::functional::HelmholtzEnergyFunctional; use super::solver::DFTSolver; +use feos_core::DensityInitialization::{Liquid, Vapor}; use feos_core::{ - Contributions, DensityInitialization, FeosError, FeosResult, ReferenceSystem, SolverOptions, - State, StateBuilder, + Composition, Contributions, DensityInitialization, FeosError, FeosResult, ReferenceSystem, + SolverOptions, State, }; -use nalgebra::{DMatrix, DVector}; +use nalgebra::{DMatrix, DVector, Dyn}; use ndarray::{Array1, Array2, Dimension, Ix1, Ix3, RemoveAxis}; use quantity::{Energy, MolarEnergy, Moles, Pressure, Temperature}; use std::iter; @@ -53,12 +54,12 @@ where } /// Calculate an adsorption isotherm (starting at low pressure) - pub fn adsorption_isotherm>( + pub fn adsorption_isotherm, X: Composition + Clone>( functional: &F, temperature: Temperature, pressure: &Pressure>, pore: &S, - molefracs: &Option>, + composition: X, solver: Option<&DFTSolver>, ) -> FeosResult> { Self::isotherm( @@ -66,19 +67,19 @@ where temperature, pressure, pore, - molefracs, + composition, DensityInitialization::Vapor, solver, ) } /// Calculate an desorption isotherm (starting at high pressure) - pub fn desorption_isotherm>( + pub fn desorption_isotherm, X: Composition + Clone>( functional: &F, temperature: Temperature, pressure: &Pressure>, pore: &S, - molefracs: &Option>, + composition: X, solver: Option<&DFTSolver>, ) -> FeosResult> { let pressure = pressure.into_iter().rev().collect(); @@ -87,7 +88,7 @@ where temperature, &pressure, pore, - molefracs, + composition, DensityInitialization::Liquid, solver, )?; @@ -98,12 +99,12 @@ where } /// Calculate an equilibrium isotherm - pub fn equilibrium_isotherm>( + pub fn equilibrium_isotherm, X: Composition + Clone>( functional: &F, temperature: Temperature, pressure: &Pressure>, pore: &S, - molefracs: &Option>, + composition: X, solver: Option<&DFTSolver>, ) -> FeosResult> { let (p_min, p_max) = (pressure.get(0), pressure.get(pressure.len() - 1)); @@ -113,7 +114,7 @@ where p_min, p_max, pore, - molefracs, + composition.clone(), solver, SolverOptions::default(), ); @@ -132,7 +133,7 @@ where temperature, &p_ads, pore, - molefracs, + composition.clone(), solver, )? .profiles; @@ -141,7 +142,7 @@ where temperature, &p_des, pore, - molefracs, + composition, solver, )? .profiles; @@ -155,7 +156,7 @@ where temperature, pressure, pore, - molefracs, + composition.clone(), solver, )?; let desorption = Self::desorption_isotherm( @@ -163,7 +164,7 @@ where temperature, pressure, pore, - molefracs, + composition, solver, )?; let omega_a = adsorption.grand_potential(); @@ -181,25 +182,24 @@ where } } - fn isotherm>( + fn isotherm, X: Composition + Clone>( functional: &F, temperature: Temperature, pressure: &Pressure>, pore: &S, - molefracs: &Option>, + composition: X, density_initialization: DensityInitialization, solver: Option<&DFTSolver>, ) -> FeosResult> { - let x = functional.validate_molefracs(molefracs)?; let mut profiles: Vec>> = Vec::with_capacity(pressure.len()); // On the first iteration, initialize the density profile according to the direction // and calculate the external potential once. - let mut bulk = State::new_xpt( + let mut bulk = State::new_npt( functional, temperature, pressure.get(0), - &x, + composition.clone(), Some(density_initialization), )?; if functional.components() > 1 && !bulk.is_stable(SolverOptions::default())? { @@ -215,11 +215,13 @@ where let mut old_density = Some(&profile.density); for i in 0..pressure.len() { - let mut bulk = StateBuilder::new(functional) - .temperature(temperature) - .pressure(pressure.get(i)) - .molefracs(&x) - .build()?; + let mut bulk = State::new_npt( + functional, + temperature, + pressure.get(i), + composition.clone(), + None, + )?; if functional.components() > 1 && !bulk.is_stable(SolverOptions::default())? { bulk = bulk .tp_flash(None, SolverOptions::default(), None)? @@ -243,37 +245,21 @@ where /// Calculate the phase transition from an empty to a filled pore. #[expect(clippy::too_many_arguments)] - pub fn phase_equilibrium>( + pub fn phase_equilibrium, X: Composition + Clone>( functional: &F, temperature: Temperature, p_min: Pressure, p_max: Pressure, pore: &S, - molefracs: &Option>, + composition: X, solver: Option<&DFTSolver>, options: SolverOptions, ) -> FeosResult> { - let x = functional.validate_molefracs(molefracs)?; - + let x = composition; // calculate density profiles for the minimum and maximum pressure - let vapor_bulk = StateBuilder::new(functional) - .temperature(temperature) - .pressure(p_min) - .molefracs(&x) - .vapor() - .build()?; - let bulk_init = StateBuilder::new(functional) - .temperature(temperature) - .pressure(p_max) - .molefracs(&x) - .liquid() - .build()?; - let liquid_bulk = StateBuilder::new(functional) - .temperature(temperature) - .pressure(p_max) - .molefracs(&x) - .vapor() - .build()?; + let vapor_bulk = State::new_npt(functional, temperature, p_min, x.clone(), Some(Vapor))?; + let bulk_init = State::new_npt(functional, temperature, p_max, x.clone(), Some(Liquid))?; + let liquid_bulk = State::new_npt(functional, temperature, p_max, x.clone(), Some(Vapor))?; let mut vapor = pore.initialize(&vapor_bulk, None, None)?.solve(solver)?; let mut liquid = pore.initialize(&bulk_init, None, None)?.solve(solver)?; @@ -287,18 +273,13 @@ where / (n_dp_drho_v / vapor_bulk.density - n_dp_drho_l / liquid_bulk.density); // update filled pore with limited step size - let mut bulk = StateBuilder::new(functional) - .temperature(temperature) - .pressure(p_max) - .molefracs(&x) - .vapor() - .build()?; + let mut bulk = State::new_npt(functional, temperature, p_max, x.clone(), Some(Vapor))?; let rho0 = liquid_bulk.density; let steps = (10.0 * (rho - rho0) / rho0).into_value().abs().ceil() as usize; let delta_rho = (rho - rho0) / steps as f64; for i in 1..=steps { let rho_i = rho0 + i as f64 * delta_rho; - bulk = State::new_intensive(functional, temperature, rho_i, &x)?; + bulk = State::new(functional, temperature, rho_i, x.clone())?; liquid = liquid.update_bulk(&bulk).solve(solver)?; } @@ -324,7 +305,7 @@ where rho += delta_rho; // update bulk phase - bulk = State::new_intensive(functional, temperature, rho, &x)?; + bulk = State::new(functional, temperature, rho, x.clone())?; } Err(FeosError::NotConverged( "Adsorption::phase_equilibrium".into(), diff --git a/crates/feos-dft/src/adsorption/pore.rs b/crates/feos-dft/src/adsorption/pore.rs index faf19e62d..ddf2f7479 100644 --- a/crates/feos-dft/src/adsorption/pore.rs +++ b/crates/feos-dft/src/adsorption/pore.rs @@ -6,9 +6,7 @@ use crate::functional_contribution::FunctionalContribution; use crate::geometry::{Axis, Geometry, Grid}; use crate::profile::{DFTProfile, MAX_POTENTIAL}; use crate::solver::DFTSolver; -use feos_core::{ - Contributions, FeosResult, ReferenceSystem, ResidualDyn, State, StateBuilder, StateHD, -}; +use feos_core::{Contributions, FeosResult, ReferenceSystem, ResidualDyn, State, StateHD}; use nalgebra::{DVector, dvector}; use ndarray::prelude::*; use ndarray::{Axis as Axis_nd, RemoveAxis}; @@ -69,10 +67,7 @@ pub trait PoreSpecification { where D::Larger: Dimension, { - let bulk = StateBuilder::new(&&Helium) - .temperature(298.0 * KELVIN) - .density(Density::from_reduced(1.0)) - .build()?; + let bulk = State::new_pure(&&Helium, 298.0 * KELVIN, Density::from_reduced(1.0))?; let pore = self.initialize(&bulk, None, None)?; let pot = Dimensionless::from_reduced( pore.profile diff --git a/crates/feos-dft/src/interface/mod.rs b/crates/feos-dft/src/interface/mod.rs index 6697535d0..1e9f02b99 100644 --- a/crates/feos-dft/src/interface/mod.rs +++ b/crates/feos-dft/src/interface/mod.rs @@ -84,8 +84,8 @@ impl PlanarInterface { let reduced_temperature = (vle.vapor().temperature / critical_temperature).into_value(); profile.profile.density = Density::from_shape_fn(profile.profile.density.raw_dim(), |(i, z)| { - let rho_v = profile.vle.vapor().partial_density.get(indices[i]); - let rho_l = profile.vle.liquid().partial_density.get(indices[i]); + let rho_v = profile.vle.vapor().partial_density().get(indices[i]); + let rho_l = profile.vle.liquid().partial_density().get(indices[i]); 0.5 * (rho_l - rho_v) * (sign * (profile.profile.grid.grids()[0][z] - z0) / 3.0 * (2.4728 - 2.3625 * reduced_temperature)) @@ -343,8 +343,9 @@ fn interp_symmetric( radius: Length, ) -> FeosResult>> { let reduced_density = Array2::from_shape_fn(rho_pdgt.raw_dim(), |(i, j)| { - ((rho_pdgt.get((i, j)) - vle_pdgt.vapor().partial_density.get(i)) - / (vle_pdgt.liquid().partial_density.get(i) - vle_pdgt.vapor().partial_density.get(i))) + ((rho_pdgt.get((i, j)) - vle_pdgt.vapor().partial_density().get(i)) + / (vle_pdgt.liquid().partial_density().get(i) + - vle_pdgt.vapor().partial_density().get(i))) .into_value() - 0.5 }); @@ -371,8 +372,8 @@ fn interp_symmetric( reduced_density.raw_dim(), |(i, j)| { reduced_density[(i, j)] - * (vle.liquid().partial_density.get(i) - vle.vapor().partial_density.get(i)) - + vle.vapor().partial_density.get(i) + * (vle.liquid().partial_density().get(i) - vle.vapor().partial_density().get(i)) + + vle.vapor().partial_density().get(i) }, )) } diff --git a/crates/feos-dft/src/pdgt.rs b/crates/feos-dft/src/pdgt.rs index 47ad359eb..eb7cfe0f2 100644 --- a/crates/feos-dft/src/pdgt.rs +++ b/crates/feos-dft/src/pdgt.rs @@ -185,7 +185,7 @@ pub trait PdgtFunctionalProperties: HelmholtzEnergyFunctional { let mu_res = vle.vapor().residual_chemical_potential(); for i in 0..self.components() { let rhoi = density.index_axis(Axis(0), i).to_owned(); - let rhoi_b = vle.vapor().partial_density.get(i); + let rhoi_b = vle.vapor().partial_density().get(i); let mui_res = mu_res.get(i); let kt = RGAS * vle.vapor().temperature; delta_omega += @@ -198,8 +198,8 @@ pub trait PdgtFunctionalProperties: HelmholtzEnergyFunctional { let drho = gradient( &density, -dx, - &vle.liquid().partial_density, - &vle.vapor().partial_density, + &vle.liquid().partial_density(), + &vle.vapor().partial_density(), ); // calculate integrand diff --git a/crates/feos-dft/src/profile/mod.rs b/crates/feos-dft/src/profile/mod.rs index edac0e7b2..941a3818b 100644 --- a/crates/feos-dft/src/profile/mod.rs +++ b/crates/feos-dft/src/profile/mod.rs @@ -232,7 +232,7 @@ where let mut bonds = bulk.eos.bond_integrals(t, &exp_dfdrho, convolver.as_ref()); bonds *= &exp_dfdrho; let mut density = Array::zeros(external_potential.raw_dim()); - let bulk_density = bulk.partial_density.to_reduced(); + let bulk_density = bulk.partial_density().into_reduced(); for (s, &c) in bulk.eos.component_index().iter().enumerate() { density.index_axis_mut(Axis_nd(0), s).assign( &(bonds.index_axis(Axis_nd(0), s).map(|is| is.min(1.0)) * bulk_density[c]), @@ -373,7 +373,7 @@ where pub fn residual(&self, log: bool) -> FeosResult<(Array, Array1, f64)> { // Read from profile let density = self.density.to_reduced(); - let partial_density = self.bulk.partial_density.to_reduced(); + let partial_density = self.bulk.partial_density().into_reduced(); let bulk_density = self .bulk .eos @@ -484,7 +484,7 @@ where // Read from profile let component_index = self.bulk.eos.component_index().into_owned(); let mut density = self.density.to_reduced(); - let partial_density = self.bulk.partial_density.to_reduced(); + let partial_density = self.bulk.partial_density().into_reduced(); let mut bulk_density = component_index .iter() .map(|&i| partial_density[i]) @@ -495,13 +495,12 @@ where // Update profile self.density = Density::from_reduced(density); - let volume = Volume::from_reduced(1.0); - let mut moles = self.bulk.moles.clone(); + let mut partial_density = self.bulk.partial_density(); bulk_density .into_iter() .enumerate() - .for_each(|(i, r)| moles.set(component_index[i], Density::from_reduced(r) * volume)); - self.bulk = State::new_nvt(&self.bulk.eos, self.bulk.temperature, volume, &moles)?; + .for_each(|(i, r)| partial_density.set(component_index[i], Density::from_reduced(r))); + self.bulk = State::new_density(&self.bulk.eos, self.bulk.temperature, partial_density)?; Ok(()) } diff --git a/crates/feos-dft/src/profile/properties.rs b/crates/feos-dft/src/profile/properties.rs index 225967409..6095fb172 100644 --- a/crates/feos-dft/src/profile/properties.rs +++ b/crates/feos-dft/src/profile/properties.rs @@ -298,7 +298,7 @@ where { fn density_derivative(&self, lhs: &Array) -> FeosResult> { let rho = self.density.to_reduced(); - let partial_density = self.bulk.partial_density.to_reduced(); + let partial_density = self.bulk.partial_density().into_reduced(); let rho_bulk = self .bulk .eos @@ -402,7 +402,7 @@ where dfdrho += &(&self.external_potential * t).mapv(|v| Dual64::from(v) / t_dual); // calculate bulk functional derivative - let partial_density = self.bulk.partial_density.to_reduced(); + let partial_density = self.bulk.partial_density().into_reduced(); let rho_bulk: Array1<_> = self .bulk .eos diff --git a/crates/feos-dft/src/solvation/pair_correlation.rs b/crates/feos-dft/src/solvation/pair_correlation.rs index d901887cb..25849af4f 100644 --- a/crates/feos-dft/src/solvation/pair_correlation.rs +++ b/crates/feos-dft/src/solvation/pair_correlation.rs @@ -59,12 +59,10 @@ impl PairCorrelation { self.profile.solve(solver, debug)?; // calculate pair correlation function + let partial_density = self.profile.bulk.partial_density(); self.pair_correlation_function = Some(Array::from_shape_fn( self.profile.density.raw_dim(), - |(i, j)| { - (self.profile.density.get((i, j)) / self.profile.bulk.partial_density.get(i)) - .into_value() - }, + |(i, j)| (self.profile.density.get((i, j)) / partial_density.get(i)).into_value(), )); // calculate self solvation free energy diff --git a/crates/feos/Cargo.toml b/crates/feos/Cargo.toml index 28e6e6109..43c429a76 100644 --- a/crates/feos/Cargo.toml +++ b/crates/feos/Cargo.toml @@ -37,7 +37,7 @@ quantity = { workspace = true, features = ["approx"] } criterion = { workspace = true } [features] -default = [] +default = ["pcsaft"] dft = ["feos-dft", "petgraph", "ndarray", "num-dual/ndarray", "feos-derive"] association = [] pcsaft = ["association"] diff --git a/crates/feos/benches/contributions.rs b/crates/feos/benches/contributions.rs index 3b2e8b367..3e74bf68c 100644 --- a/crates/feos/benches/contributions.rs +++ b/crates/feos/benches/contributions.rs @@ -74,7 +74,7 @@ fn pcsaft(c: &mut Criterion) { State::new_npt(&&eos, t, p, &moles, Some(DensityInitialization::Liquid)).unwrap(); let temperature = Dual64::from(state.temperature.into_reduced()).derivative(); let molar_volume = Dual::from(1.0 / state.density.into_reduced()); - let moles = state.moles.to_reduced().map(Dual::from); + let moles = state.moles().to_reduced().map(Dual::from); // let state_hd = state.derive1(Derivative::DT); let name1 = comp1.identifier.name.as_deref().unwrap(); let name2 = comp2.identifier.name.as_deref().unwrap(); diff --git a/crates/feos/benches/dft_pore.rs b/crates/feos/benches/dft_pore.rs index 667c5ff12..98dcba4e2 100644 --- a/crates/feos/benches/dft_pore.rs +++ b/crates/feos/benches/dft_pore.rs @@ -2,7 +2,7 @@ //! in pores at different conditions. use criterion::{Criterion, criterion_group, criterion_main}; use feos::core::parameter::IdentifierOption; -use feos::core::{PhaseEquilibrium, State, StateBuilder}; +use feos::core::{PhaseEquilibrium, State}; use feos::dft::adsorption::{ExternalPotential, Pore1D, PoreSpecification}; use feos::dft::{DFTSolver, Geometry}; use feos::gc_pcsaft::{GcPcSaftFunctional, GcPcSaftParameters}; @@ -80,11 +80,8 @@ fn pcsaft(c: &mut Criterion) { group.bench_function("butane_pentane_liquid", |b| { b.iter(|| pore.initialize(bulk, None, None).unwrap().solve(None)) }); - let bulk = StateBuilder::new(&func) - .temperature(300.0 * KELVIN) - .partial_density(&(&vle.vapor().partial_density * 0.2)) - .build() - .unwrap(); + let bulk = + State::new_density(&func, 300.0 * KELVIN, vle.vapor().partial_density() * 0.2).unwrap(); group.bench_function("butane_pentane_vapor", |b| { b.iter(|| pore.initialize(&bulk, None, None).unwrap().solve(None)) }); diff --git a/crates/feos/benches/dual_numbers.rs b/crates/feos/benches/dual_numbers.rs index 454c6828f..c3c1bd07f 100644 --- a/crates/feos/benches/dual_numbers.rs +++ b/crates/feos/benches/dual_numbers.rs @@ -21,9 +21,9 @@ use quantity::*; fn state_pcsaft(n: usize, eos: &PcSaft) -> State<&PcSaft> { let moles = DVector::from_element(n, 1.0 / n as f64) * 10.0 * MOL; let molefracs = (&moles / moles.sum()).into_value(); - let cp = State::critical_point(&eos, Some(&molefracs), None, None, Default::default()).unwrap(); + let cp = State::critical_point(&eos, molefracs, None, None, Default::default()).unwrap(); let temperature = 0.8 * cp.temperature; - State::new_nvt(&eos, temperature, cp.volume, &moles).unwrap() + State::new_nvt(&eos, temperature, cp.volume(), moles).unwrap() } /// Residual Helmholtz energy given an equation of state and a StateHD. @@ -152,11 +152,11 @@ enum Derivative { /// Creates a [StateHD] cloning temperature, volume and moles. fn derive0(state: &State) -> StateHD { - let total_moles = state.total_moles.into_reduced(); + let total_moles = state.total_moles().into_reduced(); StateHD::new( state.temperature.into_reduced(), - state.volume.into_reduced() / total_moles, - &(state.moles.to_reduced() / total_moles), + state.volume().into_reduced() / total_moles, + &(state.moles().to_reduced() / total_moles), ) } diff --git a/crates/feos/benches/dual_numbers_saftvrmie.rs b/crates/feos/benches/dual_numbers_saftvrmie.rs index 0d6329974..e7409e307 100644 --- a/crates/feos/benches/dual_numbers_saftvrmie.rs +++ b/crates/feos/benches/dual_numbers_saftvrmie.rs @@ -18,9 +18,9 @@ use quantity::*; /// - molefracs (or moles) for equimolar mixture. fn state_saftvrmie(n: usize, eos: &SaftVRMie) -> State<&SaftVRMie> { let molefracs = DVector::from_element(n, 1.0 / n as f64); - let cp = State::critical_point(&eos, Some(&molefracs), None, None, Default::default()).unwrap(); + let cp = State::critical_point(&eos, &molefracs, None, None, Default::default()).unwrap(); let temperature = 0.8 * cp.temperature; - State::new_nvt(&eos, temperature, cp.volume, &(molefracs * 10. * MOL)).unwrap() + State::new_nvt(&eos, temperature, cp.volume(), &(molefracs * 10. * MOL)).unwrap() } /// Residual Helmholtz energy given an equation of state and a StateHD. @@ -100,11 +100,11 @@ enum Derivative { /// Creates a [StateHD] cloning temperature, volume and moles. fn derive0(state: &State) -> StateHD { - let total_moles = state.total_moles.into_reduced(); + let total_moles = state.total_moles().into_reduced(); StateHD::new( state.temperature.into_reduced(), - state.volume.into_reduced() / total_moles, - &(state.moles.to_reduced() / total_moles), + state.volume().into_reduced() / total_moles, + &(state.moles().to_reduced() / total_moles), ) } diff --git a/crates/feos/benches/state_creation.rs b/crates/feos/benches/state_creation.rs index 6110abd68..3aa08778e 100644 --- a/crates/feos/benches/state_creation.rs +++ b/crates/feos/benches/state_creation.rs @@ -22,7 +22,7 @@ fn npt( } /// Evaluate critical point constructor -fn critical_point((eos, n): (&E, Option<&DVector>)) { +fn critical_point((eos, n): (&E, &DVector)) { State::critical_point(eos, n, None, None, Default::default()).unwrap(); } @@ -69,7 +69,7 @@ fn bench_states(c: &mut Criterion, group_name: &str, eos: &E) { let ncomponents = eos.components(); let x = DVector::from_element(ncomponents, 1.0 / ncomponents as f64); let n = &x * 100.0 * MOL; - let crit = State::critical_point(eos, Some(&x), None, None, Default::default()).unwrap(); + let crit = State::critical_point(eos, &x, None, None, Default::default()).unwrap(); let vle = if ncomponents == 1 { PhaseEquilibrium::pure(eos, crit.temperature * 0.95, None, Default::default()).unwrap() } else { @@ -77,7 +77,7 @@ fn bench_states(c: &mut Criterion, group_name: &str, eos: &E) { eos, crit.temperature, crit.pressure(Contributions::Total) * 0.95, - &crit.moles, + &crit.moles(), None, Default::default(), None, @@ -108,9 +108,7 @@ fn bench_states(c: &mut Criterion, group_name: &str, eos: &E) { )) }) }); - group.bench_function("critical_point", |b| { - b.iter(|| critical_point((eos, Some(&x)))) - }); + group.bench_function("critical_point", |b| b.iter(|| critical_point((eos, &x)))); if ncomponents == 2 { group.bench_function("critical_point_binary_t", |b| { b.iter(|| critical_point_binary((eos, crit.temperature))) diff --git a/crates/feos/src/epcsaft/eos/mod.rs b/crates/feos/src/epcsaft/eos/mod.rs index 03a2713d0..2edba6f69 100644 --- a/crates/feos/src/epcsaft/eos/mod.rs +++ b/crates/feos/src/epcsaft/eos/mod.rs @@ -181,7 +181,7 @@ mod tests { let v = 1e-3 * METER.powi::<3>(); let n = dvector![1.0] * MOL; let s = State::new_nvt(&&e, t, v, &n).unwrap(); - let p_ig = s.total_moles * RGAS * t / v; + let p_ig = s.total_moles() * RGAS * t / v; assert_relative_eq!(s.pressure(Contributions::IdealGas), p_ig, epsilon = 1e-10); assert_relative_eq!( s.pressure(Contributions::IdealGas) + s.pressure(Contributions::Residual), @@ -197,7 +197,7 @@ mod tests { let v = 1e-3 * METER.powi::<3>(); let n = dvector![1.0] * MOL; let s = State::new_nvt(&&e, t, v, &n).unwrap(); - let p_ig = s.total_moles * RGAS * t / v; + let p_ig = s.total_moles() * RGAS * t / v; assert_relative_eq!(s.pressure(Contributions::IdealGas), p_ig, epsilon = 1e-10); assert_relative_eq!( s.pressure(Contributions::IdealGas) + s.pressure(Contributions::Residual), @@ -269,7 +269,7 @@ mod tests { fn critical_point() { let e = ElectrolytePcSaft::new(propane_parameters()).unwrap(); let t = 300.0 * KELVIN; - let cp = State::critical_point(&&e, None, Some(t), None, Default::default()); + let cp = State::critical_point(&&e, (), Some(t), None, Default::default()); if let Ok(v) = cp { assert_relative_eq!(v.temperature, 375.1244078318015 * KELVIN, epsilon = 1e-8) } diff --git a/crates/feos/src/gc_pcsaft/eos/ad.rs b/crates/feos/src/gc_pcsaft/eos/ad.rs index 252f9a061..72b251eaa 100644 --- a/crates/feos/src/gc_pcsaft/eos/ad.rs +++ b/crates/feos/src/gc_pcsaft/eos/ad.rs @@ -510,7 +510,7 @@ pub mod test { let eos_ad = GcPcSaftAD(params); let moles = vector![1.3] * KILO * MOL; - let state = State::new_nvt(&eos_ad, temperature, volume, &moles)?; + let state = State::new_nvt(&eos_ad, temperature, volume, moles)?; let a_ad = state.residual_molar_helmholtz_energy(); let mu_ad = state.residual_chemical_potential(); let p_ad = state.pressure(Total); diff --git a/crates/feos/src/ideal_gas/dippr.rs b/crates/feos/src/ideal_gas/dippr.rs index 7553e376c..89824b07a 100644 --- a/crates/feos/src/ideal_gas/dippr.rs +++ b/crates/feos/src/ideal_gas/dippr.rs @@ -156,7 +156,7 @@ impl IdealGas for Dippr { mod tests { use approx::assert_relative_eq; use feos_core::parameter::{Identifier, PureRecord}; - use feos_core::{Contributions, EquationOfState, StateBuilder}; + use feos_core::{Contributions, EquationOfState, State}; use num_dual::first_derivative; use quantity::*; @@ -173,11 +173,7 @@ mod tests { let eos = EquationOfState::ideal_gas(dippr.clone()); let temperature = 300.0 * KELVIN; let volume = METER.powi::<3>(); - let state = StateBuilder::new(&&eos) - .temperature(temperature) - .volume(volume) - .total_moles(MOL) - .build()?; + let state = State::new_nvt(&&eos, temperature, volume, MOL)?; let t = temperature.convert_to(KELVIN); let c_p_direct = record.model_record.c_p(t); @@ -217,11 +213,7 @@ mod tests { let eos = EquationOfState::ideal_gas(dippr.clone()); let temperature = 300.0 * KELVIN; let volume = METER.powi::<3>(); - let state = StateBuilder::new(&&eos) - .temperature(temperature) - .volume(volume) - .total_moles(MOL) - .build()?; + let state = State::new_nvt(&&eos, temperature, volume, MOL)?; let t = temperature.convert_to(KELVIN); let c_p_direct = record.model_record.c_p(t); @@ -263,11 +255,7 @@ mod tests { let eos = EquationOfState::ideal_gas(dippr.clone()); let temperature = 20.0 * KELVIN; let volume = METER.powi::<3>(); - let state = StateBuilder::new(&&eos) - .temperature(temperature) - .volume(volume) - .total_moles(MOL) - .build()?; + let state = State::new_nvt(&&eos, temperature, volume, MOL)?; let t = temperature.convert_to(KELVIN); let c_p_direct = record.model_record.c_p(t); diff --git a/crates/feos/src/ideal_gas/joback.rs b/crates/feos/src/ideal_gas/joback.rs index 742c6ef2e..8f2279d5b 100644 --- a/crates/feos/src/ideal_gas/joback.rs +++ b/crates/feos/src/ideal_gas/joback.rs @@ -172,10 +172,8 @@ const KB: f64 = 1.38064852e-23; #[cfg(test)] mod tests { use approx::assert_relative_eq; - use feos_core::{ - Contributions, EquationOfState, State, StateBuilder, - parameter::{ChemicalRecord, GroupCount, Identifier, PureRecord, SegmentRecord}, - }; + use feos_core::parameter::{ChemicalRecord, GroupCount, Identifier, PureRecord, SegmentRecord}; + use feos_core::{Contributions, EquationOfState, State}; use nalgebra::dvector; use quantity::*; use std::collections::HashMap; @@ -295,11 +293,7 @@ mod tests { let temperature = 300.0 * KELVIN; let volume = METER.powi::<3>(); let moles = &dvector![1.0, 3.0] * MOL; - let state = StateBuilder::new(&&eos) - .temperature(temperature) - .volume(volume) - .moles(&moles) - .build()?; + let state = State::new_nvt(&&eos, temperature, volume, moles)?; println!( "{} {}", Joback::molar_isobaric_heat_capacity(&joback, temperature, &state.molefracs)?, diff --git a/crates/feos/src/lib.rs b/crates/feos/src/lib.rs index eab5ae128..0fbdffa71 100644 --- a/crates/feos/src/lib.rs +++ b/crates/feos/src/lib.rs @@ -23,7 +23,7 @@ //! let saft = PcSaft::new(parameters); //! //! // Define thermodynamic conditions. -//! let critical_point = State::critical_point(&&saft, Some(&dvector![1.0]), None, None, Default::default())?; +//! let critical_point = State::critical_point(&&saft, dvector![1.0], None, None, Default::default())?; //! //! // Compute properties. //! let p = critical_point.pressure(Contributions::Total); diff --git a/crates/feos/src/multiparameter/mod.rs b/crates/feos/src/multiparameter/mod.rs index 68bbd98b8..b4a6c4f9f 100644 --- a/crates/feos/src/multiparameter/mod.rs +++ b/crates/feos/src/multiparameter/mod.rs @@ -265,8 +265,13 @@ mod test { let eos = &water(); let mw = eos.molar_weight.get(0); let moles = dvector![1.8] * MOL; - let a_feos = eos.ideal_gas_helmholtz_energy(t, moles.sum() * mw / rho, &moles); - let phi_feos = (a_feos / RGAS / moles.sum() / t).into_value(); + let total_moles = moles.sum(); + let a_feos = eos.ideal_gas_helmholtz_energy( + t, + moles.sum() * mw / rho / total_moles, + &moles.convert_into(total_moles), + ); + let phi_feos = (a_feos / RGAS / t).into_value(); println!("A: {a_feos}"); println!("phi(feos): {phi_feos}"); let delta = (rho / (eos.rhoc * MOL / METER.powi::<3>() * mw)).into_value(); @@ -288,15 +293,15 @@ mod test { ..Default::default() }; let cp: State<_, Dyn, f64> = - State::critical_point(&eos, None, Some(647. * KELVIN), None, options).unwrap(); + State::critical_point(&eos, (), Some(647. * KELVIN), None, options).unwrap(); println!("{cp}"); assert_relative_eq!(cp.temperature, eos.tc * KELVIN, max_relative = 1e-13); let cp: State<_, Dyn, f64> = - State::critical_point(&eos, None, None, None, Default::default()).unwrap(); + State::critical_point(&eos, (), None, None, Default::default()).unwrap(); println!("{cp}"); assert_relative_ne!(cp.temperature, eos.tc * KELVIN, max_relative = 1e-13); let cp: State<_, Dyn, f64> = - State::critical_point(&eos, None, Some(700.0 * KELVIN), None, Default::default()) + State::critical_point(&eos, (), Some(700.0 * KELVIN), None, Default::default()) .unwrap(); println!("{cp}"); assert_relative_eq!(cp.temperature, eos.tc * KELVIN, max_relative = 1e-13) diff --git a/crates/feos/src/pcsaft/eos/mod.rs b/crates/feos/src/pcsaft/eos/mod.rs index 3a20c90e4..89bcffdb4 100644 --- a/crates/feos/src/pcsaft/eos/mod.rs +++ b/crates/feos/src/pcsaft/eos/mod.rs @@ -229,12 +229,11 @@ impl EntropyScaling for PcSaft { fn viscosity_reference( &self, temperature: Temperature, - _: Volume, - moles: &Moles>, + _: MolarVolume, + molefracs: &DVector, ) -> Viscosity { let p = &self.params; let mw = &self.parameters.molar_weight; - let x = (moles / moles.sum()).into_value(); let ce: Vec<_> = (0..self.components()) .map(|i| { let tr = (temperature / p.epsilon_k[i] / KELVIN).into_value(); @@ -247,14 +246,15 @@ impl EntropyScaling for PcSaft { for i in 0..self.components() { let denom: f64 = (0..self.components()) .map(|j| { - x[j] * (1.0 - + (ce[i] / ce[j]).into_value().sqrt() - * (mw.get(j) / mw.get(i)).powf(1.0 / 4.0)) - .powi(2) + molefracs[j] + * (1.0 + + (ce[i] / ce[j]).into_value().sqrt() + * (mw.get(j) / mw.get(i)).powf(1.0 / 4.0)) + .powi(2) / (8.0 * (1.0 + (mw.get(i) / mw.get(j)).into_value())).sqrt() }) .sum(); - ce_mix += ce[i] * x[i] / denom + ce_mix += ce[i] * molefracs[i] / denom } ce_mix } @@ -278,19 +278,19 @@ impl EntropyScaling for PcSaft { fn diffusion_reference( &self, temperature: Temperature, - volume: Volume, - moles: &Moles>, + molar_volume: MolarVolume, + _: &DVector, ) -> Diffusivity { if self.components() != 1 { panic!("Diffusion coefficients in PC-SAFT are only implemented for pure components!"); } let p = &self.params; let mw = &self.parameters.molar_weight; - let density = moles.sum() / volume; let res: Vec<_> = (0..self.components()) .map(|i| { let tr = (temperature / p.epsilon_k[i] / KELVIN).into_value(); - 3.0 / 8.0 / (p.sigma[i] * ANGSTROM).powi::<2>() / omega11(tr) / (density * NAV) + 3.0 / 8.0 / (p.sigma[i] * ANGSTROM).powi::<2>() / omega11(tr) + * (molar_volume / NAV) * (temperature * RGAS / PI / mw.get(i) / p.m[i]).sqrt() }) .collect(); @@ -321,8 +321,8 @@ impl EntropyScaling for PcSaft { fn thermal_conductivity_reference( &self, temperature: Temperature, - volume: Volume, - moles: &Moles>, + molar_volume: MolarVolume, + molefracs: &DVector, ) -> ThermalConductivity { if self.components() != 1 { panic!("Thermal conductivity in PC-SAFT is only implemented for pure components!"); @@ -331,9 +331,9 @@ impl EntropyScaling for PcSaft { let mws = self.molar_weight(); let (_, s_res) = first_derivative( partial2( - |t, &v, n| -self.residual_helmholtz_energy_unit(t, v, n) / n.sum(), - &volume, - moles, + |t, &v, n| -self.residual_molar_helmholtz_energy(t, v, n), + &molar_volume, + molefracs, ), temperature, ); @@ -399,8 +399,8 @@ mod tests { let t = 200.0 * KELVIN; let v = 1e-3 * METER.powi::<3>(); let n = dvector![1.0] * MOL; - let s = State::new_nvt(&e, t, v, &n).unwrap(); - let p_ig = s.total_moles * RGAS * t / v; + let s = State::new_nvt(&e, t, v, n).unwrap(); + let p_ig = s.total_moles() * RGAS * t / v; assert_relative_eq!(s.pressure(Contributions::IdealGas), p_ig, epsilon = 1e-10); assert_relative_eq!( s.pressure(Contributions::IdealGas) + s.pressure(Contributions::Residual), @@ -415,8 +415,8 @@ mod tests { let t = 200.0 * KELVIN; let v = 1e-3 * METER.powi::<3>(); let n = dvector![1.0] * MOL; - let s = State::new_nvt(&e, t, v, &n).unwrap(); - let p_ig = s.total_moles * RGAS * t / v; + let s = State::new_nvt(&e, t, v, n).unwrap(); + let p_ig = s.total_moles() * RGAS * t / v; assert_relative_eq!(s.pressure(Contributions::IdealGas), p_ig, epsilon = 1e-10); assert_relative_eq!( s.pressure(Contributions::IdealGas) + s.pressure(Contributions::Residual), @@ -463,7 +463,7 @@ mod tests { let t = 300.0 * KELVIN; let p = BAR; let m = dvector![1.5] * MOL; - let s = State::new_npt(&e, t, p, &m, None); + let s = State::new_npt(&e, t, p, m, None); let p_calc = if let Ok(state) = s { state.pressure(Contributions::Total) } else { @@ -490,7 +490,7 @@ mod tests { fn critical_point() { let e = &propane_parameters(); let t = 300.0 * KELVIN; - let cp = State::critical_point(&e, None, Some(t), None, Default::default()); + let cp = State::critical_point(&e, (), Some(t), None, Default::default()); if let Ok(v) = cp { assert_relative_eq!(v.temperature, 375.1244078318015 * KELVIN, epsilon = 1e-8) } @@ -507,9 +507,9 @@ mod tests { let m1m = dvector![2.0, 0.0] * MOL; let m2m = dvector![0.0, 2.0] * MOL; let s1 = State::new_nvt(&e1, t, v, &m1).unwrap(); - let s2 = State::new_nvt(&e2, t, v, &m1).unwrap(); - let s1m = State::new_nvt(&e12, t, v, &m1m).unwrap(); - let s2m = State::new_nvt(&e12, t, v, &m2m).unwrap(); + let s2 = State::new_nvt(&e2, t, v, m1).unwrap(); + let s1m = State::new_nvt(&e12, t, v, m1m).unwrap(); + let s2m = State::new_nvt(&e12, t, v, m2m).unwrap(); assert_relative_eq!( s1.pressure(Contributions::Total), s1m.pressure(Contributions::Total), @@ -528,7 +528,7 @@ mod tests { let t = 300.0 * KELVIN; let p = BAR; let n = dvector![1.0] * MOL; - let s = State::new_npt(&e, t, p, &n, None)?; + let s = State::new_npt(&e, t, p, n, None)?; assert_relative_eq!( s.viscosity(), 0.00797 * MILLI * PASCAL * SECOND, @@ -536,7 +536,7 @@ mod tests { ); assert_relative_eq!( s.ln_viscosity_reduced(), - (s.viscosity() / e.viscosity_reference(s.temperature, s.volume, &s.moles)) + (s.viscosity() / e.viscosity_reference(s.temperature, s.molar_volume, &s.molefracs)) .into_value() .ln(), epsilon = 1e-15 @@ -553,15 +553,15 @@ mod tests { let t = 303.15 * KELVIN; let p = 500.0 * BAR; let n = dvector![0.25, 0.75] * MOL; - let viscosity_mix = State::new_npt(&e, t, p, &n, None)?.viscosity(); + let viscosity_mix = State::new_npt(&e, t, p, n, None)?.viscosity(); let viscosity_paper = 0.68298 * MILLI * PASCAL * SECOND; assert_relative_eq!(viscosity_paper, viscosity_mix, epsilon = 1e-8); // Make sure pure substance case is recovered let n_pseudo_mix = dvector![1.0, 0.0] * MOL; - let viscosity_pseudo_mix = State::new_npt(&e, t, p, &n_pseudo_mix, None)?.viscosity(); + let viscosity_pseudo_mix = State::new_npt(&e, t, p, n_pseudo_mix, None)?.viscosity(); let n_nonane = dvector![1.0] * MOL; - let viscosity_nonane = State::new_npt(&nonane, t, p, &n_nonane, None)?.viscosity(); + let viscosity_nonane = State::new_npt(&nonane, t, p, n_nonane, None)?.viscosity(); assert_relative_eq!(viscosity_pseudo_mix, viscosity_nonane, epsilon = 1e-15); Ok(()) } @@ -572,7 +572,7 @@ mod tests { let t = 300.0 * KELVIN; let p = BAR; let n = dvector![1.0] * MOL; - let s = State::new_npt(&e, t, p, &n, None)?; + let s = State::new_npt(&e, t, p, n, None)?; assert_relative_eq!( s.diffusion(), 0.01505 * (CENTI * METER).powi::<2>() / SECOND, @@ -580,7 +580,7 @@ mod tests { ); assert_relative_eq!( s.ln_diffusion_reduced(), - (s.diffusion() / e.diffusion_reference(s.temperature, s.volume, &s.moles)) + (s.diffusion() / e.diffusion_reference(s.temperature, s.molar_volume, &s.molefracs)) .into_value() .ln(), epsilon = 1e-15 @@ -685,6 +685,50 @@ mod tests_parameter_fit { Ok(()) } + #[test] + fn test_boiling_temperature_derivatives_fit() -> FeosResult<()> { + let pcsaft = pcsaft_non_assoc(); + let pcsaft_ad = pcsaft.named_derivatives(["m", "sigma", "epsilon_k"]); + let pressure = BAR; + let t = pcsaft_ad.boiling_temperature(pressure)?; + let t = t.convert_into(KELVIN); + let (t, grad) = (t.re, t.eps.unwrap_generic(U3, U1)); + + println!("{t:.5}"); + println!("{grad:.5?}"); + + let (t_check, _) = PhaseEquilibrium::pure_p( + &pcsaft_ad, + Pressure::from_inner(&pressure), + None, + Default::default(), + )?; + let t_check = t_check.convert_into(KELVIN); + let (t_check, grad_check) = (t_check.re, t_check.eps.unwrap_generic(U3, U1)); + println!("{t_check:.5}"); + println!("{grad_check:.5?}"); + assert_relative_eq!(t, t_check, max_relative = 1e-15); + assert_relative_eq!(grad, grad_check, max_relative = 1e-15); + + for (i, par) in ["m", "sigma", "epsilon_k"].into_iter().enumerate() { + let mut params = pcsaft.0; + let h = params[i] * 1e-8; + params[i] += h; + let pcsaft_h = PcSaftPure(params); + let (t_h, _) = PhaseEquilibrium::pure_p(&pcsaft_h, pressure, None, Default::default())?; + let dt_h = (t_h.convert_into(KELVIN) - t) / h; + let dt = grad[i]; + println!( + "{par:12}: {:11.5} {:11.5} {:.3e}", + dt_h, + dt, + ((dt_h - dt) / dt).abs() + ); + assert_relative_eq!(dt, dt_h, max_relative = 1e-6); + } + Ok(()) + } + #[test] fn test_equilibrium_liquid_density_derivatives_fit() -> FeosResult<()> { let pcsaft = pcsaft_non_assoc(); @@ -743,14 +787,9 @@ mod tests_parameter_fit { let h = params[i] * 1e-7; params[i] += h; let pcsaft_h = PcSaftPure(params); - let rho_h = State::new_xpt( - &pcsaft_h, - temperature, - pressure, - &vector![1.0], - Some(Liquid), - )? - .density; + let rho_h = + State::new_npt(&pcsaft_h, temperature, pressure, vector![1.0], Some(Liquid))? + .density; let drho_h = (rho_h.convert_into(MOL / LITER) - rho) / h; let drho = grad[i]; println!( @@ -954,10 +993,7 @@ mod tests_parameter_fit { ..Default::default() }, )?; - let beta = vle - .vapor() - .total_moles - .convert_into(vle.vapor().total_moles + vle.liquid().total_moles); + let beta = vle.1[0]; let (beta, [[grad]]) = (beta.re, beta.eps.unwrap_generic(U1, U1).data.0); println!("{beta:.5}"); @@ -977,10 +1013,7 @@ mod tests_parameter_fit { ..Default::default() }, )?; - let beta_h = vle - .vapor() - .total_moles - .convert_into(vle.vapor().total_moles + vle.liquid().total_moles); + let beta_h = vle.1[0]; let dbeta_h = (beta_h - beta) / h; println!( "k_ij: {:11.5} {:11.5} {:.3e}", diff --git a/crates/feos/src/pcsaft/eos/pcsaft_binary.rs b/crates/feos/src/pcsaft/eos/pcsaft_binary.rs index 3eecc9db4..71b15ae88 100644 --- a/crates/feos/src/pcsaft/eos/pcsaft_binary.rs +++ b/crates/feos/src/pcsaft/eos/pcsaft_binary.rs @@ -557,7 +557,7 @@ pub mod test { let h_feos = state.residual_molar_enthalpy(); let moles = vector![1.3, 2.5] * KILO * MOL; - let state = State::new_nvt(&pcsaft, temperature, volume, &moles)?; + let state = State::new_nvt(&pcsaft, temperature, volume, moles)?; let a_ad = state.residual_molar_helmholtz_energy(); let mu_ad = state.residual_chemical_potential(); let p_ad = state.pressure(Total); diff --git a/crates/feos/src/pcsaft/eos/pcsaft_pure.rs b/crates/feos/src/pcsaft/eos/pcsaft_pure.rs index baa5bc10f..f9ea2983f 100644 --- a/crates/feos/src/pcsaft/eos/pcsaft_pure.rs +++ b/crates/feos/src/pcsaft/eos/pcsaft_pure.rs @@ -228,9 +228,9 @@ impl ParametersAD<1> for PcSaftPure { #[cfg(test)] pub mod test { - use crate::pcsaft::PcSaftRecord; use super::super::{PcSaft, PcSaftAssociationRecord, PcSaftParameters}; use super::*; + use crate::pcsaft::PcSaftRecord; use approx::assert_relative_eq; use feos_core::parameter::{AssociationRecord, PureRecord}; use feos_core::{Contributions::Total, FeosResult, State}; @@ -278,7 +278,7 @@ pub mod test { let h_feos = state.residual_molar_enthalpy(); let moles = vector![1.3] * KILO * MOL; - let state = State::new_nvt(&pcsaft, temperature, volume, &moles)?; + let state = State::new_nvt(&pcsaft, temperature, volume, moles)?; let a_ad = state.residual_molar_helmholtz_energy(); let mu_ad = state.residual_chemical_potential(); let p_ad = state.pressure(Total); diff --git a/crates/feos/src/pets/eos/mod.rs b/crates/feos/src/pets/eos/mod.rs index f1d1d76cb..e26b94960 100644 --- a/crates/feos/src/pets/eos/mod.rs +++ b/crates/feos/src/pets/eos/mod.rs @@ -152,7 +152,7 @@ mod tests { let v = 1e-3 * METER.powi::<3>(); let n = dvector![1.0] * MOL; let s = State::new_nvt(&e, t, v, &n).unwrap(); - let p_ig = s.total_moles * RGAS * t / v; + let p_ig = s.total_moles() * RGAS * t / v; assert_relative_eq!(s.pressure(Contributions::IdealGas), p_ig, epsilon = 1e-10); assert_relative_eq!( s.pressure(Contributions::IdealGas) + s.pressure(Contributions::Residual), diff --git a/crates/feos/src/uvtheory/eos/mod.rs b/crates/feos/src/uvtheory/eos/mod.rs index 9fa963ae2..8a660db81 100644 --- a/crates/feos/src/uvtheory/eos/mod.rs +++ b/crates/feos/src/uvtheory/eos/mod.rs @@ -246,8 +246,9 @@ mod test { // EoS let eos_wca = &UVTheory::new(parameters); let state_wca = State::new_nvt(&eos_wca, t_x, volume, &moles).unwrap(); - let a_wca = (state_wca.residual_helmholtz_energy() / (RGAS * t_x * state_wca.total_moles)) - .into_value(); + let a_wca = (state_wca.residual_helmholtz_energy() + / (RGAS * t_x * state_wca.total_moles())) + .into_value(); assert_relative_eq!(a_wca, -0.597791038364405, max_relative = 1e-5); Ok(()) diff --git a/crates/feos/tests/gc_pcsaft/binary.rs b/crates/feos/tests/gc_pcsaft/binary.rs index e9b9ed434..ce627ec8c 100644 --- a/crates/feos/tests/gc_pcsaft/binary.rs +++ b/crates/feos/tests/gc_pcsaft/binary.rs @@ -30,9 +30,9 @@ fn test_binary() -> FeosResult<()> { #[cfg(feature = "dft")] let func = &GcPcSaftFunctional::new(parameters_func); let molefracs = dvector![0.5, 0.5]; - let cp = State::critical_point(&eos, Some(&molefracs), None, None, Default::default())?; + let cp = State::critical_point(&eos, &molefracs, None, None, Default::default())?; #[cfg(feature = "dft")] - let cp_func = State::critical_point(&func, Some(&molefracs), None, None, Default::default())?; + let cp_func = State::critical_point(&func, molefracs, None, None, Default::default())?; println!("{}", cp.temperature); #[cfg(feature = "dft")] println!("{}", cp_func.temperature); diff --git a/crates/feos/tests/gc_pcsaft/dft.rs b/crates/feos/tests/gc_pcsaft/dft.rs index a9813c56e..ce5d103cb 100644 --- a/crates/feos/tests/gc_pcsaft/dft.rs +++ b/crates/feos/tests/gc_pcsaft/dft.rs @@ -3,7 +3,7 @@ use approx::assert_relative_eq; use feos::gc_pcsaft::{GcPcSaft, GcPcSaftFunctional, GcPcSaftParameters}; use feos_core::parameter::{ChemicalRecord, Identifier, IdentifierOption, SegmentRecord}; -use feos_core::{PhaseEquilibrium, State, StateBuilder, Verbosity}; +use feos_core::{PhaseEquilibrium, State, Verbosity}; use feos_dft::adsorption::{ExternalPotential, Pore1D, PoreSpecification}; use feos_dft::interface::PlanarInterface; use feos_dft::{DFTSolver, Geometry}; @@ -158,7 +158,7 @@ fn test_dft() -> Result<(), Box> { let t = 200.0 * KELVIN; let w = 150.0 * ANGSTROM; let points = 2048; - let tc = State::critical_point(&&func, None, None, None, Default::default())?.temperature; + let tc = State::critical_point(&&func, (), None, None, Default::default())?.temperature; let vle = PhaseEquilibrium::pure(&&func, t, None, Default::default())?; let profile = PlanarInterface::from_tanh(&vle, points, w, tc, false).solve(None)?; println!( @@ -216,10 +216,7 @@ fn test_dft_assoc() -> Result<(), Box> { let solver = DFTSolver::new(Some(Verbosity::Iter)) .picard_iteration(None, None, Some(1e-5), Some(0.05)) .anderson_mixing(None, None, None, None, None); - let bulk = StateBuilder::new(&func) - .temperature(t) - .pressure(5.0 * BAR) - .build()?; + let bulk = State::new_npt(&func, t, 5.0 * BAR, (), None)?; Pore1D::new( Geometry::Cartesian, 20.0 * ANGSTROM, @@ -252,7 +249,7 @@ fn test_dft_newton() -> Result<(), Box> { let t = 200.0 * KELVIN; let w = 150.0 * ANGSTROM; let points = 512; - let tc = State::critical_point(&&func, None, None, None, Default::default())?.temperature; + let tc = State::critical_point(&&func, (), None, None, Default::default())?.temperature; let vle = PhaseEquilibrium::pure(&&func, t, None, Default::default())?; let solver = DFTSolver::new(Some(Verbosity::Iter)) .picard_iteration(None, Some(10), None, None) diff --git a/crates/feos/tests/pcsaft/critical_point.rs b/crates/feos/tests/pcsaft/critical_point.rs index 9183e5e5d..f0f435b13 100644 --- a/crates/feos/tests/pcsaft/critical_point.rs +++ b/crates/feos/tests/pcsaft/critical_point.rs @@ -17,7 +17,7 @@ fn test_critical_point_pure() -> Result<(), Box> { )?; let saft = PcSaft::new(params); let t = 300.0 * KELVIN; - let cp = State::critical_point(&&saft, None, Some(t), None, Default::default())?; + let cp = State::critical_point(&&saft, (), Some(t), None, Default::default())?; assert_relative_eq!(cp.temperature, 375.12441 * KELVIN, max_relative = 1e-8); assert_relative_eq!( cp.density, @@ -38,7 +38,7 @@ fn test_critical_point_mix() -> Result<(), Box> { let saft = PcSaft::new(params); let t = 300.0 * KELVIN; let molefracs = dvector![0.5, 0.5]; - let cp = State::critical_point(&&saft, Some(&molefracs), Some(t), None, Default::default())?; + let cp = State::critical_point(&&saft, molefracs, Some(t), None, Default::default())?; assert_relative_eq!(cp.temperature, 407.93481 * KELVIN, max_relative = 1e-8); assert_relative_eq!( cp.density, @@ -64,10 +64,10 @@ fn test_critical_point_limits() -> Result<(), Box> { let cp_pure = State::critical_point_pure(&saft, None, None, options)?; println!("{} {}", cp_pure[0], cp_pure[1]); let molefracs = dvector![0.0, 1.0]; - let cp_2 = State::critical_point(&saft, Some(&molefracs), None, None, options)?; + let cp_2 = State::critical_point(&saft, &molefracs, None, None, options)?; println!("{}", cp_2); let molefracs = dvector![1.0, 0.0]; - let cp_1 = State::critical_point(&saft, Some(&molefracs), None, None, options)?; + let cp_1 = State::critical_point(&saft, &molefracs, None, None, options)?; println!("{}", cp_1); assert_eq!(cp_pure[0].temperature, cp_1.temperature); assert_eq!(cp_pure[0].density, cp_1.density); diff --git a/crates/feos/tests/pcsaft/dft.rs b/crates/feos/tests/pcsaft/dft.rs index e8d057e47..b04a40aa1 100644 --- a/crates/feos/tests/pcsaft/dft.rs +++ b/crates/feos/tests/pcsaft/dft.rs @@ -103,7 +103,7 @@ fn test_dft_propane() -> Result<(), Box> { let t = 200.0 * KELVIN; let w = 150.0 * ANGSTROM; let points = 2048; - let tc = State::critical_point(&&func_pure, None, None, None, Default::default())?.temperature; + let tc = State::critical_point(&&func_pure, (), None, None, Default::default())?.temperature; let vle_pure = PhaseEquilibrium::pure(&&func_pure, t, None, Default::default())?; let vle_full = PhaseEquilibrium::pure(&&func_full, t, None, Default::default())?; let vle_full_vec = PhaseEquilibrium::pure(&&func_full_vec, t, None, Default::default())?; @@ -213,7 +213,7 @@ fn test_dft_propane_newton() -> Result<(), Box> { let t = 200.0 * KELVIN; let w = 150.0 * ANGSTROM; let points = 512; - let tc = State::critical_point(&&func, None, None, None, Default::default())?.temperature; + let tc = State::critical_point(&&func, (), None, None, Default::default())?.temperature; let vle = PhaseEquilibrium::pure(&&func, t, None, Default::default())?; let solver = DFTSolver::new(Some(Verbosity::Iter)).newton(None, None, None, None); PlanarInterface::from_tanh(&vle, points, w, tc, false).solve(Some(&solver))?; @@ -234,7 +234,7 @@ fn test_dft_water() -> Result<(), Box> { let t = 400.0 * KELVIN; let w = 120.0 * ANGSTROM; let points = 2048; - let tc = State::critical_point(&&func_pure, None, None, None, Default::default())?.temperature; + let tc = State::critical_point(&&func_pure, (), None, None, Default::default())?.temperature; let vle_pure = PhaseEquilibrium::pure(&&func_pure, t, None, Default::default())?; let vle_full_vec = PhaseEquilibrium::pure(&&func_full_vec, t, None, Default::default())?; let profile_pure = PlanarInterface::from_tanh(&vle_pure, points, w, tc, false).solve(None)?; @@ -334,33 +334,39 @@ fn test_entropy_bulk_values() -> Result<(), Box> { println!("\nResidual:\n{s_res:?}"); println!( "liquid: {:?}, vapor: {:?}", - profile.vle.liquid().entropy(Contributions::Residual) / profile.vle.liquid().volume, - profile.vle.vapor().entropy(Contributions::Residual) / profile.vle.vapor().volume + profile.vle.liquid().molar_entropy(Contributions::Residual) + / profile.vle.liquid().molar_volume, + profile.vle.vapor().molar_entropy(Contributions::Residual) + / profile.vle.vapor().molar_volume ); println!("\nTotal:\n{s_tot:?}"); println!( "liquid: {:?}, vapor: {:?}", - profile.vle.liquid().entropy(Contributions::Total) / profile.vle.liquid().volume, - profile.vle.vapor().entropy(Contributions::Total) / profile.vle.vapor().volume + profile.vle.liquid().molar_entropy(Contributions::Total) + / profile.vle.liquid().molar_volume, + profile.vle.vapor().molar_entropy(Contributions::Total) / profile.vle.vapor().molar_volume ); assert_relative_eq!( s_res.get(0), - profile.vle.liquid().entropy(Contributions::Residual) / profile.vle.liquid().volume, + profile.vle.liquid().molar_entropy(Contributions::Residual) + / profile.vle.liquid().molar_volume, max_relative = 1e-8, ); assert_relative_eq!( s_res.get(2047), - profile.vle.vapor().entropy(Contributions::Residual) / profile.vle.vapor().volume, + profile.vle.vapor().molar_entropy(Contributions::Residual) + / profile.vle.vapor().molar_volume, max_relative = 1e-8, ); assert_relative_eq!( s_tot.get(0), - profile.vle.liquid().entropy(Contributions::Total) / profile.vle.liquid().volume, + profile.vle.liquid().molar_entropy(Contributions::Total) + / profile.vle.liquid().molar_volume, max_relative = 1e-8, ); assert_relative_eq!( s_tot.get(2047), - profile.vle.vapor().entropy(Contributions::Total) / profile.vle.vapor().volume, + profile.vle.vapor().molar_entropy(Contributions::Total) / profile.vle.vapor().molar_volume, max_relative = 1e-8, ); Ok(()) diff --git a/crates/feos/tests/pcsaft/properties.rs b/crates/feos/tests/pcsaft/properties.rs index 3b05c86e1..489760e5e 100644 --- a/crates/feos/tests/pcsaft/properties.rs +++ b/crates/feos/tests/pcsaft/properties.rs @@ -1,7 +1,7 @@ use approx::assert_relative_eq; use feos::pcsaft::{PcSaft, PcSaftParameters}; use feos_core::parameter::IdentifierOption; -use feos_core::{Residual, StateBuilder}; +use feos_core::{DensityInitialization::Vapor, Residual, State}; use nalgebra::dvector; use quantity::*; use std::error::Error; @@ -18,18 +18,8 @@ fn test_dln_phi_dp() -> Result<(), Box> { let t = 300.0 * KELVIN; let p = BAR; let h = 1e-1 * PASCAL; - let s = StateBuilder::new(&&saft) - .temperature(t) - .pressure(p) - .molefracs(&dvector![0.5, 0.5]) - .vapor() - .build()?; - let sh = StateBuilder::new(&&saft) - .temperature(t) - .pressure(p + h) - .molefracs(&dvector![0.5, 0.5]) - .vapor() - .build()?; + let s = State::new_npt(&&saft, t, p, dvector![0.5, 0.5], Some(Vapor))?; + let sh = State::new_npt(&&saft, t, p + h, dvector![0.5, 0.5], Some(Vapor))?; let ln_phi = s.ln_phi()[0]; let ln_phi_h = sh.ln_phi()[0]; @@ -48,7 +38,7 @@ fn test_virial_is_not_nan() -> Result<(), Box> { IdentifierOption::Name, )?; let saft = &PcSaft::new(params); - let virial_b = saft.second_virial_coefficient(300.0 * KELVIN, &None); + let virial_b = saft.second_virial_coefficient(300.0 * KELVIN, ())?; assert!(!virial_b.is_nan()); Ok(()) } diff --git a/crates/feos/tests/pcsaft/state_creation_mixture.rs b/crates/feos/tests/pcsaft/state_creation_mixture.rs index b7b93c73a..46e0442b8 100644 --- a/crates/feos/tests/pcsaft/state_creation_mixture.rs +++ b/crates/feos/tests/pcsaft/state_creation_mixture.rs @@ -2,7 +2,7 @@ use approx::assert_relative_eq; use feos::ideal_gas::{Joback, JobackParameters}; use feos::pcsaft::{PcSaft, PcSaftParameters}; use feos_core::parameter::IdentifierOption; -use feos_core::{Contributions, EquationOfState, FeosResult, StateBuilder}; +use feos_core::{Contributions, EquationOfState, FeosResult, State}; use nalgebra::dvector; use quantity::*; use std::error::Error; @@ -31,17 +31,9 @@ fn pressure_entropy_molefracs() -> Result<(), Box> { let pressure = BAR; let temperature = 300.0 * KELVIN; let x = dvector![0.3, 0.7]; - let state = StateBuilder::new(&&eos) - .temperature(temperature) - .pressure(pressure) - .molefracs(&x) - .build()?; + let state = State::new_npt(&&eos, temperature, pressure, &x, None)?; let molar_entropy = state.molar_entropy(Contributions::Total); - let state = StateBuilder::new(&&eos) - .pressure(pressure) - .molar_entropy(molar_entropy) - .molefracs(&x) - .build()?; + let state = State::new_nps(&&eos, pressure, molar_entropy, x, None, None)?; assert_relative_eq!( state.molar_entropy(Contributions::Total), molar_entropy, @@ -63,13 +55,8 @@ fn volume_temperature_molefracs() -> Result<(), Box> { let volume = 1.5e-3 * METER.powi::<3>(); let moles = MOL; let x = dvector![0.3, 0.7]; - let state = StateBuilder::new(&&saft) - .temperature(temperature) - .volume(volume) - .total_moles(moles) - .molefracs(&x) - .build()?; - assert_relative_eq!(state.volume, volume, max_relative = 1e-10); + let state = State::new_nvt(&&saft, temperature, volume, (x, moles))?; + assert_relative_eq!(state.volume(), volume, max_relative = 1e-10); Ok(()) } @@ -80,15 +67,9 @@ fn temperature_partial_density() -> Result<(), Box> { let x = dvector![0.3, 0.7]; let partial_density = x.clone() * MOL / METER.powi::<3>(); let density = partial_density.sum(); - let state = StateBuilder::new(&&saft) - .temperature(temperature) - .partial_density(&partial_density) - .build()?; + let state = State::new_density(&&saft, temperature, partial_density)?; assert_relative_eq!(x, state.molefracs, max_relative = 1e-10); assert_relative_eq!(density, state.density, max_relative = 1e-10); - // Zip::from(&state.partial_density.to_reduced(reference)) - // .and(&partial_density.into_value()?) - // .for_each(|&r1, &r2| assert_relative_eq!(r1, r2, max_relative = 1e-10)); Ok(()) } @@ -98,11 +79,7 @@ fn temperature_density_molefracs() -> Result<(), Box> { let temperature = 300.0 * KELVIN; let x = dvector![0.3, 0.7]; let density = MOL / METER.powi::<3>(); - let state = StateBuilder::new(&&saft) - .temperature(temperature) - .density(density) - .molefracs(&x) - .build()?; + let state = State::new(&&saft, temperature, density, &x)?; state .molefracs .iter() @@ -118,11 +95,7 @@ fn temperature_pressure_molefracs() -> Result<(), Box> { let temperature = 300.0 * KELVIN; let pressure = BAR; let x = dvector![0.3, 0.7]; - let state = StateBuilder::new(&&saft) - .temperature(temperature) - .pressure(pressure) - .molefracs(&x) - .build()?; + let state = State::new_npt(&&saft, temperature, pressure, &x, None)?; state .molefracs .iter() diff --git a/crates/feos/tests/pcsaft/state_creation_pure.rs b/crates/feos/tests/pcsaft/state_creation_pure.rs index d5bc9769e..cda35b75f 100644 --- a/crates/feos/tests/pcsaft/state_creation_pure.rs +++ b/crates/feos/tests/pcsaft/state_creation_pure.rs @@ -1,12 +1,10 @@ use approx::assert_relative_eq; use feos::ideal_gas::{Joback, JobackParameters}; use feos::pcsaft::{PcSaft, PcSaftParameters}; +use feos_core::DensityInitialization::{InitialDensity, Liquid, Vapor}; use feos_core::parameter::IdentifierOption; -use feos_core::{ - Contributions, EquationOfState, FeosResult, PhaseEquilibrium, State, StateBuilder, Total, -}; +use feos_core::{Contributions, EquationOfState, FeosResult, PhaseEquilibrium, State, Total}; use quantity::*; -use std::error::Error; fn propane_parameters() -> FeosResult<(PcSaftParameters, Vec)> { let saft = PcSaftParameters::from_json( @@ -25,77 +23,59 @@ fn propane_parameters() -> FeosResult<(PcSaftParameters, Vec)> { } #[test] -fn temperature_volume() -> Result<(), Box> { +fn temperature_volume() -> FeosResult<()> { let saft = PcSaft::new(propane_parameters()?.0); let temperature = 300.0 * KELVIN; let volume = 1.5e-3 * METER.powi::<3>(); let moles = MOL; - let state = StateBuilder::new(&&saft) - .temperature(temperature) - .volume(volume) - .total_moles(moles) - .build()?; - assert_relative_eq!(state.volume, volume, max_relative = 1e-10); + let state = State::new_nvt(&&saft, temperature, volume, moles)?; + assert_relative_eq!(state.volume(), volume, max_relative = 1e-10); Ok(()) } #[test] -fn temperature_density() -> Result<(), Box> { +fn temperature_density() -> FeosResult<()> { let saft = PcSaft::new(propane_parameters()?.0); let temperature = 300.0 * KELVIN; let density = MOL / METER.powi::<3>(); - let state = StateBuilder::new(&&saft) - .temperature(temperature) - .density(density) - .build()?; + let state = State::new_pure(&&saft, temperature, density)?; assert_relative_eq!(state.density, density, max_relative = 1e-10); Ok(()) } #[test] -fn temperature_total_moles_volume() -> Result<(), Box> { +fn temperature_total_moles_volume() -> FeosResult<()> { let saft = PcSaft::new(propane_parameters()?.0); let temperature = 300.0 * KELVIN; let total_moles = MOL; let volume = METER.powi::<3>(); - let state = StateBuilder::new(&&saft) - .temperature(temperature) - .volume(volume) - .total_moles(total_moles) - .build()?; - assert_relative_eq!(state.volume, volume, max_relative = 1e-10); - assert_relative_eq!(state.total_moles, total_moles, max_relative = 1e-10); + let state = State::new_nvt(&&saft, temperature, volume, total_moles)?; + assert_relative_eq!(state.volume(), volume, max_relative = 1e-10); + assert_relative_eq!(state.total_moles(), total_moles, max_relative = 1e-10); Ok(()) } #[test] -fn temperature_total_moles_density() -> Result<(), Box> { +fn temperature_total_moles_density() -> FeosResult<()> { let saft = PcSaft::new(propane_parameters()?.0); let temperature = 300.0 * KELVIN; let total_moles = MOL; let density = MOL / METER.powi::<3>(); - let state = StateBuilder::new(&&saft) - .temperature(temperature) - .density(density) - .total_moles(total_moles) - .build()?; + let state = State::new_pure(&&saft, temperature, density)?.set_total_moles(total_moles); assert_relative_eq!(state.density, density, max_relative = 1e-10); - assert_relative_eq!(state.total_moles, total_moles, max_relative = 1e-10); - assert_relative_eq!(state.volume, total_moles / density, max_relative = 1e-10); + assert_relative_eq!(state.total_moles(), total_moles, max_relative = 1e-10); + assert_relative_eq!(state.volume(), total_moles / density, max_relative = 1e-10); Ok(()) } // Pressure constructors #[test] -fn pressure_temperature() -> Result<(), Box> { +fn pressure_temperature() -> FeosResult<()> { let saft = PcSaft::new(propane_parameters()?.0); let pressure = BAR; let temperature = 300.0 * KELVIN; - let state = StateBuilder::new(&&saft) - .temperature(temperature) - .pressure(pressure) - .build()?; + let state = State::new_npt(&&saft, temperature, pressure, (), None)?; assert_relative_eq!( state.pressure(Contributions::Total), pressure, @@ -105,15 +85,11 @@ fn pressure_temperature() -> Result<(), Box> { } #[test] -fn pressure_temperature_phase() -> Result<(), Box> { +fn pressure_temperature_phase() -> FeosResult<()> { let saft = PcSaft::new(propane_parameters()?.0); let pressure = BAR; let temperature = 300.0 * KELVIN; - let state = StateBuilder::new(&&saft) - .temperature(temperature) - .pressure(pressure) - .liquid() - .build()?; + let state = State::new_npt(&&saft, temperature, pressure, (), Some(Liquid))?; assert_relative_eq!( state.pressure(Contributions::Total), pressure, @@ -123,15 +99,12 @@ fn pressure_temperature_phase() -> Result<(), Box> { } #[test] -fn pressure_temperature_initial_density() -> Result<(), Box> { +fn pressure_temperature_initial_density() -> FeosResult<()> { let saft = PcSaft::new(propane_parameters()?.0); let pressure = BAR; let temperature = 300.0 * KELVIN; - let state = StateBuilder::new(&&saft) - .temperature(temperature) - .pressure(pressure) - .initial_density(MOL / METER.powi::<3>()) - .build()?; + let init = Some(InitialDensity(MOL / METER.powi::<3>())); + let state = State::new_npt(&&saft, temperature, pressure, (), init)?; assert_relative_eq!( state.pressure(Contributions::Total), pressure, @@ -141,17 +114,13 @@ fn pressure_temperature_initial_density() -> Result<(), Box> { } #[test] -fn pressure_enthalpy_vapor() -> Result<(), Box> { +fn pressure_enthalpy_vapor() -> FeosResult<()> { let (saft_params, joback) = propane_parameters()?; let saft = PcSaft::new(saft_params); let eos = EquationOfState::new(joback, saft); let pressure = 0.3 * BAR; let molar_enthalpy = 2000.0 * JOULE / MOL; - let state = StateBuilder::new(&&eos) - .pressure(pressure) - .molar_enthalpy(molar_enthalpy) - .vapor() - .build()?; + let state = State::new_nph(&&eos, pressure, molar_enthalpy, (), Some(Vapor), None)?; assert_relative_eq!( state.molar_enthalpy(Contributions::Total), molar_enthalpy, @@ -163,11 +132,7 @@ fn pressure_enthalpy_vapor() -> Result<(), Box> { max_relative = 1e-10 ); - let state = StateBuilder::new(&&eos) - .volume(state.volume) - .temperature(state.temperature) - .moles(&state.moles) - .build()?; + let state = State::new(&&eos, state.temperature, state.density, state.molefracs)?; assert_relative_eq!( state.molar_enthalpy(Contributions::Total), molar_enthalpy, @@ -182,24 +147,22 @@ fn pressure_enthalpy_vapor() -> Result<(), Box> { } #[test] -fn density_internal_energy() -> Result<(), Box> { +fn density_internal_energy() -> FeosResult<()> { let (saft_params, joback) = propane_parameters()?; let saft = PcSaft::new(saft_params); let eos = EquationOfState::new(joback, saft); let pressure = 5.0 * BAR; let temperature = 315.0 * KELVIN; let total_moles = 2.5 * MOL; - let state = StateBuilder::new(&&eos) - .pressure(pressure) - .temperature(temperature) - .total_moles(total_moles) - .build()?; + let state = State::new_npt(&&eos, temperature, pressure, total_moles, None)?; let molar_internal_energy = state.molar_internal_energy(Contributions::Total); - let state_nvu = StateBuilder::new(&&eos) - .volume(state.volume) - .molar_internal_energy(molar_internal_energy) - .total_moles(total_moles) - .build()?; + let state_nvu = State::new_nvu( + &&eos, + state.volume(), + molar_internal_energy, + total_moles, + None, + )?; assert_relative_eq!( molar_internal_energy, state_nvu.molar_internal_energy(Contributions::Total), @@ -211,19 +174,21 @@ fn density_internal_energy() -> Result<(), Box> { } #[test] -fn pressure_enthalpy_total_moles_vapor() -> Result<(), Box> { +fn pressure_enthalpy_total_moles_vapor() -> FeosResult<()> { let (saft_params, joback) = propane_parameters()?; let saft = PcSaft::new(saft_params); let eos = EquationOfState::new(joback, saft); let pressure = 0.3 * BAR; let molar_enthalpy = 2000.0 * JOULE / MOL; let total_moles = 2.5 * MOL; - let state = StateBuilder::new(&&eos) - .pressure(pressure) - .molar_enthalpy(molar_enthalpy) - .total_moles(total_moles) - .vapor() - .build()?; + let state = State::new_nph( + &&eos, + pressure, + molar_enthalpy, + total_moles, + Some(Vapor), + None, + )?; assert_relative_eq!( state.molar_enthalpy(Contributions::Total), molar_enthalpy, @@ -235,11 +200,12 @@ fn pressure_enthalpy_total_moles_vapor() -> Result<(), Box> { max_relative = 1e-10 ); - let state = StateBuilder::new(&&eos) - .volume(state.volume) - .temperature(state.temperature) - .total_moles(state.total_moles) - .build()?; + let state = State::new_nvt( + &&eos, + state.temperature, + state.volume(), + state.total_moles(), + )?; assert_relative_eq!( state.molar_enthalpy(Contributions::Total), molar_enthalpy, @@ -254,17 +220,13 @@ fn pressure_enthalpy_total_moles_vapor() -> Result<(), Box> { } #[test] -fn pressure_entropy_vapor() -> Result<(), Box> { +fn pressure_entropy_vapor() -> FeosResult<()> { let (saft_params, joback) = propane_parameters()?; let saft = PcSaft::new(saft_params); let eos = EquationOfState::new(joback, saft); let pressure = 0.3 * BAR; let molar_entropy = -2.0 * JOULE / MOL / KELVIN; - let state = StateBuilder::new(&&eos) - .pressure(pressure) - .molar_entropy(molar_entropy) - .vapor() - .build()?; + let state = State::new_nps(&&eos, pressure, molar_entropy, (), Some(Vapor), None)?; assert_relative_eq!( state.molar_entropy(Contributions::Total), molar_entropy, @@ -276,11 +238,7 @@ fn pressure_entropy_vapor() -> Result<(), Box> { max_relative = 1e-10 ); - let state = StateBuilder::new(&&eos) - .volume(state.volume) - .temperature(state.temperature) - .moles(&state.moles) - .build()?; + let state = State::new(&&eos, state.temperature, state.density, state.molefracs)?; assert_relative_eq!( state.molar_entropy(Contributions::Total), molar_entropy, @@ -295,24 +253,20 @@ fn pressure_entropy_vapor() -> Result<(), Box> { } #[test] -fn temperature_entropy_vapor() -> Result<(), Box> { +fn temperature_entropy_vapor() -> FeosResult<()> { let (saft_params, joback) = propane_parameters()?; let saft = PcSaft::new(saft_params); let eos = EquationOfState::new(joback, saft); let pressure = 3.0 * BAR; let temperature = 315.15 * KELVIN; let total_moles = 3.0 * MOL; - let state = StateBuilder::new(&&eos) - .temperature(temperature) - .pressure(pressure) - .total_moles(total_moles) - .build()?; + let state = State::new_npt(&&eos, temperature, pressure, total_moles, None)?; let s = State::new_nts( &&eos, temperature, state.molar_entropy(Contributions::Total), - &state.moles, + state.moles(), None, )?; assert_relative_eq!( @@ -354,7 +308,7 @@ fn assert_multiple_states( } #[test] -fn test_consistency() -> Result<(), Box> { +fn test_consistency() -> FeosResult<()> { let (saft_params, joback) = propane_parameters()?; let saft = PcSaft::new(saft_params); let eos = EquationOfState::new(joback, saft); @@ -362,10 +316,7 @@ fn test_consistency() -> Result<(), Box> { let pressures = [1.0 * BAR, 2.0 * BAR, 3.0 * BAR]; for (&temperature, &pressure) in temperatures.iter().zip(pressures.iter()) { - let state = StateBuilder::new(&&eos) - .pressure(pressure) - .temperature(temperature) - .build()?; + let state = State::new_npt(&&eos, temperature, pressure, (), None)?; assert_relative_eq!( state.pressure(Contributions::Total), pressure, @@ -379,49 +330,26 @@ fn test_consistency() -> Result<(), Box> { let molar_entropy = state.molar_entropy(Contributions::Total); let density = state.density; - let state_tv = StateBuilder::new(&&eos) - .temperature(temperature) - .density(density) - .build()?; + let state_tv = State::new_pure(&&eos, temperature, density)?; let vle = PhaseEquilibrium::pure(&&eos, temperature, None, Default::default()); let eos = &eos; - let builder = if let Ok(ps) = vle { + let phase = if let Ok(ps) = vle { let p_sat = ps.liquid().pressure(Contributions::Total); - if pressure > p_sat { - StateBuilder::new(&eos).liquid() - } else { - StateBuilder::new(&eos).vapor() - } + if pressure > p_sat { Liquid } else { Vapor } } else { - StateBuilder::new(&eos).vapor() + Vapor }; - let state_ts = builder - .clone() - .temperature(temperature) - .molar_entropy(molar_entropy) - .build()?; + let state_ts = State::new_nts(&eos, temperature, molar_entropy, (), Some(phase))?; - let state_ps = builder - .clone() - .pressure(pressure) - .molar_entropy(molar_entropy) - .build()?; + let state_ps = State::new_nps(&eos, pressure, molar_entropy, (), Some(phase), None)?; dbg!("ph"); - let state_ph = builder - .clone() - .pressure(pressure) - .molar_enthalpy(molar_enthalpy) - .build()?; + let state_ph = State::new_nph(&eos, pressure, molar_enthalpy, (), Some(phase), None)?; dbg!("th"); - let state_th = builder - .clone() - .temperature(temperature) - .molar_enthalpy(molar_enthalpy) - .build()?; + let state_th = State::new_nth(&eos, temperature, molar_enthalpy, (), Some(phase))?; dbg!("assertions"); assert_multiple_states( diff --git a/crates/feos/tests/saftvrmie/critical_properties.rs b/crates/feos/tests/saftvrmie/critical_properties.rs index c78eb110a..1f3334ea3 100644 --- a/crates/feos/tests/saftvrmie/critical_properties.rs +++ b/crates/feos/tests/saftvrmie/critical_properties.rs @@ -52,7 +52,7 @@ fn critical_properties_pure() { let option = SolverOptions::default(); let p = parameters.remove(name).unwrap(); let eos = SaftVRMie::new(p); - let cp = State::critical_point(&&eos, None, t0, None, option).unwrap(); + let cp = State::critical_point(&&eos, (), t0, None, option).unwrap(); assert_relative_eq!(cp.temperature, data.0, max_relative = 2e-3); assert_relative_eq!( cp.pressure(feos_core::Contributions::Total), diff --git a/py-feos/src/ad/mod.rs b/py-feos/src/ad/mod.rs index 1286f6c14..75958d7ca 100644 --- a/py-feos/src/ad/mod.rs +++ b/py-feos/src/ad/mod.rs @@ -50,13 +50,39 @@ type GradResult<'py> = ( #[pyfunction] pub fn vapor_pressure_derivatives<'py>( model: PyEquationOfStateAD, - parameter_names: Bound<'py, PyAny>, + parameter_names: &Bound<'py, PyAny>, parameters: PyReadonlyArray2, input: PyReadonlyArray2, ) -> GradResult<'py> { _vapor_pressure_derivatives(model, parameter_names, parameters, input) } +/// Calculate boiling temperatures and derivatives w.r.t. model parameters. +/// +/// Parameters +/// ---------- +/// model: EquationOfStateAD +/// The equation of state to use. +/// parameter_names: List[string] +/// The name of the parameters for which derivatives are calculated. +/// parameters: np.ndarray[float] +/// The parameters for every data point. +/// input: np.ndarray[float] +/// The pressure (in Pa) for every data point. +/// +/// Returns +/// ------- +/// (np.ndarray[float], np.ndarray[float], np.ndarray[bool]): The boiling temperature (in K), gradients, and convergence status. +#[pyfunction] +pub fn boiling_temperature_derivatives<'py>( + model: PyEquationOfStateAD, + parameter_names: Bound<'py, PyAny>, + parameters: PyReadonlyArray2, + input: PyReadonlyArray2, +) -> GradResult<'py> { + _boiling_temperature_derivatives(model, parameter_names, parameters, input) +} + /// Calculate liquid densities and derivatives w.r.t. model parameters. /// /// Parameters @@ -76,7 +102,7 @@ pub fn vapor_pressure_derivatives<'py>( #[pyfunction] pub fn liquid_density_derivatives<'py>( model: PyEquationOfStateAD, - parameter_names: Bound<'py, PyAny>, + parameter_names: &Bound<'py, PyAny>, parameters: PyReadonlyArray2, input: PyReadonlyArray2, ) -> GradResult<'py> { @@ -102,7 +128,7 @@ pub fn liquid_density_derivatives<'py>( #[pyfunction] pub fn equilibrium_liquid_density_derivatives<'py>( model: PyEquationOfStateAD, - parameter_names: Bound<'py, PyAny>, + parameter_names: &Bound<'py, PyAny>, parameters: PyReadonlyArray2, input: PyReadonlyArray2, ) -> GradResult<'py> { @@ -129,7 +155,7 @@ pub fn equilibrium_liquid_density_derivatives<'py>( #[pyfunction] pub fn bubble_point_pressure_derivatives<'py>( model: PyEquationOfStateAD, - parameter_names: Bound<'py, PyAny>, + parameter_names: &Bound<'py, PyAny>, parameters: PyReadonlyArray2, input: PyReadonlyArray2, ) -> GradResult<'py> { @@ -156,7 +182,7 @@ pub fn bubble_point_pressure_derivatives<'py>( #[pyfunction] pub fn dew_point_pressure_derivatives<'py>( model: PyEquationOfStateAD, - parameter_names: Bound<'py, PyAny>, + parameter_names: &Bound<'py, PyAny>, parameters: PyReadonlyArray2, input: PyReadonlyArray2, ) -> GradResult<'py> { @@ -169,7 +195,7 @@ macro_rules! expand_models { #[pyfunction] fn [<_ $prop _derivatives>]<'py>( model: PyEquationOfStateAD, - parameter_names: Bound<'py, PyAny>, + parameter_names: &Bound<'py, PyAny>, parameters: PyReadonlyArray2, input: PyReadonlyArray2, ) -> GradResult<'py> { @@ -194,7 +220,7 @@ macro_rules! impl_evaluate_gradients { expand_models!($enum, $prop, $($model: $type),*); paste!( fn $prop<'py, R: ParametersAD<$n>>( - parameter_names: Bound<'py, PyAny>, + parameter_names: &Bound<'py, PyAny>, parameters: PyReadonlyArray2, input: PyReadonlyArray2, ) -> ( @@ -222,7 +248,7 @@ macro_rules! impl_evaluate_gradients { impl_evaluate_gradients!( pure, - [vapor_pressure, liquid_density, equilibrium_liquid_density], + [vapor_pressure, boiling_temperature, liquid_density, equilibrium_liquid_density], {PcSaftNonAssoc: PcSaftPure, PcSaftFull: PcSaftPure} ); diff --git a/py-feos/src/dft/adsorption/mod.rs b/py-feos/src/dft/adsorption/mod.rs index 1f006be0a..1c5040003 100644 --- a/py-feos/src/dft/adsorption/mod.rs +++ b/py-feos/src/dft/adsorption/mod.rs @@ -1,9 +1,9 @@ use super::PyDFTSolver; -use crate::eos::{parse_molefracs, PyEquationOfState}; +use crate::PyVerbosity; +use crate::eos::{Compositions, PyEquationOfState}; use crate::error::PyFeosError; use crate::ideal_gas::IdealGasModel; use crate::residual::ResidualModel; -use crate::PyVerbosity; use feos_core::EquationOfState; use feos_dft::adsorption::{Adsorption, Adsorption1D, Adsorption3D}; use nalgebra::DMatrix; @@ -45,8 +45,8 @@ macro_rules! impl_adsorption_isotherm { /// The pressures for which the profiles are calculated. /// pore : Pore /// The pore parameters. - /// molefracs: numpy.ndarray[float], optional - /// For a mixture, the molefracs of the bulk system. + /// composition : float | SINumber | numpy.ndarray[float] | SIArray1 | list[float], optional + /// The composition of the mixture. /// solver: DFTSolver, optional /// Custom solver options. /// @@ -55,14 +55,14 @@ macro_rules! impl_adsorption_isotherm { /// Adsorption /// #[staticmethod] - #[pyo3(text_signature = "(functional, temperature, pressure, pore, molefracs=None, solver=None)")] - #[pyo3(signature = (functional, temperature, pressure, pore, molefracs=None, solver=None))] + #[pyo3(text_signature = "(functional, temperature, pressure, pore, composition=None, solver=None)")] + #[pyo3(signature = (functional, temperature, pressure, pore, composition=None, solver=None))] fn adsorption_isotherm( functional: &PyEquationOfState, temperature: Temperature, pressure: Pressure>, pore: &$py_pore, - molefracs: Option>, + composition: Option<&Bound<'_, PyAny>>, solver: Option, ) -> PyResult { Ok(Self(Adsorption::adsorption_isotherm( @@ -70,7 +70,7 @@ macro_rules! impl_adsorption_isotherm { temperature, &pressure, &pore.0, - &parse_molefracs(molefracs), + Compositions::try_from(composition)?, solver.map(|s| s.0).as_ref(), ).map_err(PyFeosError::from)?)) } @@ -89,8 +89,8 @@ macro_rules! impl_adsorption_isotherm { /// The pressures for which the profiles are calculated. /// pore : Pore /// The pore parameters. - /// molefracs: numpy.ndarray[float], optional - /// For a mixture, the molefracs of the bulk system. + /// composition : float | SINumber | numpy.ndarray[float] | SIArray1 | list[float], optional + /// The composition of the mixture. /// solver: DFTSolver, optional /// Custom solver options. /// @@ -99,14 +99,14 @@ macro_rules! impl_adsorption_isotherm { /// Adsorption /// #[staticmethod] - #[pyo3(text_signature = "(functional, temperature, pressure, pore, molefracs=None, solver=None)")] - #[pyo3(signature = (functional, temperature, pressure, pore, molefracs=None, solver=None))] + #[pyo3(text_signature = "(functional, temperature, pressure, pore, composition=None, solver=None)")] + #[pyo3(signature = (functional, temperature, pressure, pore, composition=None, solver=None))] fn desorption_isotherm( functional: &PyEquationOfState, temperature: Temperature, pressure: Pressure>, pore: &$py_pore, - molefracs: Option>, + composition: Option<&Bound<'_, PyAny>>, solver: Option, ) -> PyResult { Ok(Self(Adsorption::desorption_isotherm( @@ -114,7 +114,7 @@ macro_rules! impl_adsorption_isotherm { temperature, &pressure, &pore.0, - &parse_molefracs(molefracs), + Compositions::try_from(composition)?, solver.map(|s| s.0).as_ref(), ).map_err(PyFeosError::from)?)) } @@ -136,8 +136,8 @@ macro_rules! impl_adsorption_isotherm { /// The pressures for which the profiles are calculated. /// pore : Pore /// The pore parameters. - /// molefracs: numpy.ndarray[float], optional - /// For a mixture, the molefracs of the bulk system. + /// composition : float | SINumber | numpy.ndarray[float] | SIArray1 | list[float], optional + /// The composition of the mixture. /// solver: DFTSolver, optional /// Custom solver options. /// @@ -146,14 +146,14 @@ macro_rules! impl_adsorption_isotherm { /// Adsorption /// #[staticmethod] - #[pyo3(text_signature = "(functional, temperature, pressure, pore, molefracs=None, solver=None)")] - #[pyo3(signature = (functional, temperature, pressure, pore, molefracs=None, solver=None))] + #[pyo3(text_signature = "(functional, temperature, pressure, pore, composition=None, solver=None)")] + #[pyo3(signature = (functional, temperature, pressure, pore, composition=None, solver=None))] fn equilibrium_isotherm( functional: &PyEquationOfState, temperature: Temperature, pressure: Pressure>, pore: &$py_pore, - molefracs: Option>, + composition: Option<&Bound<'_, PyAny>>, solver: Option, ) -> PyResult { Ok(Self(Adsorption::equilibrium_isotherm( @@ -161,7 +161,7 @@ macro_rules! impl_adsorption_isotherm { temperature, &pressure, &pore.0, - &parse_molefracs(molefracs), + Compositions::try_from(composition)?, solver.map(|s| s.0).as_ref(), ).map_err(PyFeosError::from)?)) } @@ -180,8 +180,8 @@ macro_rules! impl_adsorption_isotherm { /// A suitable upper limit for the pressure. /// pore : Pore /// The pore parameters. - /// molefracs: numpy.ndarray[float], optional - /// For a mixture, the molefracs of the bulk system. + /// composition : float | SINumber | numpy.ndarray[float] | SIArray1 | list[float], optional + /// The composition of the mixture. /// solver: DFTSolver, optional /// Custom solver options. /// max_iter : int, optional @@ -196,8 +196,8 @@ macro_rules! impl_adsorption_isotherm { /// Adsorption /// #[staticmethod] - #[pyo3(text_signature = "(functional, temperature, p_min, p_max, pore, molefracs=None, solver=None, max_iter=None, tol=None, verbosity=None)")] - #[pyo3(signature = (functional, temperature, p_min, p_max, pore, molefracs=None, solver=None, max_iter=None, tol=None, verbosity=None))] + #[pyo3(text_signature = "(functional, temperature, p_min, p_max, pore, composition=None, solver=None, max_iter=None, tol=None, verbosity=None)")] + #[pyo3(signature = (functional, temperature, p_min, p_max, pore, composition=None, solver=None, max_iter=None, tol=None, verbosity=None))] #[expect(clippy::too_many_arguments)] fn phase_equilibrium( functional: &PyEquationOfState, @@ -205,7 +205,7 @@ macro_rules! impl_adsorption_isotherm { p_min: Pressure, p_max: Pressure, pore: &$py_pore, - molefracs: Option>, + composition: Option<&Bound<'_, PyAny>>, solver: Option, max_iter: Option, tol: Option, @@ -217,7 +217,7 @@ macro_rules! impl_adsorption_isotherm { p_min, p_max, &pore.0, - &parse_molefracs(molefracs), + Compositions::try_from(composition)?, solver.map(|s| s.0).as_ref(), (max_iter, tol, verbosity.map(|v| v.into())).into(), ).map_err(PyFeosError::from)?)) diff --git a/py-feos/src/eos/mod.rs b/py-feos/src/eos/mod.rs index 1147ac020..cd50de42f 100644 --- a/py-feos/src/eos/mod.rs +++ b/py-feos/src/eos/mod.rs @@ -3,9 +3,9 @@ use crate::ideal_gas::IdealGasModel; use crate::residual::ResidualModel; use feos_core::*; use indexmap::IndexMap; -use nalgebra::{DVector, DVectorView, Dyn}; +use nalgebra::{DVector, DVectorView, Dyn, U1}; use numpy::{PyArray1, PyReadonlyArray1, ToPyArray}; -use pyo3::prelude::*; +use pyo3::{exceptions::PyValueError, prelude::*}; use quantity::*; use std::ops::Div; use std::sync::Arc; @@ -58,17 +58,17 @@ impl PyEquationOfState { /// /// Parameters /// ---------- - /// molefracs : np.ndarray[float], optional + /// composition : float | SINumber | numpy.ndarray[float] | SIArray1 | list[float], optional /// The composition of the mixture. /// /// Returns /// ------- /// SINumber - #[pyo3(text_signature = "(molefracs=None)", signature = (molefracs=None))] - fn max_density(&self, molefracs: Option>) -> PyResult { + #[pyo3(text_signature = "(composition=None)", signature = (composition=None))] + fn max_density(&self, composition: Option<&Bound<'_, PyAny>>) -> PyResult { Ok(self .0 - .max_density(&parse_molefracs(molefracs)) + .max_density(Compositions::try_from(composition)?) .map_err(PyFeosError::from)?) } @@ -78,20 +78,22 @@ impl PyEquationOfState { /// ---------- /// temperature : SINumber /// The temperature for which B should be computed. - /// molefracs : np.ndarray[float], optional + /// composition : float | SINumber | numpy.ndarray[float] | SIArray1 | list[float], optional /// The composition of the mixture. /// /// Returns /// ------- /// SINumber - #[pyo3(text_signature = "(temperature, molefracs=None)", signature = (temperature, molefracs=None))] + #[pyo3(text_signature = "(temperature, composition=None)", signature = (temperature, composition=None))] fn second_virial_coefficient( &self, temperature: Temperature, - molefracs: Option>, - ) -> Quot { - self.0 - .second_virial_coefficient(temperature, &parse_molefracs(molefracs)) + composition: Option<&Bound<'_, PyAny>>, + ) -> PyResult> { + Ok(self + .0 + .second_virial_coefficient(temperature, Compositions::try_from(composition)?) + .map_err(PyFeosError::from)?) } /// Calculate the third Virial coefficient C(T,x). @@ -100,20 +102,22 @@ impl PyEquationOfState { /// ---------- /// temperature : SINumber /// The temperature for which C should be computed. - /// molefracs : np.ndarray[float], optional + /// composition : float | SINumber | numpy.ndarray[float] | SIArray1 | list[float], optional /// The composition of the mixture. /// /// Returns /// ------- /// SINumber - #[pyo3(text_signature = "(temperature, molefracs=None)", signature = (temperature, molefracs=None))] + #[pyo3(text_signature = "(temperature, composition=None)", signature = (temperature, composition=None))] fn third_virial_coefficient( &self, temperature: Temperature, - molefracs: Option>, - ) -> Quot, Density> { - self.0 - .third_virial_coefficient(temperature, &parse_molefracs(molefracs)) + composition: Option<&Bound<'_, PyAny>>, + ) -> PyResult, Density>> { + Ok(self + .0 + .third_virial_coefficient(temperature, Compositions::try_from(composition)?) + .map_err(PyFeosError::from)?) } /// Calculate the derivative of the second Virial coefficient B(T,x) @@ -123,22 +127,25 @@ impl PyEquationOfState { /// ---------- /// temperature : SINumber /// The temperature for which B' should be computed. - /// molefracs : np.ndarray[float], optional + /// composition : float | SINumber | numpy.ndarray[float] | SIArray1 | list[float], optional /// The composition of the mixture. /// /// Returns /// ------- /// SINumber - #[pyo3(text_signature = "(temperature, molefracs=None)", signature = (temperature, molefracs=None))] + #[pyo3(text_signature = "(temperature, composition=None)", signature = (temperature, composition=None))] fn second_virial_coefficient_temperature_derivative( &self, temperature: Temperature, - molefracs: Option>, - ) -> Quot, Temperature> { - self.0.second_virial_coefficient_temperature_derivative( - temperature, - &parse_molefracs(molefracs), - ) + composition: Option<&Bound<'_, PyAny>>, + ) -> PyResult, Temperature>> { + Ok(self + .0 + .second_virial_coefficient_temperature_derivative( + temperature, + Compositions::try_from(composition)?, + ) + .map_err(PyFeosError::from)?) } /// Calculate the derivative of the third Virial coefficient C(T,x) @@ -148,22 +155,26 @@ impl PyEquationOfState { /// ---------- /// temperature : SINumber /// The temperature for which C' should be computed. - /// molefracs : np.ndarray[float], optional + /// composition : float | SINumber | numpy.ndarray[float] | SIArray1 | list[float], optional /// The composition of the mixture. /// /// Returns /// ------- /// SINumber - #[pyo3(text_signature = "(temperature, molefracs=None)", signature = (temperature, molefracs=None))] + #[expect(clippy::type_complexity)] + #[pyo3(text_signature = "(temperature, composition=None)", signature = (temperature, composition=None))] fn third_virial_coefficient_temperature_derivative( &self, temperature: Temperature, - molefracs: Option>, - ) -> Quot, Density>, Temperature> { - self.0.third_virial_coefficient_temperature_derivative( - temperature, - &parse_molefracs(molefracs), - ) + composition: Option<&Bound<'_, PyAny>>, + ) -> PyResult, Density>, Temperature>> { + Ok(self + .0 + .third_virial_coefficient_temperature_derivative( + temperature, + Compositions::try_from(composition)?, + ) + .map_err(PyFeosError::from)?) } } @@ -192,36 +203,64 @@ pub(crate) fn parse_molefracs(molefracs: Option>) -> Optio }) } -// impl_state_entropy_scaling!(EquationOfState, ResidualModel>, PyEquationOfState); -// impl_phase_equilibrium!(EquationOfState, ResidualModel>, PyEquationOfState); - -// #[cfg(feature = "estimator")] -// impl_estimator!(EquationOfState, ResidualModel>, PyEquationOfState); -// #[cfg(all(feature = "estimator", feature = "pcsaft"))] -// impl_estimator_entropy_scaling!(EquationOfState, ResidualModel>, PyEquationOfState); - -// #[pymodule] -// pub fn eos(m: &Bound<'_, PyModule>) -> PyResult<()> { -// m.add_class::()?; -// m.add_class::()?; - -// m.add_class::()?; -// m.add_class::()?; -// m.add_class::()?; -// m.add_class::()?; -// m.add_class::()?; +#[derive(Clone)] +pub enum Compositions { + None, + Scalar(f64), + TotalMoles(Moles), + Molefracs(DVector), + Moles(Moles>), + PartialDensity(Density>), +} -// #[cfg(feature = "estimator")] -// m.add_wrapped(wrap_pymodule!(estimator_eos))?; +impl Composition for Compositions { + fn into_molefracs>( + self, + eos: &E, + ) -> FeosResult<(DVector, Option>)> { + match self { + Self::None => ().into_molefracs(eos), + Self::Scalar(x) => x.into_molefracs(eos), + Self::TotalMoles(total_moles) => total_moles.into_molefracs(eos), + Self::Molefracs(molefracs) => molefracs.into_molefracs(eos), + Self::Moles(moles) => moles.into_molefracs(eos), + Self::PartialDensity(partial_density) => partial_density.into_molefracs(eos), + } + } -// Ok(()) -// } + fn density(&self) -> Option> { + if let Self::PartialDensity(partial_density) = self { + partial_density.density() + } else { + None + } + } +} -// #[cfg(feature = "estimator")] -// #[pymodule] -// pub fn estimator_eos(m: &Bound<'_, PyModule>) -> PyResult<()> { -// m.add_class::()?; -// m.add_class::()?; -// m.add_class::()?; -// m.add_class::() -// } +impl TryFrom>> for Compositions { + type Error = PyErr; + fn try_from(composition: Option<&Bound<'_, PyAny>>) -> PyResult { + let Some(composition) = composition else { + return Ok(Compositions::None); + }; + if let Ok(x) = composition.extract::>() + && let Some(x) = x.try_as_matrix::() + { + Ok(Compositions::Molefracs(x.clone_owned())) + } else if let Ok(x) = composition.extract::>() { + Ok(Compositions::Molefracs(DVector::from_vec(x))) + } else if let Ok(x) = composition.extract::() { + Ok(Compositions::Scalar(x)) + } else if let Ok(n) = composition.extract::>>() { + Ok(Compositions::Moles(n)) + } else if let Ok(n) = composition.extract::() { + Ok(Compositions::TotalMoles(n)) + } else if let Ok(rho) = composition.extract::>>() { + Ok(Compositions::PartialDensity(rho)) + } else { + Err(PyErr::new::(format!( + "failed to parse value '{composition}' as composition." + ))) + } + } +} diff --git a/py-feos/src/lib.rs b/py-feos/src/lib.rs index b4f94441a..7ae4d269c 100644 --- a/py-feos/src/lib.rs +++ b/py-feos/src/lib.rs @@ -91,6 +91,7 @@ fn feos(m: &Bound<'_, PyModule>) -> PyResult<()> { #[cfg(feature = "ad")] { m.add_function(wrap_pyfunction!(ad::vapor_pressure_derivatives, m)?)?; + m.add_function(wrap_pyfunction!(ad::boiling_temperature_derivatives, m)?)?; m.add_function(wrap_pyfunction!(ad::liquid_density_derivatives, m)?)?; m.add_function(wrap_pyfunction!( ad::equilibrium_liquid_density_derivatives, diff --git a/py-feos/src/parameter/mod.rs b/py-feos/src/parameter/mod.rs index 115f8f79a..aabbed5ab 100644 --- a/py-feos/src/parameter/mod.rs +++ b/py-feos/src/parameter/mod.rs @@ -557,6 +557,7 @@ impl PyGcParameters { .map_err(PyFeosError::from)?) } + #[cfg(feature = "gc_pcsaft")] pub fn try_convert_heterosegmented( self, ) -> PyResult> diff --git a/py-feos/src/phase_equilibria.rs b/py-feos/src/phase_equilibria.rs index c3786e300..38d0d2d12 100644 --- a/py-feos/src/phase_equilibria.rs +++ b/py-feos/src/phase_equilibria.rs @@ -61,7 +61,7 @@ impl PyPhaseEquilibrium { #[pyo3(signature = (eos, temperature_or_pressure, initial_state=None, max_iter=None, tol=None, verbosity=None))] pub(crate) fn pure( eos: &PyEquationOfState, - temperature_or_pressure: Bound<'_, PyAny>, + temperature_or_pressure: &Bound<'_, PyAny>, initial_state: Option<&PyPhaseEquilibrium>, max_iter: Option, tol: Option, @@ -198,9 +198,9 @@ impl PyPhaseEquilibrium { #[expect(clippy::too_many_arguments)] pub(crate) fn bubble_point<'py>( eos: &PyEquationOfState, - temperature_or_pressure: Bound<'_, PyAny>, + temperature_or_pressure: &Bound<'_, PyAny>, liquid_molefracs: PyReadonlyArray1<'py, f64>, - tp_init: Option>, + tp_init: Option<&Bound<'_, PyAny>>, vapor_molefracs: Option>, max_iter_inner: Option, max_iter_outer: Option, @@ -286,9 +286,9 @@ impl PyPhaseEquilibrium { #[expect(clippy::too_many_arguments)] pub(crate) fn dew_point<'py>( eos: &PyEquationOfState, - temperature_or_pressure: Bound<'_, PyAny>, + temperature_or_pressure: &Bound<'_, PyAny>, vapor_molefracs: PyReadonlyArray1<'py, f64>, - tp_init: Option>, + tp_init: Option<&Bound<'_, PyAny>>, liquid_molefracs: Option>, max_iter_inner: Option, max_iter_outer: Option, @@ -398,7 +398,7 @@ impl PyPhaseEquilibrium { #[staticmethod] fn vle_pure_comps( eos: &PyEquationOfState, - temperature_or_pressure: Bound<'_, PyAny>, + temperature_or_pressure: &Bound<'_, PyAny>, ) -> PyResult>> { if let Ok(t) = temperature_or_pressure.extract::() { Ok(PhaseEquilibrium::vle_pure_comps(&eos.0, t) @@ -514,9 +514,9 @@ impl PyPhaseEquilibrium { #[expect(clippy::too_many_arguments)] fn heteroazeotrope( eos: &PyEquationOfState, - temperature_or_pressure: Bound<'_, PyAny>, + temperature_or_pressure: &Bound<'_, PyAny>, x_init: (f64, f64), - tp_init: Option>, + tp_init: Option<&Bound<'_, PyAny>>, max_iter: Option, tol: Option, verbosity: Option, @@ -1224,7 +1224,7 @@ impl PyPhaseDiagram { #[expect(clippy::too_many_arguments)] pub(crate) fn binary_vle( eos: &PyEquationOfState, - temperature_or_pressure: Bound<'_, PyAny>, + temperature_or_pressure: &Bound<'_, PyAny>, npoints: Option, x_lle: Option<(f64, f64)>, max_iter_inner: Option, @@ -1299,7 +1299,7 @@ impl PyPhaseDiagram { #[pyo3(signature = (eos, temperature_or_pressure, feed, min_tp, max_tp, npoints=None))] pub(crate) fn lle( eos: &PyEquationOfState, - temperature_or_pressure: Bound<'_, PyAny>, + temperature_or_pressure: &Bound<'_, PyAny>, feed: Moles>, min_tp: Bound<'_, PyAny>, max_tp: Bound<'_, PyAny>, @@ -1388,10 +1388,10 @@ impl PyPhaseDiagram { #[expect(clippy::too_many_arguments)] pub(crate) fn binary_vlle( eos: &PyEquationOfState, - temperature_or_pressure: Bound<'_, PyAny>, + temperature_or_pressure: &Bound<'_, PyAny>, x_lle: (f64, f64), - tp_lim_lle: Option>, - tp_init_vlle: Option>, + tp_lim_lle: Option<&Bound<'_, PyAny>>, + tp_init_vlle: Option<&Bound<'_, PyAny>>, npoints_vle: Option, npoints_lle: Option, max_iter_inner: Option, diff --git a/py-feos/src/state.rs b/py-feos/src/state.rs index e7b54211c..6d5ff09f8 100644 --- a/py-feos/src/state.rs +++ b/py-feos/src/state.rs @@ -1,4 +1,4 @@ -use crate::eos::parse_molefracs; +use crate::eos::{Compositions, parse_molefracs}; use crate::{ PyVerbosity, eos::PyEquationOfState, error::PyFeosError, ideal_gas::IdealGasModel, residual::ResidualModel, @@ -13,15 +13,13 @@ use pyo3::exceptions::{PyIndexError, PyValueError}; use pyo3::prelude::*; use quantity::*; use std::collections::HashMap; -use std::ops::{Deref, Div, Neg, Sub}; +use std::ops::{Deref, Div, Neg}; use std::sync::Arc; type Quot = >::Output; -type DpDn = Quantity>::Output>; type InvT = Quantity::Output>; type InvP = Quantity::Output>; -type InvM = Quantity::Output>; /// Possible contributions that can be computed. #[derive(Clone, Copy, PartialEq)] @@ -69,14 +67,8 @@ impl From for Contributions { /// Volume. /// density : SINumber, optional /// Molar density. -/// partial_density : SIArray1, optional -/// Partial molar densities. -/// total_moles : SINumber, optional -/// Total amount of substance (of a mixture). -/// moles : SIArray1, optional -/// Amount of substance for each component. -/// molefracs : numpy.ndarray[float] -/// Molar fraction of each component. +/// composition : float | SINumber | numpy.ndarray[float] | SIArray1 | list[float], optional +/// Composition of the mixture. /// pressure : SINumber, optional /// Pressure. /// molar_enthalpy : SINumber, optional @@ -110,19 +102,16 @@ pub struct PyState(pub State, ResidualMod impl PyState { #[new] #[pyo3( - text_signature = "(eos, temperature=None, volume=None, density=None, partial_density=None, total_moles=None, moles=None, molefracs=None, pressure=None, molar_enthalpy=None, molar_entropy=None, molar_internal_energy=None, density_initialization=None, initial_temperature=None)" + text_signature = "(eos, temperature=None, volume=None, density=None, composition=None, pressure=None, molar_enthalpy=None, molar_entropy=None, molar_internal_energy=None, density_initialization=None, initial_temperature=None)" )] - #[pyo3(signature = (eos, temperature=None, volume=None, density=None, partial_density=None, total_moles=None, moles=None, molefracs=None, pressure=None, molar_enthalpy=None, molar_entropy=None, molar_internal_energy=None, density_initialization=None, initial_temperature=None))] + #[pyo3(signature = (eos, temperature=None, volume=None, density=None, composition=None, pressure=None, molar_enthalpy=None, molar_entropy=None, molar_internal_energy=None, density_initialization=None, initial_temperature=None))] #[expect(clippy::too_many_arguments)] pub fn new<'py>( eos: &PyEquationOfState, temperature: Option, volume: Option, density: Option, - partial_density: Option>>, - total_moles: Option, - moles: Option>>, - molefracs: Option>, + composition: Option<&Bound<'py, PyAny>>, pressure: Option, molar_enthalpy: Option, molar_entropy: Option, @@ -130,7 +119,7 @@ impl PyState { density_initialization: Option<&Bound<'py, PyAny>>, initial_temperature: Option, ) -> PyResult { - let x = parse_molefracs(molefracs); + let composition = Compositions::try_from(composition)?; let density_init = if let Some(di) = density_initialization { if let Ok(d) = di.extract::().as_deref() { match d { @@ -147,25 +136,23 @@ impl PyState { } } else { Ok(None) - }; - let s = State::new_full( - &eos.0, - temperature, - volume, - density, - partial_density.as_ref(), - total_moles, - moles.as_ref(), - x.as_ref(), - pressure, - molar_enthalpy, - molar_entropy, - molar_internal_energy, - density_init?, - initial_temperature, - ) - .map_err(PyFeosError::from)?; - Ok(Self(s)) + }?; + Ok(Self( + State::build_full( + &eos.0, + temperature, + volume, + density, + composition, + pressure, + molar_enthalpy, + molar_entropy, + molar_internal_energy, + density_init, + initial_temperature, + ) + .map_err(PyFeosError::from)?, + )) } /// Return a list of thermodynamic state at critical conditions @@ -233,12 +220,12 @@ impl PyState { /// State : State at critical conditions. #[staticmethod] #[pyo3( - text_signature = "(eos, molefracs=None, initial_temperature=None, initial_density=None, max_iter=None, tol=None, verbosity=None)" + text_signature = "(eos, composition=None, initial_temperature=None, initial_density=None, max_iter=None, tol=None, verbosity=None)" )] - #[pyo3(signature = (eos, molefracs=None, initial_temperature=None, initial_density=None, max_iter=None, tol=None, verbosity=None))] + #[pyo3(signature = (eos, composition=None, initial_temperature=None, initial_density=None, max_iter=None, tol=None, verbosity=None))] fn critical_point<'py>( eos: &PyEquationOfState, - molefracs: Option>, + composition: Option<&Bound<'py, PyAny>>, initial_temperature: Option, initial_density: Option, max_iter: Option, @@ -248,7 +235,7 @@ impl PyState { Ok(PyState( State::critical_point( &eos.0, - parse_molefracs(molefracs).as_ref(), + Compositions::try_from(composition)?, initial_temperature, initial_density, (max_iter, tol, verbosity.map(|v| v.into())).into(), @@ -287,7 +274,7 @@ impl PyState { #[expect(clippy::too_many_arguments)] fn critical_point_binary( eos: &PyEquationOfState, - temperature_or_pressure: Bound<'_, PyAny>, + temperature_or_pressure: &Bound<'_, PyAny>, initial_temperature: Option, initial_molefracs: Option<[f64; 2]>, initial_density: Option, @@ -353,13 +340,13 @@ impl PyState { /// (State, State): Spinodal states. #[staticmethod] #[pyo3( - text_signature = "(eos, temperature, molefracs=None, max_iter=None, tol=None, verbosity=None)" + text_signature = "(eos, temperature, composition=None, max_iter=None, tol=None, verbosity=None)" )] - #[pyo3(signature = (eos, temperature, molefracs=None, max_iter=None, tol=None, verbosity=None))] + #[pyo3(signature = (eos, temperature, composition=None, max_iter=None, tol=None, verbosity=None))] fn spinodal<'py>( eos: &PyEquationOfState, temperature: Temperature, - molefracs: Option>, + composition: Option<&Bound<'py, PyAny>>, max_iter: Option, tol: Option, verbosity: Option, @@ -367,7 +354,7 @@ impl PyState { let [state1, state2] = State::spinodal( &eos.0, temperature, - parse_molefracs(molefracs).as_ref(), + Compositions::try_from(composition)?, (max_iter, tol, verbosity.map(|v| v.into())).into(), ) .map_err(PyFeosError::from)?; @@ -476,7 +463,7 @@ impl PyState { self.0.compressibility(contributions.into()) } - /// Return partial derivative of pressure w.r.t. volume. + /// Return partial derivative of pressure w.r.t. molar volume. /// /// Parameters /// ---------- @@ -488,7 +475,7 @@ impl PyState { /// ------- /// SINumber #[pyo3(signature = (contributions=PyContributions::Total), text_signature = "($self, contributions)")] - fn dp_dv(&self, contributions: PyContributions) -> Quot { + fn dp_dv(&self, contributions: PyContributions) -> Quot { self.0.dp_dv(contributions.into()) } @@ -536,11 +523,11 @@ impl PyState { /// ------- /// SIArray1 #[pyo3(signature = (contributions=PyContributions::Total), text_signature = "($self, contributions)")] - fn dp_dni(&self, contributions: PyContributions) -> DpDn> { - self.0.dp_dni(contributions.into()) + fn dp_dni(&self, contributions: PyContributions) -> Pressure> { + self.0.n_dp_dni(contributions.into()) } - /// Return second partial derivative of pressure w.r.t. volume. + /// Return second partial derivative of pressure w.r.t. molar volume. /// /// Parameters /// ---------- @@ -552,7 +539,10 @@ impl PyState { /// ------- /// SINumber #[pyo3(signature = (contributions=PyContributions::Total), text_signature = "($self, contributions)")] - fn d2p_dv2(&self, contributions: PyContributions) -> Quot, Volume> { + fn d2p_dv2( + &self, + contributions: PyContributions, + ) -> Quot, MolarVolume> { self.0.d2p_dv2(contributions.into()) } @@ -652,8 +642,8 @@ impl PyState { /// ------- /// SIArray2 #[pyo3(signature = (contributions=PyContributions::Total), text_signature = "($self, contributions)")] - fn dmu_dni(&self, contributions: PyContributions) -> Quot>, Moles> { - self.0.dmu_dni(contributions.into()) + fn n_dmu_dni(&self, contributions: PyContributions) -> MolarEnergy> { + self.0.n_dmu_dni(contributions.into()) } /// Return logarithmic fugacity coefficient. @@ -771,8 +761,8 @@ impl PyState { /// Returns /// ------- /// SIArray2 - fn dln_phi_dnj(&self) -> InvM> { - self.0.dln_phi_dnj() + fn n_dln_phi_dnj<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray2> { + self.0.n_dln_phi_dnj().to_pyarray(py) } /// Return thermodynamic factor. @@ -848,7 +838,7 @@ impl PyState { self.0.entropy(contributions.into()) } - /// Return derivative of entropy with respect to temperature. + /// Return derivative of molar entropy with respect to temperature. /// /// Parameters /// ---------- @@ -860,7 +850,7 @@ impl PyState { /// ------- /// SINumber #[pyo3(signature = (contributions=PyContributions::Total), text_signature = "($self, contributions)")] - fn ds_dt(&self, contributions: PyContributions) -> Quot { + fn ds_dt(&self, contributions: PyContributions) -> Quot { self.0.ds_dt(contributions.into()) } @@ -1357,7 +1347,7 @@ impl PyState { #[getter] fn get_total_moles(&self) -> Moles { - self.0.total_moles + self.0.total_moles() } #[getter] @@ -1367,7 +1357,7 @@ impl PyState { #[getter] fn get_volume(&self) -> Volume { - self.0.volume + self.0.volume() } #[getter] @@ -1377,12 +1367,12 @@ impl PyState { #[getter] fn get_moles(&self) -> Moles> { - self.0.moles.clone() + self.0.moles() } #[getter] fn get_partial_density(&self) -> Density> { - self.0.partial_density.clone() + self.0.partial_density() } #[getter]