started impl for yaeos interface

g-bauer · g-bauer · commit df3088554f0c · 2024-07-11T18:27:16.000+02:00
diff --git 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
@@ -1,6 +1,6 @@
 #include <stdio.h>
 #include <stdint.h>
-#include <feos.h>
+#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
@@ -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<feos_density_initialization_t> 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<feos_contributions_t> 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<EquationOfState<IdealGasModel, ResidualModel>>);
 
 #[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<PureRecord<PcSaftRecord>> =
-                serde_json::from_value(v["residual_substance_parameters"].clone()).unwrap();    
+                serde_json::from_value(v["residual_substance_parameters"].clone()).unwrap();
             let binary_records: Vec<BinaryRecord<Identifier, PcSaftBinaryRecord>> =
-                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 = <PcSaftParameters as Parameter>::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<EquationOfState<IdealGasModel, ResidualModel>>);
 
 #[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,64 +137,53 @@ 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)))
 }
 
 #[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
