Merge branch 'main' into dev

.github/workflows/CI.yml

@@ -67,7 +67,7 @@ jobs:
 
     - name: Upload coverage reports to Codecov
       if: ${{ env.FC == 'gfortran' }}
-      uses: codecov/codecov-action@v4.0.1
+      uses: codecov/codecov-action@v4.4.0
       with:
         token: ${{ secrets.CODECOV_TOKEN }}
         slug: ipqa-research/yaeos

README.md

@@ -44,6 +44,8 @@ files with [FORD](https://github.com/Fortran-FOSS-Programmers/ford). On the
 other hand, the Python API documentation can be generated by processing the
 source files with [Sphinx](https://github.com/sphinx-doc/sphinx).
 
+## Developers
+This library is currently maintained by the research group of Prof. Martín Cismondi-Duarte at [IPQA (UNC-CONICET)](https://ipqa.unc.edu.ar/en/)
 
 ## Available models
 - CubicEoS
@@ -110,7 +112,7 @@ with:
  fpm run --example <example name here>
 ```
 
-Not providing any example will show all the possible examples that can be run.
+Not providing any examples will show all the possible examples that can be run.
 
 # How to install/run it
 
@@ -188,9 +190,3 @@ A complete implementation of the PR76 Equation of State can me found in
 ## Tapenade-based autodiff
 It is also possible to differentiate with `tapenade`. Examples can be seen
 in the documentation pages or in [The tools directory](tools/tapenade_diff/)
-
-# Documentation
-The latest API documentation for the `main` branch can be found
-[here](https://ipqa-research.github.io/yaeos). This was generated from the source
-code using [FORD](https://github.com/Fortran-FOSS-Programmers/ford). We're
-working in extending it more.

doc/page/usage/eos/cubics/index.md

@@ -2,13 +2,73 @@
 title: Cubics
 ---
 
-** Table of contents**
-[TOC]
+All our Cubic Equations of State are implemented based on the generic Cubic
+Equation:
+
+\[
+    P = \frac{RT}{V-b} - \frac{a_c\alpha(T_r)}{(V + \delta_1 b)(V - \delta_2 b)}
+\]
 
 ## SoaveRedlichKwong
+[[SoaveRedlichKwong]]
+
+Using the critical constants setup the parameters to use the 
+SoaveRedlichKwong Equation of State
+
+- \[\alpha(T_r) = (1 + k (1 - \sqrt{T_r}))^2\]
+- \[k = 0.48 + 1.574 \omega - 0.175 \omega^2 \]
+- \[a_c = 0.427480  R^2 * T_c^2/P_c\]
+- \[b = 0.086640  R T_c/P_c\]
+- \[\delta_1 = 1\]
+- \[\delta_2 = 0\]
+
+There is also the optional posibility to include the k_{ij} and l_{ij}
+matrices. Using by default Classic Van der Waals mixing rules. For more information
+about mixing rules look at [Mixing Rules](mixing.md)
 
 ## PengRobinson76
+[[PengRobinson76]]
+
+- \[\alpha(T_r) = (1 + k (1 - \sqrt{T_r}))^2\]
+- \[k = 0.37464 + 1.54226 * \omega - 0.26993 \omega^2 \]
+- \[a_c = 0.45723553  R^2 T_c^2 / P_c\]
+- \[b = 0.07779607r  R T_c/P_c\]
+- \[\delta_1 = 1 + \sqrt{2}\]
+- \[\delta_2 = 1 - \sqrt{2}\]
 
 ## PengRobinson78
+[[PengRobinson78]]
+
+- \[\alpha(T_r) = (1 + k (1 - \sqrt{T_r}))^2\]
+- \[k = 0.37464 + 1.54226 \omega - 0.26992 \omega^2  \text{ where } \omega <=0.491\]
+- \[k = 0.37464 + 1.48503 \omega - 0.16442 \omega^2  + 0.016666 \omega^3 \text{ where } \omega > 0.491\]
+- \[a_c = 0.45723553  R^2 T_c^2 / P_c\]
+- \[b = 0.07779607r  R T_c/P_c\]
+- \[\delta_1 = 1 + \sqrt{2}\]
+- \[\delta_2 = 1 - \sqrt{2}\]
+
+## RKPR
+[[RKPR]]
+
+The RKPR EoS extends the classical formulation of Cubic Equations of State by
+freeing the parameter \(\delta_1\) and setting 
+\(\delta_2 = \frac{1+\delta_1}{1-\delta_1}\).
+This extra degree provides extra ways of implementing the equation in
+comparison of other Cubic EoS (like PR and SRK) which are limited to
+definition of their critical constants.
+
+Besides that extra parameter, the RKRR includes another \(\alpha\)
+function:
+
+\[
+ \alpha(T_r) = \left(\frac{3}{2+T_r}\right)^k
+\]
+
+In this implementation we take the simplest form which correlates
+the extra parameter to the critical compressibility factor \(Z_c\) and
+the \(k\) parameter of the \(\alpha\) function to \(Z_c\) and \(\omega\):
+
+\[\delta_1 = d_1 + d_2 (d_3 - Z_c)^d_4 + d_5 (d_3 - Z_c) ^ d_6\]
+\[k = (A_1  Z_c + A_0)\omega^2 + (B_1 Z_c + B_0)\omega + (C_1 Z_c + C_0)\]
 
-## RKPR
+It is also possible to include the parameters as optional arguments.

doc/page/usage/eos/cubics/mixing.md

@@ -62,8 +62,21 @@ subroutine my_aij_implementation(self, ai, daidt, daidt2, aij, daijdt, daijdt2)
 end subroutine
 ```
 
-### Constant \(k_{ij}\)
-
 ## \(G^E\) Models Mixing Rules
+It is possible to mix the attractive parameter of Cubic Equations with an 
+excess Gibbs-based model.
+This can be useful for cases of polar molecules and/or systems that have been 
+fitted to \(G^E\) models.
+
+### Michelsen's Modified Huron-Vidal Mixing Rules
+This mixing rule is based on the aproximate zero-pressure limit 
+of a cubic equation of state. At the aproximate zero-pressure limit the
+attractive parameter can be expressed as:
+
+\[
+\frac{D}{RTB}(n, T) = \sum_i n_i \frac{a_i(T)}{b_i} + \frac{1}{q}
+\left(\frac{G^E(n, T)}{RT} + \sum_i n_i \ln \frac{B}{nb_i} \right)
+\]
 
-### Michelsen's Modified Huron-Vidal Mixing Rules
+Where \(q\) is a weak function of temperature. In the case of `MHV`
+and simplicity it is considered that depends on the model used.

doc/page/usage/newmodels/autodiff.md

@@ -56,61 +56,61 @@ module hyperdual_pr76
 
 contains
 
-   type(PR76) function setup(tc, pc, w, kij, lij) result(self)
-      !! Function to obtain a defined PR76 model with setted up parameters
-      !! as function of Tc, Pc, and w
-      real(pr) :: tc(:)
-      real(pr) :: pc(:)
-      real(pr) :: w(:)
-      real(pr) :: kij(:, :)
-      real(pr) :: lij(:, :)
-
-      self%components%tc = tc
-      self%components%pc = pc
-      self%components%w = w
-
-      self%ac = 0.45723553_pr * R**2 * tc**2 / pc
-      self%b = 0.07779607_pr * R * tc/pc
-      self%k = 0.37464_pr + 1.54226_pr * w - 0.26993_pr * w**2
-
-      self%kij = kij
-      self%lij = lij
-   end function setup
-
-   function arfun(self, n, v, t) result(ar)
-      !! Residual Helmholtz calculation for a generic cubic with
-      !! quadratic mixing rules.
-      use hyperdual_mod
-      class(PR76) :: self
-      type(hyperdual), intent(in) :: n(:), v, t
-      type(hyperdual) :: ar
-
-      type(hyperdual) :: amix, a(size(n)), ai(size(n)), n2(size(n))
-      type(hyperdual) :: bmix
-      type(hyperdual) :: b_v, nij
-
-      integer :: i, j
-
-      ! Associate allows us to keep the later expressions simple.
-      associate(&
-         pc => self%components%pc, ac => self%ac, b => self%b, k => self%k,&
-         kij => self%kij, lij => self%lij, tc => self%components%tc &
-         )
-
-         ! Soave alpha function
-         a = 1.0_pr + k * (1.0_pr - sqrt(t/tc))
-         a = ac * a ** 2
-         ai = sqrt(a)
-
-         ! Quadratic Mixing Rule
-         amix = 0.0_pr
-         bmix = 0.0_pr
-
-         do i=1,size(n)-1
-            do j=i+1,size(n)
-               nij = n(i) * n(j)
-               amix = amix + 2 * nij * (ai(i) * ai(j)) * (1 - kij(i, j))
-               bmix = bmix + nij * (b(i) + b(j)) * (1 - lij(i, j))
+    type(PR76) function setup(tc, pc, w, kij, lij) result(self)
+        !! Function to obtain a defined PR76 model with setted up parameters
+        !! as function of Tc, Pc, and w
+        real(pr) :: tc(:)
+        real(pr) :: pc(:)
+        real(pr) :: w(:)
+        real(pr) :: kij(:, :)
+        real(pr) :: lij(:, :)
+
+        self%composition%tc = tc
+        self%composition%pc = pc
+        self%composition%w = w
+
+        self%ac = 0.45723553_pr * R**2 * tc**2 / pc
+        self%b = 0.07779607_pr * R * tc_in/pc_in
+        self%k = 0.37464_pr + 1.54226_pr * w - 0.26993_pr * w**2
+
+        self%kij = kij
+        self%lij = lij
+    end function
+
+    function arfun(self, n, v, t) result(ar)
+        !! Residual Helmholtz calculation for a generic cubic with
+        !! quadratic mixing rules.
+        class(PR76) :: self
+        type(hyperdual), intent(in) :: n(:), v, t
+        type(hyperdual) :: ar
+    
+        type(hyperdual) :: amix, a(size(n)), ai(size(n)), n2(size(n))
+        type(hyperdual) :: bmix
+        type(hyperdual) :: b_v, nij
+
+        integer :: i, j
+
+        ! Associate allows us to keep the later expressions simple.
+        associate(&
+            pc => self%composition%pc, ac => self%ac, b => self%b, k => self%k,&
+            kij => self%kij, lij => self%lij, tc => self%compostion%tc & 
+            )
+
+            ! Soave alpha function
+            a = 1.0_pr + k * (1.0_pr - sqrt(t/tc))
+            a = ac * a ** 2
+            ai = sqrt(a)
+
+            ! Quadratic Mixing Rule
+            amix = 0.0_pr
+            bmix = 0.0_pr
+
+            do i=1,size(n)-1
+                do j=i+1,size(n)
+                    nij = n(i) * n(j)
+                    amix = amix + 2 * nij * (ai(i) * ai(j)) * (1 - kij(i, j))
+                    bmix = bmix + nij * (b(i) + b(j)) * (1 - lij(i, j))
+                end do
             end do
          end do