From a1665f1d234a0b699803be3f9662d59739ea1cbe Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 3 Jul 2024 14:44:31 +0100 Subject: [PATCH 1/7] Refactor plasma_profiles.py and profiles.py to use scipy.special.gamma instead of maths_library.gamfun --- process/plasma_profiles.py | 8 +++----- process/profiles.py | 13 ++++++------- 2 files changed, 9 insertions(+), 12 deletions(-) diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py index cae753ed08..2d83703d12 100644 --- a/process/plasma_profiles.py +++ b/process/plasma_profiles.py @@ -6,7 +6,6 @@ from process.fortran import ( constants, divertor_variables, - maths_library, physics_variables, ) @@ -96,14 +95,13 @@ def parabolic_paramterisation(self): ) # Line averaged electron density (IPDG89) - # 0.5*gamfun(0.5) = 0.5*sqrt(pi) = 0.886227 - + physics_variables.dnla = ( physics_variables.dene * (1.0 + physics_variables.alphan) * 0.886227 - * maths_library.gamfun(physics_variables.alphan + 1.0) - / maths_library.gamfun(physics_variables.alphan + 1.5e0) + * sp.special.gamma(physics_variables.alphan + 1.0) + / sp.special.gamma(physics_variables.alphan + 1.5) ) # Density-weighted temperatures diff --git a/process/profiles.py b/process/profiles.py index 58d2efbd63..37484629a8 100644 --- a/process/profiles.py +++ b/process/profiles.py @@ -1,9 +1,10 @@ import numpy as np import logging from scipy import integrate +import scipy as sp from abc import ABC, abstractmethod -from process.fortran import maths_library, physics_variables, error_handling +from process.fortran import physics_variables, error_handling logger = logging.getLogger(__name__) # Logging handler for console output @@ -261,19 +262,17 @@ def tcore(rhopedt, teped, tesep, tav, alphat, tbeta): :return: core temperature :rtype: numpy.array """ - # For integer values of alphat, the limit of - # gamfun(-alphat)*sin(pi*alphat) needs to be calculated directly gamfac = ( - maths_library.gamfun(1 + alphat + 2 / tbeta) - / maths_library.gamfun((2 + tbeta) / tbeta) + sp.special.gamma(1 + alphat + 2 / tbeta) + / sp.special.gamma((2 + tbeta) / tbeta) / rhopedt**2 ) if abs(alphat - np.around(alphat)) <= 1e-7: - gamfac = -gamfac / maths_library.gamfun(1 + alphat) + gamfac = -gamfac / sp.special.gamma(1 + alphat) else: gamfac = ( - gamfac * maths_library.gamfun(-alphat) * np.sin(np.pi * alphat) / np.pi + gamfac * sp.special.gamma(-alphat) * np.sin(np.pi * alphat) / np.pi ) # Calculate core temperature From 34219f9b1d5ce3735d6a99b0f7cfddb5ab0aa760 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 3 Jul 2024 15:24:38 +0100 Subject: [PATCH 2/7] Refactor maths_library.f90 to remove gamfun function and use scipy.special.gamma instead --- source/fortran/maths_library.f90 | 55 -------------------------------- 1 file changed, 55 deletions(-) diff --git a/source/fortran/maths_library.f90 b/source/fortran/maths_library.f90 index 57ab63b2f2..5bd4770ad2 100644 --- a/source/fortran/maths_library.f90 +++ b/source/fortran/maths_library.f90 @@ -336,61 +336,6 @@ function binomial(n,k) result(coefficient) end function binomial ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - recursive function gamfun(x) result(gamma) - - !! Calculates the gamma function for arbitrary real x - !! author: P J Knight, CCFE, Culham Science Centre - !! x : input real : gamma function argument - !! This routine evaluates the gamma function, using an - !! asymptotic expansion based on Stirling's approximation. - !! http://en.wikipedia.org/wiki/Gamma_function - !! T&M/PKNIGHT/LOGBOOK24, p.5 - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - implicit none - - ! Arguments - - real(dp), intent(in) :: x - real(dp) :: gamma - - ! Local variables - - real(dp), parameter :: sqtwopi = 2.5066282746310005D0 - real(dp), parameter :: c1 = 8.3333333333333333D-2 ! 1/12 - real(dp), parameter :: c2 = 3.4722222222222222D-3 ! 1/288 - real(dp), parameter :: c3 = 2.6813271604938272D-3 ! 139/51840 - real(dp), parameter :: c4 = 2.2947209362139918D-4 ! 571/2488320 - real(dp) :: summ, denom - integer :: i,n - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - if (x > 1.0D0) then - - summ = 1.0D0 + c1/x + c2/x**2 - c3/x**3 - c4/x**4 - gamma = exp(-x) * x**(x-0.5D0) * sqtwopi * summ - - else - - ! Use recurrence formula to shift the argument to >1 - ! gamma(x) = gamma(x+n) / (x*(x+1)*(x+2)*...*(x+n-1)) - ! where n is chosen to make x+n > 1 - - n = int(-x) + 2 - denom = x - do i = 1,n-1 - denom = denom*(x+i) - end do - gamma = gamfun(x+n)/denom - - end if - - end function gamfun - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - function binarysearch(length, array, value, delta) ! Given an array and a value, returns the index of the element that ! is closest to, but less than, the given value. From 2612326ef4193eef1bd65874e7aae12f63621006 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 3 Jul 2024 15:28:06 +0100 Subject: [PATCH 3/7] Black format --- process/plasma_profiles.py | 2 +- process/profiles.py | 4 +--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py index 2d83703d12..d69f96724d 100644 --- a/process/plasma_profiles.py +++ b/process/plasma_profiles.py @@ -95,7 +95,7 @@ def parabolic_paramterisation(self): ) # Line averaged electron density (IPDG89) - + physics_variables.dnla = ( physics_variables.dene * (1.0 + physics_variables.alphan) diff --git a/process/profiles.py b/process/profiles.py index 37484629a8..3ad8184e94 100644 --- a/process/profiles.py +++ b/process/profiles.py @@ -271,9 +271,7 @@ def tcore(rhopedt, teped, tesep, tav, alphat, tbeta): if abs(alphat - np.around(alphat)) <= 1e-7: gamfac = -gamfac / sp.special.gamma(1 + alphat) else: - gamfac = ( - gamfac * sp.special.gamma(-alphat) * np.sin(np.pi * alphat) / np.pi - ) + gamfac = gamfac * sp.special.gamma(-alphat) * np.sin(np.pi * alphat) / np.pi # Calculate core temperature From d01a326dd85ecf27dad4b7a1db0379e8476c5c41 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 4 Jul 2024 10:35:25 +0100 Subject: [PATCH 4/7] Update unit tests --- tests/unit/test_plasma_profiles.py | 44 +++++++++++++++--------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/tests/unit/test_plasma_profiles.py b/tests/unit/test_plasma_profiles.py index c839b8fabd..e428d63562 100644 --- a/tests/unit/test_plasma_profiles.py +++ b/tests/unit/test_plasma_profiles.py @@ -86,15 +86,15 @@ class TProfileParam(NamedTuple): alphat=1.45, teped=5.5, expected_tprofile=[ - 18.84880504, - 18.73828506, - 18.40251429, - 17.82801477, - 16.98907408, - 15.84091405, - 14.3037138, - 12.21867652, - 9.17694539, + 18.85, + 18.739472621498333, + 18.403679368327712, + 17.829141392239833, + 16.990144534125456, + 15.841907634203302, + 14.30460446612375, + 12.219427594826554, + 9.177492824141694, 1.0, ], ), @@ -118,7 +118,7 @@ def test_tcore(): tbeta = 2.0 assert tprofile.tcore(rhopedt, tped, tsep, tav, alphat, tbeta) == pytest.approx( - 28.089723663920328 + 28.09093632260765 ) @@ -246,15 +246,15 @@ class PlasmaProfilesParam(NamedTuple): gradient_length_te=0.0, rminor=2.9264516129032256, expected_prn1=0.45623618297262686, - expected_ten=14.521871327399182, - expected_tin=14.521871327399182, + expected_ten=14.52233022043558, + expected_tin=14.52233022043558, expected_alphap=2.4500000000000002, - expected_te0=27.369013322953624, - expected_p0=868071.46874220832, - expected_pcoef=1.1110842637642833, + expected_te0=27.370104119511087, + expected_p0=868106.0658743214, + expected_pcoef=1.111119374172577, expected_ni0=9.210720071916929e19, expected_ne0=1.0585658890823703e20, - expected_ti0=27.369013322953624, + expected_ti0=27.370104119511087, expected_dnla=8.8687354645836431e19, expected_ti=13.07, ), @@ -291,15 +291,15 @@ class PlasmaProfilesParam(NamedTuple): gradient_length_te=0.0, rminor=2.9264516129032256, expected_prn1=0.45623618297262686, - expected_ten=14.521871327399182, - expected_tin=14.521871327399182, + expected_ten=14.52233022043558, + expected_tin=14.52233022043558, expected_alphap=2.4500000000000002, - expected_te0=27.369013322953624, - expected_p0=868071.46874220832, - expected_pcoef=1.1110842637642833, + expected_te0=27.370104119511087, + expected_p0=868106.0658743214, + expected_pcoef=1.111119374172577, expected_ni0=9.210720071916929e19, expected_ne0=1.0585658890823703e20, - expected_ti0=27.369013322953624, + expected_ti0=27.370104119511087, expected_dnla=8.8687354645836431e19, expected_ti=13.07, ), From 735b6143152960cd538caf2b8d1e581e963f59d8 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 4 Jul 2024 10:49:13 +0100 Subject: [PATCH 5/7] Revert "Refactor maths_library.f90 to remove gamfun function and use scipy.special.gamma instead" This reverts commit 34219f9b1d5ce3735d6a99b0f7cfddb5ab0aa760. --- source/fortran/maths_library.f90 | 55 ++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/source/fortran/maths_library.f90 b/source/fortran/maths_library.f90 index 5bd4770ad2..57ab63b2f2 100644 --- a/source/fortran/maths_library.f90 +++ b/source/fortran/maths_library.f90 @@ -336,6 +336,61 @@ function binomial(n,k) result(coefficient) end function binomial ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + recursive function gamfun(x) result(gamma) + + !! Calculates the gamma function for arbitrary real x + !! author: P J Knight, CCFE, Culham Science Centre + !! x : input real : gamma function argument + !! This routine evaluates the gamma function, using an + !! asymptotic expansion based on Stirling's approximation. + !! http://en.wikipedia.org/wiki/Gamma_function + !! T&M/PKNIGHT/LOGBOOK24, p.5 + ! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + implicit none + + ! Arguments + + real(dp), intent(in) :: x + real(dp) :: gamma + + ! Local variables + + real(dp), parameter :: sqtwopi = 2.5066282746310005D0 + real(dp), parameter :: c1 = 8.3333333333333333D-2 ! 1/12 + real(dp), parameter :: c2 = 3.4722222222222222D-3 ! 1/288 + real(dp), parameter :: c3 = 2.6813271604938272D-3 ! 139/51840 + real(dp), parameter :: c4 = 2.2947209362139918D-4 ! 571/2488320 + real(dp) :: summ, denom + integer :: i,n + + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + if (x > 1.0D0) then + + summ = 1.0D0 + c1/x + c2/x**2 - c3/x**3 - c4/x**4 + gamma = exp(-x) * x**(x-0.5D0) * sqtwopi * summ + + else + + ! Use recurrence formula to shift the argument to >1 + ! gamma(x) = gamma(x+n) / (x*(x+1)*(x+2)*...*(x+n-1)) + ! where n is chosen to make x+n > 1 + + n = int(-x) + 2 + denom = x + do i = 1,n-1 + denom = denom*(x+i) + end do + gamma = gamfun(x+n)/denom + + end if + + end function gamfun + + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + function binarysearch(length, array, value, delta) ! Given an array and a value, returns the index of the element that ! is closest to, but less than, the given value. From 25512e57febb83586211f2f903acbdd4dc58bda2 Mon Sep 17 00:00:00 2001 From: Christopher Ashe- UKAEA <91618944+chris-ashe@users.noreply.github.com> Date: Fri, 5 Jul 2024 15:15:39 +0100 Subject: [PATCH 6/7] Update process/profiles.py Co-authored-by: Timothy <75321887+timothy-nunn@users.noreply.github.com> --- process/profiles.py | 1 - 1 file changed, 1 deletion(-) diff --git a/process/profiles.py b/process/profiles.py index 3ad8184e94..2d8585d4b0 100644 --- a/process/profiles.py +++ b/process/profiles.py @@ -1,6 +1,5 @@ import numpy as np import logging -from scipy import integrate import scipy as sp from abc import ABC, abstractmethod From 92ab6f59ee7b9f22f4224d184d160157ad55d47f Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 8 Jul 2024 09:30:25 +0100 Subject: [PATCH 7/7] chore: Update scipy import in profiles.py --- process/profiles.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/process/profiles.py b/process/profiles.py index 2d8585d4b0..37355877b2 100644 --- a/process/profiles.py +++ b/process/profiles.py @@ -49,7 +49,7 @@ def integrate_profile_y(self): """ Integrate profile_y values using scipy.integrate.simpson() function. """ - self.profile_integ = integrate.simpson( + self.profile_integ = sp.integrate.simpson( self.profile_y, x=self.profile_x, dx=self.profile_dx )