From 9a7a317420d48c6e340637e3f5c68f679ac8e43b Mon Sep 17 00:00:00 2001 From: fedebenelli Date: Sun, 29 Dec 2024 23:35:57 -0300 Subject: [PATCH 1/4] Hyperdual API for Ge models --- src/adiff/autodiff_api/gemodel_adiff_api.f90 | 162 +++++++++++++++++++ src/adiff/autodiff_api/tapenade_ge_api.f90 | 13 ++ test/test_ge_hyperdual.f90 | 105 ++++++++++++ 3 files changed, 280 insertions(+) create mode 100644 src/adiff/autodiff_api/gemodel_adiff_api.f90 create mode 100644 test/test_ge_hyperdual.f90 diff --git a/src/adiff/autodiff_api/gemodel_adiff_api.f90 b/src/adiff/autodiff_api/gemodel_adiff_api.f90 new file mode 100644 index 000000000..b246fd1ce --- /dev/null +++ b/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 diff --git a/src/adiff/autodiff_api/tapenade_ge_api.f90 b/src/adiff/autodiff_api/tapenade_ge_api.f90 index 44bdbf086..077f60d1e 100644 --- a/src/adiff/autodiff_api/tapenade_ge_api.f90 +++ b/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 diff --git a/test/test_ge_hyperdual.f90 b/test/test_ge_hyperdual.f90 new file mode 100644 index 000000000..7e0189b71 --- /dev/null +++ b/test/test_ge_hyperdual.f90 @@ -0,0 +1,105 @@ +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(-alpha * 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(:,:) + real(pr) :: alpha + 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%alpha * 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(2), T + + real(pr) :: Geog, Geimp + real(pr) :: GeogT, GeimpT + real(pr) :: GeogT2, GeimpT2 + real(pr) :: Geogn(2), Geimpn(2) + real(pr) :: GeognT(2), GeimpNt(2) + real(pr) :: Geogn2(2,2), Geimpn2(2,2) + + og = binary_NRTL_tape() + + imp%A = og%a + imp%B = og%b + imp%alpha = og%C(1,2) + + n = [0.2, 0.8] + 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 \ No newline at end of file From a726056c5b9b77adf32a341bb79a1c98a8348ae6 Mon Sep 17 00:00:00 2001 From: fedebenelli Date: Sun, 29 Dec 2024 23:42:48 -0300 Subject: [PATCH 2/4] version --- python/pyproject.toml | 2 +- python/yaeos/__init__.py | 5 ++++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/python/pyproject.toml b/python/pyproject.toml index 067daea11..a4d3d277b 100644 --- a/python/pyproject.toml +++ b/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" diff --git a/python/yaeos/__init__.py b/python/yaeos/__init__.py index eb38b8c90..39ac634fe 100644 --- a/python/yaeos/__init__.py +++ b/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") From 3807ac122058d660e7f85917b6c9525cc5b5b844 Mon Sep 17 00:00:00 2001 From: fedebenelli Date: Sun, 29 Dec 2024 23:45:25 -0300 Subject: [PATCH 3/4] test title --- test/test_tx.f90 | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/test/test_tx.f90 b/test/test_tx.f90 index a5b5dae7e..33a3adbb1 100644 --- a/test/test_tx.f90 +++ b/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") From 62b1af94d35bc28c62c3f99596a09d84d1ede054 Mon Sep 17 00:00:00 2001 From: fedebenelli Date: Sun, 29 Dec 2024 23:51:30 -0300 Subject: [PATCH 4/4] better coverage --- test/test_ge_hyperdual.f90 | 38 +++++++++++++++++++++++++++----------- 1 file changed, 27 insertions(+), 11 deletions(-) diff --git a/test/test_ge_hyperdual.f90 b/test/test_ge_hyperdual.f90 index 7e0189b71..95db5d9b5 100644 --- a/test/test_ge_hyperdual.f90 +++ b/test/test_ge_hyperdual.f90 @@ -4,7 +4,7 @@ module nrtl_hd !! The NRTL model is represented by the equations: !! !! \[tau = A + B/T\] - !! \[G = exp(-alpha * tau)\] + !! \[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 @@ -13,8 +13,7 @@ module nrtl_hd type, extends(GeModelAdiff) :: NRTLHD !! NRTL model with automatic differentiation. - real(pr), allocatable :: A(:,:), B(:,:) - real(pr) :: alpha + real(pr), allocatable :: A(:,:), B(:,:), C(:,:) contains procedure :: Ge => Ge end type NRTLHD @@ -33,7 +32,7 @@ function ge(self, n, t) integer :: i, j tau = self%a(:, :) + self%b(:, :)/T - G = exp(-self%alpha * tau) + G = exp(-self%C * tau) ge = 0._pr do i=1,size(n) @@ -52,22 +51,39 @@ program test type(NRTL) :: og type(NRTLHD) :: imp - real(pr) :: n(2), T + real(pr) :: n(3), T real(pr) :: Geog, Geimp real(pr) :: GeogT, GeimpT real(pr) :: GeogT2, GeimpT2 - real(pr) :: Geogn(2), Geimpn(2) - real(pr) :: GeognT(2), GeimpNt(2) - real(pr) :: Geogn2(2,2), Geimpn2(2,2) + real(pr) :: Geogn(3), Geimpn(3) + real(pr) :: GeognT(3), GeimpNt(3) + real(pr) :: Geogn2(3,3), Geimpn2(3,3) - og = binary_NRTL_tape() + 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%alpha = og%C(1,2) + imp%C = og%C - n = [0.2, 0.8] + n = [0.2, 0.8, 0.1] T = 273 write(*, *) test_title("NRTL model with automatic differentiation")