From df3088554f0c1cae5b5b0b25c60fdef56462dc68 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Thu, 11 Jul 2024 18:27:16 +0200 Subject: [PATCH] started impl for yaeos interface --- Cargo.toml | 1 + src/c/test.c | 13 +++-- src/lib.rs | 109 ++++++++++++++++++++++++---------- src/residual.rs | 151 ++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 237 insertions(+), 37 deletions(-) create mode 100644 src/residual.rs diff --git a/Cargo.toml b/Cargo.toml index 1531072..7e25033 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -16,6 +16,7 @@ serde_json = "1.0" libc = "0.2" ndarray = "0.15" num-dual = "0.8" +typenum = "1.17.0" [build-dependencies] cbindgen = { version = "0.26", default-features = false } diff --git a/src/c/test.c b/src/c/test.c index 5465e6c..9ba4faa 100644 --- a/src/c/test.c +++ b/src/c/test.c @@ -1,6 +1,6 @@ #include #include -#include +#include "feos.h" // typedef struct feos_eos feos_eos_t; // extern feos_eos_t *feos_eos_from_json(const char *str); // extern void feos_eos_free(feos_eos_t *); @@ -109,25 +109,28 @@ int main(void) // feos_equation_of_state_t *eos = feos_eos_from_json(str2); // state - double temperature = 200.0; // K - double pressure = 15.0; // bar + double temperature = 300.0; // K + double pressure = 1.0; // bar double moles[] = {1.0}; // mol for each component size_t n = sizeof(moles) / sizeof(size_t); // number of components + feos_density_initialization_t density_init = {initial_density, 0.01}; // {tag, density in mol/meter³} + feos_state_t *state = feos_state_new_npt( eos, // equation of state temperature, // temperature in K pressure, // pressure in bar moles, // number of mols n, // number of components - "stable" // which phase to compute (stable, liquid, vapor) + density_init // which density to use as starting point ); - double pressure_calculated = feos_state_pressure(state, 2); + double pressure_calculated = feos_state_pressure(state, Total); double density = feos_state_density(state); printf("State stable? %s\n", feos_state_is_stable(state) ? "true" : "false"); printf("pressure: %f bar\n", pressure); + printf("pressure_calculated: %f bar\n", pressure_calculated); printf("density: %f mol/m3\n", density); feos_state_free(state); diff --git a/src/lib.rs b/src/lib.rs index 9f66879..394dc4e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,23 +1,78 @@ use feos::pcsaft::{PcSaft, PcSaftBinaryRecord, PcSaftParameters, PcSaftRecord}; use feos::{ideal_gas::IdealGasModel, ResidualModel}; use feos_core::parameter::{BinaryRecord, Identifier, IdentifierOption, Parameter, PureRecord}; -use feos_core::si::{BAR, KELVIN, KILOGRAM, METER, MOL}; -use feos_core::{Contributions, DensityInitialization, EquationOfState, SolverOptions, State}; +use feos_core::si::{Pressure, ANGSTROM, BAR, KELVIN, KILOGRAM, METER, MOL, PASCAL}; +use feos_core::{ + Contributions, DensityInitialization, Derivative, EquationOfState, SolverOptions, State, StateHD +}; use libc::{c_char, size_t}; use ndarray::Array1; +use num_dual::HyperDual64; use serde_json::Value; use std::ffi::CStr; use std::slice; use std::sync::Arc; +use typenum::P3; pub mod pcsaft; +pub mod residual; + +#[repr(C)] +#[derive(Clone, Copy, Debug)] +#[allow(non_camel_case_types)] +pub enum feos_density_initialization_t { + /// Calculate a vapor phase by initializing using the ideal gas. + vapor, + /// Calculate a liquid phase by using the `max_density`. + liquid, + /// Use the given density as initial value. + initial_density(f64), + /// Calculate the most stable phase by calculating both a vapor and a liquid + /// and return the one with the lower molar Gibbs energy. + none, +} + +impl From for DensityInitialization { + fn from(value: feos_density_initialization_t) -> Self { + match value { + feos_density_initialization_t::vapor => Self::Vapor, + feos_density_initialization_t::liquid => Self::Liquid, + feos_density_initialization_t::initial_density(rho) => { + Self::InitialDensity(rho.clone() * MOL / (METER * METER * METER)) + } + feos_density_initialization_t::none => Self::None, + } + } +} + +/// Possible contributions that can be computed. +#[repr(C)] +#[derive(Clone, Copy)] +pub enum feos_contributions_t { + /// Only compute the ideal gas contribution + IdealGas, + /// Only compute the difference between the total and the ideal gas contribution + Residual, + /// Compute ideal gas and residual contributions + Total, +} + +impl From for Contributions { + fn from(value: feos_contributions_t) -> Self { + match value { + feos_contributions_t::IdealGas => Contributions::IdealGas, + feos_contributions_t::Residual => Contributions::Residual, + feos_contributions_t::Total => Contributions::Total, + } + } +} #[allow(non_camel_case_types)] pub struct feos_equation_of_state_t(Arc>); #[no_mangle] #[allow(non_camel_case_types)] -pub extern fn feos_eos_from_json(json: *const c_char) -> *mut feos_equation_of_state_t { +pub extern "C" fn feos_eos_from_json(json: *const c_char) -> *mut feos_equation_of_state_t { let c_str = unsafe { assert!(!json.is_null()); CStr::from_ptr(json) @@ -34,9 +89,10 @@ pub extern fn feos_eos_from_json(json: *const c_char) -> *mut feos_equation_of_s match residual_model_name.to_lowercase().as_str() { "pc-saft" | "pcsaft" => { let pure_records: Vec> = - serde_json::from_value(v["residual_substance_parameters"].clone()).unwrap(); + serde_json::from_value(v["residual_substance_parameters"].clone()).unwrap(); let binary_records: Vec> = - serde_json::from_value(v["residual_binary_parameters"].clone()).unwrap_or(Vec::new()); + serde_json::from_value(v["residual_binary_parameters"].clone()) + .unwrap_or(Vec::new()); let binary_matrix = ::binary_matrix_from_records( &pure_records, &binary_records, @@ -57,7 +113,7 @@ pub extern fn feos_eos_from_json(json: *const c_char) -> *mut feos_equation_of_s #[no_mangle] #[allow(non_camel_case_types)] -pub extern fn feos_eos_free(ptr: *mut feos_equation_of_state_t) { +pub extern "C" fn feos_eos_free(ptr: *mut feos_equation_of_state_t) { if !ptr.is_null() { unsafe { drop(Box::from_raw(ptr)); @@ -71,7 +127,7 @@ pub struct feos_state_t(State>); #[no_mangle] #[allow(non_camel_case_types)] -pub extern fn feos_state_free(ptr: *mut feos_state_t) { +pub extern "C" fn feos_state_free(ptr: *mut feos_state_t) { if !ptr.is_null() { unsafe { drop(Box::from_raw(ptr)); @@ -81,38 +137,32 @@ pub extern fn feos_state_free(ptr: *mut feos_state_t) { #[no_mangle] #[allow(non_camel_case_types)] -pub extern fn feos_state_new_npt( +pub extern "C" fn feos_state_new_npt( eos_ptr: *const feos_equation_of_state_t, t_k: f64, p_bar: f64, moles: *const f64, len: size_t, - phase: *const c_char, + density_initialization: feos_density_initialization_t, ) -> *mut feos_state_t { let eos = unsafe { assert!(!eos_ptr.is_null()); &*eos_ptr }; - let phase_str = unsafe { - assert!(!phase.is_null()); - CStr::from_ptr(phase) - }; + // let initialization = unsafe { + // assert!(!density_initialization.is_null()); + // &*density_initialization + // }; let n = unsafe { assert!(!moles.is_null()); Array1::from_iter(slice::from_raw_parts(moles, len).iter().cloned()) * MOL }; - - let density_initialization = match phase_str.to_str().unwrap() { - "liquid" => DensityInitialization::Liquid, - "vapor" => DensityInitialization::Vapor, - _ => DensityInitialization::None, - }; let state = State::new_npt( &eos.0, t_k * KELVIN, p_bar * BAR, &n, - density_initialization, + density_initialization.into(), ) .unwrap(); Box::into_raw(Box::new(feos_state_t(state))) @@ -120,25 +170,20 @@ pub extern fn feos_state_new_npt( #[no_mangle] #[allow(non_camel_case_types)] -pub extern fn feos_state_pressure( +pub extern "C" fn feos_state_pressure( state_ptr: *const feos_state_t, - contributions: size_t, + contributions: feos_contributions_t, ) -> f64 { let state = unsafe { assert!(!state_ptr.is_null()); &*state_ptr }; - let c = match contributions { - 0 => Contributions::IdealGas, - 1 => Contributions::Residual, - _ => Contributions::Total, - }; - state.0.pressure(c).convert_into(BAR) + state.0.pressure(contributions.into()).convert_into(BAR) } #[no_mangle] #[allow(non_camel_case_types)] -pub extern fn feos_state_density(state_ptr: *const feos_state_t) -> f64 { +pub extern "C" fn feos_state_density(state_ptr: *const feos_state_t) -> f64 { let state = unsafe { assert!(!state_ptr.is_null()); &*state_ptr @@ -147,12 +192,12 @@ pub extern fn feos_state_density(state_ptr: *const feos_state_t) -> f64 { } /// Mass density of the state. -/// +/// /// @params state_ptr Pointer to the state. /// @returns mass density in units of kg/m3 #[no_mangle] #[allow(non_camel_case_types)] -pub extern fn feos_state_mass_density(state_ptr: *const feos_state_t) -> f64 { +pub extern "C" fn feos_state_mass_density(state_ptr: *const feos_state_t) -> f64 { let state = unsafe { assert!(!state_ptr.is_null()); &*state_ptr @@ -170,7 +215,7 @@ pub extern fn feos_state_mass_density(state_ptr: *const feos_state_t) -> f64 { /// @returns True if the state is stable, false otherwise. #[no_mangle] #[allow(non_camel_case_types)] -pub extern fn feos_state_is_stable(state_ptr: *const feos_state_t) -> bool { +pub extern "C" fn feos_state_is_stable(state_ptr: *const feos_state_t) -> bool { let state = unsafe { assert!(!state_ptr.is_null()); &*state_ptr diff --git a/src/residual.rs b/src/residual.rs new file mode 100644 index 0000000..7d1604c --- /dev/null +++ b/src/residual.rs @@ -0,0 +1,151 @@ +use feos::pcsaft::{PcSaft, PcSaftBinaryRecord, PcSaftParameters, PcSaftRecord}; +use feos::ResidualModel; +use feos_core::parameter::{BinaryRecord, Identifier, IdentifierOption, Parameter, PureRecord}; +use feos_core::si::{ANGSTROM, KELVIN, MOL}; +use feos_core::{Derivative, Residual, State}; +use libc::{c_char, size_t}; +use ndarray::{Array1, ArrayViewMut}; +use num_dual::Dual2_64; +use serde_json::Value; +use std::ffi::CStr; +use std::slice; +use std::sync::Arc; +use typenum::P3; + +#[allow(non_camel_case_types)] +pub struct feos_residual_t(Arc); + +#[no_mangle] +#[allow(non_camel_case_types)] +pub extern "C" fn feos_residual_from_json(json: *const c_char) -> *mut feos_residual_t { + let c_str = unsafe { + assert!(!json.is_null()); + CStr::from_ptr(json) + }; + let r_str = c_str.to_str().unwrap(); + + let v: Value = serde_json::from_str(r_str).unwrap(); + let residual_model_name = match &v["residual_model"] { + Value::Null => return std::ptr::null_mut(), + Value::String(s) => s, + _ => return std::ptr::null_mut(), + }; + + match residual_model_name.to_lowercase().as_str() { + "pc-saft" | "pcsaft" => { + let pure_records: Vec> = + serde_json::from_value(v["residual_substance_parameters"].clone()).unwrap(); + let binary_records: Vec> = + serde_json::from_value(v["residual_binary_parameters"].clone()) + .unwrap_or(Vec::new()); + let binary_matrix = ::binary_matrix_from_records( + &pure_records, + &binary_records, + IdentifierOption::Name, + ); + let parameters = + Arc::new(PcSaftParameters::from_records(pure_records, binary_matrix).unwrap()); + Box::into_raw(Box::new(feos_residual_t(Arc::new(ResidualModel::PcSaft( + PcSaft::new(parameters), + ))))) + } + _ => std::ptr::null_mut(), + } +} + +#[no_mangle] +#[allow(non_camel_case_types)] +pub extern "C" fn feos_residual_free(ptr: *mut feos_residual_t) { + if !ptr.is_null() { + unsafe { + drop(Box::from_raw(ptr)); + } + } +} + +#[no_mangle] +#[allow(non_camel_case_types)] +pub extern "C" fn feos_residual_derivatives( + eos_ptr: *mut feos_residual_t, + n: *const f64, + v: f64, + t: f64, + ar: &mut f64, + arv: &mut f64, + art: &mut f64, + art2: &mut f64, + artv: &mut f64, + arv2: &mut f64, + arn: *mut f64, + arvn: *mut f64, + artn: *mut f64, + arn2: *mut f64, + len: size_t, +) { + let eos = unsafe { + assert!(!eos_ptr.is_null()); + &*eos_ptr + }; + let moles = unsafe { + assert!(!n.is_null()); + Array1::from_iter(slice::from_raw_parts(n, len).iter().cloned()) * MOL + }; + let mut arn_ = unsafe { + assert!(!arn.is_null()); + ArrayViewMut::from_shape_ptr(len, arn) + }; + let mut arvn_ = unsafe { + assert!(!arn.is_null()); + ArrayViewMut::from_shape_ptr(len, arvn) + }; + let mut artn_ = unsafe { + assert!(!arn.is_null()); + ArrayViewMut::from_shape_ptr(len, artn) + }; + let mut arnn_ = unsafe { + assert!(!arn.is_null()); + ArrayViewMut::from_shape_ptr((len, len), arn2) + }; + + let state = State::new_nvt(&eos.0, t * KELVIN, v * ANGSTROM.powi::(), &moles).unwrap(); + + // Derivatives + let new_state = state.derive2_mixed(Derivative::DT, Derivative::DV); + let a = &eos.0.evaluate_residual(&new_state) * new_state.temperature; + *ar = a.re; + *art = a.eps1; + *arv = a.eps2; + *artv = a.eps1eps2; + + let new_state = state.derive2(Derivative::DV); + let a = &eos.0.evaluate_residual(&new_state) * new_state.temperature; + *arv2 = a.v2; + + let new_state = state.derive2(Derivative::DV); + let a = &eos.0.evaluate_residual(&new_state) * new_state.temperature; + *art2 = a.v2; + + for i in 0..len { + for j in i + 1..len { + let new_state = state.derive2_mixed(Derivative::DN(i), Derivative::DN(j)); + let a = &eos.0.evaluate_residual(&new_state) * new_state.temperature; + arnn_[[i, j]] = a.eps1eps2; + arnn_[[j, i]] = a.eps1eps2; + } + } + for i in 0..len { + let new_state = state.derive2(Derivative::DN(i)); + let Dual2_64 { re: _, v1, v2, .. } = + &eos.0.evaluate_residual(&new_state) * new_state.temperature; + arn_[i] = v1; + arnn_[[i, i]] = v2; + + let new_state = state.derive2_mixed(Derivative::DN(i), Derivative::DV); + let a = &eos.0.evaluate_residual(&new_state) * new_state.temperature; + arvn_[i] = a.eps1eps2; + + let new_state = state.derive2_mixed(Derivative::DN(i), Derivative::DT); + let a = &eos.0.evaluate_residual(&new_state) * new_state.temperature; + artn_[i] = a.eps1eps2; + } +}