module salpha_equilibrium_m !! Module containing the implementation of an salpha like equilibrium !! type. !! N.b. this equilibrium is axisymmetric - phi, when required, is not used !! !!************************************************************************* !! !! Note: This equilibrium is currently implemented in an inconsitent way, !! i.e. the analytical q-profile does not match the calculated one. !! To be consistent the q-profile needs to stay as it is while the !! magnetic field has to be b_x = -rho*sin(theta)/(x*rho*sqrt(1-rho**2)) !! and b_y = rho*cos(theta)/(x*rho*sqrt(1-rho**2)) !! A consistent equilibrium of this form with a smaller set of parameters !! is implemented in circular_toroidal_equilibrium_m.f90 !! !!************************************************************************* !! use precision_m, only: FP, FP_EPS use screen_io_m, only: get_stdout use comm_handling_m, only: is_master use descriptors_m, only: DISTRICT_WALL, DISTRICT_SOL, & DISTRICT_CLOSED, DISTRICT_CORE use constants_m, only: PI, TWO_PI use equilibrium_m, only: equilibrium_t use params_equi_salpha_m, only: read_params_salpha, get_salpha_rhomin, & get_salpha_rhomax, get_salpha_minor_r, get_salpha_B_ref, & get_salpha_L_ref, get_salpha_q_ref, get_salpha_shear, & get_salpha_dtheta_limiter, get_salpha_rho_limiter, & get_salpha_theta_limiter use error_handling_m, only: handle_error, error_info_t use status_codes_m, only: PARALLAX_ERR_EQUILIBRIUM implicit none private type, public, extends(equilibrium_t) :: salpha_t !! Type implementing an salpha like equilibrium. The equilibrium has !! concentric flux surfaces in a toroidal geometry. Shear as well as !! the q-factor can be specified. To get the correct poloidal flux !! in Wb = V*s the reference magnetic field strength has to be provided !! in Tesla and the reference length scale has to be provided in !! meters. private real(FP) :: q_ref !! q value at the core real(FP) :: shear !! Magnetic shear parameter real(FP) :: minor_r !! Minor radius of the tokamak normalized to L_ref real(FP) :: L_ref !! Reference length scale real(FP) :: B_ref !! Reference magnetic field real(FP) :: theta_limiter !! Poloidal position of limiter real(FP) :: dtheta_limiter !! Poloidal width of limiter real(FP), private :: thlim1 real(FP), private :: thlim2 !! Poloidal positions of limiter real(FP), public :: rho_limiter !! Radial position (start) of limiter contains procedure, public :: init procedure, public :: display procedure, public :: debug procedure, public :: is_axisymmetric procedure, public :: rho procedure, public :: psi procedure, public :: bx procedure, public :: by procedure, public :: btor procedure, public :: jacobian procedure, public :: epol procedure, public :: erad procedure, public :: district procedure, public :: in_vessel procedure, public :: mag_axis_loc procedure, public :: qval procedure, public, nopass :: theta procedure, public :: x_point_psi procedure, public :: o_point_psi end type salpha_t contains subroutine init(self,filename) !! Reads input parameters and initializes the equilibrium class(salpha_t), intent(inout) :: self character(len = *), intent(in), optional :: filename !! Filename, where parameters are read from real(FP), parameter :: boxfac = 1.25_FP ! Factor by which rectangular box is increased w.r.t. rhomax if (present(filename)) then call read_params_salpha(filename) endif self%x0 = 1.0_FP self%y0 = 0.0_FP self%rhomin = get_salpha_rhomin() self%rhomax = get_salpha_rhomax() self%minor_r = get_salpha_minor_r() self%xmin = 1.0_FP - self%minor_r * self%rhomax * boxfac self%xmax = 1.0_FP + self%minor_r * self%rhomax * boxfac self%ymin = -self%minor_r * self%rhomax * boxfac self%ymax = self%minor_r * self%rhomax * boxfac self%q_ref = get_salpha_q_ref() self%shear = get_salpha_shear() self%L_ref = get_salpha_L_ref() self%B_ref = get_salpha_B_ref() self%theta_limiter = get_salpha_theta_limiter() self%dtheta_limiter = get_salpha_dtheta_limiter() self%rho_limiter = get_salpha_rho_limiter() ! Compute poloidal position of limiter plates self%thlim1 = self%theta_limiter + self%dtheta_limiter self%thlim1 = mod(self%thlim1, 2.0_FP * pi) self%thlim2 = self%theta_limiter - self%dtheta_limiter self%thlim2 = mod(self%thlim2, 2.0_FP * pi) if(self%shear < 0.0_FP) then call handle_error("Salpha equilibrium shear < 0!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__, & additional_info=& error_info_t("shear:", [self%shear])) endif if(self%minor_r < 0.0_FP) then call handle_error("Salpha equilibrium minor_r < 0!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__, & additional_info=& error_info_t("minor_r:", [self%minor_r])) endif if(self%q_ref < 0.0_FP) then call handle_error("Salpha equilibrium q_ref < 0!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__, & additional_info=& error_info_t("q_ref:", [self%q_ref])) endif ! If no limiter set, set rho_limiter to a large number if (self%dtheta_limiter <= 0.0_FP) then self%rho_limiter = 10.0_FP * self%rhomax endif ! Limiter width greater than pi will cause errors in the chi function, ! thus we raise an error in that case. if (self%dtheta_limiter >= pi) then call handle_error("Salpha equilibrium limiter width greater & &than pi!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__, & additional_info=& error_info_t("dtheta:", [self%dtheta_limiter])) endif self%initialized = .true. end subroutine init subroutine display(self) !! Prints to console information about the equilibrium class(salpha_t), intent(in) :: self write(get_stdout(), 99) & self%rhomin, self%rhomax, self%rho_limiter, & self%theta_limiter, self%dtheta_limiter 99 FORMAT( / & 1X,'S-ALPHA equilibrium:'/, & 1X,'rhomin = ',ES14.6E3 /, & 1X,'rhomax = ',ES14.6E3 /, & 1X,'rho_limiter = ',ES14.6E3 /, & 1X,'theta_limiter = ',ES14.6E3 /, & 1X,'dtheta_limiter = ',ES14.6E3 / & ) end subroutine display subroutine debug(self) !! Prints to console information about the equilibrium class(salpha_t), intent(in) :: self write(get_stdout(), 100) & self%x0, & self%y0, & self%rhomin, & self%rhomax, & self%xmin, & self%xmax, & self%ymin, & self%ymax, & self%q_ref, & self%shear, & self%minor_r, & self%L_ref, & self%B_ref, & self%theta_limiter, & self%dtheta_limiter, & self%thlim1, & self%thlim2, & self%rho_limiter 100 FORMAT( / & 1X,'S-Alpha equilibrium:'/, & 1X,'x0 = ',ES14.6E3 /, & 1X,'y0 = ',ES14.6E3 /, & 1X,'rhomin = ',ES14.6E3 /, & 1X,'rhomax = ',ES14.6E3 /, & 1X,'xmin = ',ES14.6E3 /, & 1X,'xmax = ',ES14.6E3 /, & 1X,'ymin = ',ES14.6E3 /, & 1X,'ymax = ',ES14.6E3 /, & 1X,'q_ref = ',ES14.6E3 /, & 1X,'shear = ',ES14.6E3 /, & 1X,'minor_r = ',ES14.6E3 /, & 1X,'L_ref = ',ES14.6E3 /, & 1X,'B_ref = ',ES14.6E3 /, & 1X,'theta_limiter = ',ES14.6E3 /, & 1X,'dtheta_limiter = ',ES14.6E3 /, & 1X,'thlim1 = ',ES14.6E3 /, & 1X,'thlim2 = ',ES14.6E3 /, & 1X,'rho_limiter = ',ES14.6E3 / & ) end subroutine debug logical function is_axisymmetric(self) class(salpha_t), intent(in) :: self is_axisymmetric = .true. end function real(FP) function rho(self, x, y, phi) !! Returns the flux surface label class(salpha_t),intent(in) :: self real(FP), intent(in) :: x, y, phi rho = sqrt((x - 1.0_FP)**2 + y**2) / self%minor_r ! NOTE: For simplicity we define rho as the radial coordinate in ! polar coordinates. As the flux surfaces are circular and ! concentric this is a valid label for the flux surfaces. Below ! is an alternative, more complicated flux surface label using ! the poloidal flux function. ! rho = sqrt((self%psi(x, y) - self%o_point_psi()) & ! / (self%x_point_psi() - self%o_point_psi())) end function rho real(FP) function psi(self, x, y) !! Returns the poloidal flux function in SI units class(salpha_t),intent(in) :: self real(FP), intent(in) :: x, y real(FP) :: r r = sqrt((x - 1.0_FP)**2 + y**2) associate(a => self%q_ref, b => self%shear, c => self%minor_r) ! This formula comes from analytically calculating psi if (abs(self%shear) > FP_EPS) then ! Only for non-zero shear psi = TWO_PI*self%L_ref**2 * self%B_ref * & (- c**2*Log((a*c**2)/(b + a*c**2)) & + c**2*Log((a*c**2 + b*r**2)/(a*c**2)))/(2.0_FP*b) else ! Limit of psi for shear going to zero psi = TWO_PI*self%L_ref**2 * self%B_ref * & (r**2 + 1.0_FP) / (2.0_FP * a) end if end associate end function real(FP) function x_point_psi(self) !! Return the value of psi at last closed flux surface class(salpha_t),intent(in) :: self x_point_psi = self%psi(1.0_FP + self%minor_r, 0.0_FP) end function real(FP) function o_point_psi(self) !! Return the value of psi at magnetic axis class(salpha_t),intent(in) :: self o_point_psi = self%psi(1.0_FP, 0.0_FP) end function real(FP) function bx(self, x, y, phi) !! Returns magnetic field component bx normalised to B_0 (on axis) class(salpha_t),intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP) :: r r = sqrt((x - 1.0_FP)**2 + y**2) bx = -r / x * 1.0_FP/ self%qval(r/self%minor_r) * sin(self%theta(x, y)) end function bx real(FP) function by(self, x, y, phi) !! Return the magnetic field component by normalised to B_0 (on axis) class(salpha_t),intent(in) :: self !! Instnace of the type real(FP), intent(in) :: x, y, phi real(FP) :: r r = sqrt((x - 1.0_FP)**2 + y**2) by = r / x * 1.0_FP/ self%qval(r/self%minor_r) * cos(self%theta(x, y)) end function by real(FP) function btor(self, x, y, phi) !! Returns the magnetic field component btor normalised to B_0 (on axis) class(salpha_t),intent(in) :: self real(FP), intent(in) :: x, y, phi btor = 1.0_FP / x end function btor real(FP) function jacobian(self, x, y, phi) !! Returns the jacobian of the geometry. The geometry is cylindrical !! and thus J = x class(salpha_t),intent(in) :: self real(FP), intent(in) :: x, y, phi jacobian = x end function jacobian subroutine epol(self, x, y, phi, epolx, epoly) !! Calculates the unit vector along poloidal (along flux surface) !! direction class(salpha_t),intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP), intent(out) :: epolx !! Unit vector component in x-direction real(FP), intent(out) :: epoly !! Unit vector component in y-direction epolx = -y/sqrt((x - 1.0_FP)**2+y**2) epoly = (x - 1.0_FP)/sqrt((x - 1.0_FP)**2+y**2) end subroutine epol subroutine erad(self, x, y, phi, eradx, erady) !! Calculates unit vector along radial (across flux surface) direction class(salpha_t),intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP), intent(out) :: eradx !! Unit vector component in x-direction real(FP), intent(out) :: erady !! Unit vector component in y-direction eradx = (x - 1.0_FP)/sqrt((x - 1.0_FP)**2+y**2) erady = y/sqrt((x - 1.0_FP)**2+y**2) end subroutine erad integer function district(self, x, y, phi) !! Returns in which district point (x, y) is (see module descriptors_m) class(salpha_t),intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP) :: rhot rhot = self%rho(x, y, phi) if (rhot > self%rhomax - FP_EPS) then district = DISTRICT_WALL return endif if (rhot < self%rhomin + FP_EPS) then if(self%rhomin > FP_EPS) then ! if rhomin == 0 we do not have a core boundary district = DISTRICT_CORE return endif endif if (rhot <= self%rho_limiter) then district = DISTRICT_CLOSED return else district = DISTRICT_SOL return endif end function district real(FP) function qval(self, rho_pt) !! Returns the safety factor class(salpha_t),intent(in) :: self real(FP), intent(in) :: rho_pt !! Flux surface label qval = self%q_ref + self%shear * rho_pt**2 end function qval real(FP) function theta(x, y) !! Returns the geometric poloidal angle real(FP), intent(in) :: x, y theta = atan2(y, x - 1.0_FP) if(theta < 0) theta = theta + TWO_PI end function theta ! TODO: Fix radial penalisation ! THIS IS THE WORKING VERSION OF in_vessel FOR CIRCULAR. ! THE (uncommented) VERSION IS A HACK TO GET RADIAL PENALISATION WORKING! ! logical function in_vessel(self, x, y) ! !! Returns a boolean of whether a given point is within the reactor ! !! vessel (inside the bounding surface of ! !! the first-wall and divertor) or outside this bounding surface. ! class(circular_t),intent(in) :: self ! !! circular equilibrium ! real(FP), intent(in) :: x ! !! x-coordinate ! real(FP), intent(in) :: y ! !! y-coordinate ! real(FP) :: rho, theta ! rho = self%rho(x, y) ! theta = self%theta(x, y) ! if (rho < self%rho_limiter) then ! ! Core/open-fieldline region ! in_vessel = .true. ! elseif(rho > self%rhomax) then ! ! First-wall ! in_vessel = .false. ! elseif((theta1_centred() < 0.0_FP) .and. (theta2_centred() > 0.0_FP)) then ! ! In limiter ! in_vessel = .false. ! else ! ! In SOL ! in_vessel = .true. ! endif ! contains ! real(FP) function theta1_centred() ! ! poloidal distance to limiter 1 ! theta1_centred = mod(theta - self%thlim1 + two_pi, two_pi) ! if (theta1_centred > pi - self%dtheta_limiter) then ! theta1_centred = theta1_centred - two_pi ! endif ! end ! real(FP) function theta2_centred() ! ! poloidal distance to limiter 2 ! theta2_centred = mod(theta - self%thlim2 + two_pi, two_pi) ! if (theta2_centred > pi + self%dtheta_limiter) then ! theta2_centred = theta2_centred - two_pi ! endif ! end ! end function in_vessel logical function in_vessel(self, x, y, phi) !! Returns a boolean of whether a given point is within theta1 and !! theta2 (angular extent of limiter) class(salpha_t),intent(in) :: self real(FP), intent(in) :: x !! x-coordinate real(FP), intent(in) :: y !! y-coordinate real(FP), intent(in) :: phi !! toroidal coordinate (unused) real(FP) :: rho_pt, theta_pt ! TODO: This routine has been copied from circular_t and has not been ! tested or adapted to the salpha geometry rho_pt = self%rho(x, y, phi) theta_pt = self%theta(x, y) ! Count the open-fieldline region near the limiter as 'out of vessel', ! then apply radial step. ! if (rho_t < self%rho_limiter) then ! ! Core/open-fieldline region ! in_vessel = .true. if(rho_pt > self%rhomax) then ! First-wall in_vessel = .false. elseif((theta1_centred() < 0.0_FP) .and. (theta2_centred() > 0.0_FP)) then ! In limiter in_vessel = .false. else ! In SOL in_vessel = .true. endif contains real(FP) function theta1_centred() ! poloidal distance to limiter 1 theta1_centred = mod(theta_pt - self%thlim1 + two_pi, two_pi) if (theta1_centred > pi - self%dtheta_limiter) then theta1_centred = theta1_centred - two_pi endif end real(FP) function theta2_centred() ! poloidal distance to limiter 2 theta2_centred = mod(theta_pt - self%thlim2 + two_pi, two_pi) if (theta2_centred > pi + self%dtheta_limiter) then theta2_centred = theta2_centred - two_pi endif end end function in_vessel subroutine mag_axis_loc(self, phi, axis_x, axis_y) !! Returns the coordinates of magnetic axis class(salpha_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: phi !! Toroidal angle real(FP), intent(out) :: axis_x !! x-coordinate of the magnetic axis real(FP), intent(out) :: axis_y !! y-coordinate of the magnetic axis axis_x = self%x0 axis_y = self%y0 end subroutine mag_axis_loc end module salpha_equilibrium_m