merged with origional upstream and solved differences

feos-org · JelleLagerweij · May 16, 2025 · May 16, 2025 · May 16, 2025 · May 17, 2025
commit 07bf0238eea15212bd9f2150b750bd5e70ef5b00
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
@@ -15,7 +15,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        crate: [feos-core, feos-dft]
+        crate: [feos-core, feos-ad, feos-dft]
     steps:
       - uses: actions/checkout@v4
       - name: Build
@@ -50,3 +50,15 @@ jobs:
         run: cargo build --release --features "${{ matrix.model }} dft"
       - name: Run tests
         run: cargo test --release --features "${{ matrix.model }} dft"
+
+  test_models_ad:
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        model: [pcsaft, gc_pcsaft]
+
+    steps:
+      - uses: actions/checkout@v4
+      - name: Run tests
+        run: cargo test --release -p feos-ad --features ${{ matrix.model }}
diff --git a/Cargo.toml b/Cargo.toml
@@ -43,6 +43,7 @@ libm = "0.2"
 gauss-quad = "0.2"
 approx = "0.5"
 criterion = "0.5"
+arrayvec = "0.7"
 
 feos-core = { version = "0.8", path = "crates/feos-core" }
 feos-dft = { version = "0.8", path = "crates/feos-dft" }

diff --git a/crates/feos-ad/CHANGELOG.md b/crates/feos-ad/CHANGELOG.md
@@ -0,0 +1,34 @@
+# Changelog
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
+and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+
+## [Unreleased]
+
+## [0.2.3] - 2025-05-28
+### Changed
+- Simplify the `Eigen` trait for better readable trait bounds. [#5](https://github.com/feos-org/feos-ad/pull/5)
+
+## [0.2.2] - 2025-05-28
+### Fixed
+- Export `Eigen` to be able to calculate critical points for pure components and binary mixtures generically. [#4](https://github.com/feos-org/feos-ad/pull/4)
+
+## [0.2.1] - 2025-04-14
+### Added
+- Added `StateAD::molar_isochoric_heat_capacity` and `StateAD::molar_isobaric_heat_capacity`. [#3](https://github.com/feos-org/feos-ad/pull/3)
+
+## [0.2.0] - 2025-01-27
+### Changed
+- Made `PcSaftBinary` generic for associating/non-associating systems. [#2](https://github.com/feos-org/feos-ad/pull/2)
+
+### Removed
+- Removed `ChemicalRecord`. [#2](https://github.com/feos-org/feos-ad/pull/2)
+
+## [0.1.1] - 2025-01-21
+### Added
+- Implemented PC-SAFT for a binary mixture. [#1](https://github.com/feos-org/feos-ad/pull/1)
+
+## [0.1.0] - 2025-01-08
+### Added
+- Initial release
diff --git a/crates/feos-ad/Cargo.toml b/crates/feos-ad/Cargo.toml
@@ -0,0 +1,29 @@
+[package]
+name = "feos-ad"
+edition.workspace = true
+version = "0.2.3"
+authors = ["Philipp Rehner <prehner@ethz.ch"]
+homepage.workspace = true
+license.workspace = true
+repository.workspace = true
+keywords.workspace = true
+categories.workspace = true
+description = "FeOs-AD - Implicit automatic differentiation of equations of state and phase equilibria."
+
+[dependencies]
+num-dual = { workspace = true }
+quantity = { workspace = true, features = ["num-dual"] }
+nalgebra = { workspace = true }
+feos-core = { workspace = true }
+ndarray = { workspace = true }
+
+[dev-dependencies]
+feos = { workspace = true }
+approx = { workspace = true }
+quantity = { workspace = true, features = ["num-dual", "approx"] }
+
+[features]
+default = []
+pcsaft = ["feos/pcsaft"]
+gc_pcsaft = ["feos/gc_pcsaft"]
+parameter_fit = []
diff --git a/crates/feos-ad/README.md b/crates/feos-ad/README.md
@@ -0,0 +1,33 @@
+# FeOs-AD - Implicit automatic differentiation of equations of state and phase equilibria
+
+[![crate](https://img.shields.io/crates/v/feos-ad.svg)](https://crates.io/crates/feos-ad)
+[![documentation](https://docs.rs/feos-ad/badge.svg)](https://docs.rs/feos-ad)
+
+The `FeOs-AD` crate builds on the implementation of phase equilibrium calculations in `FeOs` to provide implicit automatic differentiation of properties and phase equilibria based on Helmholtz energy equations of state. Derivatives can be determined for any inputs, like temperature or pressure, but also model parameters.
+
+Derivatives are calculated using forward automatic differentiation with generalized (hyper-) dual numbers from the [`num-dual`](https://github.com/itt-ustutt/num-dual) crate.
+
+## Contents
+For now, the most important properties and phase equilibria are implemented:
+- **State construction**
+    - from temperature and pressure
+    - from pressure and entropy
+    - from pressure and enthalpy
+    - critical points
+- **Phase equilibria**
+    - bubble points
+    - dew points
+    - tp flash
+- **Properties**
+    - pressure, molar entropy, molar enthalpy
+
+The following **models** are implemented:
+- **PC-SAFT** for pure components and binary mixtures
+- heterosegmented **gc-PC-SAFT** for pure components and mixtures
+- The **Joback & Reid** GC method for ideal gas heat capacities
+
+## Installation
+Just add the dependency to your `Cargo.toml`
+```toml
+feos-ad = "0.2"
+```
diff --git a/crates/feos-ad/license-apache b/crates/feos-ad/license-apache
@@ -0,0 +1 @@
+../../license-apache
diff --git a/crates/feos-ad/license-mit b/crates/feos-ad/license-mit
@@ -0,0 +1 @@
+../../license-mit
diff --git a/crates/feos-ad/src/core/mod.rs b/crates/feos-ad/src/core/mod.rs
@@ -0,0 +1,102 @@
+use feos_core::{Components, IdealGas, Residual, StateHD};
+use nalgebra::{Const, SVector, U1};
+use ndarray::{arr1, Array1, ScalarOperand};
+use num_dual::{Derivative, DualNum, DualVec};
+use std::sync::Arc;
+
+mod phase_equilibria;
+mod residual;
+mod state;
+mod total;
+pub use phase_equilibria::PhaseEquilibriumAD;
+pub use residual::{ParametersAD, ResidualHelmholtzEnergy};
+pub use state::{Eigen, StateAD};
+pub use total::{EquationOfStateAD, IdealGasAD, TotalHelmholtzEnergy};
+
+/// Used internally to implement the [Residual] and [IdealGas] traits from FeOs.
+pub struct FeOsWrapper<E, const N: usize>(E);
+
+impl<R: ParametersAD, const N: usize> Components for FeOsWrapper<R, N> {
+    fn components(&self) -> usize {
+        N
+    }
+
+    fn subset(&self, _: &[usize]) -> Self {
+        panic!("Calculating properties of subsets of models is not possible with AD.")
+    }
+}
+
+impl<R: ResidualHelmholtzEnergy<N>, const N: usize> Residual for FeOsWrapper<R, N> {
+    fn compute_max_density(&self, moles: &Array1<f64>) -> f64 {
+        let total_moles = moles.sum();
+        let molefracs = SVector::from_fn(|i, _| moles[i] / total_moles);
+        self.0.compute_max_density(&molefracs)
+    }
+
+    fn residual_helmholtz_energy_contributions<D: DualNum<f64> + Copy + ScalarOperand>(
+        &self,
+        state: &StateHD<D>,
+    ) -> Vec<(String, D)> {
+        let temperature = state.temperature;
+        let volume = state.volume;
+        let density = SVector::from_column_slice(state.partial_density.as_slice().unwrap());
+        let parameters = self.0.params();
+        let a = R::residual_helmholtz_energy_density(&parameters, temperature, &density) * volume
+            / temperature;
+        vec![(R::RESIDUAL.into(), a)]
+    }
+}
+
+impl<E: TotalHelmholtzEnergy<N>, const N: usize> IdealGas for FeOsWrapper<E, N> {
+    fn ln_lambda3<D: DualNum<f64> + Copy>(&self, temperature: D) -> Array1<D> {
+        let parameters = self.0.params();
+        arr1(&E::ln_lambda3(&parameters, temperature).data.0[0])
+    }
+
+    fn ideal_gas_model(&self) -> String {
+        E::IDEAL_GAS.into()
+    }
+}
+
+/// Struct that stores the reference to the equation of state together with the (possibly) dual parameters.
+pub struct HelmholtzEnergyWrapper<E: ParametersAD, D: DualNum<f64> + Copy, const N: usize> {
+    pub eos: Arc<FeOsWrapper<E, N>>,
+    pub parameters: E::Parameters<D>,
+}
+
+impl<E: ParametersAD, const N: usize> HelmholtzEnergyWrapper<E, f64, N> {
+    /// Manually set the parameters and their derivatives.
+    pub fn derivatives<D: DualNum<f64> + Copy>(
+        &self,
+        parameters: E::Parameters<D>,
+    ) -> HelmholtzEnergyWrapper<E, D, N> {
+        HelmholtzEnergyWrapper {
+            eos: self.eos.clone(),
+            parameters,
+        }
+    }
+}
+
+/// Models for which derivatives with respect to individual parameters can be calculated.
+pub trait NamedParameters: ParametersAD {
+    /// Return a mutable reference to the parameter named by `index` from the parameter set.
+    fn index_parameters_mut<'a, D: DualNum<f64> + Copy>(
+        parameters: &'a mut Self::Parameters<D>,
+        index: &str,
+    ) -> &'a mut D;
+}
+
+impl<E: NamedParameters, const N: usize> HelmholtzEnergyWrapper<E, f64, N> {
+    /// Initialize the parameters to calculate their derivatives.
+    pub fn named_derivatives<const P: usize>(
+        &self,
+        parameters: [&str; P],
+    ) -> HelmholtzEnergyWrapper<E, DualVec<f64, f64, Const<P>>, N> {
+        let mut params: E::Parameters<DualVec<f64, f64, Const<P>>> = self.eos.0.params();
+        for (i, p) in parameters.into_iter().enumerate() {
+            E::index_parameters_mut(&mut params, p).eps =
+                Derivative::derivative_generic(Const::<P>, U1, i)
+        }
+        self.derivatives(params)
+    }
+}