module circular_toroidal_equilibrium_m !! Implementation of a analytical toroidal equilibrium with nested circular !! fluxsurfaces. For an extensive reference we refer to: !! https://youjunhu.github.io/research_notes/tokamak_equilibrium_htlatex/ !! tokamak_equilibrium.html#x1-9400016 use precision_m, only: FP, FP_EPS use constants_m, only: TWO_PI use equilibrium_m, only: equilibrium_t use kisslinger_m, only : kisslinger_t use params_equi_circular_toroidal_m, only: read_params_circular_toroidal,& get_circular_toroidal_rhomin, get_circular_toroidal_rhomax, & get_circular_toroidal_q_0, get_circular_toroidal_q_quad_param, & get_circular_toroidal_hel_amp, get_circular_toroidal_hel_m, & get_circular_toroidal_hel_n, get_circular_toroidal_hel_rho, & get_circular_toroidal_hel_sigma, read_circtor_bnd_params, & write_circtor_bnd_params, get_circular_toroidal_kiss_boundary_on,& nbnd_seg, bnd_paths, bnd_symms use descriptors_m, only: DISTRICT_CORE, DISTRICT_CLOSED, DISTRICT_WALL use error_handling_m, only: handle_error, error_info_t use status_codes_m, only: PARALLAX_ERR_PARAMETERS,PARALLAX_ERR_EQUILIBRIUM, & PARALLAX_WRN_GENERAL use screen_io_m, only: get_stdout use euclidean_geo_m, only : intersection_lines implicit none private type, public, extends(equilibrium_t) :: circular_toroidal_t !! Type implementing a consistent analytical toroidal !! equilibrium with nested circular fluxsurfaces private real(FP) :: q_0 ! Value of q(0) on the magnetic axis real(FP) :: q_quad_param ! Quadratic parameter of the q-profile parabola real(FP) :: hel_amp ! Amplitude of the helical perturbation integer :: hel_m ! Poloidal mode number of the helical perturbation integer :: hel_n ! Toroidal mode number of the helical perturbation real(FP) :: hel_rho ! Radial position of the Gaussian radial envelope function real(FP) :: hel_sigma ! Variance of the Gaussian radial envelope function logical :: kiss_boundary_on ! True if some boundary description in Kisslinger format is applied type(kisslinger_t), private, allocatable, dimension(:) :: kiss_bnd ! Description of boundary in Kisslinger's format 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 :: theta procedure, private :: eta procedure, private :: radenv procedure, private :: dradenv_drho procedure, private :: deta_drho procedure, private :: deta_dtheta end type circular_toroidal_t contains subroutine init(self, filename) !! Initialises the toroidal equilibrium class(circular_toroidal_t), intent(inout) :: self !! Instance of class character(len=*), intent(in), optional :: filename integer :: ibnd if (present(filename)) then call read_params_circular_toroidal(filename) endif self%x0 = 1.0_FP self%y0 = 0.0_FP self%rhomin = get_circular_toroidal_rhomin() self%rhomax = get_circular_toroidal_rhomax() self%xmin = FP_EPS self%xmax = 2.0_FP - FP_EPS self%ymin = -1.0_FP + FP_EPS self%ymax = 1.0_FP - FP_EPS self%q_0 = get_circular_toroidal_q_0() self%q_quad_param = get_circular_toroidal_q_quad_param() self%hel_amp = get_circular_toroidal_hel_amp() self%hel_m = get_circular_toroidal_hel_m() self%hel_n = get_circular_toroidal_hel_n() self%hel_rho = get_circular_toroidal_hel_rho() self%hel_sigma = get_circular_toroidal_hel_sigma() self%kiss_boundary_on = get_circular_toroidal_kiss_boundary_on() if (self%rhomin > self%rhomax) then call handle_error("Circular toroidal equilibrium: rhomin can't be & larger than rhomax!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__, & additional_info = & error_info_t("rhomin, rhomax:", & r_info=[self%rhomin, self%rhomax])) elseif (self%rhomax >= 1.0_FP) then call handle_error("Circular toroidal equilibrium: rhomax can't be & larger than 1!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__, & additional_info = & error_info_t("rhomax:", r_info=[self%rhomax])) endif ! Initialise description of boundaries if (present(filename)) then call read_circtor_bnd_params(filename) endif if(self%kiss_boundary_on) then allocate(self%kiss_bnd(nbnd_seg)) do ibnd = 1, nbnd_seg call self%kiss_bnd(ibnd)%init(bnd_paths(ibnd), & self%x0, self%y0, bnd_symms(ibnd)) enddo endif self%initialized = .true. end subroutine subroutine display(self) !! Prints to console information about the equilibrium and the boundary class(circular_toroidal_t), intent(in) :: self !! Instance of class integer :: ibnd write(get_stdout(), 99) & self%rhomin, self%rhomax, self%q_0, self%q_quad_param, & self%hel_amp, self%hel_m, self%hel_n, self%hel_rho, self%hel_sigma 99 FORMAT(1X,'CIRCULAR TOROIDAL equilibrium:'/, & 1X,'rhomin = ',ES14.6E3 /, & 1X,'rhomax = ',ES14.6E3 /, & 1X,'q_0 = ',ES14.6E3 /, & 1X,'q_quad_param = ',ES14.6E3 /, & 1X,'hel_amp = ',ES14.6E3 /, & 1X,'hel_m = ',I8 /, & 1X,'hel_n = ',I8 /, & 1X,'hel_rho = ',ES14.6E3 /, & 1X,'hel_sigma = ',ES14.6E3 / ) if(self%kiss_boundary_on) then do ibnd = 1, nbnd_seg call self%kiss_bnd(ibnd)%display() enddo endif end subroutine subroutine debug(self) !! Prints to console extended debug information about the equilibrium class(circular_toroidal_t), intent(in) :: self !! Instance of class write(get_stdout(), 100) & self%x0, self%y0, self%rhomin, self%rhomax, self%xmin, self%xmax, & self%ymin, self%ymax, self%q_0, self%q_quad_param, & self%hel_amp, self%hel_m, self%hel_n, self%hel_rho, self%hel_sigma 100 FORMAT(1X,'CIRCULAR TOROIDAL 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_0 = ',ES14.6E3 /, & 1X,'q_quad_param = ',ES14.6E3 /, & 1X,'hel_amp = ',ES14.6E3 /, & 1X,'hel_m = ',I8 /, & 1X,'hel_n = ',I8 /, & 1X,'hel_rho = ',ES14.6E3 /, & 1X,'hel_sigma = ',ES14.6E3 / ) end subroutine logical function is_axisymmetric(self) !! Returns if the equilibrium is axisymmetric or not class(circular_toroidal_t), intent(in) :: self !! Instance of class if (self%hel_amp > FP_EPS) then is_axisymmetric = .false. else is_axisymmetric = .true. endif end function real(FP) function rho(self, x, y, phi) !! Returns the radial coordinate, i.e. distance to torus axis class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) rho = sqrt((x - 1.0_FP)**2 + y**2) end function real(FP) function theta(self, x, y, phi) !! Returns the polar angle in polar coordinates class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) theta = atan2(y, (x - 1.0_FP)) if (theta < 0.0_FP) theta = theta + TWO_PI end function real(FP) function qval(self, rho_pt) !! Returns the safety factor class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: rho_pt !! rho-coordinate (along minor radius) qval = self%q_0 + self%q_quad_param * rho_pt**2 end function real(FP) function bx(self, x, y, phi) !! Returns the x-component of the magnetic field class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) real(FP) :: rho_pt, theta_pt, qval_pt, eta, radenv, dradenv_drho real(FP) :: deta_dy, deta_drho, deta_dtheta, drho_dy, dtheta_dy, amp integer :: mode_m, mode_n rho_pt = self%rho(x, y, phi) theta_pt = self%theta(x, y, phi) qval_pt = self%qval(rho_pt) bx = - rho_pt * sin(theta_pt) / (x * qval_pt * sqrt(1.0_FP - rho_pt**2)) ! Before applying the perturbation, check if we are on the axis, ! since the perturbed bx ~ 1/rho_pt if (rho_pt < FP_EPS) return ! Set the values for m and n and the amplitude of the ! helical perturbation mode_m = self%hel_m mode_n = self%hel_n amp = self%hel_amp ! Calculate the straight field line angle eta eta = self%eta(rho_pt, theta_pt) ! Calculate the envelope function and its derivative radenv = self%radenv(rho_pt) dradenv_drho = self%dradenv_drho(rho_pt) ! Calculate the derivative of eta with respect to y and all ! derivatives that are necessary for that deta_drho = self%deta_drho(rho_pt, theta_pt) deta_dtheta = self%deta_dtheta(rho_pt, theta_pt) drho_dy = y / rho_pt dtheta_dy = (x - 1.0_FP) / rho_pt**2 deta_dy = deta_drho * drho_dy + deta_dtheta * dtheta_dy ! Add the helical perturbation to the x-component of the magnetic field bx = bx + amp / x * & (-dradenv_drho * drho_dy * sin(mode_m * eta - mode_n * phi) & - mode_m * radenv * deta_dy * cos(mode_m * eta - mode_n * phi)) end function real(FP) function by(self, x, y, phi) !! Returns the y-component of the magnetic field class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) real(FP) :: rho_pt, theta_pt, qval_pt , eta, radenv, dradenv_drho real(FP) :: deta_dx, deta_drho, deta_dtheta, drho_dx, dtheta_dx, amp integer :: mode_m, mode_n rho_pt = self%rho(x, y, phi) theta_pt = self%theta(x, y, phi) qval_pt = self%qval(rho_pt) by = rho_pt * cos(theta_pt) / (x * qval_pt * sqrt(1.0_FP - rho_pt**2)) ! Before applying the perturbation, check if we are on the axis, ! since the perturbed by ~ 1/rho_pt if (rho_pt < FP_EPS) return ! Set the values for m and n and the amplitude of the ! helical perturbation mode_m = self%hel_m mode_n = self%hel_n amp = self%hel_amp ! Calculate the straight field line angle eta eta = self%eta(rho_pt, theta_pt) ! Calculate the envelope function and its derivative radenv = self%radenv(rho_pt) dradenv_drho = self%dradenv_drho(rho_pt) ! Calculate the derivative of eta with respect to y and all ! derivatives that are necessary for that deta_drho = self%deta_drho(rho_pt, theta_pt) deta_dtheta = self%deta_dtheta(rho_pt, theta_pt) drho_dx = (x - 1.0_FP) / rho_pt dtheta_dx = -y / rho_pt**2 deta_dx = deta_drho * drho_dx + deta_dtheta * dtheta_dx ! Add the helical perturbation to the y-component of the magnetic field by = by + amp / x * & (dradenv_drho * drho_dx * sin(mode_m * eta - mode_n * phi) & + mode_m * radenv * deta_dx * cos(mode_m * eta - mode_n * phi)) end function real(FP) function btor(self, x, y, phi) !! Returns the toroidal component btor of the magnetic field class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) btor = 1.0_FP / x end function real(FP) function jacobian(self, x, y, phi) !! Returns the Jacobian of the geometry class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) jacobian = x end function subroutine epol(self, x, y, phi, epolx, epoly) !! Calculates the the x- and y-component of the !! poloidal unit vector epol (along flux surfaces) class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) real(FP), intent(out) :: epolx, epoly !! Components of poloidal unit vector erad real(FP) :: rho_pt rho_pt = self%rho(x, y, phi) ! Since we have to divide by rho_loc, we have to make sure that we ! do not divide by zero on the magnetic axis. For this case we set ! epolx and eploy a default value, which ensures |epol| = 1 if (rho_pt < FP_EPS) then epolx = -1.0_FP / sqrt(2.0_FP) epoly = 1.0_FP / sqrt(2.0_FP) else epolx = -y / rho_pt epoly = (x - 1.0_FP) / rho_pt endif end subroutine subroutine erad(self, x, y, phi, eradx, erady) !! Calculates the the x- and y-component of the !! radial unit vector erad (across flux surfaces) class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) real(FP), intent(out) :: eradx, erady !! Components of radial unit vector erad real(FP) :: rho_pt rho_pt = self%rho(x, y, phi) ! Since we have to divide by rho_loc, we have to make sure that we ! do not divide by zero on the magnetic axis. For this case we set ! eradx and erady a default value, which ensures |erad| = 1 and that ! erad is orthogonal to epol if (rho_pt < FP_EPS) then eradx = 1.0_FP / sqrt(2.0_FP) erady = 1.0_FP / sqrt(2.0_FP) else eradx = (x - 1.0_FP) / rho_pt erady = y / rho_pt endif end subroutine integer function district(self, x, y, phi) !! Returns in which district point (x, y) is (see module descriptors_m) class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) if (self%rho(x, y, phi) < self%rhomin) then district = DISTRICT_CORE elseif (self%rho(x, y, phi) > self%rhomax) then district = DISTRICT_WALL else district = DISTRICT_CLOSED endif end function logical function in_vessel(self, x, y, phi) !! Returns whether a given point is within the vessel or not class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) integer :: npoly, ibnd, info, ipoly real(FP), allocatable, dimension(:) :: xpoly, ypoly real(FP), dimension(2) :: pa, qa, pb, qb real(FP):: ta, tb, xi, yi in_vessel = .true. do ibnd = 1, nbnd_seg ! Running through the boundary elements call self%kiss_bnd(ibnd)%polygon_at_tor_xsection( & phi, npoly, xpoly, ypoly) pa(1) = self%x0 pa(2) = self%y0 qa(1) = x qa(2) = y do ipoly = 1, npoly-1 pb(1) = xpoly(ipoly) pb(2) = ypoly(ipoly) qb(1) = xpoly(ipoly+1) qb(2) = ypoly(ipoly+1) call intersection_lines(pa, qa, pb, qb, xi, yi, info, & ta, tb) if (info == 0) then ! Checking if the intersection is on the 2nd line if (ta > 0.0_FP .and. ta < 1.0_FP) then in_vessel = .false. endif ! Checking if the intersection is on the 1st line if (tb > 1.0_FP .or. tb < 0.0_FP) then in_vessel = .true. endif if(.not. in_vessel) then exit endif else if (info == 1 .or. info == -2) then call handle_error(& "Line intersection detected degeneracy", & PARALLAX_WRN_GENERAL, __LINE__, __FILE__) endif enddo if (allocated(xpoly)) deallocate(xpoly) if (allocated(ypoly)) deallocate(ypoly) if(.not. in_vessel) then exit endif enddo end function subroutine mag_axis_loc(self, phi, axis_x, axis_y) !! Calculated the coordinates of the magnetic axis class(circular_toroidal_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 real(FP) function eta(self, rho_pt, theta_pt) !! Returns the straight field line angle eta in the range [0,2*pi] class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: rho_pt !! rho-coordinate (along minor radius) real(FP), intent(in) :: theta_pt !! Geometric poloidal angle eta = -2.0_FP * atan((rho_pt - 1.0_FP) * tan(theta_pt / 2.0_FP) / & sqrt(1.0_FP - rho_pt**2)) if (eta < 0.0_FP) eta = eta + TWO_PI end function real(FP) function radenv(self, rho_pt) !! Returns the radial envelope, which is a Gaussian at position !! hel_rho and with the variance hel_sigma class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: rho_pt !! rho-coordinate (along minor radius) radenv = exp(-(rho_pt - self%hel_rho)**2 / (2.0_FP * self%hel_sigma**2)) end function real(FP) function dradenv_drho(self, rho_pt) !! Returns the derivative of the radial envelope with respect to rho class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: rho_pt !! rho-coordinate (along minor radius) dradenv_drho = -(rho_pt - self%hel_rho) / self%hel_sigma**2 * & exp(-(rho_pt - self%hel_rho)**2 / (2.0_FP * self%hel_sigma**2)) end function real(FP) function deta_drho(self, rho_pt, theta_pt) !! Returns the derivative of eta with respect to rho class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: rho_pt !! rho-coordinate (along minor radius) real(FP), intent(in) :: theta_pt !! Geometric poloidal angle deta_drho = -sin(theta_pt) / (sqrt(1.0_FP - rho_pt**2) * & (1.0_FP + rho_pt * cos(theta_pt))) end function real(FP) function deta_dtheta(self, rho_pt, theta_pt) !! Returns the derivative of eta with respect to theta class(circular_toroidal_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: rho_pt !! rho-coordinate (along minor radius) real(FP), intent(in) :: theta_pt !! Geometric poloidal angle deta_dtheta = sqrt(1.0_FP - rho_pt**2) / & (1.0_FP + rho_pt * cos(theta_pt)) end function end module circular_toroidal_equilibrium_m