module circular_equilibrium_m !! Definition of circular equilibrium, based on prescribed q-profile !! N.b. this equilibrium is axisymmetric - phi, when required, is not used use precision_m, only: FP use screen_io_m, only: get_stdout use descriptors_m, only: DISTRICT_WALL, DISTRICT_SOL, & DISTRICT_CLOSED, DISTRICT_CORE use constants_m, only: pi, two_pi use elementary_functions_m, only: step_tanh use equilibrium_m, only: equilibrium_t use params_equi_circular_m, only: read_params_circular, & get_circular_rhomin, get_circular_rhomax, get_circular_qtype, & get_circular_rhoq_ref, get_circular_q_ref, get_circular_shear, & get_circular_dtheta_limiter, get_circular_rho_limiter, & get_circular_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) :: circular_t private integer :: qtype ! indicator for type of q-profile ! 0 => Linear q profile ! 1 => Inverse (1/x hyperbolic) q-profile real(FP) :: rhoq_ref ! reference flux surface for definition of q-profile real(FP) :: q_ref ! q value at reference real(FP) :: shear ! magnetic shear 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 :: 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 end type circular_t contains subroutine init(self, filename) !! Initialises the circular equilibrium class(circular_t), intent(inout) :: self !! Instance of the type character(len=*), intent(in), optional :: filename !! Filename, where parameters are read from ! NOTE: This number was chosen for backwards compatibility ! with the cmod-testsuite real(FP), parameter :: boxfac = 0.4_FP/0.3685_FP if (present(filename)) then call read_params_circular(filename) endif self%x0 = 0.0_FP self%y0 = 0.0_FP self%rhomin = get_circular_rhomin() self%rhomax = get_circular_rhomax() self%xmin = -boxfac * self%rhomax self%xmax = boxfac * self%rhomax self%ymin = -boxfac * self%rhomax self%ymax = boxfac * self%rhomax self%qtype = get_circular_qtype() self%rhoq_ref = get_circular_rhoq_ref() self%q_ref = get_circular_q_ref() self%shear = get_circular_shear() self%theta_limiter = get_circular_theta_limiter() self%dtheta_limiter = get_circular_dtheta_limiter() self%rho_limiter = get_circular_rho_limiter() ! Compute poloidal position of limiter plates self%thlim1 = self%theta_limiter + self%dtheta_limiter self%thlim1 = mod(self%thlim1, two_pi) self%thlim2 = self%theta_limiter - self%dtheta_limiter self%thlim2 = mod(self%thlim2, two_pi) ! 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("Circular 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(circular_t), intent(in) :: self write(get_stdout(), 99) & self%rhomin, self%rhomax, self%rho_limiter, & self%theta_limiter, self%dtheta_limiter 99 FORMAT( / & 1X,'CIRCULAR 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(circular_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%qtype, self%rhoq_ref, self%q_ref, self%shear, & self%theta_limiter, self%dtheta_limiter, self%thlim1, self%thlim2, & self%rho_limiter 100 FORMAT( / & 1X,'Circular 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,'qtype = ',I4 /, & 1X,'rhoq_ref = ',ES14.6E3 /, & 1X,'q_ref = ',ES14.6E3 /, & 1X,'shear = ',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(circular_t), intent(in) :: self is_axisymmetric = .true. end function real(FP) function rho(self, x, y, phi) !! flux surface label class(circular_t),intent(in) :: self real(FP), intent(in) :: x, y, phi rho = sqrt(x*x+y*y) end function rho real(FP) function bx(self, x, y, phi) !! magnetic field component bx normalised to B_0 (on axis) class(circular_t),intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP) :: rhot rhot = self%rho(x, y, phi) bx = -1.0_FP / self%qval(rhot) * rhot*sin(self%theta(x, y)) end function bx real(FP) function by(self, x, y, phi) !! magnetic field component by normalised to B_0 (on axis) class(circular_t),intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP) :: rhot rhot = self%rho(x, y, phi) by = 1.0_FP / self%qval(rhot) * rhot*cos(self%theta(x, y)) end function by real(FP) function btor(self, x, y, phi) !! magnetic field component btor normalised to B_0 (on axis) class(circular_t),intent(in) :: self real(FP), intent(in) :: x, y, phi btor = 1.0_FP end function btor real(FP) function jacobian(self, x, y, phi) !! Jacobian of geometry (J=1, i.e. infinite aspect ratio assumption) class(circular_t),intent(in) :: self real(FP), intent(in) :: x, y, phi jacobian = 1.0_FP end function jacobian subroutine epol(self, x, y, phi, epolx, epoly) !! unit vector along poloidal (along flux surface) direction class(circular_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**2+y**2) epoly = x/sqrt(x**2+y**2) end subroutine epol subroutine erad(self, x, y, phi, eradx, erady) !! unit vector along radial (across flux surface) direction class(circular_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/sqrt(x**2+y**2) erady = y/sqrt(x**2+y**2) end subroutine erad integer function district(self, x, y, phi) !! returns in which district point (x, y, phi) is (see module descriptors_m) class(circular_t),intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP) :: rhot rhot = self%rho(x, y, phi) if (rhot > self%rhomax) then district = DISTRICT_WALL return endif if (rhot < self%rhomin) then district = DISTRICT_CORE return 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) !! safety factor class(circular_t),intent(in) :: self real(FP), intent(in) :: rho_pt !! flux label if (self%qtype == 0) then qval = self%q_ref + self%shear * (rho_pt - self%rhoq_ref) elseif (self%qtype == 1) then qval = self%q_ref / (1.0_FP - self%shear * (rho_pt - self%rhoq_ref)) else call handle_error("Safety factor type for circular equilibrium & ¬ valid!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__, & additional_info=& error_info_t("qtype:", [self%qtype])) endif end function qval real(FP) function theta(x, y) !! geometric poloidal angle real(FP), intent(in) :: x, y if (x > 0.0_FP) then theta = atan(y/x) elseif (x < 0.0_FP) then theta = atan(y/x)+pi else if (y > 0.0_FP) then theta = pi/2.0_FP elseif (y < 0.0_FP) then theta = pi*3.0_FP/2.0_FP else theta = 0.0_FP endif endif theta = mod(theta, two_pi) if (theta < 0.0_FP) then theta = theta+two_pi endif end function theta logical function in_vessel(self, x, y, phi) !! Returns a boolean of whether a given point is within limiter region class(circular_t),intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP) :: rho_pt, theta_pt rho_pt = self%rho(x, y, phi) theta_pt = self%theta(x, y) in_vessel = .true. if((theta1_centred() < 0.0_FP) .and. (theta2_centred() > 0.0_FP)) then ! In limiter if (rho_pt >=self%rho_limiter) then in_vessel = .false. endif 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(circular_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 circular_equilibrium_m