type Other: fmt::Display + TemperatureOrPressure;

const IDENTIFIER: &'static str;

fn temperature_pressure(&self, tp_init: Self::Other) -> (Temperature, Pressure);

fn from_state<E: Residual>(state: &State<E>) -> Self::Other;

fn linspace(

&self,

start: Self::Other,

end: Self::Other,

n: usize,

) -> (Temperature<Array1<f64>>, Pressure<Array1<f64>>);

fn bubble_dew_point<E: Residual>(

eos: &Arc<E>,

tp_spec: Self,

tp_init: Option<Self::Other>,

molefracs_spec: &Array1<f64>,

molefracs_init: Option<&Array1<f64>>,

bubble: bool,

options: (SolverOptions, SolverOptions),

) -> EosResult<PhaseEquilibrium<E, 2>>;

fn adjust_t_p<E: Residual>(

tp_spec: Self,

var: &mut Self::Other,

state1: &mut State<E>,

state2: &mut State<E>,

verbosity: Verbosity,

) -> EosResult<f64>;

fn newton_step<E: Residual>(

tp_spec: Self,

var: &mut Self::Other,

state1: &mut State<E>,

state2: &mut State<E>,

verbosity: Verbosity,

) -> EosResult<f64>;

type Other = Pressure;

const IDENTIFIER: &'static str = "temperature";

fn temperature_pressure(&self, tp_init: Self::Other) -> (Temperature, Pressure) {

(*self, tp_init)

}

fn from_state<E: Residual>(state: &State<E>) -> Self::Other {

state.pressure(Contributions::Total)

}

fn linspace(

&self,

start: Pressure,

end: Pressure,

n: usize,

) -> (Temperature<Array1<f64>>, Pressure<Array1<f64>>) {

(

Temperature::linspace(*self, *self, n),

Pressure::linspace(start, end, n),

)

}

fn bubble_dew_point<E: Residual>(

eos: &Arc<E>,

temperature: Self,

p_init: Option<Pressure>,

molefracs_spec: &Array1<f64>,

molefracs_init: Option<&Array1<f64>>,

bubble: bool,

options: (SolverOptions, SolverOptions),

) -> EosResult<PhaseEquilibrium<E, 2>> {

// First use given initial pressure if applicable

if let Some(p) = p_init {

return PhaseEquilibrium::iterate_bubble_dew(

eos,

temperature,

p,

molefracs_spec,

molefracs_init,

bubble,

options,

);

}

// Next try to initialize with an ideal gas assumption

let vle =

PhaseEquilibrium::starting_pressure_ideal_gas(eos, temperature, molefracs_spec, bubble)

.and_then(|(p, x)| {

PhaseEquilibrium::iterate_bubble_dew(

eos,

temperature,

p,

molefracs_spec,

molefracs_init.or(Some(&x)),

bubble,

options,

)

});

// Finally use the spinodal to initialize the calculation

vle.or_else(|_| {

PhaseEquilibrium::iterate_bubble_dew(

eos,

temperature,

PhaseEquilibrium::starting_pressure_spinodal(eos, temperature, molefracs_spec)?,

molefracs_spec,

molefracs_init,

bubble,

options,

)

})

}

fn adjust_t_p<E: Residual>(

temperature: Temperature,

pressure: &mut Pressure,

state1: &mut State<E>,

state2: &mut State<E>,

verbosity: Verbosity,

) -> EosResult<f64> {

// calculate K = phi_1/phi_2 = x_2/x_1

let ln_phi_1 = state1.ln_phi();

let ln_phi_2 = state2.ln_phi();

let k = (&ln_phi_1 - &ln_phi_2).mapv(f64::exp);

// calculate residual

let f = (&state1.molefracs * &k).sum() - 1.0;

// Derivative w.r.t. ln(pressure)

let ln_phi_1_dp = state1.dln_phi_dp();

let ln_phi_2_dp = state2.dln_phi_dp();

let df =

(((ln_phi_1_dp - ln_phi_2_dp) * *pressure).into_value() * &state1.molefracs * &k).sum();

let mut lnpstep = -f / df;

// catch too big p-steps

lnpstep = lnpstep.clamp(-MAX_LNPSTEP, MAX_LNPSTEP);

// Update p

*pressure *= lnpstep.exp();

// update states with new temperature/pressure

adjust_states(temperature, *pressure, state1, state2, None)?;

// log

log_iter!(

verbosity,

"{:14} | {:<14.8e} | {:12.8} | {:.8}",

"",

f.abs(),

pressure,

state2.molefracs

);

Ok(f.abs())

}

fn newton_step<E: Residual>(

_: Temperature,

pressure: &mut Pressure,

state1: &mut State<E>,

state2: &mut State<E>,

verbosity: Verbosity,

) -> EosResult<f64> {

let dmu_drho_1 = (state1.dmu_dni(Contributions::Total) * state1.volume)

.to_reduced()

.dot(&state1.molefracs);

let dmu_drho_2 = (state2.dmu_dni(Contributions::Total) * state2.volume).to_reduced();

let dp_drho_1 = (state1.dp_dni(Contributions::Total) * state1.volume)

.to_reduced()

.dot(&state1.molefracs);

let dp_drho_2 = (state2.dp_dni(Contributions::Total) * state2.volume).to_reduced();

let mu_1_res = state1.residual_chemical_potential().to_reduced();

let mu_2_res = state2.residual_chemical_potential().to_reduced();

let p_1 = state1.pressure(Contributions::Total).to_reduced();

let p_2 = state2.pressure(Contributions::Total).to_reduced();

// calculate residual

let dmu_ig = (RGAS * state1.temperature).to_reduced()

* (&state1.partial_density / &state2.partial_density)

.into_value()

.mapv(f64::ln);

let res = concatenate![Axis(0), mu_1_res - mu_2_res + dmu_ig, arr1(&[p_1 - p_2])];

let error = norm(&res);

// calculate Jacobian

let jacobian = concatenate![

Axis(1),

concatenate![Axis(0), -dmu_drho_2, -dp_drho_2.insert_axis(Axis(0))],

concatenate![

Axis(0),

dmu_drho_1.insert_axis(Axis(1)),

arr2(&[[dp_drho_1]])

]

];

// calculate Newton step

let dx = LU::new(jacobian)?.solve(&res);

// apply Newton step

let rho_l1 = state1.density - Density::from_reduced(dx[dx.len() - 1]);

let rho_l2 =

state2.partial_density.clone() - Density::from_reduced(dx.slice(s![0..-1]).to_owned());

// update states

*state1 = StateBuilder::new(&state1.eos)

.temperature(state1.temperature)

.density(rho_l1)

.molefracs(&state1.molefracs)

.build()?;

*state2 = StateBuilder::new(&state2.eos)

.temperature(state2.temperature)

.partial_density(&rho_l2)

.build()?;

*pressure = state1.pressure(Contributions::Total);

log_iter!(

verbosity,

"{:<14.8e} | {:14} | {:12.8} | {:.8} NEWTON",

error,

"",

pressure,

state2.molefracs

);

Ok(error)

}

type Other = Temperature;

const IDENTIFIER: &'static str = "pressure";

fn temperature_pressure(&self, tp_init: Self::Other) -> (Temperature, Pressure) {

(tp_init, *self)

}

fn from_state<E: Residual>(state: &State<E>) -> Self::Other {

state.temperature

}

fn linspace(

&self,

start: Temperature,

end: Temperature,

n: usize,

) -> (Temperature<Array1<f64>>, Pressure<Array1<f64>>) {

(

Temperature::linspace(start, end, n),

Pressure::linspace(*self, *self, n),

)

}

fn bubble_dew_point<E: Residual>(

eos: &Arc<E>,

pressure: Self,

t_init: Option<Temperature>,

molefracs_spec: &Array1<f64>,

molefracs_init: Option<&Array1<f64>>,

bubble: bool,

options: (SolverOptions, SolverOptions),

) -> EosResult<PhaseEquilibrium<E, 2>> {

let temperature = t_init.expect("An initial temperature is required for the calculation of bubble/dew points at given pressure!");

PhaseEquilibrium::iterate_bubble_dew(

eos,

pressure,

temperature,

molefracs_spec,

molefracs_init,

bubble,

options,

)

}

fn adjust_t_p<E: Residual>(

pressure: Pressure,

temperature: &mut Temperature,

state1: &mut State<E>,

state2: &mut State<E>,

verbosity: Verbosity,

) -> EosResult<f64> {

// calculate K = phi_1/phi_2 = x_2/x_1

let ln_phi_1 = state1.ln_phi();

let ln_phi_2 = state2.ln_phi();

let k = (&ln_phi_1 - &ln_phi_2).mapv(f64::exp);

// calculate residual

let f = (&state1.molefracs * &k).sum() - 1.0;

// Derivative w.r.t. temperature

let ln_phi_1_dt = state1.dln_phi_dt();

let ln_phi_2_dt = state2.dln_phi_dt();

let df = ((ln_phi_1_dt - ln_phi_2_dt) * Dimensionless::from(&state1.molefracs * &k)).sum();

let mut tstep = -f / df;

// catch too big t-steps

if tstep < -Temperature::from_reduced(MAX_TSTEP) {

tstep = -Temperature::from_reduced(MAX_TSTEP);

} else if tstep > Temperature::from_reduced(MAX_TSTEP) {

tstep = Temperature::from_reduced(MAX_TSTEP);

}

// Update t

*temperature += tstep;

// update states with new temperature

adjust_states(*temperature, pressure, state1, state2, None)?;

// log

log_iter!(

verbosity,

"{:14} | {:<14.8e} | {:12.8} | {:.8}",

"",

f.abs(),

temperature,

state2.molefracs

);

Ok(f.abs())

}

fn newton_step<E: Residual>(

pressure: Pressure,

temperature: &mut Temperature,

state1: &mut State<E>,

state2: &mut State<E>,

verbosity: Verbosity,

) -> EosResult<f64> {

let dmu_drho_1 = (state1.dmu_dni(Contributions::Total) * state1.volume)

.to_reduced()

.dot(&state1.molefracs);

let dmu_drho_2 = (state2.dmu_dni(Contributions::Total) * state2.volume).to_reduced();

let dmu_res_dt_1 = state1.dmu_res_dt().to_reduced();

let dmu_res_dt_2 = state2.dmu_res_dt().to_reduced();

let dp_drho_1 = (state1.dp_dni(Contributions::Total) * state1.volume)

.to_reduced()

.dot(&state1.molefracs);

let dp_dt_1 = state1.dp_dt(Contributions::Total).to_reduced();

let dp_dt_2 = state2.dp_dt(Contributions::Total).to_reduced();

let dp_drho_2 = (state2.dp_dni(Contributions::Total) * state2.volume).to_reduced();

let mu_1_res = state1.residual_chemical_potential().to_reduced();

let mu_2_res = state2.residual_chemical_potential().to_reduced();

let p_1 = state1.pressure(Contributions::Total).to_reduced();

let p_2 = state2.pressure(Contributions::Total).to_reduced();

let p = pressure.to_reduced();

// calculate residual

let delta_dmu_ig_dt = (&state1.partial_density / &state2.partial_density)

.into_value()

.mapv(f64::ln);

let delta_mu_ig = (RGAS * state1.temperature).to_reduced() * &delta_dmu_ig_dt;

let res = concatenate![

Axis(0),

mu_1_res - mu_2_res + delta_mu_ig,

arr1(&[p_1 - p]),

arr1(&[p_2 - p])

];

let error = norm(&res);

// calculate Jacobian

let jacobian = concatenate![

Axis(1),

concatenate![

Axis(0),

-dmu_drho_2,

Array2::zeros((1, res.len() - 2)),

dp_drho_2.insert_axis(Axis(0))

],

concatenate![

Axis(0),

dmu_drho_1.insert_axis(Axis(1)),

arr2(&[[dp_drho_1], [0.0]])

],

concatenate![

Axis(0),

(dmu_res_dt_1 - dmu_res_dt_2 + delta_dmu_ig_dt).insert_axis(Axis(1)),

arr2(&[[dp_dt_1], [dp_dt_2]])

]

];

// calculate Newton step

let dx = LU::new(jacobian)?.solve(&res);

// apply Newton step

let rho_l1 = state1.density - Density::from_reduced(dx[dx.len() - 2]);

let rho_l2 =

state2.partial_density.clone() - Density::from_reduced(dx.slice(s![0..-2]).to_owned());

let t = state1.temperature - Temperature::from_reduced(dx[dx.len() - 1]);

// update states

*state1 = StateBuilder::new(&state1.eos)

.temperature(t)

.density(rho_l1)

.molefracs(&state1.molefracs)

.build()?;

*state2 = StateBuilder::new(&state2.eos)

.temperature(t)

.partial_density(&rho_l2)

.build()?;

*temperature = t;

log_iter!(

verbosity,

"{:<14.8e} | {:14} | {:12.8} | {:.8} NEWTON",

error,

"",

temperature,

state2.molefracs

);

Ok(error)

}

/// Calculate a phase equilibrium for a given temperature

/// or pressure and composition of the liquid phase.

pub fn bubble_point<TP: TemperatureOrPressure>(

eos: &Arc<E>,

temperature_or_pressure: TP,

liquid_molefracs: &Array1<f64>,

tp_init: Option<TP::Other>,

vapor_molefracs: Option<&Array1<f64>>,

options: (SolverOptions, SolverOptions),

) -> EosResult<Self> {

Self::bubble_dew_point(

eos,

temperature_or_pressure,

tp_init,

liquid_molefracs,

vapor_molefracs,

true,

options,

)

}

/// Calculate a phase equilibrium for a given temperature

/// or pressure and composition of the vapor phase.

pub fn dew_point<TP: TemperatureOrPressure>(

eos: &Arc<E>,

temperature_or_pressure: TP,

vapor_molefracs: &Array1<f64>,

tp_init: Option<TP::Other>,

liquid_molefracs: Option<&Array1<f64>>,

options: (SolverOptions, SolverOptions),

) -> EosResult<Self> {

Self::bubble_dew_point(

eos,

temperature_or_pressure,

tp_init,

vapor_molefracs,

liquid_molefracs,

false,

options,

)

}

pub(super) fn bubble_dew_point<TP: TemperatureOrPressure>(

eos: &Arc<E>,

tp_spec: TP,

tp_init: Option<TP::Other>,

molefracs_spec: &Array1<f64>,

molefracs_init: Option<&Array1<f64>>,

bubble: bool,

options: (SolverOptions, SolverOptions),

) -> EosResult<Self> {

TP::bubble_dew_point(

eos,

tp_spec,

tp_init,

molefracs_spec,

molefracs_init,

bubble,

options,

)

}

fn iterate_bubble_dew<TP: TemperatureOrPressure>(

eos: &Arc<E>,

tp_spec: TP,

tp_init: TP::Other,

molefracs_spec: &Array1<f64>,

molefracs_init: Option<&Array1<f64>>,

bubble: bool,

options: (SolverOptions, SolverOptions),

) -> EosResult<Self> {

let (t, p) = tp_spec.temperature_pressure(tp_init);

let [state1, state2] = if bubble {

starting_x2_bubble(eos, t, p, molefracs_spec, molefracs_init)

} else {

starting_x2_dew(eos, t, p, molefracs_spec, molefracs_init)

}?;

bubble_dew(tp_spec, tp_init, state1, state2, bubble, options)

}

fn starting_pressure_ideal_gas(

eos: &Arc<E>,

temperature: Temperature,

molefracs_spec: &Array1<f64>,

bubble: bool,

) -> EosResult<(Pressure, Array1<f64>)> {

if bubble {

Self::starting_pressure_ideal_gas_bubble(eos, temperature, molefracs_spec)

} else {

Self::starting_pressure_ideal_gas_dew(eos, temperature, molefracs_spec)

}

}

pub(super) fn starting_pressure_ideal_gas_bubble(

eos: &Arc<E>,

temperature: Temperature,

liquid_molefracs: &Array1<f64>,

) -> EosResult<(Pressure, Array1<f64>)> {

let m = Moles::from_reduced(liquid_molefracs.to_owned());

let density = 0.75 * eos.max_density(Some(&m))?;

let liquid = State::new_nvt(eos, temperature, m.sum() / density, &m)?;

let v_l = liquid.partial_molar_volume();

let p_l = liquid.pressure(Contributions::Total);

let mu_l = liquid.residual_chemical_potential();

let p_i = (liquid_molefracs * temperature * density * RGAS)

* Dimensionless::from(

((mu_l - p_l * v_l) / (RGAS * temperature))

.into_value()

.mapv(f64::exp),

);

let p = p_i.sum();

let y = (p_i / p).into_value();

Ok((p, y))

}

fn starting_pressure_ideal_gas_dew(

eos: &Arc<E>,

temperature: Temperature,

vapor_molefracs: &Array1<f64>,

) -> EosResult<(Pressure, Array1<f64>)> {

let mut p: Option<Pressure> = None;

let mut x = vapor_molefracs.clone();

for _ in 0..5 {

let m = Moles::from_reduced(x);

let density = 0.75 * eos.max_density(Some(&m))?;

let liquid = State::new_nvt(eos, temperature, m.sum() / density, &m)?;

let v_l = liquid.partial_molar_volume();

let p_l = liquid.pressure(Contributions::Total);

let mu_l = liquid.residual_chemical_potential();

let k = vapor_molefracs

/ ((mu_l - p_l * v_l) / (RGAS * temperature))

.into_value()

.mapv(f64::exp);

let p_new = RGAS * temperature * density / k.sum();

x = &k / k.sum();

if let Some(p_old) = p {

if ((p_new - p_old) / p_old).into_value().abs() < 1e-5 {

p = Some(p_new);

break;

}

}

p = Some(p_new);

}

Ok((p.unwrap(), x))

}

pub(super) fn starting_pressure_spinodal(

eos: &Arc<E>,

temperature: Temperature,

molefracs: &Array1<f64>,

) -> EosResult<Pressure> {

let moles = Moles::from_reduced(molefracs.clone());

let [sp_v, sp_l] = State::spinodal(eos, temperature, Some(&moles), 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))

}

eos: &Arc<E>,

temperature: Temperature,

pressure: Pressure,

liquid_molefracs: &Array1<f64>,

vapor_molefracs: Option<&Array1<f64>>,

) -> EosResult<[State<E>; 2]> {

let liquid_state = State::new_npt(

eos,

temperature,

pressure,

&Moles::from_reduced(liquid_molefracs.clone()),

Liquid,

)?;

let xv = match vapor_molefracs {

Some(xv) => xv.clone(),

None => liquid_state.ln_phi().mapv(f64::exp) * liquid_molefracs,

};

let vapor_state = State::new_npt(eos, temperature, pressure, &Moles::from_reduced(xv), Vapor)?;

Ok([liquid_state, vapor_state])

eos: &Arc<E>,

temperature: Temperature,

pressure: Pressure,

vapor_molefracs: &Array1<f64>,

liquid_molefracs: Option<&Array1<f64>>,

) -> EosResult<[State<E>; 2]> {

let vapor_state = State::new_npt(

eos,

temperature,

pressure,

&Moles::from_reduced(vapor_molefracs.clone()),

Vapor,

)?;

let xl = match liquid_molefracs {

Some(xl) => xl.clone(),

None => {

let xl = vapor_state.ln_phi().mapv(f64::exp) * vapor_molefracs;

let liquid_state =

State::new_npt(eos, temperature, pressure, &Moles::from_reduced(xl), Liquid)?;

(vapor_state.ln_phi() - liquid_state.ln_phi()).mapv(f64::exp) * vapor_molefracs

}

};

let liquid_state =

State::new_npt(eos, temperature, pressure, &Moles::from_reduced(xl), Liquid)?;

Ok([vapor_state, liquid_state])

tp_spec: TP,

mut var_tp: TP::Other,

mut state1: State<E>,

mut state2: State<E>,

bubble: bool,

options: (SolverOptions, SolverOptions),

) -> EosResult<PhaseEquilibrium<E, 2>> {

let (options_inner, options_outer) = options;

// initialize variables

let mut err_out = 1.0;

let mut k_out = 0;

if PhaseEquilibrium::is_trivial_solution(&state1, &state2) {

log_iter!(options_outer.verbosity, "Trivial solution encountered!");

return Err(EosError::TrivialSolution);

}

log_iter!(

options_outer.verbosity,

"res outer loop | res inner loop | {:^16} | molefracs second phase",

TP::IDENTIFIER

);

log_iter!(options_outer.verbosity, "{:-<85}", "");

log_iter!(

options_outer.verbosity,

"{:14} | {:14} | {:12.8} | {:.8}",

"",

"",

var_tp,

state2.molefracs

);

// Outer loop for finding x2

for ko in 0..options_outer.max_iter.unwrap_or(MAX_ITER_OUTER) {

// Iso-Fugacity equation

err_out = if err_out > NEWTON_TOL {

// Inner loop for finding T or p

for _ in 0..options_inner.max_iter.unwrap_or(MAX_ITER_INNER) {

if TP::adjust_t_p(

tp_spec,

&mut var_tp,

&mut state1,

&mut state2,

options_inner.verbosity,

)? < options_inner.tol.unwrap_or(TOL_INNER)

{

break;

}

}

adjust_x2(&state1, &mut state2, options_outer.verbosity)

} else {

TP::newton_step(

tp_spec,

&mut var_tp,

&mut state1,

&mut state2,

options_outer.verbosity,

)

}?;

if PhaseEquilibrium::is_trivial_solution(&state1, &state2) {

log_iter!(options_outer.verbosity, "Trivial solution encountered!");

return Err(EosError::TrivialSolution);

}

if err_out < options_outer.tol.unwrap_or(TOL_OUTER) {

k_out = ko + 1;

break;

}

}

if err_out < options_outer.tol.unwrap_or(TOL_OUTER) {

log_result!(

options_outer.verbosity,

"Bubble/dew point: calculation converged in {} step(s)\n",

k_out

);

if bubble {

Ok(PhaseEquilibrium([state2, state1]))

} else {

Ok(PhaseEquilibrium([state1, state2]))

}

} else {

// not converged, return EosError

Err(EosError::NotConverged(String::from("bubble-dew-iteration")))

}

temperature: Temperature,

pressure: Pressure,

state1: &mut State<E>,

state2: &mut State<E>,

moles_state2: Option<&Moles<Array1<f64>>>,

) -> EosResult<()> {

*state1 = State::new_npt(

&state1.eos,

temperature,

pressure,

&state1.moles,

InitialDensity(state1.density),

)?;

*state2 = State::new_npt(

&state2.eos,

temperature,

pressure,

moles_state2.unwrap_or(&state2.moles),

InitialDensity(state2.density),

)?;

Ok(())

state1: &State<E>,

state2: &mut State<E>,

verbosity: Verbosity,

) -> EosResult<f64> {

let x1 = &state1.molefracs;

let ln_phi_1 = state1.ln_phi();

let ln_phi_2 = state2.ln_phi();

let k = (ln_phi_1 - ln_phi_2).mapv(f64::exp);

let err_out = (&k * x1 / &state2.molefracs - 1.0).mapv(f64::abs).sum();

let x2 = (x1 * &k) / (&k * x1).sum();

log_iter!(verbosity, "{:<14.8e} | {:14} | {:16} |", err_out, "", "");

*state2 = State::new_npt(

&state2.eos,

state2.temperature,

state2.pressure(Contributions::Total),

&Moles::from_reduced(x2),

InitialDensity(state2.density),

)?;

Ok(err_out)