module cerfons_equilibrium_m !! Module containing the implementation of an analytic tokamak equilibrium as !! described by Cerfon and Freidberg (Phys. Plasmas 2010 https://doi.org/10.1063/1.3328818) !! N.b. this equilibrium is axisymmetric - phi, when required, is not used 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 analytic_divertor_equilibrium_m, only: & analytic_divertor_equilibrium_t, max_polygon_N_pts use descriptors_m use elementary_functions_m, only : smoothstep => step_hermite implicit none ! Poloidal field is multiplied with tanh step, setting it to smoothly zero ! in weird regions real(FP), parameter, private :: & y_bpolzero = -2.0_FP, & wy_bpolzero = 5.0E-3_FP type, public, extends(analytic_divertor_equilibrium_t) :: cerfons_t real(FP) :: a, rx2, zx2, normpol, normtor real(FP) :: c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12 contains procedure, public :: init procedure, public :: display procedure, public :: debug procedure, public :: psi procedure, public :: bx procedure, public :: by procedure, public :: btor end type cerfons_t contains subroutine init(self, filename) ! Initialisation function for cerfons equilibrium class(cerfons_t), intent(inout) :: self character(len=*), intent(in), optional :: filename integer :: io_error 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_cerfons_params / & 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 ! parameters as used in Held paper self%a = 0.0_FP self%c1 = 0.0735011444550040364542354505071_FP self%c2 = -0.0866241743631724894540195833599_FP self%c3 = -0.146393154340110145566533995171_FP self%c4 = -0.0763123710053627236722106917166_FP self%c5 = 0.0903179011379421292953276540381_FP self%c6 = -0.0915754123901870980044268057017_FP self%c7 = -0.00389228297983755811356246513352_FP self%c8 = 0.0427189122507639110238025727235_FP self%c9 = 0.227554564600278024358552580518_FP self%c10 = -0.130472413601776209642913874538_FP self%c11 = -0.0300697410847693691591393359449_FP self%c12 = 0.00421267189210391228482454852824_FP self%R0 = 1.073539190422562_FP self%Z0 = 0.029233201183845_FP self%RX = 0.733701975249365_FP self%ZX = -0.763711941107498 self%rx2 = 10.0_FP self%zx2 = 1.0_FP self%O_point_psi = -0.058149455020455_FP self%X_point_psi = 0.0_FP self%normpol = 1.0_FP self%axis_Btor = self%normpol self%normtor = 0.0_FP xmin = 0.4_FP xmax = 1.5_FP ymin = -1.0_FP ymax = 0.8_FP rhomin = 0.5_FP rhomax = 1.1_FP rhomin_privflux = 0.95_FP district_definition = 'flux' divertor_polygon_N_pts = 8 divertor_polygon_X_pts(1:divertor_polygon_N_pts) = & (/+0.470_FP, +0.000_FP, +0.00_FP, +10.0_FP, +10.00_FP, +1.085_FP, +0.724_FP, +0.470_FP/) divertor_polygon_Y_pts(1:divertor_polygon_N_pts) = & (/-0.812_FP, -0.812_FP, -10.0_FP, -10.0_FP, -0.744_FP, -0.744_FP, -0.998_FP, -0.812_FP/) exclusion_polygon_N_pts = 8 exclusion_polygon_X_pts(1:exclusion_polygon_N_pts) = & (/+0.470_FP, +0.000_FP, +0.00_FP, +10.0_FP, +10.00_FP, +1.085_FP, +0.724_FP, +0.470_FP/) exclusion_polygon_Y_pts(1:exclusion_polygon_N_pts) = & (/-1.0120_FP, -1.0120_FP, -10.2000_FP, -10.2000_FP, -0.9440_FP, -0.9440_FP, -1.1980_FP, -1.0120_FP/) self%invert_divertor_polygon = .true. self%invert_exclusion_polygon = .true. 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_cerfons_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%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) class(cerfons_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,"CERFONS 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(cerfons_t), intent(in) :: self write(get_stdout(), 100) & 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 A.J.Cerfon et al, Physics of Plasmas 17,& & p. 032502, (2010), section IX, up-down asymmetric formulation ' , /, & 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(cerfons_t), intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP) :: R, Z, logr R = x*self%r0 Z = self%r0*(y+self%z0/self%r0) !.... ! This is an ugly hack that prevents a float invalid in in the field line tracing ! The initial step size in dop853 is in some cases too large with conditional trace and ! the trace jumps too far. Eventually the first step is anyway corrected and the result is fine if (R<=0.0) then logr = -100.0_FP else logr = log(R) endif ! .... psi = R**4/8.0_FP+self%a*(R**2*logr/2.0_FP-R**4/8.0_FP) & + self%c1 * ( 1.0_FP) & + self%c2 * ( R**2) & + self%c3 * ( Z**2-R**2*logr ) & + self%c4 * ( R**4-4.0_FP*R**2*Z**2 ) & + self%c5 * ( 2._FP*Z**4-9._FP*Z**2*R**2+(3._FP*R**4-12._FP*R**2*Z**2)*logr ) & + self%c6 * ( R**6-12.0_FP*R**4*Z**2+8.0_FP*R**2*Z**4 ) & + self%c7 * ( 8.0_FP*Z**6-140.0_FP*Z**4*R**2+75.0_FP*Z**2*R**4 & + (-15.0_FP*R**6+180.0_FP*R**4*Z**2-120.0_FP*R**2*Z**4) * logr ) & + self%c8 * (Z) & + self%c9 * (R**2*Z ) & + self%c10* (Z**3-3.0_FP*R**2*Z*logr ) & + self%c11* (3.0_FP*R**4*Z-4.0_FP*R**2*Z**3) & + self%c12* (8.0_FP*Z**5-45.0_FP*R**4*Z & + (-80.0_FP*R**2*Z**3+60.0_FP*R**4*Z)*logr ) end function psi real(FP) function bx(self, x, y, phi) real(FP), intent(in) :: x, y, phi class(cerfons_t), intent(in) :: self real(FP) :: R, Z, logr R = x*self%r0 Z = self%r0*(y+self%z0/self%r0) !.... ugly hack see psi if (R<=0.0) then logr = -100.0_FP else logr = log(R) endif bx = 1/R*( & 0.0_FP & + self%c1 * ( 0.0_FP) & + self%c2 * ( 0.0_FP) & + self%c3 * ( -2.0_FP*Z ) & + self%c4 * ( 8.0_FP*R**2*Z ) & + self%c5 * ( -8.0_FP*Z**3 + 18.0_FP*R**2*Z+24.0_FP*R**2*Z*logr ) & + self%c6 * ( 24.0_FP*R**4*Z-32.0_FP*R**2*Z**3 ) & + self%c7 * ( -48.0_FP*Z**5+560.0_FP*R**2*Z**3-150.0_FP*R**4*Z & + (-360.0_FP*R**4*Z+480.0_FP*R**2*Z**3)*logr ) & + self%c8 * ( -1.0_FP ) & + self%c9 * ( -R**2) & + self%c10* ( -3.0_FP*Z**2+3.0_FP*R**2*logr) & + self%c11* ( -3.0_FP*R**4+12.0_FP*R**2*Z**2) & + self%c12* ( -40.0_FP*Z**4+45.0_FP*R**4 & + (240.0_FP*R**2*Z**2-60.0_FP*R**4)*logr ) ) bx = self%normpol*bx ! set bx to zero below y_bpolzero (region out of interest) bx = bx * smoothstep(y_bpolzero,y,wy_bpolzero) end function bx real(FP) function by(self, x, y, phi) real(FP), intent(in) :: x, y, phi class(cerfons_t), intent(in) :: self real(FP) :: R, Z, logr R = x*self%r0 Z = self%r0*(y+self%z0/self%r0) !.... ugly hack see psi if (R<=0.0) then logr = -100.0_FP else logr = log(R) endif by = 1/R*( & R**3/2.0_FP +self%a*(R*logr + R/2.0_FP - R**3/2.0_FP) & + self%c1 * ( 0.0_FP ) & + self%c2 * ( 2.0_FP*R ) & + self%c3 * ( -2.0_FP*R*logr-R ) & + self%c4 * ( 4.0_FP*R**3-8.0_FP*R*Z**2 ) & + self%c5 * ( -30.0_FP*R*Z**2+3.0_FP*R**3+(12.0_FP*R**3-24.0_FP*R*Z**2)*logr) & + self%c6 * ( 6.0_FP*R**5-48.0_FP*R**3*Z**2+16.0_FP*R*Z**4) & + self%c7 * ( -400.0_FP*R*Z**4+480.0_FP*R**3*Z**2-15.0_FP*R**5 & + (-90.0_FP*R**5+720.0_FP*R**3*Z**2-240.0_FP*R*Z**4)*logr ) & + self%c8 * ( 0.0_FP) & + self%c9 * ( 2.0_FP*R*Z) & + self%c10* ( -6.0_FP*R*Z*logr-3.0_FP*R*Z) & + self%c11* ( 12.0_FP*R**3*Z-8.0_FP*R*Z**3) & + self%c12* ( -120.0_FP*R**3*Z -80.0_FP*R*Z**3 & + (-160.0_FP*R*Z**3+240.0_FP*R**3*Z)*logr ) ) by = self%normpol*by ! set by to zero below y_bpolzero (region out of interest) by = by * smoothstep(y_bpolzero,y,wy_bpolzero) end function by real(FP) function btor(self, x, y, phi) real(FP), intent(in) :: x, y, phi class(cerfons_t), intent(in) :: self real(FP) :: R, Z, psi_pt R = x*self%r0 Z = self%r0*(y+self%z0/self%r0) psi_pt = self%psi(x, y, phi) btor = 1.0_FP/R*sqrt(1.0_FP-2.0_FP*self%normtor*psi_pt) end function btor end module cerfons_equilibrium_m