Merge pull request #134 from ipqa-research/dev_autodiff_ge

Hyperdual API for Ge models

python/pyproject.toml

@@ -14,7 +14,7 @@ name = "yaeos"
 description = "Thermodynamic modelling with Equation of State"
 readme = "README.md"
 dependencies = ["numpy"]
-dynamic = ["version"]
+version = "2.3.0"
 requires-python = ">=3.10"
 
 

python/yaeos/__init__.py

@@ -4,7 +4,7 @@
 relevant constants, procedures and objects to have better access to them.
 """
 
-__version__ = "2.3.0"
+import importlib.metadata
 
 import yaeos.constants as constants
 from yaeos.lib import yaeos_c
@@ -38,3 +38,6 @@
     "MHV",
     "HV",
 ]
+
+
+__version__ = importlib.metadata.version("yaeos")

src/adiff/autodiff_api/gemodel_adiff_api.f90

@@ -0,0 +1,162 @@
+module yaeos__adiff_hyperdual_ge_api
+   !! Module that contains the automatic differentiation logic for an Ge model.
+   !!
+   !! All that is needed to define an Ge model that uses automatic
+   !! differentiation with hyperdual numbers is to define a new derived type
+   !! that overloads the method to the Ge function that you want to use.
+   use yaeos__constants, only: pr
+   use yaeos__models_ge, only: GeModel
+   use hyperdual_mod
+
+   implicit none
+
+   type, abstract, extends(GeModel) :: GeModelAdiff
+   contains
+      procedure(hyperdual_ge), deferred :: Ge
+      procedure :: excess_gibbs => excess_gibbs
+   end type GeModelAdiff
+
+   abstract interface
+      type(hyperdual) function hyperdual_Ge(self, n, t)
+         import hyperdual, GeModelAdiff
+         class(GeModelAdiff) :: self
+         type(hyperdual), intent(in) :: n(:)
+         type(hyperdual), intent(in) :: t
+      end function hyperdual_Ge
+   end interface
+contains
+
+   subroutine excess_gibbs(self, n, T, Ge, GeT, GeT2, Gen, GeTn, Gen2)
+      class(GeModelAdiff), intent(in) :: self
+      real(pr), intent(in) :: n(:)
+      real(pr), intent(in) :: t
+      real(pr), optional, intent(out) :: Ge, GeT, GeT2
+      real(pr), optional, dimension(size(n)), intent(out) :: Gen, GeTn
+      real(pr), optional, intent(out) :: Gen2(size(n), size(n))
+
+      type(hyperdual) :: d_t, d_n(size(n))
+      type(hyperdual) :: d_Ge
+
+      real(pr) :: dGe(size(n)+1, size(n)+1)
+      integer :: nc
+
+      logical :: any_deriv
+
+      any_deriv = .false.
+
+      nc = size(n)
+
+      if (present(GeT) .or. present(GeT2)) then
+         any_deriv = .true.
+         if (present(GeT2)) then
+            call get_dgedt2
+         end if
+
+         if (.not. (present(GeT2) .and. .not. present(GeTn))) then
+            call get_dgedt
+         end if
+      end if
+
+      if (present(GeTn)) then
+         any_deriv = .true.
+         call get_dgedtn
+      end if
+
+
+      if (present(Gen) .or. present(Gen2)) then
+         any_deriv = .true.
+         if (present(Gen2)) then
+            call get_dgedn2
+         else
+            call get_dgedn
+         end if
+      end if
+
+      if (present(Ge)) then
+         if (.not. any_deriv) then
+            call reset_vars
+            d_ge = self%Ge(d_n, d_t)
+         end if
+         Ge = d_ge%f0
+      end if
+
+   contains
+
+      subroutine get_dgedn()
+         integer :: i, j
+
+         do i=2, size(n), 2
+            call reset_vars
+            d_n(i-1)%f1 = 1
+            d_n(i )%f2 = 1
+
+            d_Ge = self%Ge(d_n, d_t)
+
+            Gen(i-1) = d_Ge%f1
+            Gen(i) = d_Ge%f2
+         end do
+
+         if (mod(size(n), 2) /= 0) then
+            call reset_vars
+            d_n(size(n))%f1 = 1
+            d_Ge = self%Ge(d_n, d_t)
+            Gen(size(n)) = d_Ge%f1
+         end if
+
+      end subroutine get_dgedn
+
+      subroutine get_dgedn2()
+         integer :: i, j
+
+         do i=1,size(n)
+            do j=i,size(n)
+               call reset_vars
+               d_n(i)%f1 = 1
+               d_n(j)%f2 = 1
+
+               d_Ge = self%Ge(d_n, d_t)
+
+               if(present(Gen)) Gen(i) = d_Ge%f1
+               Gen2(i, j) = d_Ge%f12
+               Gen2(j, i) = d_Ge%f12
+            end do
+         end do
+      end subroutine get_dgedn2
+
+      subroutine get_dgedtn()
+         integer :: i
+
+         do i=1,size(n)
+            call reset_vars
+            d_n(i)%f1 = 1
+            d_t%f2 = 1
+            d_Ge = self%Ge(d_n, d_t)
+            if (present(Gen)) Gen(i) = d_Ge%f1
+            if (present(GeT)) GeT = d_Ge%f2
+            GeTn(i) = d_Ge%f12
+         end do
+      end subroutine get_dgedtn
+
+      subroutine get_dgedt()
+         call reset_vars
+         d_t%f1 = 1
+         d_Ge = self%Ge(d_n, d_t)
+         GeT = d_Ge%f1
+      end subroutine get_dgedt
+
+      subroutine get_dgedt2()
+         call reset_vars
+         d_t%f1 = 1
+         d_t%f2 = 1
+         d_Ge = self%Ge(d_n, d_t)
+         if (present(GeT)) GeT = d_Ge%f1
+         if (present(GeT2)) GeT2 = d_Ge%f12
+      end subroutine get_dgedt2
+
+      subroutine reset_vars()
+         d_n = n
+         d_t = t
+      end subroutine reset_vars
+   end subroutine excess_gibbs
+
+end module yaeos__adiff_hyperdual_ge_api

src/adiff/autodiff_api/tapenade_ge_api.f90

@@ -118,6 +118,7 @@ subroutine excess_gibbs(&
             Gen = nb
             if (present(GeT)) GeT = tb
       else
+         if (present(GeT)) GeT = get_dGedT()
       end if
 
       if (present(GeTn)) GeTn = get_GenT()
@@ -159,6 +160,18 @@ function get_dGedT2()
          get_dGedT2 = gedd
       end function
 
+      function get_dGedT()
+         real(pr) :: get_dGedT
+         call reset_vars
+
+         td = 1
+
+         call self%ge_d(&
+            n, nd, t, td, ge, ged &
+            )
+         get_dGedT = ged
+      end function
+
       function get_GenT()
          real(pr) :: get_GenT(size(n))
          call reset_vars

test/test_ge_hyperdual.f90

@@ -0,0 +1,121 @@
+module nrtl_hd
+   !! Example implementation of the NRTL model using automatic differentiation.
+   !!
+   !! The NRTL model is represented by the equations:
+   !!
+   !! \[tau = A + B/T\]
+   !! \[G = exp(-C * tau)\]
+   !! \[Ge = R * T * sum(n * sum(n * tau * G) / sum(n * G))\]
+   use yaeos
+   use yaeos__adiff_hyperdual_ge_api, only: GeModelAdiff
+   use hyperdual_mod
+   implicit none
+
+   type, extends(GeModelAdiff) :: NRTLHD
+      !! NRTL model with automatic differentiation.
+      real(pr), allocatable :: A(:,:), B(:,:), C(:,:)
+   contains
+      procedure :: Ge => Ge
+   end type NRTLHD
+
+contains
+   function ge(self, n, t)
+      class(NRTLHD) :: self
+      type(hyperdual), intent(in) :: n(:)
+      type(hyperdual), intent(in) :: t
+      type(hyperdual) :: ge
+
+      type(hyperdual) ::tau(size(n), size(n)), G(size(n), size(n))
+
+
+      real(pr) :: down
+      integer :: i, j
+
+      tau = self%a(:, :) + self%b(:, :)/T
+      G = exp(-self%C * tau)
+
+      ge = 0._pr
+      do i=1,size(n)
+         ge = ge + n(i) * sum(n(:) * tau(:, i) * G(:, i)) / sum(n(:) *  g(:, i))
+      end do
+      ge = R * T * ge
+   end function ge
+end module nrtl_hd
+
+program test
+   use yaeos
+   use nrtl_hd
+   use fixtures_models, only: binary_NRTL_tape
+   use testing_aux, only: test_title, assert
+
+   type(NRTL) :: og
+   type(NRTLHD) :: imp
+
+   real(pr) :: n(3), T
+
+   real(pr) :: Geog, Geimp
+   real(pr) :: GeogT, GeimpT
+   real(pr) :: GeogT2, GeimpT2
+   real(pr) :: Geogn(3), Geimpn(3)
+   real(pr) :: GeognT(3), GeimpNt(3)
+   real(pr) :: Geogn2(3,3), Geimpn2(3,3)
+
+   real(pr) :: A(3,3), B(3,3), C(3,3)
+
+   A = 0; B=0; C=0   
+
+   A(1, 2) = 0.1; A(1, 3) = 0.2
+   A(2, 1) = 0.3; A(2, 3) = 0.4
+   A(3, 1) = 0.5; A(3, 2) = 0.6
+
+   B(1, 2) = 300; B(1, 3) = 500
+   B(2, 1) = 700; B(2, 3) = 900
+   B(3, 1) = 1100; B(3, 2) = 1300
+
+   C(1, 2) = 0.7; C(1, 3) = 0.8; C(2, 3) = 1.0
+   C(2, 1) = 1.1; C(3, 1) = 1.2; C(3, 2) = 1.3
+
+   og%a = A
+   og%b = B
+   og%C = C
+
+   imp%A = og%a
+   imp%B = og%b
+   imp%C = og%C
+
+   n = [0.2, 0.8, 0.1]
+   T = 273
+
+   write(*, *) test_title("NRTL model with automatic differentiation")
+
+   call og%excess_gibbs(n, T, Geog)
+   call imp%excess_gibbs(n, T, Geimp)
+
+   call assert(abs(Geog - Geimp) < 1e-5, "Ge")
+
+   call og%excess_gibbs(n, T, Geog, GeT=GeogT)
+   call imp%excess_gibbs(n, T, Geimp, GeT=GeimpT)
+
+   call assert(abs(GeogT - GeimpT) < 1e-5, "GeT")
+
+   call og%excess_gibbs(n, T, Ge=Geog, GeT2=GeogT2)
+   call imp%excess_gibbs(n, T, Ge=Geimp, GeT2=GeimpT2)
+
+   call assert(abs(GeogT2 - GeimpT2) < 1e-5, "GeT2")
+
+   call og%excess_gibbs(n, T, Ge=Geog, GeTn=GeognT)
+   call imp%excess_gibbs(n, T, Ge=Geimp, GeTn=GeimpnT)
+
+   call assert(all(abs(GeognT - GeimpnT) < 1e-5), "GeTn")
+
+   call og%excess_gibbs(n, T, Geog, Gen=Geogn)
+   call imp%excess_gibbs(n, T, Geimp, Gen=Geimpn)
+
+   call assert(all(abs(Geogn - Geimpn) < 1e-5), "Gen")
+
+   call og%excess_gibbs(n, T, Geog, Gen2=Geogn2)
+   call imp%excess_gibbs(n, T, Geog, Gen2=Geimpn2)
+
+   call assert(all(abs(Geogn2 - Geimpn2) < 1e-5), "Gen2")
+
+end program

test/test_tx.f90

@@ -18,12 +18,9 @@ program main
    type(Substance) :: sus(nc)
    type(PurePsat) :: vp1, vp2
 
-   integer :: i, tf_brute, tf_criti, tf_envel
-
-   open(newunit=tf_brute, file="test_tx_brute")
-   open(newunit=tf_criti, file="test_tx_criti")
-   open(newunit=tf_envel, file="test_tx_envel")
+   integer :: i
 
+   write(*, *) test_title("TX Envelope 2ph")
 
    forsus_dir = "build/dependencies/forsus/" // forsus_default_dir
    sus(1) = Substance("methane")