diff --git a/CMakeLists.txt b/CMakeLists.txt index 9988cdd44b..a585af5e32 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -79,7 +79,6 @@ LIST(APPEND PROCESS_SRCS pfcoil.f90 reinke_module.f90 sctfcoil.f90 - plasma_profiles.f90 final_module.f90 cost_variables.f90 divertor_variables.f90 diff --git a/process/current_drive.py b/process/current_drive.py index 87d6e12b67..c39a039571 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -1,20 +1,22 @@ import numpy as np +from process.plasma_profiles import PlasmaProfile + from process.fortran import ( heat_transport_variables, current_drive_variables, physics_variables, cost_variables, constants, - profiles_module, process_output as po, error_handling as eh, ) class CurrentDrive: - def __init__(self): + def __init__(self, plasma_profile: PlasmaProfile): self.outfile = constants.nout + self.plasma_profile = plasma_profile def cudriv(self, output: bool): """Routine to calculate the current drive power requirements @@ -1375,7 +1377,7 @@ def cullhy(self): # Local density, temperature, toroidal field at this minor radius - dlocal = 1.0e-19 * profiles_module.nprofile( + dlocal = 1.0e-19 * self.plasma_profile.neprofile.calculate_profile_y( rratio, physics_variables.rhopedn, physics_variables.ne0, @@ -1383,7 +1385,7 @@ def cullhy(self): physics_variables.nesep, physics_variables.alphan, ) - tlocal = profiles_module.tprofile( + tlocal = self.plasma_profile.teprofile.calculate_profile_y( rratio, physics_variables.rhopedt, physics_variables.te0, @@ -1439,7 +1441,7 @@ def culecd(self): rrr = 1.0e0 / 3.0e0 # Temperature - tlocal = profiles_module.tprofile( + tlocal = self.plasma_profile.teprofile.calculate_profile_y( rrr, physics_variables.rhopedt, physics_variables.te0, @@ -1450,7 +1452,7 @@ def culecd(self): ) # Density (10**20 m**-3) - dlocal = 1.0e-20 * profiles_module.nprofile( + dlocal = 1.0e-20 * self.plasma_profile.neprofile.calculate_profile_y( rrr, physics_variables.rhopedn, physics_variables.ne0, @@ -1829,7 +1831,7 @@ def lheval(self, drfind, rratio): routine lhrad. AEA FUS 172: Physics Assessment for the European Reactor Study """ - dlocal = 1.0e-19 * profiles_module.nprofile( + dlocal = 1.0e-19 * self.plasma_profile.neprofile.calculate_profile_y( rratio, physics_variables.rhopedn, physics_variables.ne0, @@ -1840,7 +1842,7 @@ def lheval(self, drfind, rratio): # Local electron temperature - tlocal = profiles_module.tprofile( + tlocal = self.plasma_profile.teprofile.calculate_profile_y( rratio, physics_variables.rhopedt, physics_variables.te0, diff --git a/process/main.py b/process/main.py index df8d035b91..67e5bb77cb 100644 --- a/process/main.py +++ b/process/main.py @@ -616,7 +616,7 @@ def __init__(self): self.fw = Fw() self.blanket_library = BlanketLibrary(fw=self.fw) self.ccfe_hcpb = CCFE_HCPB(blanket_library=self.blanket_library) - self.current_drive = CurrentDrive() + self.current_drive = CurrentDrive(plasma_profile=self.plasma_profile) self.physics = Physics( plasma_profile=self.plasma_profile, current_drive=self.current_drive ) diff --git a/source/fortran/plasma_profiles.f90 b/source/fortran/plasma_profiles.f90 deleted file mode 100644 index 5b6bb71fc9..0000000000 --- a/source/fortran/plasma_profiles.f90 +++ /dev/null @@ -1,438 +0,0 @@ -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -module profiles_module - - !! Density and temperature profiles - !! author: P J Knight, CCFE, Culham Science Centre - !! N/A - !! This module contains routines that give the density and temperature - !! profile quantities - !! T&M/PKNIGHT/LOGBOOK24, pp.4-7 - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -#ifndef dp - use, intrinsic :: iso_fortran_env, only: dp=>real64 -#endif - private - public :: plasma_profiles, ncore, nprofile, tcore, tprofile - -contains - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine plasma_profiles - - !! Calculates density and temperature profile quantities - !! author: P J Knight, CCFE, Culham Science Centre - !! None - !! This subroutine initialises the density and temperature - !! profile averages and peak values, given the main - !! parameters describing these profiles. - !! T&M/PKNIGHT/LOGBOOK24, pp.4-7 - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use constants, only: echarge - use divertor_variables, only: prn1 - use maths_library, only: gamfun, sumup3 - use physics_variables, only: rhopedt, ten, tin, alphap, tbeta, te0, p0, & - nesep, tesep, pcoef, ipedestal, ni0, ne0, ti0, tratio, dnla, alphat, & - dnitot, neped, ti, rhopedn, dene, teped, alphan, te, rho_ne_max, & - rho_te_max, gradient_length_ne, gradient_length_te, rminor - implicit none - - ! Arguments - - ! Local variables - - integer, parameter :: nrho = 501 - integer :: irho - real(dp) :: drho, rho, integ1, integ2, dens, temp, te_max, ne_max, dndrho_max, dtdrho_max - real(dp), dimension(nrho) :: arg1, arg2, arg3 - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Volume-averaged ion temperature - ! (input value used directly if tratio=0.0) - - if (tratio > 0.0D0) ti = tratio * te - - if (ipedestal == 0) then - - ! Reset pedestal values to agree with original parabolic profiles - - rhopedt = 1.0D0 ; rhopedn = 1.0D0 - teped = 0.0D0 ; tesep = 0.0D0 - neped = 0.0D0 ; nesep = 0.0D0 - tbeta = 2.0D0 - - ! Profile factor; ratio of density-weighted to volume-averaged - ! temperature - - pcoef = (1.0D0 + alphan)*(1.0D0 + alphat)/(1.0D0+alphan+alphat) - - ! Line averaged electron density (IPDG89) - ! 0.5*gamfun(0.5) = 0.5*sqrt(pi) = 0.886227 - - dnla = dene*(1.0D0+alphan) * 0.886227D0 * gamfun(alphan+1.0D0) / & - gamfun(alphan+1.5D0) - - ! Density-weighted temperatures - - ten = te * pcoef - tin = ti * pcoef - - ! Central values for temperature (keV) and density (m**-3) - - te0 = te * (1.0D0+alphat) - ti0 = ti * (1.0D0+alphat) - - ne0 = dene * (1.0D0+alphan) - ni0 = dnitot * (1.0D0+alphan) - - else - - ! The following reproduces the above results within sensible - ! tolerances if rhopedt = rhopedn = 1.0, teped = tesep = neped - ! = nesep = 0.0, and tbeta = 2.0 - - ! Central values for temperature (keV) and density (m**-3) - - te0 = tcore(rhopedt,teped,tesep,te,alphat,tbeta) - ti0 = ti/te * te0 - - ne0 = ncore(rhopedn,neped,nesep,dene,alphan) - ni0 = dnitot/dene * ne0 - - ! Perform integrations to calculate ratio of density-weighted - ! to volume-averaged temperature, etc. - ! Density-weighted temperature = integral(n.T dV) / integral(n dV) - ! which is approximately equal to the ratio - ! integral(rho.n(rho).T(rho) drho) / integral(rho.n(rho) drho) - - drho = 1.0D0/(nrho-1) - do irho = 1,nrho - rho = (irho-1.0D0)/(nrho-1) - dens = nprofile(rho,rhopedn,ne0,neped,nesep,alphan) - temp = tprofile(rho,rhopedt,te0,teped,tesep,alphat,tbeta) - arg1(irho) = rho*dens*temp - arg2(irho) = rho*dens - arg3(irho) = dens - end do - call sumup3(drho,arg1,integ1,nrho) - call sumup3(drho,arg2,integ2,nrho) - - ! Density-weighted temperatures - - ten = integ1/integ2 - tin = ti/te * ten - - ! Profile factor; ratio of density-weighted to volume-averaged - ! temperature - - pcoef = ten / te - - ! Line-averaged electron density - ! = integral(n(rho).drho) - - call sumup3(drho,arg3,dnla,nrho) - - ! Scrape-off density / volume averaged density - ! (Input value is used if ipedestal = 0) - - prn1 = max(0.01D0, nesep/dene) ! preventing division by zero later - - end if - - ! Central pressure (Pa), from ideal gas law : p = nkT - - p0 = (ne0*te0 + ni0*ti0) * 1.0D3 * echarge - - ! Pressure profile index (N.B. no pedestal effects included here) - ! N.B. p0 is NOT equal to

* (1 + alphap), but p(rho) = n(rho)*T(rho) - ! and

= .T_n where <...> denotes volume-averages and T_n is the - ! density-weighted temperature - - alphap = alphan + alphat - - ! The gradient information for ipedestal = 0: - ! All formulas can be obtained from the analytical parametric form of the ipedestal profiles - ! rho_max is obtained by equalling the second derivative to zero e.g. - - if (ipedestal == 0) then - if(alphat > 1.0d0) then - ! Rho (normalized radius), where temperature derivative is largest - rho_te_max = 1.0d0/sqrt(-1.0d0 +2.0d0 * alphat) - dtdrho_max = -2.0d0**alphat*(-1.0d0 + alphat)**(-1.0d0 + alphat)*alphat*(-1.0d0 + & - 2.0d0 * alphat)**(0.5d0 - alphat)*te0 - te_max = te0*(1d0 - rho_te_max**2)**alphat - elseif (alphat .le. 1.0d0 .and. alphaT > 0.0d0) then - ! This makes the profiles very 'boxy' - ! The gradient diverges here at the edge so define some 'wrong' value of 0.9 - ! to approximate the gradient - rho_te_max = 0.9d0 - dtdrho_max = -2.0d0 * alphat * rho_te_max*(1.d0-rho_te_max**2)**(-1.0d0+alphat)*te0 - te_max = te0*(1d0 - rho_te_max**2)**alphat - else - print *, "ERROR: alphat is negative!" - call exit(1) - end if - - ! Same for density - if(alphan > 1.0d0) then - rho_ne_max = 1.0d0/sqrt(-1.0d0 +2.0d0 * alphan) - dndrho_max = -2.0d0**alphan*(-1.0d0 + alphan)**(-1.0d0 + alphan)*alphan*(-1.0d0 + & - 2.0d0 * alphan)**(0.5d0 - alphan)*ne0 - ne_max = ne0*(1d0 - rho_ne_max**2)**alphan - elseif (alphan .le. 1.0d0 .and. alphan > 0.0d0) then - ! This makes the profiles very 'boxy' - ! The gradient diverges here at the edge so define some 'wrong' value of 0.9 - ! to approximate the gradient - rho_ne_max = 0.9d0 - dndrho_max = -2.0d0 * alphan * rho_ne_max*(1.d0-rho_ne_max**2)**(-1.0d0+alphan)*ne0 - ne_max = ne0*(1d0 - rho_ne_max**2)**alphan - else - print *, "ERROR: alphan is negative!" - call exit(1) - end if - - ! set normalized gradient length - ! te at rho_te_max - gradient_length_te = -dtdrho_max * rminor*rho_te_max/te_max - ! same for density: - gradient_length_ne = -dndrho_max * rminor*rho_ne_max/ne_max - - end if - - - - end subroutine plasma_profiles - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function tcore(rhopedt, tped, tsep, tav, alphat, tbeta) - - !! Central temperature for a pedestal profile - !! author: R Kemp, CCFE, Culham Science Centre - !! author: H Lux, CCFE, Culham Science Centre - !! author: P J Knight, CCFE, Culham Science Centre - !! rhopedt : input real : normalised minor radius pedestal position - !! tped : input real : pedestal temperature (keV) - !! tsep : input real : separatrix temperature (keV) - !! tav : input real : volume average temperature (keV) - !! alphat : input real : temperature peaking parameter - !! tbeta : input real : second temperature exponent - !! This routine calculates the core temperature (keV) - !! of a pedestalised profile. - !! J.Johner, Fusion Science and Technology 59 (2011), pp 308-349 - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use constants, only: pi - use maths_library, only: gamfun - implicit none - - real(dp) :: tcore - - ! Arguments - - real(dp), intent(in) :: rhopedt, tped, tsep, tav, alphat, tbeta - - ! Local variables - - real(dp), parameter :: numacc = 1.0D-7 - real(dp) :: gamfac - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! For integer values of alphat, the limit of - ! gamfun(-alphat)*sin(pi*alphat) needs to be calculated directly - - gamfac = gamfun(1.0D0 + alphat + 2.0D0/tbeta) / & - gamfun((2.0D0 + tbeta)/tbeta) / rhopedT**2 - - if (abs(alphat-nint(alphat)) <= numacc) then - gamfac = -gamfac / gamfun(1.0D0 + alphat) - else - gamfac = gamfac * gamfun(-alphat)*sin(pi*alphat)/pi - end if - - ! Calculate core temperature - - tcore = tped + gamfac * ( tped*rhopedT**2 - tav + & - (1.0D0 - rhopedT)/3.0D0 * & - ( (1.0D0 + 2.0D0*rhopedT)*tped + (2.0D0 + rhopedT)*tsep ) ) - - end function tcore - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function tprofile(rho, rhopedt, t0, tped, tsep, alphat, tbeta) - - !! Implementation of HELIOS-type temperature pedestal profile - !! author: R Kemp, CCFE, Culham Science Centre - !! author: H Lux, CCFE, Culham Science Centre - !! author: P J Knight, CCFE, Culham Science Centre - !! rho : input real : normalised minor radius - !! rhopedt : input real : normalised minor radius pedestal position - !! t0 : input real : central temperature (keV) - !! tped : input real : pedestal temperature (keV) - !! tsep : input real : separatrix temperature (keV) - !! alphat : input real : temperature peaking parameter - !! tbeta : input real : second temperature exponent - !! trho : output real : T(rho) (keV) - !! This routine calculates the temperature at a normalised minor - !! radius position rho for a pedestalised profile. - !!

If ipedestal = 0 the original parabolic - !! profile form is used instead. - !! J.Johner, Fusion Science and Technology 59 (2011), pp 308-349 - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use error_handling, only: fdiags, report_error - use physics_variables, only: ipedestal - implicit none - - real(dp) :: tprofile - - ! Arguments - - real(dp), intent(in) :: rho, rhopedt, t0, tped, tsep, alphat, tbeta - - ! Local variables - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - if (ipedestal == 0) then - tprofile = t0 * (1.0D0 - rho**2)**alphat - return - end if - - ! Error trap; shouldn't happen unless volume-averaged temperature has - ! been allowed to drop below tped. This may happen during a HYBRD case, - ! but should have been prevented for optimisation runs. - - if (t0 < tped) then - fdiags(1) = tped ; fdiags(2) = t0 - call report_error(148) - end if - - if (rho <= rhopedt) then - tprofile = tped + (t0 - tped) * (1.0D0 - (rho/rhopedT)**tbeta)**alphat - else - tprofile = tsep + (tped - tsep) * (1.0D0 - rho)/(1.0D0 - rhopedt) - end if - - end function tprofile - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function ncore(rhopedn, nped, nsep, nav, alphan) - - !! Central density of a pedestal profile - !! author: R Kemp, CCFE, Culham Science Centre - !! author: H Lux, CCFE, Culham Science Centre - !! author: P J Knight, CCFE, Culham Science Centre - !! rhopedn : input real : normalised minor radius pedestal position - !! nped : input real : pedestal density (/m3) - !! nsep : input real : separatrix density (/m3) - !! nav : input real : volume average density (/m3) - !! alphan : input real : density peaking parameter - !! This routine calculates the central density - !! of a pedestalised profile. - !! J.Johner, Fusion Science and Technology 59 (2011), pp 308-349 - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use error_handling, only: report_error - implicit none - - real(dp) :: ncore - - ! Arguments - - real(dp), intent(in) :: rhopedn, nped, nsep, nav, alphan - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ncore = 1.0D0 / (3.0D0*rhopedn**2) * ( & - 3.0D0*nav*(1.0D0 + alphan) & - + nsep*(1.0D0 + alphan)*(-2.0D0 + rhopedn + rhopedn**2) & - - nped*( (1.0D0 + alphan)*(1.0D0 + rhopedn) + & - (alphan - 2.0D0)*rhopedn**2 ) ) - - if (ncore < 1.0D-6) then - ! Prevent ncore from going negative (and terminating the optimisation) by - ! kludging to small positive value. Allows solver to continue and - ! hopefully be constrained away from this point (e.g. constraint 81, ne0 > neped) - ncore = 1.0D-6 - ! write(*,*) 'Error in ncore: negative central density' - ! write(*,*) 'nped =', nped, ' nsep =', nsep - ! write(*,*) 'nav =', nav, ' ncore =', ncore - ! call report_error(212) - end if - - end function ncore - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function nprofile(rho, rhopedn, n0, nped, nsep, alphan) - - !! Implementation of HELIOS-type density pedestal profile - !! author: R Kemp, CCFE, Culham Science Centre - !! author: H Lux, CCFE, Culham Science Centre - !! author: P J Knight, CCFE, Culham Science Centre - !! rho : input real : normalised minor radius - !! rhopedn : input real : normalised minor radius pedestal position - !! n0 : input real : central density (/m3) - !! nped : input real : pedestal density (/m3) - !! nsep : input real : separatrix density (/m3) - !! alphan : input real : density peaking parameter - !! This routine calculates the density at a normalised minor - !! radius position rho for a pedestalised profile. - !!

If ipedestal = 0 the original parabolic - !! profile form is used instead. - !! J.Johner, Fusion Science and Technology 59 (2011), pp 308-349 - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use error_handling, only: fdiags, report_error - use physics_variables, only: ipedestal - implicit none - - real(dp) :: nprofile - - ! Arguments - - real(dp), intent(in) :: rho, rhopedn, n0, nped, nsep, alphan - - ! Local variables - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - if (ipedestal == 0) then - nprofile = n0 * (1.0D0 - rho**2)**alphan - return - end if - - ! Error trap; shouldn't happen unless volume-averaged density has - ! been allowed to drop below nped. This may happen during a HYBRD case, - ! but should have been prevented for optimisation runs. - - ! Input checks - - if (n0 < nped) then - fdiags(1) = nped ; fdiags(2) = n0 - call report_error(153) - end if - - if (rho <= rhopedn) then - nprofile = nped + (n0 - nped) * (1.0D0 - (rho/rhopedn)**2)**alphan - else - nprofile = nsep + (nped - nsep) * (1.0D0 - rho)/(1.0D0 - rhopedn) - end if - - end function nprofile - -end module profiles_module diff --git a/tests/unit/test_current_drive.py b/tests/unit/test_current_drive.py index d7b88d011f..9a88a102fc 100644 --- a/tests/unit/test_current_drive.py +++ b/tests/unit/test_current_drive.py @@ -9,6 +9,7 @@ heat_transport_variables, ) from process.current_drive import CurrentDrive +from process.plasma_profiles import PlasmaProfile @pytest.fixture @@ -18,7 +19,7 @@ def current_drive(): :returns current_drive: initialised CurrentDrive object :rtype: process.current_drive.CurrentDrive """ - return CurrentDrive() + return CurrentDrive(PlasmaProfile()) class CudrivParam(NamedTuple): diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 402a67322f..ca43a7cffd 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -36,7 +36,7 @@ def physics(): :returns: initialised Physics object :rtype: process.physics.Physics """ - return Physics(PlasmaProfile(), CurrentDrive()) + return Physics(PlasmaProfile(), CurrentDrive(PlasmaProfile())) def test_beta_poloidal(): diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 06ffcb6ed9..207f58e28d 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -42,8 +42,8 @@ def stellarator(): Power(), PlasmaProfile(), CCFE_HCPB(BlanketLibrary(Fw())), - CurrentDrive(), - Physics(PlasmaProfile(), CurrentDrive()), + CurrentDrive(PlasmaProfile()), + Physics(PlasmaProfile(), CurrentDrive(PlasmaProfile())), )