Repository navigation
Description
Dear all,
I am a PhD in the Engineering Thermodynamics group at the TU Delft. We mostly work on molecular simulations with Monte Carlo and molecular dynamics simulations in Fortran and C++. In my free time, I am learning Rust (I just finished coding my own LJ Monte Carlo program for fun) and came across this package because an office mate was using FeOs. He is using FeOs (from the python side) for phase equilibrium predictions of multi-component mixtures, where he computes among others spinodal lines.
Besides FeOs he is using Refprop and the GERG-2008 EoS (doi: 10.1021/je300655b) and asked me if it would be possible to implement this EoS as well into FeOs. I looked into the Peng-Robinson implementation (python example in the tutorial and the Rust code in the feos-core crate) and in the half hour I spend on it, it seems that FeOs requires Helmholtz free energy functions for correct predictions of its state. This is similar to GERG-2008, although that is a multivariable empirical relation instead of something as elegant and simple as the already implemented PR EoS. I propose the following process, as I have mediocre experience in Rust:
- I implement and test the cubic equations of states first in python to check if I have proper understanding of what I am doing.
- I implement additional cubic equation of states first to get some Rust practice. I went through the code and expect that I would not have too many problems in implementing the VdW and SRK EoSs for single and multi-component fluids.
- When that is implemented and merged successfully in the main FeOs branch, I will start on implementing empirical Helmholtz free energy models such as GERG-2008 and maybe IAPWS95.
For this, I would need the following from the FeOs maintainers:
- Some instructions on structure, styling, and testing guidelines in FeOs (as to integrate it well).
- Are there other dependencies and requirements besides those mentioned on the prerequisite page in the documentation or included in the Cargo.toml files? (maybe there is some C++ Makefile trauma speaking here :) .)
- Besides the: feos/crates/feos-core/src/cubic.rs, which other files should be changed to produce all required python bindings and such?
- And connected to that. Why are the cubic EoS rust files in the src root directory and not in the equation_of_state subdirectory like the ideal gas law?
- While I have some ideas already on how VdW and SRK can be implemented, I expect more issues with the GERG model. Do you foresee pitfalls for such empirical models, and can you give some advises on how to approach it?
I believe that FeOs would get the following gains from this proposal:
- Having extensive, fast cubic EoS support might attract CFD crates to use FeOs for the prediction of fluid properties during simulations. Compressible CFD simulations require good and most of all fast EoS predictions for density, heat capacity and viscosity. I believe that FeOs can be attractive for such programs when there are more cubic equations of state to choose from.
- When GERG-2008 is implemented, FeOs can claim that it can compute fluid properties compliant with the ISO 20765-2/3 standard. Scientifically, I do not see too much value in this, however, it might make FeOs suitable for usage in industry.
I cannot work on this project professionally, as my PhD is about computing transport properties in complex aqueous salt mixtures and I do not intend to use FeOs for my own projects. However, I think I can learn aplenty from this side project (especially regarding traits) and I am willing to spend some evenings on it.
Would you guys be interested in such extensions to FeOs?
Kind regards,
Jelle Lagerweij
Activity
g-bauer commentedon May 15, 2025
Hi Jelle,
thank you very much for the issue and your interest in FeOs. I'd happily guide you through the steps needed to implement a new model. You can also ask here at any point or join our discord server (link is in the repo-readme).
The general procedure could look as follows:
Helmholtz energy
In FeOs, residual equations of state (without the ideal gas part) are implemented by defining the Helmholtz energy as function of the thermodynamic conditions. These conditions are not encoded as regular floating point numbers but as generic (hyper-) dual numbers. In essence, these numbers enable us to write a single function (Helmholtz energy) and yield all exact partial derivatives without an analytical implementation. So all you have to care about is how to implement the model - properties are calculated via automatic differentiation.
Because we need to formulate the model as dual numbers, you'll find that the relevant functions in the
Residualtrait, which you'll have to implement, are parameterized over theDualNumtrait.DualNumimplements all relevant arithmetic operations so you can use dual numbers almost like regular floats (there are some quirks to get used to).So in a nutshell: Implementing a new equation of state means implementing the
Residualtrait (and theIdealGastrait if your model needs a specific ideal gas model).Your questions
feoscrate (crates/feosdirectory). Each model has its own module infeosand is gated behind a cargo-feature. We userustfmtto format everything andclippyas linter. Tests are usually added inside the module where the functions are defined (see existing models for examples).Cargo.toml:Dfeos-core. New cubic models will be in thefeoscrate.Awesome, that's how I started with Rust as well :D Our molecular simulation lectures use a Rust code with Python bindings nowadays (an updated version of this code).
g-bauer commentedon May 15, 2025
I started adding cubic equations of state (in a generic way) some time ago. The branch is still here, although we since changed the project's structure a bit (see here). Maybe you can use that as a starting point.
Also, when working with a new model in FeOs, configure your IDE/editor so that only your model is active (disable all other feos-features) and use
cargo build(develop profile, not release or even lto). That way, the crates compile way faster.prehner commentedon May 16, 2025
Thank you @JelleLagerweij for the interest in contributing to feos and thank you @g-bauer for the detailed tutorial.
JelleLagerweij commentedon May 16, 2025
Dear Gernot,
Thanks for you quick reply. I see that your impl_cubics branch is indeed quite different (not only in project structure, but it also has a more general cubic eos setup compared to the cubic.rs file in the main branch). To start, I will do the following:
If that work out, then I will see how we can continue with the empirical models. Which, I expect, shall be their own module within the feos crate.
Then regarding sargas, I had some fun with reading through it. We of course have our own molecular simulation course, for which I have been the TA a few times already. We just give the students the book of Frenkel and Smit and tell them you have 3 weeks for coding your own MC simulation for a LJ system. Then give them 3 weeks for their own MD code (with Nosé-Hoover as thermostat) and lastly 3 weeks for transport properties post-processing script to get diffusion and viscosity. They can use the programming language of their choice (all choose python). We noticed that students struggle with the freedom and do not know how to structure their code. Additionally, they have problems with computational speed (#python_and_for_loops), especially for when they need to perform longer simulations to get the transport properties.
The way you structured Sargas might actually avoid these issues, as the structure is already there (included tests for their assignments) and the Rust code will probably perform quite a bit better. We will probably keep our own assignments, but sargas gives us some ideas on how we could make the student experience a bit better.
g-bauer commentedon May 16, 2025
Correct, new models should be in
crates/feos/src/model_name.Feel free to take what I've been building, although some renaming (e.g. of the mixing rule) might be needed before merging. Refactoring things like that is quite easy though. If you want to start with less generic implementations - that's also fine. Even starting with the empirical models wouldn't be necessarily harder to do.
Regarding teaching molecular simulations - that's certainly a fascinating approach. But not keep this issue focused, let's move this to a different conversation. :)
JelleLagerweij commentedon May 16, 2025
Dear all,
I cloned the main branch of your repository and I am currently refactoring the cubic parts of the impl_cubic into the new structure. The current progress can be found in https://github.com/JelleLagerweij/feos_cubics (or in your fork list of course :) ).
I started by creating the
mod.rsand theparameters.rsfiles. I had some difficulties with a few imports:feos_core::{EoSresult}andfeos_core::parameter::{ParameterError}. The last one I found, you put all errors together in anerror.rsfile which matches the error case using enums. Therefore, I changed theparameters.rsfile to usefeos_core::{FeosError, FeosResult}.However, I had more issues with understanding
feos_core::{EoSresult}. The impl_cubic implementation for the pub fn new isEosResult<Self>. However, I noticed that the other EoSs return justSelf. I guess that the difference comes from the generality that is intended for the cubic EoSs behavior (that all take critical point inputs ect.). What type shall be returned from the new() function in cubic?Examples:
In crates/feos/src/cubic/mod.rs
In for example crates/feos/src/saftvrmie/eos/mod.rs
g-bauer commentedon May 16, 2025
Yeah we changed that recently. Now there only exists a single error type,
FeosError.There is no hard rule for the return type of the
newfunction. If it makes sense to return aResult, because something can go wrong, that's fine. In this case returning aResultmakes sense because thealpha.validatefunction is fallible.JelleLagerweij commentedon Jun 2, 2025
Dear all,
I have made little progress regarding expanding the cubic equations of state support. I had a few busy weeks. However, in my free time, I tried to fiddle around a bit more with feos::cubic, mostly regarding the testing of the PengRobinson1976 which fails.
I compared with Clapeyron.jl for a logical liquid state point (where we should deviate strongly from the ideal gas behavior)
It returns the following: P = 1.5476632110008812e6 Pa
This is very close (3.4458935260772705e-8 Pa difference) to the PengRobinson1976 as implemented here.
And quite different from the comparison state we are using (feos-core::cubic::PengRobinson is about 8047.94 Pa off compared to the the feos:cubic::PengRobinson1976). So I guess that the PengRobinson in
feos-core::cubic somewhere has an issue. Any ideas where to start when investigating this? Or is dropping the comparision with feos-core::cubic fine if we replace the test with appropriate comparisons with Clapeyron.jl for now?
Besides that, I intended to use
State<Cubic>instead ofStateHD<f64>for testing our implementation (as that gives inputs to which we can control the units of the input, and it interacts more similar with the EoS compared with real usecases so we could use it as an example as well. I could use it easily to compute the pressure comparison.However, if I wanted to get the Helmholtz energy (
a = helmholtz_energy) I get a trait error:This is fully logical, because in feos-core/src/tate/properties.rs we require
impl<E: Residual + IdealGas> State<E> {for the entire code block. This is not required for the pressure function in the residual_properties.rs file. Could you explain to me how to make Cubic satisfy feos_core::IdealGas?cjsisco commentedon Jun 2, 2025
@JelleLagerweij My first attempt at coding in Rust was to expand the cubic EOS to include the temperature-dependent volume translation:
c = c0 + c1*(T - 288.15). I have put it on the back-burner because I have not yet had enough time or expertise to tackle testing in Rust. If it's not too much trouble, I would appreciate it if you could add this to your cubic EOS implementation.This was my attempt at generalizing the cubic to calculate
a, b, kappadepending on the specified version. There are undoubtedly better ways to do it that allow more flexibility with the alpha-function. It seems that @g-bauer was on his way to implementing that in his version.Note that the RKS19 and PR19 are based on the work from Pina-Martinez et al: 10.1016/j.fluid.2018.12.007.
prehner commentedon Jun 2, 2025
Hi, in this case, you just want to use
state.residual_helmholtz_energyorstate.residual_molar_helmholtz_energy, which are independent of the ideal gas contribution and probably the values you want to compare to the results from other implementations.Regarding the deviations in
feos-core::cubic, I don't have the overview over the cubics right now to answer that, maybe @g-bauer has a spontaneous idea, otherwise we can investigate it more thoroughly.In that context also thanks to @cjsisco for the suggestion! Sounds like a good addition.
g-bauer commentedon Jun 2, 2025
As Philipp wrote, since we only want to compare the residual models, we don't need the ideal gas model. For completeness, though, you can use the
EquationOfStatestruct from core - but then you have to specify the ideal gas model:Again, not necessary here, but might come into play when you think about GERG.
Regarding the differences between implementations: The difference stems from the two
omegavalues. In our core implementation, we usewhereas in the generic cubic implementation, we follow Privat et al., which allows to generate PR and SRK (and possibly others) from a single expression:
which yields
If you simply replace the calculated omega's with the values of the core implementation, you'll find that the values match very well. E.g. for the Helmholtz energy:
So as summary - it is probably more sensible to compare against Clapeyron.jl than against our core implementation. ;)
g-bauer commentedon Jun 2, 2025
@cjsisco Thanks for the message. All alpha functions are implemented in this file (might be a good idea to take a very close look at the used constants). A future optimization could be to extract the
m(yourkappa) and only compute it once.And yes - also including volume translation is a good idea. A good place would be to put it into the
CubicOptionsthat also contain the alpha function and the mixing rule. Once everything is properly implemented and tested, it'd make sense to take a close look at the structures and see if we can simplify it a bit to make it more readable.I appreciate all the input.
cjsisco commentedon Jun 2, 2025
@g-bauer According to Jaubert et al. (10.1016/j.fluid.2016.03.012), the Helmholtz energy is unaffected by volume translation (even the temperature-dependent form), so it seems natural to include a
volume_shift_c0andvolume_shift_c1in the default parameter set and default their values=0. Unless there are other volume translation forms than the linear Péneloux form, then I don't foresee the user needing to specify anything inCubicOptions.If they don't provide
volume_shift_*values, then use the defaults. If they do, then apply the simple Péneloux form.JelleLagerweij commentedon Jun 3, 2025
Ok perfect,
I already noticed that I need to get some more hands-on experience before implementing my own EOS fully, so such smaller changes are welcome to me :).
To implement all the above correctly, I propose to restructure the cubic file tree a bit. This would than allow anyone to pick 1 out of every submodule to combine them till their own EoS. I could implement a default mixing rule and a default no-translation method in their specific subtree, and we can use the
CubicOptionsto make people switch between them. If as @cjsisco states that volume translations does not affect the residual Helmholtz energy (I went to a conference where Jaubert was talking about this as well), we could also couple it less strictly to the cubic parts of FeOs. But I believe that from a user perspective, it is more logical to combine all in one.CubicOptionseven if this isn't strictly required.Is this file structure ok?
PS: I will still think about implementing the testing in such a way that we do not depend on a magical external value which is computed through Clapeyron. For the general non-mixed PR and (S)RK eos I can probably compute the pressure directly from their original expressions and compare that to the pressure computed from the implemented residual Helmholtz energy. I will come back to this :).
g-bauer commentedon Jun 3, 2025
Sounds good. I have no strong opinion about having each alpha-function, mixing rule and translation method in its own file.
If you have a good idea how to circumvent this, cool. But it is also common practice to use values from other libraries for comparison. We could, for example, use two libraries such as Clapeyron and thermopack as a first check and then use their values (if they agree) to test our output. That would be totally fine.
cjsisco commentedon Jun 3, 2025
@g-bauer What would it look like to repurpose elements of SAFT for the cubic in
feos?For example, the CPA model of Kontogeorgis is the same association term originally proposed by Walter Chapman for the SAFT model and simply added to the basic repulsive and attractive Helmholtz terms of the cubic. A few years ago, my group hybridized the chain and association terms from the original SAFT into the cubic to make what we call CPCA (cubic+chain+association, 10.1021/acs.iecr.3c02774). This involved a little bit of a change to the repulsive and attractive Helmholtz terms because of the introduction of the chain length.
How much modification do you think would be required to make generic chain and association files that could be used by other models?
g-bauer commentedon Jun 3, 2025
@cjsisco We already have association as separate module here and reuse it across the different implementations. Now, we'd have to take a look at the specific model to see if and how to extend what we already have.
The code is probably not so easy to read if you're not used to generics in Rust. Importantly, one would need to implement the AssociationStrength trait if the current implementation does not suffice.
The chain term on the other hand, we added in the SAFT model's modules because it is a compact implementation. Generalizing contributions often comes at the cost of more complexity because we'd probably have to introduce a new trait. It is always possible to copy the code into the cubic model as a first step, though.