adjust dft

feos-org · g-bauer · Feb 26, 2024 · Feb 12, 2024 · Feb 12, 2024 · Feb 13, 2024
commit f45378d830af40d9a1590f01c6634516dec73af0
diff --git a/feos-core/src/cubic.rs b/feos-core/src/cubic.rs
@@ -8,7 +8,7 @@ use crate::equation_of_state::{Components, Residual};
 use crate::parameter::{Identifier, Parameter, ParameterError, PureRecord};
 use crate::si::{MolarWeight, GRAM, MOL};
 use crate::state::StateHD;
-use ndarray::{Array1, Array2};
+use ndarray::{Array1, Array2, ScalarOperand};
 use num_dual::DualNum;
 use serde::{Deserialize, Serialize};
 use std::f64::consts::SQRT_2;
@@ -146,16 +146,6 @@ impl Parameter for PengRobinsonParameters {
     }
 }
 
-struct PengRobinsonContribution {
-    parameters: Arc<PengRobinsonParameters>,
-}
-
-impl fmt::Display for PengRobinsonContribution {
-    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
-        write!(f, "Peng Robinson")
-    }
-}
-
 /// A simple version of the Peng-Robinson equation of state.
 pub struct PengRobinson {
     /// Parameters
@@ -215,7 +205,7 @@ impl Residual for PengRobinson {
                 * ((v + b * n * (1.0 + SQRT_2)) / (v + b * n * (1.0 - SQRT_2))).ln())
     }
 
-    fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy>(
+    fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy + ScalarOperand>(
         &self,
         state: &StateHD<D>,
     ) -> Vec<(String, D)> {

diff --git a/feos-core/src/equation_of_state/ideal_gas.rs b/feos-core/src/equation_of_state/ideal_gas.rs
@@ -2,12 +2,13 @@ use super::Components;
 use crate::StateHD;
 use ndarray::Array1;
 use num_dual::DualNum;
-use std::fmt::Display;
 
 /// Ideal gas Helmholtz energy contribution.
-pub trait IdealGas: Components + Sync + Send + Display {
+pub trait IdealGas: Components + Sync + Send {
     fn ln_lambda3<D: DualNum<f64> + Copy>(&self, temperature: D) -> Array1<D>;
 
+    fn ideal_gas_model(&self) -> String;
+
     /// Evaluate the ideal gas Helmholtz energy contribution for a given state.
     ///
     /// In some cases it could be advantageous to overwrite this

diff --git a/feos-core/src/equation_of_state/mod.rs b/feos-core/src/equation_of_state/mod.rs
@@ -2,8 +2,8 @@ use crate::{
     si::{Diffusivity, MolarWeight, Moles, Temperature, ThermalConductivity, Viscosity, Volume},
     EosResult,
 };
-use ndarray::Array1;
-use std::{fmt::Display, sync::Arc};
+use ndarray::{Array1, ScalarOperand};
+use std::sync::Arc;
 
 mod ideal_gas;
 mod residual;
@@ -52,12 +52,6 @@ impl<I: IdealGas> EquationOfState<I, NoResidual> {
     }
 }
 
-impl<I: IdealGas, R> Display for EquationOfState<I, R> {
-    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-        write!(f, "{}", self.ideal_gas)
-    }
-}
-
 impl<I: Components, R: Components> Components for EquationOfState<I, R> {
     fn components(&self) -> usize {
         assert_eq!(
@@ -80,21 +74,18 @@ impl<I: IdealGas, R: Components + Sync + Send> IdealGas for EquationOfState<I, R
     fn ln_lambda3<D: num_dual::DualNum<f64> + Copy>(&self, temperature: D) -> Array1<D> {
         self.ideal_gas.ln_lambda3(temperature)
     }
+
+    fn ideal_gas_model(&self) -> String {
+        self.ideal_gas.ideal_gas_model()
+    }
 }
 
 impl<I: IdealGas, R: Residual> Residual for EquationOfState<I, R> {
     fn compute_max_density(&self, moles: &Array1<f64>) -> f64 {
         self.residual.compute_max_density(moles)
     }
 
-    fn residual_helmholtz_energy<D: num_dual::DualNum<f64> + Copy>(
-        &self,
-        state: &crate::StateHD<D>,
-    ) -> D {
-        self.residual.residual_helmholtz_energy(state)
-    }
-
-    fn residual_helmholtz_energy_contributions<D: num_dual::DualNum<f64> + Copy>(
+    fn residual_helmholtz_energy_contributions<D: num_dual::DualNum<f64> + Copy + ScalarOperand>(
         &self,
         state: &crate::StateHD<D>,
     ) -> Vec<(String, D)> {

diff --git a/feos-core/src/equation_of_state/residual.rs b/feos-core/src/equation_of_state/residual.rs
@@ -3,6 +3,7 @@ use crate::si::*;
 use crate::StateHD;
 use crate::{EosError, EosResult};
 use ndarray::prelude::*;
+use ndarray::ScalarOperand;
 use num_dual::*;
 use num_traits::{One, Zero};
 use std::ops::Div;
@@ -24,13 +25,16 @@ pub trait Residual: Components + Send + Sync {
 
     /// Evaluate the reduced Helmholtz energy of each individual contribution
     /// and return them together with a string representation of the contribution.
-    fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy>(
+    fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy + ScalarOperand>(
         &self,
         state: &StateHD<D>,
     ) -> Vec<(String, D)>;
 
     /// Evaluate the residual reduced Helmholtz energy $\beta A^\mathrm{res}$.
-    fn residual_helmholtz_energy<D: DualNum<f64> + Copy>(&self, state: &StateHD<D>) -> D {
+    fn residual_helmholtz_energy<D: DualNum<f64> + Copy + ScalarOperand>(
+        &self,
+        state: &StateHD<D>,
+    ) -> D {
         self.residual_helmholtz_energy_contributions(state)
             .iter()
             .fold(D::zero(), |acc, (_, a)| acc + a)

diff --git a/feos-core/src/lib.rs b/feos-core/src/lib.rs
@@ -126,18 +126,11 @@ mod tests {
     use crate::EosResult;
     use crate::StateBuilder;
     use approx::*;
-    use std::fmt::Display;
     use std::sync::Arc;
 
     // Only to be able to instantiate an `EquationOfState`
     struct NoIdealGas;
 
-    impl Display for NoIdealGas {
-        fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-            write!(f, "NoIdealGas")
-        }
-    }
-
     impl Components for NoIdealGas {
         fn components(&self) -> usize {
             1
@@ -155,6 +148,10 @@ mod tests {
         ) -> ndarray::prelude::Array1<D> {
             unreachable!()
         }
+
+        fn ideal_gas_model(&self) -> String {
+            "NoIdealGas".into()
+        }
     }
 
     fn pure_record_vec() -> Vec<PureRecord<PengRobinsonRecord>> {

diff --git a/feos-core/src/parameter/mod.rs b/feos-core/src/parameter/mod.rs
@@ -84,7 +84,7 @@ where
     /// `pure_records`, the `Default` implementation of Self::Binary is used.
     #[allow(clippy::expect_fun_call)]
     fn binary_matrix_from_records(
-        pure_records: &Vec<PureRecord<Self::Pure>>,
+        pure_records: &[PureRecord<Self::Pure>],
         binary_records: &[BinaryRecord<Identifier, Self::Binary>],
         identifier_option: IdentifierOption,
     ) -> Option<Array2<Self::Binary>> {

diff --git a/feos-core/src/state/properties.rs b/feos-core/src/state/properties.rs
@@ -257,7 +257,7 @@ impl<E: Residual + IdealGas> State<E> {
         let contributions = self.eos.residual_helmholtz_energy_contributions(&new_state);
         let mut res = Vec::with_capacity(contributions.len() + 1);
         res.push((
-            self.eos.to_string(),
+            self.eos.ideal_gas_model(),
             MolarEnergy::from_reduced(
                 (self.eos.ideal_gas_helmholtz_energy(&new_state) * new_state.temperature).eps,
             ),

diff --git a/feos-core/tests/parameters.rs b/feos-core/tests/parameters.rs
@@ -82,7 +82,7 @@ fn from_records() -> Result<(), ParameterError> {
             }
         ]
         "#;
-    let pure_records = serde_json::from_str(pr_json).expect("Unable to parse json.");
+    let pure_records: Vec<_> = serde_json::from_str(pr_json).expect("Unable to parse json.");
     let binary_records: Vec<_> = serde_json::from_str(br_json).expect("Unable to parse json.");
     let binary_matrix = MyParameter::binary_matrix_from_records(
         &pure_records,
@@ -201,7 +201,7 @@ fn from_records_missing_binary() -> Result<(), ParameterError> {
             }
         ]
         "#;
-    let pure_records = serde_json::from_str(pr_json).expect("Unable to parse json.");
+    let pure_records: Vec<_> = serde_json::from_str(pr_json).expect("Unable to parse json.");
     let binary_records: Vec<_> = serde_json::from_str(br_json).expect("Unable to parse json.");
     let binary_matrix = MyParameter::binary_matrix_from_records(
         &pure_records,
@@ -266,7 +266,7 @@ fn from_records_correct_binary_order() -> Result<(), ParameterError> {
             }
         ]
         "#;
-    let pure_records = serde_json::from_str(pr_json).expect("Unable to parse json.");
+    let pure_records: Vec<_> = serde_json::from_str(pr_json).expect("Unable to parse json.");
     let binary_records: Vec<_> = serde_json::from_str(br_json).expect("Unable to parse json.");
     let binary_matrix = MyParameter::binary_matrix_from_records(
         &pure_records,

diff --git a/feos-derive/src/dft.rs b/feos-derive/src/dft.rs
@@ -95,7 +95,7 @@ fn impl_helmholtz_energy_functional(
     let contributions = variants.iter().map(|v| {
         let name = &v.ident;
         quote! {
-            Self::#name(functional) => functional.contributions()
+            Self::#name(functional) => Box::new(functional.contributions().map(FunctionalContributionVariant::from))
         }
     });
 
@@ -121,6 +121,7 @@ fn impl_helmholtz_energy_functional(
 
     Ok(quote! {
         impl HelmholtzEnergyFunctional for FunctionalVariant {
+            type Contribution = FunctionalContributionVariant;
             fn molecule_shape(&self) -> MoleculeShape {
                 match self {
                     #(#molecule_shape,)*
@@ -131,7 +132,7 @@ fn impl_helmholtz_energy_functional(
                     #(#compute_max_density,)*
                 }
             }
-            fn contributions(&self) -> &[Box<dyn FunctionalContribution>] {
+            fn contributions(&self) -> Box<dyn Iterator<Item = FunctionalContributionVariant>> {
                 match self {
                     #(#contributions,)*
                 }

diff --git a/feos-derive/src/functional_contribution.rs b/feos-derive/src/functional_contribution.rs
@@ -0,0 +1,115 @@
+use quote::quote;
+use syn::{DeriveInput, Ident};
+
+pub(crate) fn expand_functional_contribution(
+    input: DeriveInput,
+) -> syn::Result<proc_macro2::TokenStream> {
+    let ident = input.ident;
+    let variants = match input.data {
+        syn::Data::Enum(syn::DataEnum { ref variants, .. }) => variants,
+        _ => panic!("this derive macro only works on enums"),
+    };
+
+    let functional_contribution = impl_functional_contribution(&ident, variants);
+    let display = impl_display(&ident, variants);
+    let from = impl_from(&ident, variants);
+    Ok(quote! {
+        #functional_contribution
+        #display
+        #from
+    })
+}
+
+fn impl_functional_contribution(
+    ident: &Ident,
+    variants: &syn::punctuated::Punctuated<syn::Variant, syn::token::Comma>,
+) -> proc_macro2::TokenStream {
+    let weight_functions = variants.iter().map(|v| {
+        let name = &v.ident;
+        quote! {
+            Self::#name(functional_contribution) => functional_contribution.weight_functions(temperature)
+        }
+    });
+    let weight_functions_pdgt = variants.iter().map(|v| {
+        let name = &v.ident;
+        quote! {
+            Self::#name(functional_contribution) => functional_contribution.weight_functions_pdgt(temperature)
+        }
+    });
+    let helmholtz_energy_density = variants.iter().map(|v| {
+        let name = &v.ident;
+        quote! {
+            Self::#name(functional_contribution) => functional_contribution.helmholtz_energy_density(temperature, weighted_densities)
+        }
+    });
+
+    quote! {
+        impl FunctionalContribution for #ident {
+            fn weight_functions<N: DualNum<f64> + Copy+ScalarOperand>(&self, temperature: N) -> WeightFunctionInfo<N> {
+                match self {
+                    #(#weight_functions,)*
+                }
+            }
+            fn weight_functions_pdgt<N: DualNum<f64> + Copy+ScalarOperand>(&self, temperature: N) -> WeightFunctionInfo<N> {
+                match self {
+                    #(#weight_functions_pdgt,)*
+                }
+            }
+            fn helmholtz_energy_density<N: DualNum<f64> + Copy+ScalarOperand>(
+                &self,
+                temperature: N,
+                weighted_densities: ArrayView2<N>,
+            ) -> EosResult<Array1<N>> {
+                match self {
+                    #(#helmholtz_energy_density,)*
+                }
+            }
+        }
+    }
+}
+
+fn impl_display(
+    ident: &Ident,
+    variants: &syn::punctuated::Punctuated<syn::Variant, syn::token::Comma>,
+) -> proc_macro2::TokenStream {
+    let fmt = variants.iter().map(|v| {
+        let name = &v.ident;
+        quote! {
+            Self::#name(functional_contribution) => functional_contribution.fmt(f)
+        }
+    });
+
+    quote! {
+        impl std::fmt::Display for #ident {
+            fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+                match self {
+                    #(#fmt,)*
+                }
+            }
+        }
+    }
+}
+
+fn impl_from(
+    ident: &Ident,
+    variants: &syn::punctuated::Punctuated<syn::Variant, syn::token::Comma>,
+) -> proc_macro2::TokenStream {
+    let from = variants.iter().map(|v| {
+        let name = &v.ident;
+        let syn::Fields::Unnamed(syn::FieldsUnnamed { unnamed, .. }) = &v.fields else {
+            panic!("All variants must be tuple structs!")
+        };
+        let inner = &unnamed.first().unwrap().ty;
+        quote! {
+            impl From<#inner> for #ident {
+                fn from(variant: #inner) -> Self {
+                    Self::#name(variant)
+                }
+            }
+        }
+    });
+
+    quote! {
+        #(#from)*
+    }
+}