Skip to content

Common association implementation for different SAFT models #239

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
prehner merged 7 commits into from
Apr 19, 2024
Merged

Common association implementation for different SAFT models #239

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
70cc640
Common association implementation for different SAFT models
prehner Apr 17, 2024
a3a51eb
more generic association implementation
prehner Apr 17, 2024
7ef046f
major refactor
prehner Apr 17, 2024
f9cda7e
Require kappa_ab and epsilon_k_ab
prehner Apr 18, 2024
d6a9229
fix epcsaft
prehner Apr 19, 2024
c4c7057
fix wheels
prehner Apr 19, 2024
7f4b58a
fix docs and remove mu an q from epcsaft
prehner Apr 19, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
Common association implementation for different SAFT models
  • Loading branch information
@prehner
prehner committed Apr 19, 2024
commit 70cc640e4b23a68013a43b490116a29ef2d188ab
90 changes: 52 additions & 38 deletions src/association/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -134,16 +134,15 @@ pub struct AssociationParameters {
sites_a: Array1<AssociationSite>,
sites_b: Array1<AssociationSite>,
sites_c: Array1<AssociationSite>,
pub sigma3_kappa_ab: Array2<f64>,
pub sigma3_kappa_cc: Array2<f64>,
pub kappa_ab: Array2<f64>,
pub kappa_cc: Array2<f64>,
pub epsilon_k_ab: Array2<f64>,
pub epsilon_k_cc: Array2<f64>,
}

impl AssociationParameters {
pub fn new(
records: &[Vec<AssociationRecord>],
sigma: &Array1<f64>,
binary_records: &[((usize, usize), BinaryAssociationRecord)],
component_index: Option<&Array1<usize>>,
) -> Self {
Expand Down Expand Up @@ -201,14 +200,11 @@ impl AssociationParameters {
.map(|(i, site)| ((site.assoc_comp, site.site_index), i))
.collect();

let mut sigma3_kappa_ab =
Array2::from_shape_fn([sites_a.len(), sites_b.len()], |(i, j)| {
(sigma[sites_a[i].assoc_comp] * sigma[sites_b[j].assoc_comp]).powf(1.5)
* (sites_a[i].kappa_ab * sites_b[j].kappa_ab).sqrt()
});
let mut sigma3_kappa_cc = Array2::from_shape_fn([sites_c.len(); 2], |(i, j)| {
(sigma[sites_c[i].assoc_comp] * sigma[sites_c[j].assoc_comp]).powf(1.5)
* (sites_c[i].kappa_ab * sites_c[j].kappa_ab).sqrt()
let mut kappa_ab = Array2::from_shape_fn([sites_a.len(), sites_b.len()], |(i, j)| {
(sites_a[i].kappa_ab * sites_b[j].kappa_ab).sqrt()
});
let mut kappa_cc = Array2::from_shape_fn([sites_c.len(); 2], |(i, j)| {
(sites_c[i].kappa_ab * sites_c[j].kappa_ab).sqrt()
});
let mut epsilon_k_ab = Array2::from_shape_fn([sites_a.len(), sites_b.len()], |(i, j)| {
0.5 * (sites_a[i].epsilon_k_ab + sites_b[j].epsilon_k_ab)
Expand All @@ -224,15 +220,15 @@ impl AssociationParameters {
epsilon_k_ab[[*x, *y]] = epsilon_k_aibj;
}
if let Some(kappa_aibj) = record.kappa_ab {
sigma3_kappa_ab[[*x, *y]] = (sigma[i] * sigma[j]).powf(1.5) * kappa_aibj;
kappa_ab[[*x, *y]] = kappa_aibj;
}
}
if let (Some(y), Some(x)) = (indices_b.get(&(i, a)), indices_a.get(&(j, b))) {
if let Some(epsilon_k_aibj) = record.epsilon_k_ab {
epsilon_k_ab[[*x, *y]] = epsilon_k_aibj;
}
if let Some(kappa_aibj) = record.kappa_ab {
sigma3_kappa_ab[[*x, *y]] = (sigma[i] * sigma[j]).powf(1.5) * kappa_aibj;
kappa_ab[[*x, *y]] = kappa_aibj;
}
}
if let (Some(x), Some(y)) = (indices_c.get(&(i, a)), indices_c.get(&(j, b))) {
Expand All @@ -241,8 +237,8 @@ impl AssociationParameters {
epsilon_k_cc[[*y, *x]] = epsilon_k_aibj;
}
if let Some(kappa_aibj) = record.kappa_ab {
sigma3_kappa_cc[[*x, *y]] = (sigma[i] * sigma[j]).powf(1.5) * kappa_aibj;
sigma3_kappa_cc[[*y, *x]] = (sigma[i] * sigma[j]).powf(1.5) * kappa_aibj;
kappa_cc[[*x, *y]] = kappa_aibj;
kappa_cc[[*y, *x]] = kappa_aibj;
}
}
}
Expand All @@ -254,8 +250,8 @@ impl AssociationParameters {
sites_a: Array1::from_vec(sites_a),
sites_b: Array1::from_vec(sites_b),
sites_c: Array1::from_vec(sites_c),
sigma3_kappa_ab,
sigma3_kappa_cc,
kappa_ab,
kappa_cc,
epsilon_k_ab,
epsilon_k_cc,
}
Expand Down Expand Up @@ -306,26 +302,33 @@ impl<P: HardSphereProperties> Association<P> {
fn association_strength<D: DualNum<f64> + Copy>(
&self,
temperature: D,
sigma: &Array1<f64>,
diameter: &Array1<D>,
n2: D,
n3i: D,
xi: D,
) -> [Array2<D>; 2] {
let p = &self.association_parameters;
let delta_ab = Array2::from_shape_fn([p.sites_a.len(), p.sites_b.len()], |(i, j)| {
let si = sigma[p.sites_a[i].assoc_comp];
let sj = sigma[p.sites_b[j].assoc_comp];
let di = diameter[p.sites_a[i].assoc_comp];
let dj = diameter[p.sites_b[j].assoc_comp];
let k = di * dj / (di + dj) * (n2 * n3i);
n3i * (k * xi * (k / 18.0 + 0.5) + 1.0)
* p.sigma3_kappa_ab[(i, j)]
* p.kappa_ab[(i, j)]
* (si * sj).powf(1.5)
* (temperature.recip() * p.epsilon_k_ab[(i, j)]).exp_m1()
});
let delta_cc = Array2::from_shape_fn([p.sites_c.len(); 2], |(i, j)| {
let si = sigma[p.sites_c[i].assoc_comp];
let sj = sigma[p.sites_c[j].assoc_comp];
let di = diameter[p.sites_c[i].assoc_comp];
let dj = diameter[p.sites_c[j].assoc_comp];
let k = di * dj / (di + dj) * (n2 * n3i);
n3i * (k * xi * (k / 18.0 + 0.5) + 1.0)
* p.sigma3_kappa_cc[(i, j)]
* p.kappa_cc[(i, j)]
* (si * sj).powf(1.5)
* (temperature.recip() * p.epsilon_k_cc[(i, j)]).exp_m1()
});
[delta_ab, delta_cc]
Expand All @@ -337,6 +340,7 @@ impl<P: HardSphereProperties> Association<P> {
pub fn helmholtz_energy<D: DualNum<f64> + Copy>(
&self,
state: &StateHD<D>,
sigma: &Array1<f64>,
diameter: &Array1<D>,
) -> D {
let p: &P = &self.parameters;
Expand All @@ -349,7 +353,7 @@ impl<P: HardSphereProperties> Association<P> {

// association strength
let [delta_ab, delta_cc] =
self.association_strength(state.temperature, diameter, n2, n3i, D::one());
self.association_strength(state.temperature, sigma, diameter, n2, n3i, D::one());

match (
a.sites_a.len() * a.sites_b.len(),
Expand Down Expand Up @@ -576,7 +580,6 @@ mod tests {
AssociationRecord::new(0.3, 2000., 0.0, 1.0, 0.0),
];
let records = [comp1, comp2, comp3, comp4];
let sigma = arr1(&[3.0, 3.0, 3.0, 3.0]);
let binary = [
(
(0, 1),
Expand All @@ -591,7 +594,7 @@ mod tests {
BinaryAssociationRecord::new(Some(3.5), Some(3333.), Some([0, 1])),
),
];
let assoc = AssociationParameters::new(&records, &sigma, &binary, None);
let assoc = AssociationParameters::new(&records, &binary, None);
println!("{}", assoc.epsilon_k_ab);
let epsilon_k_ab = arr2(&[
[2500., 1234., 3140., 2250.],
Expand All @@ -606,19 +609,18 @@ mod tests {
let comp1 = vec![AssociationRecord::new(0.1, 2500., 1.0, 1.0, 0.0)];
let comp2 = vec![AssociationRecord::new(0.1, -500., 0.0, 1.0, 0.0)];
let comp3 = vec![AssociationRecord::new(0.0, 0.0, 0.0, 1.0, 0.0)];
let sigma = arr1(&[3.0, 3.5]);
let binary = [(
(0, 1),
BinaryAssociationRecord::new(Some(0.1), Some(1000.), None),
)];
let assoc1 = AssociationParameters::new(&[comp1.clone(), comp2], &sigma, &[], None);
let assoc2 = AssociationParameters::new(&[comp1, comp3], &sigma, &binary, None);
let assoc1 = AssociationParameters::new(&[comp1.clone(), comp2], &[], None);
let assoc2 = AssociationParameters::new(&[comp1, comp3], &binary, None);
println!("{}", assoc1.epsilon_k_ab);
println!("{}", assoc2.epsilon_k_ab);
assert_eq!(assoc1.epsilon_k_ab, assoc2.epsilon_k_ab);
println!("{}", assoc1.sigma3_kappa_ab);
println!("{}", assoc2.sigma3_kappa_ab);
assert_eq!(assoc1.sigma3_kappa_ab, assoc2.sigma3_kappa_ab);
println!("{}", assoc1.kappa_ab);
println!("{}", assoc2.kappa_ab);
assert_eq!(assoc1.kappa_ab, assoc2.kappa_ab);
}
}

Expand All @@ -640,7 +642,7 @@ mod tests_pcsaft {
let n = 1.23;
let s = StateHD::new(t, v, arr1(&[n]));
let d = params.hs_diameter(t);
let a_rust = assoc.helmholtz_energy(&s, &d) / n;
let a_rust = assoc.helmholtz_energy(&s, &params.sigma, &d) / n;
assert_relative_eq!(a_rust, -4.229878997054543, epsilon = 1e-10);
}

Expand All @@ -653,7 +655,7 @@ mod tests_pcsaft {
let n = 1.23;
let s = StateHD::new(t, v, arr1(&[n]));
let d = params.hs_diameter(t);
let a_rust = assoc.helmholtz_energy(&s, &d) / n;
let a_rust = assoc.helmholtz_energy(&s, &params.sigma, &d) / n;
assert_relative_eq!(a_rust, -4.229878997054543, epsilon = 1e-10);
}

Expand All @@ -673,8 +675,8 @@ mod tests_pcsaft {
let n = 1.23;
let s = StateHD::new(t, v, arr1(&[n]));
let d = params.hs_diameter(t);
let a_assoc = assoc.helmholtz_energy(&s, &d) / n;
let a_cross_assoc = cross_assoc.helmholtz_energy(&s, &d) / n;
let a_assoc = assoc.helmholtz_energy(&s, &params.sigma, &d) / n;
let a_cross_assoc = cross_assoc.helmholtz_energy(&s, &params.sigma, &d) / n;
assert_relative_eq!(a_assoc, a_cross_assoc, epsilon = 1e-10);
Ok(())
}
Expand Down Expand Up @@ -704,8 +706,12 @@ mod tests_gc_pcsaft {
arr1(&[Dual64::from_re(moles)]),
);
let diameter = params.hs_diameter(state.temperature);
let pressure =
Pressure::from_reduced(-contrib.helmholtz_energy(&state, &diameter).eps * temperature);
let pressure = Pressure::from_reduced(
-contrib
.helmholtz_energy(&state, &params.sigma, &diameter)
.eps
* temperature,
);
assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10);
}

Expand All @@ -722,8 +728,12 @@ mod tests_gc_pcsaft {
arr1(&[Dual64::from_re(moles)]),
);
let diameter = params.hs_diameter(state.temperature);
let pressure =
Pressure::from_reduced(-contrib.helmholtz_energy(&state, &diameter).eps * temperature);
let pressure = Pressure::from_reduced(
-contrib
.helmholtz_energy(&state, &params.sigma, &diameter)
.eps
* temperature,
);
assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10);
}

Expand All @@ -740,8 +750,12 @@ mod tests_gc_pcsaft {
moles.mapv(Dual64::from_re),
);
let diameter = params.hs_diameter(state.temperature);
let pressure =
Pressure::from_reduced(-contrib.helmholtz_energy(&state, &diameter).eps * temperature);
let pressure = Pressure::from_reduced(
-contrib
.helmholtz_energy(&state, &params.sigma, &diameter)
.eps
* temperature,
);
assert_relative_eq!(pressure, -26.105606376765632 * PASCAL, max_relative = 1e-10);
}
}
26 changes: 19 additions & 7 deletions src/gc_pcsaft/eos/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,7 @@ impl Residual for GcPcSaft {
if let Some(association) = self.association.as_ref() {
v.push((
association.to_string(),
association.helmholtz_energy(state, &d),
association.helmholtz_energy(state, &self.parameters.sigma, &d),
))
}
v
Expand Down Expand Up @@ -202,8 +202,12 @@ mod test {
arr1(&[Dual64::from_re(moles)]),
);
let diameter = parameters.hs_diameter(state.temperature);
let pressure =
Pressure::from_reduced(-contrib.helmholtz_energy(&state, &diameter).eps * temperature);
let pressure = Pressure::from_reduced(
-contrib
.helmholtz_energy(&state, &parameters.sigma, &diameter)
.eps
* temperature,
);
assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10);
}

Expand All @@ -221,8 +225,12 @@ mod test {
arr1(&[Dual64::from_re(moles)]),
);
let diameter = parameters.hs_diameter(state.temperature);
let pressure =
Pressure::from_reduced(-contrib.helmholtz_energy(&state, &diameter).eps * temperature);
let pressure = Pressure::from_reduced(
-contrib
.helmholtz_energy(&state, &parameters.sigma, &diameter)
.eps
* temperature,
);
assert_relative_eq!(pressure, -3.6819598891967344 * PASCAL, max_relative = 1e-10);
}

Expand All @@ -239,8 +247,12 @@ mod test {
moles.mapv(Dual64::from_re),
);
let diameter = parameters.hs_diameter(state.temperature);
let pressure =
Pressure::from_reduced(-contrib.helmholtz_energy(&state, &diameter).eps * temperature);
let pressure = Pressure::from_reduced(
-contrib
.helmholtz_energy(&state, &parameters.sigma, &diameter)
.eps
* temperature,
);
assert_relative_eq!(pressure, -26.105606376765632 * PASCAL, max_relative = 1e-10);
}
}
2 changes: 1 addition & 1 deletion src/gc_pcsaft/eos/parameter.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -224,7 +224,7 @@ impl ParameterHetero for GcPcSaftEosParameters {
let sigma = Array1::from_vec(sigma);
let component_index = Array1::from_vec(component_index);
let association =
AssociationParameters::new(&association_records, &sigma, &[], Some(&component_index));
AssociationParameters::new(&association_records, &[], Some(&component_index));

Ok(Self {
molarweight,
Expand Down
6 changes: 3 additions & 3 deletions src/pcsaft/eos/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -173,7 +173,7 @@ impl Residual for PcSaft {
if let Some(association) = self.association.as_ref() {
v.push((
association.to_string(),
association.helmholtz_energy(state, &d),
association.helmholtz_energy(state, &self.parameters.sigma, &d),
))
}
v
Expand Down Expand Up @@ -449,7 +449,7 @@ mod tests {
let n = 1.23;
let s = StateHD::new(t, v, arr1(&[n]));
let d = parameters.hs_diameter(t);
let a_rust = assoc.helmholtz_energy(&s, &d) / n;
let a_rust = assoc.helmholtz_energy(&s, &parameters.sigma, &d) / n;
assert_relative_eq!(a_rust, -4.229878997054543, epsilon = 1e-10);
}

Expand All @@ -463,7 +463,7 @@ mod tests {
let n = 1.23;
let s = StateHD::new(t, v, arr1(&[n]));
let d = parameters.hs_diameter(t);
let a_rust = assoc.helmholtz_energy(&s, &d) / n;
let a_rust = assoc.helmholtz_energy(&s, &parameters.sigma, &d) / n;
assert_relative_eq!(a_rust, -4.229878997054543, epsilon = 1e-10);
}

Expand Down
2 changes: 1 addition & 1 deletion src/pcsaft/parameters.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -429,7 +429,7 @@ impl Parameter for PcSaftParameters {
})
.collect();
let association =
AssociationParameters::new(&association_records, &sigma, &binary_association, None);
AssociationParameters::new(&association_records, &binary_association, None);

let k_ij = binary_records.as_ref().map(|br| br.map(|br| br.k_ij));
let mut sigma_ij = Array::zeros((n, n));
Expand Down