Added ideal gas implementation to user_defined module

feos-org · g-bauer · Feb 26, 2024 · Feb 12, 2024 · Feb 12, 2024 · Feb 13, 2024
commit 2a793e6e83b27d3726dbf358060e408642d5e9b7
diff --git a/feos-core/src/python/user_defined.rs b/feos-core/src/python/user_defined.rs
@@ -11,8 +11,6 @@ use std::any::Any;
 use std::convert::TryInto;
 use std::fmt;
 
-struct PyHelmholtzEnergy(Py<PyAny>);
-
 pub struct PyIdealGas(Py<PyAny>);
 
 impl PyIdealGas {
@@ -61,13 +59,45 @@ impl Components for PyIdealGas {
     }
 }
 
-impl IdealGas for PyIdealGas {
-    fn ideal_gas_name(&self) -> String {
-        unimplemented!()
-    }
-    fn ln_lambda3<D: DualNum<f64> + Copy>(&self, temperature: D) -> Array1<D> {
-        unimplemented!()
-    }
+macro_rules! impl_ideal_gas {
+    ($($py_hd_id:ident, $hd_ty:ty);*) => {
+        impl IdealGas for PyIdealGas {
+            fn ideal_gas_name(&self) -> String {
+                "Python".to_string()
+            }
+
+            fn ln_lambda3<D: DualNum<f64> + Copy>(&self, temperature: D) -> Array1<D> {
+                let mut result = Array1::from_elem((self.components(),), D::zero());
+
+                $(
+                    if let Some(t) = (&temperature as &dyn Any).downcast_ref::<$hd_ty>() {
+                        let l3_any = (&mut result as &mut dyn Any).downcast_mut::<Array1<$hd_ty>>().unwrap();
+                        *l3_any = Python::with_gil(|py| {
+                            let py_result = self
+                                .0
+                                .as_ref(py)
+                                .call_method1("ln_lambda3", (<$py_hd_id>::from(t.clone()),))
+                                .unwrap();
+
+                            // f64
+                            if let Ok(r) = py_result.extract::<PyReadonlyArray1<f64>>() {
+                                r.to_owned_array()
+                                    .mapv(|ri| <$hd_ty>::from(ri))
+                            // anything but f64
+                            } else if let Ok(r) = py_result.extract::<PyReadonlyArray1<PyObject>>() {
+                                r.to_owned_array()
+                                    .mapv(|ri| <$hd_ty>::from(ri.extract::<$py_hd_id>(py).unwrap()))
+                            } else {
+                                    panic!("ln_lambda3: data type of result must be one-dimensional numpy ndarray")
+                            }
+                        });
+                        return result
+                    }
+                )*
+                panic!("ln_lambda3: input data type not understood")
+            }
+        }
+    };
 }
 
 impl fmt::Display for PyIdealGas {
@@ -141,7 +171,7 @@ impl Components for PyResidual {
     }
 }
 
-macro_rules! residual {
+macro_rules! impl_residual {
     ($($py_state_id:ident, $py_hd_id:ident, $hd_ty:ty);*) => {
         impl Residual for PyResidual {
             fn compute_max_density(&self, moles: &Array1<f64>) -> f64 {
@@ -173,28 +203,14 @@ macro_rules! residual {
                         return a
                     }
                 )*
-
-                // // Dual64
-                // if let Some(s) = (state as &dyn Any).downcast_ref::<StateHD<Dual64>>() {
-                //     let d = (&mut a as &mut dyn Any).downcast_mut::<Dual64>().unwrap();
-                //     *d = Python::with_gil(|py| {
-                //         let py_result = self
-                //             .obj
-                //             .as_ref(py)
-                //             .call_method1("helmholtz_energy", (<PyStateD>::from(s.clone()),))
-                //             .unwrap();
-                //         <Dual64>::from(py_result.extract::<PyDual64>().unwrap())
-                //     });
-                //     return a
-                // }
-                panic!("Something went wrong!.")
+                panic!("helmholtz_energy: input data type not understood")
             }
 
             fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy>(
                     &self,
                     state: &StateHD<D>,
                 ) -> Vec<(String, D)> {
-                unimplemented!()
+                vec![("Python".to_string(), self.residual_helmholtz_energy(state))]
             }
 
             fn molar_weight(&self) -> MolarWeight<Array1<f64>> {
@@ -217,12 +233,6 @@ macro_rules! residual {
     }
 }
 
-impl fmt::Display for PyHelmholtzEnergy {
-    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
-        write!(f, "Custom")
-    }
-}
-
 macro_rules! state {
     ($py_state_id:ident, $py_hd_id:ident, $hd_ty:ty) => {
         #[pyclass]
@@ -332,14 +342,11 @@ macro_rules! dual_number {
 //             }
 //         }
 //     };
-// }
 
 macro_rules! impl_dual_state_helmholtz_energy {
     ($py_state_id:ident, $py_hd_id:ident, $hd_ty:ty, $py_field_ty:ty) => {
         dual_number!($py_hd_id, $hd_ty, $py_field_ty);
         state!($py_state_id, $py_hd_id, $hd_ty);
-        // helmholtz_energy!($py_state_id, $py_hd_id, $hd_ty);
-        // de_broglie_wavelength!($py_hd_id, $hd_ty);
     };
 }
 
@@ -399,7 +406,30 @@ impl_dual_state_helmholtz_energy!(
     PyDualVec3
 );
 
-residual!(
+impl_ideal_gas!(
+    f64, f64;
+    PyDual64, Dual64;
+    PyDualDualVec3,
+    Dual<DualSVec64<3>, f64>;
+    PyHyperDual64, HyperDual64;
+     PyDual2_64, Dual2_64;
+     PyDual3_64, Dual3_64;
+     PyHyperDualDual64, HyperDual<Dual64, f64>;
+    PyHyperDualVec2,
+    HyperDual<DualSVec64<2>, f64>;
+    PyHyperDualVec3,
+    HyperDual<DualSVec64<3>, f64>;
+    PyDual2Dual64,
+    Dual2<Dual64, f64>;
+    PyDual3Dual64,
+    Dual3<Dual64, f64>;
+    PyDual3DualVec2,
+    Dual3<DualSVec64<2>, f64>;
+    PyDual3DualVec3,
+    Dual3<DualSVec64<3>, f64>
+);
+
+impl_residual!(
     PyStateF, f64, f64;
     PyStateD, PyDual64, Dual64;
     PyStateDualDualVec3,
@@ -427,4 +457,4 @@ residual!(
     PyStateD3DVec3,
     PyDual3DualVec3,
     Dual3<DualSVec64<3>, f64>
-);
+);