diff --git a/python/pyproject.toml b/python/pyproject.toml index 067daea1..a4d3d277 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 eb38b8c9..39ac634f 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") 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 00000000..b246fd1c --- /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 44bdbf08..077f60d1 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 00000000..95db5d9b --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/test_tx.f90 b/test/test_tx.f90 index a5b5dae7..33a3adbb 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")