Skip to content

Commit

Permalink
testsss
Browse files Browse the repository at this point in the history
  • Loading branch information
@fedebenelli
fedebenelli committed Mar 10, 2024
1 parent 906f4bf commit ced15ca
Show file tree
Hide file tree
Showing 9 changed files with 518 additions and 285 deletions.
19 changes: 19 additions & 0 deletions test/fixtures/auxiliar_funcions.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
module auxiliar_functions
use yaeos_constants, only: pr
contains
elemental function rel_error(x, y)
real(pr), intent(in) :: x, y
real(pr) :: rel_error

rel_error = abs(x - y)/abs(x)
end function rel_error

function allclose(x, y, atol)
real(pr), intent(in) :: x(:)
real(pr), intent(in) :: y(:)
real(pr), intent(in) :: atol

logical :: allclose
allclose = maxval(rel_error(x, y)) < atol
end function allclose
end module auxiliar_functions
33 changes: 32 additions & 1 deletion test/fixtures/fix_models.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@ type(CubicEoS) function binary_PR76() result(eos)
integer, parameter :: n=2
real(pr) :: tc(n), pc(n), w(n)
real(pr) :: kij(n, n), lij(n, n)
! z = [0.3_pr, 0.7_pr]
tc = [190._pr, 310._pr]
pc = [14._pr, 30._pr]
w = [0.001_pr, 0.03_pr]
Expand All @@ -19,5 +18,37 @@ type(CubicEoS) function binary_PR76() result(eos)

eos = PengRobinson76(tc, pc, w, kij, lij)
end function

type(CubicEoS) function binary_PR78() result(eos)
use yaeos, only: PengRobinson78
integer, parameter :: n=2
real(pr) :: tc(n), pc(n), w(n)
real(pr) :: kij(n, n), lij(n, n)
tc = [190._pr, 310._pr]
pc = [14._pr, 30._pr]
w = [0.001_pr, 0.03_pr]

kij = reshape([0., 0.1, 0.1, 0.], [n,n])
lij = kij / 2


eos = PengRobinson78(tc, pc, w, kij, lij)
end function

type(CubicEoS) function binary_SRK() result(eos)
use yaeos, only: SoaveRedlichKwong
integer, parameter :: n=2
real(pr) :: tc(n), pc(n), w(n)
real(pr) :: kij(n, n), lij(n, n)
tc = [190._pr, 310._pr]
pc = [14._pr, 30._pr]
w = [0.001_pr, 0.03_pr]

kij = reshape([0., 0.1, 0.1, 0.], [n,n])
lij = kij / 2


eos = SoaveRedlichKwong(tc, pc, w, kij, lij)
end function

end module
89 changes: 60 additions & 29 deletions test/test_alphas.f90
View file
Original file line number Diff line number Diff line change
@@ -1,40 +1,71 @@
module test_cubic_alphas
use testdrive, only: new_unittest, unittest_type, error_type, check
implicit none
use yaeos_constants, only: pr
use testdrive, only: new_unittest, unittest_type, error_type, check
use auxiliar_functions, only: allclose
implicit none

real(pr) :: absolute_tolerance = 1e-5_pr

contains
subroutine collect_suite(testsuite)
!> Collection of tests
type(unittest_type), allocatable, intent(out) :: testsuite(:)
subroutine collect_suite(testsuite)
!> Collection of tests
type(unittest_type), allocatable, intent(out) :: testsuite(:)

testsuite = [ &
new_unittest("AlphaSoave", test_alpha_soave), &
new_unittest("AlphaRKPR", test_alpha_RKPR) &
]
end subroutine collect_suite

subroutine test_alpha_soave(error)
use yaeos_constants, only: pr
use yaeos_models_ar_cubic_alphas, only: AlphaSoave
type(error_type), allocatable, intent(out) :: error

type(AlphaSoave) :: alpha

integer, parameter :: n = 2
real(pr) :: Tr(n), k(n)
real(pr) :: a(n), dadt(n), dadt2(n)

real(pr) :: aval(n) = [1.1524213446238356, 1.0594365081389596]
real(pr) :: davaldt(n) = [-0.33947331922020552, -0.14556349186104045]
real(pr) :: davaldt2(n) = [0.47434164902525683, 0.15556349186104046]

Tr = [0.4_pr, 0.5_pr]
k = [0.2_pr, 0.1_pr]

alpha%k = k
call alpha%alpha(Tr, a, dadt, dadt2)

testsuite = [ &
new_unittest("AlphaSoave", test_alpha_soave) &
]
end subroutine
call check(error, allclose(a, aval, absolute_tolerance))
call check(error, allclose(dadt, davaldt, absolute_tolerance))
call check(error, allclose(dadt2, davaldt2, absolute_tolerance))
end subroutine test_alpha_soave

subroutine test_alpha_soave(error)
use yaeos_constants, only: pr
use yaeos_models_ar_cubic_alphas, only: AlphaSoave
type(error_type), allocatable, intent(out) :: error
subroutine test_alpha_RKPR(error)
use yaeos_constants, only: pr
use yaeos_models_ar_cubic_alphas, only: AlphaRKPR
type(error_type), allocatable, intent(out) :: error

type(AlphaSoave) :: alpha
type(AlphaRKPR) :: alpha

integer, parameter :: n = 2
real(pr) :: Tr(n), k(n)
real(pr) :: a(n), dadt(n), dadt2(n)
integer, parameter :: n = 2
real(pr) :: Tr(n), k(n)
real(pr) :: a(n), dadt(n), dadt2(n)

real(pr) :: aval(n) = [1.1524213446238356, 1.0594365081389596]
real(pr) :: davaldt(n) = [-0.33947331922020552,-0.14556349186104045]
real(pr) :: davaldt2(n) = [0.47434164902525683,0.15556349186104046]
real(pr) :: aval(n) = [1.0456395525912732, 1.0183993761470242]
real(pr) :: davaldt(n) = [-8.7136629382606107E-002, -4.0735975045880973E-002]
real(pr) :: davaldt2(n) = [4.3568314691303053E-002, 1.7923829020187628E-002]

Tr = [0.4_pr, 0.5_pr]
k = [0.2_pr, 0.1_pr]
Tr = [0.4_pr, 0.5_pr]
k = [0.2_pr, 0.1_pr]

alpha%k = k
call alpha%alpha(Tr, a, dadt, dadt2)
alpha%k = k
call alpha%alpha(Tr, a, dadt, dadt2)

call check(error, maxval(abs(a - aval)) < 1e-5)
call check(error, maxval(abs(dadt - davaldt)) < 1e-5)
call check(error, maxval(abs(dadt2 - davaldt2)) < 1e-5)
end subroutine
end module
call check(error, allclose(a, aval, absolute_tolerance))
call check(error, allclose(dadt, davaldt, absolute_tolerance))
call check(error, allclose(dadt2, davaldt2, absolute_tolerance))
end subroutine test_alpha_RKPR
end module test_cubic_alphas
165 changes: 165 additions & 0 deletions test/test_cubic_implementations.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,165 @@
module test_cubic_implementations
use yaeos_constants, only: pr
use testdrive, only: new_unittest, unittest_type, error_type, check
use auxiliar_functions, only: allclose
implicit none

real(pr) :: absolute_tolerance = 1e-7_pr

contains

subroutine collect_suite(testsuite)
!> Collection of tests
type(unittest_type), allocatable, intent(out) :: testsuite(:)

testsuite = [ &
new_unittest("SRK", test_srk), &
new_unittest("PR76", test_pr76), &
new_unittest("PR78", test_pr78) &
]
end subroutine collect_suite

subroutine test_srk(error)
use yaeos_constants, only: pr
use fixtures_models, only: binary_SRK
use yaeos, only: ArModel
type(error_type), allocatable, intent(out) :: error

class(ArModel), allocatable :: eos
integer, parameter :: n = 2
real(pr) :: z(n), V, T
real(pr) :: Ar, ArV, ArV2, ArT, ArTV, ArT2
real(pr) :: Arn(n), ArVn(n), ArTn(n), Arn2(n, n)

real(pr) :: Ar_val, ArV_val, ArV2_val, ArT_val, ArTV_val, ArT2_val
real(pr) :: Arn_val(n), ArVn_val(n), ArTn_val(n), Arn2_val(n, n)

Ar_val = -9.4849231269705072
ArV_val = 9.0478810077323164
ArT_val = 3.0631941020155939E-002
ArT2_val = -1.0589478951539604E-004
ArV2_val = -17.255504598247207
ArTV_val = -3.0039878324119831E-002
Arn_val = [-14.710404803520872, -20.170975369630906]
ArVn_val = [13.488065586019152, 18.870121409429380]
ArTn_val = [5.7833039664119255E-002, 6.1888439276263030E-002]
Arn2_val(1, :) = [-11.980899399513874, -14.133993988331257]
Arn2_val(2, :) = [-14.133993988331257, -20.899890419408248]

eos = binary_SRK()
z = [0.3, 0.7]
v = 1
T = 150
call eos%residual_helmholtz( &
z, V, T, Ar=Ar, ArV=ArV, ArV2=ArV2, ArT=ArT, ArTV=ArTV, &
ArT2=ArT2, Arn=Arn, ArVn=ArVn, ArTn=ArTn, Arn2=Arn2 &
)

call check(error, allclose([Ar], [Ar_val], absolute_tolerance))
call check(error, allclose([ArV], [ArV_val], absolute_tolerance))
call check(error, allclose([ArT], [ArT_val], absolute_tolerance))
call check(error, allclose([ArTV], [ArTV_val], absolute_tolerance))
call check(error, allclose([ArV2], [ArV2_val], absolute_tolerance))
call check(error, allclose([ArT2], [ArT2_val], absolute_tolerance))

call check(error, allclose([ArVn], [ArVn_val], absolute_tolerance))
call check(error, allclose([ArTn], [ArTn_val], absolute_tolerance))
call check(error, allclose([Arn2], [Arn2_val], absolute_tolerance))
end subroutine test_srk

subroutine test_pr76(error)
use yaeos_constants, only: pr
use fixtures_models, only: binary_PR76
use yaeos, only: ArModel
type(error_type), allocatable, intent(out) :: error

class(ArModel), allocatable :: eos
integer, parameter :: n = 2
real(pr) :: z(n), V, T
real(pr) :: Ar, ArV, ArV2, ArT, ArTV, ArT2
real(pr) :: Arn(n), ArVn(n), ArTn(n), Arn2(n, n)

real(pr) :: Ar_val, ArV_val, ArV2_val, ArT_val, ArTV_val, ArT2_val
real(pr) :: Arn_val(n), ArVn_val(n), ArTn_val(n), Arn2_val(n, n)

Ar_val = -9.5079213387597061
ArV_val = 8.8348105702414230
ArT_val = 2.5288760006412853E-002
ArT2_val = -8.1263714911056052E-005
ArV2_val = -16.452712169871607
ArTV_val = -2.4354181554918298E-002
Arn_val = [-14.760083989416412, -19.878152533126190]
ArVn_val = [12.970846906902654, 17.944940224423746]
ArTn_val = [4.7299709855544367E-002, 5.0647183777961201E-002]
Arn2_val(1, :) = [-11.697767407192650, -13.516452437750393]
Arn2_val(2, :) = [-13.516452437750393, -19.842863669307611]

eos = binary_PR76()
z = [0.3, 0.7]
v = 1
T = 150
call eos%residual_helmholtz( &
z, V, T, Ar=Ar, ArV=ArV, ArV2=ArV2, ArT=ArT, ArTV=ArTV, &
ArT2=ArT2, Arn=Arn, ArVn=ArVn, ArTn=ArTn, Arn2=Arn2 &
)

call check(error, allclose([Ar], [Ar_val], absolute_tolerance))
call check(error, allclose([ArV], [ArV_val], absolute_tolerance))
call check(error, allclose([ArT], [ArT_val], absolute_tolerance))
call check(error, allclose([ArTV], [ArTV_val], absolute_tolerance))
call check(error, allclose([ArV2], [ArV2_val], absolute_tolerance))
call check(error, allclose([ArT2], [ArT2_val], absolute_tolerance))

call check(error, allclose([ArVn], [ArVn_val], absolute_tolerance))
call check(error, allclose([ArTn], [ArTn_val], absolute_tolerance))
call check(error, allclose([Arn2], [Arn2_val], absolute_tolerance))
end subroutine test_pr76

subroutine test_pr78(error)
use yaeos_constants, only: pr
use fixtures_models, only: binary_PR78
use yaeos, only: ArModel
type(error_type), allocatable, intent(out) :: error

class(ArModel), allocatable :: eos
integer, parameter :: n = 2
real(pr) :: z(n), V, T
real(pr) :: Ar, ArV, ArV2, ArT, ArTV, ArT2
real(pr) :: Arn(n), ArVn(n), ArTn(n), Arn2(n, n)

real(pr) :: Ar_val, ArV_val, ArV2_val, ArT_val, ArTV_val, ArT2_val
real(pr) :: Arn_val(n), ArVn_val(n), ArTn_val(n), Arn2_val(n, n)

Ar_val = -9.5079213387597061
ArV_val = 8.8348105702414230
ArT_val = 2.5288760006412853E-002
ArT2_val = -8.1263714911056052E-005
ArV2_val = -16.452712169871607
ArTV_val = -2.4354181554918298E-002
Arn_val = [-14.760083989416412, -19.878152533126190]
ArVn_val = [12.970846906902654, 17.944940224423746]
ArTn_val = [4.7299709855544367E-002, 5.0647183777961201E-002]
Arn2_val(1, :) = [-11.697767407192650, -13.516452437750393]
Arn2_val(2, :) = [-13.516452437750393, -19.842863669307611]

eos = binary_PR78()
z = [0.3, 0.7]
v = 1
T = 150
call eos%residual_helmholtz( &
z, V, T, Ar=Ar, ArV=ArV, ArV2=ArV2, ArT=ArT, ArTV=ArTV, &
ArT2=ArT2, Arn=Arn, ArVn=ArVn, ArTn=ArTn, Arn2=Arn2 &
)

call check(error, allclose([Ar], [Ar_val], absolute_tolerance))
call check(error, allclose([ArV], [ArV_val], absolute_tolerance))
call check(error, allclose([ArT], [ArT_val], absolute_tolerance))
call check(error, allclose([ArTV], [ArTV_val], absolute_tolerance))
call check(error, allclose([ArV2], [ArV2_val], absolute_tolerance))
call check(error, allclose([ArT2], [ArT2_val], absolute_tolerance))

call check(error, allclose([ArVn], [ArVn_val], absolute_tolerance))
call check(error, allclose([ArTn], [ArTn_val], absolute_tolerance))
call check(error, allclose([Arn2], [Arn2_val], absolute_tolerance))
end subroutine test_pr78
end module test_cubic_implementations
Loading

0 comments on commit ced15ca

Please sign in to comment.