module carthy_equilibrium_m !! Module containing the implementation of an analytic tokamak equilibrium as !! described by Carthy (Phys. Plasmas 1999 https://doi.org/10.1063/1.873630) !! N.b. this equilibrium is axisymmetric - phi, when required, is not used # ifdef __INTEL_COMPILER use ifport # endif use precision_m, only : FP use screen_io_m, only : get_stdout, get_stderr use comm_handling_m, only: is_master use constants_m, only : pi use descriptors_m use analytic_divertor_equilibrium_m, only: & analytic_divertor_equilibrium_t, max_polygon_N_pts use elementary_functions_m, only : smoothstep => step_hermite implicit none private ! select the correct bessel function based on the selected precision #ifdef DOUBLE_PREC # if defined __NVCOMPILER || defined _CRAYFTN # define XBESJ0 bessel_j0 # define XBESJ1 bessel_j1 # define XBESY0 bessel_y0 # define XBESY1 bessel_y1 # else # define XBESJ0 DBesJ0 # define XBESJ1 DBesJ1 # define XBESY0 DBesY0 # define XBESY1 DBesY1 # endif #else # if defined __NVCOMPILER || defined _CRAYFTN # define XBESJ0 bessel_j0 # define XBESJ1 bessel_j1 # define XBESY0 bessel_y0 # define XBESY1 bessel_y1 # else # define XBESJ0 BesJ0 # define XBESJ1 BesJ1 # define XBESY0 BesY0 # define XBESY1 BesY1 # endif #endif type, public, extends(analytic_divertor_equilibrium_t) :: carthy_t contains procedure, public :: init procedure, public :: display procedure, public :: debug procedure, public :: psi procedure, public :: bx procedure, public :: by procedure, public :: btor end type carthy_t contains subroutine init(self, filename) ! Initialisation function for carthy equilibrium class(carthy_t), intent(inout) :: self character(len=*), intent(in), optional :: filename ! Path to params.in file integer :: io_error real(FP) :: axis_Btor real(FP) :: xmin, xmax, ymin, ymax, rhomin, rhomax, rhomin_privflux character(len=:), allocatable :: district_definition integer :: divertor_polygon_N_pts, exclusion_polygon_N_pts real(FP), dimension(max_polygon_N_pts) :: divertor_polygon_X_pts real(FP), dimension(max_polygon_N_pts) :: divertor_polygon_Y_pts real(FP), dimension(max_polygon_N_pts) :: exclusion_polygon_X_pts real(FP), dimension(max_polygon_N_pts) :: exclusion_polygon_Y_pts integer :: iunit = 35 character(len=100) :: line namelist / equi_carthy_params / axis_Btor, & xmin, xmax, ymin, ymax, rhomin, rhomax, rhomin_privflux, & divertor_polygon_N_pts, divertor_polygon_X_pts, divertor_polygon_Y_pts, & exclusion_polygon_N_pts, exclusion_polygon_X_pts, exclusion_polygon_Y_pts, & district_definition ! Constant equilibrium parameters self%R0 = 1.7686_FP self%Z0 = -0.2309_FP self%O_point_psi = 1.931594598_FP self%X_point_psi = 0.0_FP self%RX = 0.848237_FP self%ZX = -0.378087_FP ! Default values, in case not defined in namelist axis_Btor = 5.0_FP xmin = 0.65_FP xmax = 1.25_FP ymin = -0.56_FP ymax = 0.35_FP rhomin = 0.96_FP rhomax = 1.04_FP rhomin_privflux = 0.99_FP district_definition = 'flux' divertor_polygon_N_pts = 8 divertor_polygon_X_pts(1:divertor_polygon_N_pts) = & (/ +0.745_FP, +0.831_FP, +0.966_FP, +10.00_FP, +10.0_FP, +0.00_FP, +0.000_FP, +0.745_FP /) divertor_polygon_Y_pts(1:divertor_polygon_N_pts) = & (/ -0.372_FP, -0.494_FP, -0.418_FP, -0.418_FP, +10.0_FP, +10.0_FP, -0.372_FP, -0.372_FP /) exclusion_polygon_N_pts = 8 exclusion_polygon_X_pts(1:exclusion_polygon_N_pts) = & (/ +0.680_FP, +0.831_FP, +0.966_FP, +10.00_FP, +10.0_FP, +0.00_FP, +0.000_FP, +0.680_FP /) exclusion_polygon_Y_pts(1:exclusion_polygon_N_pts) = & (/ -0.422_FP, -0.594_FP, -0.518_FP, -0.518_FP, +10.0_FP, +10.0_FP, -0.472_FP, -0.422_FP /) self%invert_divertor_polygon = .false. self%invert_exclusion_polygon = .false. if (present(filename)) then open(unit = iunit, file = filename, status = 'old', action = 'read', iostat = io_error) if (io_error /= 0) then write(get_stderr(), *) & 'error(equi%init): file', filename, 'not existent' ERROR STOP ERR_PARAMETER_FILE endif read(iunit, nml = equi_carthy_params, iostat = io_error) if (io_error == -1) then write(get_stderr(),'(A)') & 'End-of-file reached: either missing a parameter group, & &or missing a blank line at the end of file' write(get_stderr(), *) & 'error(equi%init): erroneous namelist', io_error ERROR STOP ERR_PARAMETER_FILE elseif(io_error /= 0) then backspace(iunit) read(iunit,fmt='(A)') line write(get_stderr(),'(A)') & 'Invalid line in namelist: '//trim(line) write(get_stderr(), *) & 'error(equi%init): erroneous namelist', io_error ERROR STOP ERR_PARAMETER_FILE endif close(iunit) endif self%axis_Btor = axis_Btor self%x0 = 1.0_FP self%y0 = 0.0_FP self%xmin = xmin self%xmax = xmax self%ymin = ymin self%ymax = ymax self%rhomin = rhomin self%rhomax = rhomax self%rhomin_privflux = rhomin_privflux self%district_definition = district_definition call self%make_polygon_from_params(divertor_polygon_N_pts, divertor_polygon_X_pts, divertor_polygon_Y_pts, & self%divertor_polygon) call self%make_polygon_from_params(exclusion_polygon_N_pts, exclusion_polygon_X_pts, exclusion_polygon_Y_pts, & self%exclusion_polygon) self%initialized = .true. end subroutine init subroutine display(self) ! Print to stdout information for carthy equilibrium class(carthy_t), intent(in) :: self write(get_stdout(), 99) & self%rhomin, self%rhomax, self%rhomin_privflux, & self%district_definition, & trim(self%divertor_polygon%polygon_string()), & trim(self%exclusion_polygon%polygon_string()) 99 FORMAT( & 3X,"CARTHY equilibrium", /, & 3X,'rhomin =', ES14.6E3, /, & 3X,'rhomax =', ES14.6E3, /, & 3X,'rhomin_privflux =', ES14.6E3, /, & 3X,'district_definition = ', A32, /, & 1X, '<<Divertor polygon>>'/, & 1X, A /, & 1X, '<<Exclusion polygon>>' /, & 1X, A & /) end subroutine display subroutine debug(self) class(carthy_t), intent(in) :: self write(get_stdout(), 100) & self%R0, self%Z0, self%axis_Btor, self%RX, self%ZX, & self%x0, self%y0, self%xmin, self%xmax, self%ymin, self%ymax, & self%rhomin, self%rhomax, self%rhomin_privflux, & trim(self%divertor_polygon%polygon_string()), & trim(self%exclusion_polygon%polygon_string()) 100 FORMAT( & 1X,'Information about equilibrium:' , /, & 3X,'Equilibrium from P.J. Mc Carthy, Physics of Plasmas 6, p. 3554, (1999)' , /,& 3X,'r0 =',ES14.6E3 , /, & 3X,'z0 =',ES14.6E3 , /, & 3X,'B0 =',ES14.6E3 , /, & 3X,'rx =',ES14.6E3 , /, & 3X,'zx =',ES14.6E3 , /, & 3X,'x0 =',ES14.6E3 , /, & 3X,'y0 =',ES14.6E3 , /, & 3X,'xmin =',ES14.6E3 , /, & 3X,'xmax =',ES14.6E3 , /, & 3X,'ymin =',ES14.6E3 , /, & 3X,'ymax =',ES14.6E3 , /, & 3X,'rhomin =',ES14.6E3 , /, & 3X,'rhomax =',ES14.6E3 , /, & 3X,'rhomin_privflux =',ES14.6E3 , /, & 1X, '<<Divertor polygon>>' /, & 1X, A /, & 1X, '<<Exclusion polygon>>' /, & 1X, A & / ) end subroutine debug real(FP) function psi(self, x, y, phi) class(carthy_t), intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP) :: R, Z R = x*self%R0 Z = self%R0*(y+self%Z0/self%R0) psi = 1.5614_FP & - 1.14118_FP * R**2 & - 0.145736_FP * R * XBESJ1(3.3_FP*R) & - 0.650909_FP * R * Z * XBESJ1(3.3_FP*R) & - 3.69755_FP * R * XBESY1(3.3_FP*R) & - 2.55895_FP * R * XBESJ1(2.64_FP*R) * cos(1.98_FP*Z) & + 2.00008_FP * R * XBESJ1(1.98_FP*R) * cos(2.64_FP*Z) & + 1.2129_FP * R * XBESJ1(0.33_FP*R) * cos(3.2834585424518457_FP*Z) & - 2.13022_FP * R * XBESJ1(3.8587562763149474_FP*R) * cosh(2.0_FP*Z) & + 1.32363_FP * R * XBESJ1(2.64_FP*R) * sin(1.98_FP*Z) & - 0.910325_FP * sin(3.3_FP*Z) end function psi real(FP) function bx(self, x, y, phi) real(FP), intent(in) :: x, y, phi class(carthy_t), intent(in) :: self real(FP) :: R, Z R = x*self%R0 Z = self%R0*(y+self%Z0/self%R0) bx = 1.0/(2.0_FP*pi*R)*( & - 0.650909_FP * R * XBESJ1(3.3_FP*R) & + 5.0667210_FP * R * XBESJ1(2.64_FP*R) * sin(1.98_FP*Z) & - 5.2802112_FP * R * XBESJ1(1.98_FP*R) * sin(2.64_FP*Z) & - 3.982506866_FP * R * XBESJ1(0.33_FP*R) * sin(3.2834585424518457_FP*Z) & - 4.26044_FP * R * XBESJ1(3.8587562763149474_FP*R) * sinh(2.0_FP*Z) & + 2.6207874_FP * R * XBESJ1(2.64_FP*R) * cos(1.98_FP*Z) & - 3.0040725_FP * cos(3.3_FP*Z) & ) bx = bx / self%axis_Btor end function bx real(FP) function by(self, x, y, phi) real(FP), intent(in) :: x, y, phi class(carthy_t), intent(in) :: self real(FP)::R,Z R = x*self%R0 Z = self%R0*(y+self%Z0/self%R0) by= -1.0/(2.0_FP*pi*R)*( & - 2.28236_FP * R & - 0.145736_FP * XBESJ1(3.3_FP*R) & - 0.145736_FP * R * ( 3.3_FP*XBESJ0(3.3_FP*R) - XBESJ1(3.3_FP*R)/R ) & - 0.650909_FP * Z * XBESJ1(3.3_FP*R) & - 0.650909_FP * R*Z * ( 3.3_FP*XBESJ0(3.3_FP*R) - XBESJ1(3.3_FP*R)/R ) & - 3.69755_FP * XBESY1(3.3_FP*R) & - 3.69755_FP * R * ( 3.3_FP*XBESY0(3.3_FP*R) - XBESY1(3.3_FP*R)/R ) & - 2.55895_FP * XBESJ1(2.64_FP*R) * cos(1.98_FP*Z) & - 2.55895_FP * R * ( 2.64_FP*XBESJ0(2.64_FP*R) - XBESJ1(2.64_FP*R)/R) & * cos(1.98_FP*Z) & + 2.00008_FP * XBESJ1(1.98_FP*R) * cos(2.64_FP*Z) & + 2.00008_FP * R * ( 1.98_FP*XBESJ0(1.98_FP*R) - XBESJ1(1.98_FP*R)/R) & * cos(2.64_FP*Z) & + 1.2129_FP * XBESJ1(0.33_FP*R) * cos(3.2834585424518457_FP*Z) & + 1.2129_FP * R * ( 0.33_FP*XBESJ0(0.33_FP*R) - XBESJ1(0.33_FP*R)/R) & * cos(3.2834585424518457_FP*Z) & - 2.13022_FP * XBESJ1(3.8587562763149474_FP*R) * cosh(2.0_FP*Z) & - 2.13022_FP * R * ( 3.8587562763149474_FP*XBESJ0(3.8587562763149474_FP*R) & - XBESJ1(3.8587562763149474_FP*R)/R) * cosh(2.0_FP*Z) & + 1.32363_FP * XBESJ1(2.64_FP*R) * sin(1.98_FP*Z) & + 1.32363_FP * R * (2.64_FP*XBESJ0(2.64_FP*R) & - XBESJ1(2.64_FP*R)/R) * sin(1.98_FP*Z) & ) by=by/self%axis_Btor end function by real(FP) function btor(self, x, y, phi) real(FP), intent(in) :: x, y, phi class(carthy_t), intent(in) :: self btor = 1.0_FP/x end function btor #undef XBESJ1 #undef XBESY1 end module carthy_equilibrium_m