module dommaschk_equilibrium_m !! [*] W. Dommaschk, "Representations for vacuum potentials in !! stellarators", Computer Physics Communications 40, pg. 203 (1986) !! !! Definition of Dommaschk potentials [*]. Fully 3D equilibrium use precision_m, only: FP, FP_EPS use screen_io_m, only: get_stdout use error_handling_m, only: handle_error, error_info_t use status_codes_m, only: PARALLAX_WRN_GENERAL, PARALLAX_ERR_EQUILIBRIUM, & PARALLAX_SUCCESS use constants_m, only: TWO_PI, PI use elementary_functions_m, only: factorial use fieldline_tracer_m, only: trace use descriptors_m, only: DISTRICT_CLOSED, DISTRICT_WALL, DISTRICT_OUT, & DISTRICT_CORE use polygon_m, only: closed_polygon2d_t, create_closed_polygon2d_t use params_equi_dommaschk_m, only: & read_params_dommaschk, write_params_dommaschk, & get_dommaschk_x0, get_dommaschk_y0, get_dommaschk_phi0, & get_dommaschk_xmin, get_dommaschk_xmax, & get_dommaschk_ymin, get_dommaschk_ymax, & get_dommaschk_rhomin, get_dommaschk_rhomax, & get_dommaschk_l_pol, get_dommaschk_m_tor_consecutive, & get_dommaschk_num_field_periods, & get_dommaschk_fitting_coef, get_dommaschk_bndry_from_file, & get_dommaschk_inner_bndry_file, get_dommaschk_inner_bndry_var, & get_dommaschk_outer_bndry_file, get_dommaschk_outer_bndry_var, & get_dommaschk_exclusion_from_file, & get_dommaschk_exclusion_file, get_dommaschk_exclusion_var, & get_dommaschk_rho_from_file, get_dommaschk_rho_file use equilibrium_m, only: equilibrium_t use interpolation_m, only: interpol1d implicit none private :: point_inside_polygons type, public, extends(equilibrium_t) :: dommaschk_t !! Type implementing 3D vacuum fields as described by Dommaschk [*]. !! These are determined with toroidal mode numbers m, poloidal mode !! numbers l, the number of field periods, and fitting coefficients !! which are specific to each magnetic configuration. integer, private :: l_pol !! Poloidal periodicity upper bound, see Eq. 1 in [*] integer, private :: m_tor_consecutive !! Toroidal periodicity upper bound, see Eq. 1 in [*] ! NOTE: This will be multiplied by number of field periods ! (m = m_tor_consecutive * num_field_periods) to satisfy the ! condition that m can only be integer multiples of the field ! periodicity (see Eq. 1 in [*]) integer, private :: num_field_periods !! Number of field periods real(FP), private, dimension(:,:,:), pointer :: fitting_coef !! This represents the coefficients a_{m,l}, b_{m,l}, c_{m,l}, and !! d_{m.l} from Eq. 12 of [*] !! These coefficients are specific to each magnetic configuration and !! modulate the amplitudes in phi of the D_{m,l}(R,Z) and !! N_{m,l-1}(R,Z) components !! The first dimension 1:4 gives a, b, c, or d !! The second dimension is 0:m_tor_consecutive, and the third 0:l_pol ! NOTE: These attributes are public in order to unit test them real(FP), public, dimension(:,:,:), allocatable :: CD_r_coef !! Coefficients of the R terms of the C^D_{m,k} function. See Eqs. 31 !! and 32 in [*] real(FP), public, dimension(:,:,:), allocatable :: CN_r_coef !! Coefficients of the R terms of the C^N_{m,k} function. See Eqs. 31 !! and 32 in [*] real(FP), public, dimension(:,:,:), allocatable :: C_r_power !! Powers of the R terms of the C^D_{m,k} and C^N_{m,k} functions. See !! Eqs. 31 and 32 in [*] real(FP), public, dimension(:,:,:), allocatable :: CD_lnr_coef !! Coefficients of the log(R) terms of the C^D_{m,k} function. See Eqs. !! 31 and 32 in [*] real(FP), public, dimension(:,:,:), allocatable :: CN_lnr_coef !! Coefficients of the log(R) terms of the C^N_{m,k} function. See Eqs. !! 31 and 32 in [*] real(FP), public, dimension(:,:,:), allocatable :: C_lnr_power !! Powers of the log(R) terms of the C^D_{m,k} and C^N_{m,k} functions. !! See Eqs. 31 and 32 in [*] real(FP), public, dimension(:,:,:), allocatable :: dCD_dR_r_coef !! Coefficients of the R terms of the derivative with respect to R of !! the C^D_{m,k} function real(FP), public, dimension(:,:,:), allocatable :: dCN_dR_r_coef !! Coefficients of the R terms of the derivative with respect to R of !! the C^N_{m,k} function real(FP), public, dimension(:,:,:), allocatable :: dC_dR_r_power !! Powers of the R terms of the derivatives with respect to R of the !! C^D_{m,k} and C^N_{m,k} functions. See Eqs. 31 and 32 in [*] real(FP), public, dimension(:,:,:), allocatable :: dCD_dR_lnr_coef !! Coefficients of the log(R) terms of the derivative with respect to R !! of the C^D_{m,k} function real(FP), public, dimension(:,:,:), allocatable :: dCN_dR_lnr_coef !! Coefficients of the log(R) terms of the derivative with respect to R !! of the C^N_{m,k} function real(FP), public, dimension(:,:,:), allocatable :: dC_dR_lnr_power !! Powers of the log(R) terms of the derivatives with respect to R of !! the C^D_{m,k} and C^N_{m,k} functions real(FP), public, dimension(:,:,:), allocatable :: I_mn_coef !! Coefficients of the Z terms of the I_{m,n} function. See Eqs. 3, 8, !! and 9 in [*] real(FP), public, dimension(:,:,:), allocatable :: I_mn_power !! Powers of the Z terms of the I_{m,n} function. See Eqs. 3, 8, and 9 !! in [*] real(FP), public, dimension(:,:,:), allocatable :: dI_dZ_coef !! Coefficients of the Z terms of the derivative with respect to Z of !! the I_{m,n} function real(FP), public, dimension(:,:,:), allocatable :: dI_dZ_power !! Powers of the Z terms of the derivative with respect to Z of the !! I_{m,n} function. real(FP), private :: B_norm_axis !! Normalization factor for magnetic field, calculated such that the !! toroidal field on axis is 1.0 logical, private :: bndry_from_file !! Switch whether to use box limits or polygons from file for !! defining the boundaries real(FP), dimension(:,:,:), allocatable, private :: outer_bndry_polygons !! Polygons which describe the outer boundary real(FP), dimension(:), allocatable, private :: outer_bndry_polygons_phi !! Toroidal angles of each outer boundary polygon logical, private :: use_inner_bndry !! Whether to use an inner boundary real(FP), dimension(:,:,:), allocatable, private :: inner_bndry_polygons !! Polygons which describe the inner boundary real(FP), dimension(:), allocatable, private :: inner_bndry_polygons_phi !! Toroidal angles of each inner boundary polygon logical, private :: exclusion_from_file !! Switch whether to use box limits or polygons from file for !! defining exclusion region (where DISTRICT_OUT) real(FP), dimension(:,:,:), allocatable, private :: exclusion_polygons !! Polygons which describe the domain exclusion boundaries real(FP), dimension(:), allocatable, private :: exclusion_polygons_phi !! Toroidal angles of exclusion polygons logical, private :: rho_from_file !! Whether to calculate an approximate rho from flux surface polygons !! read from file real(FP), dimension(:), allocatable, private :: rho_polygons_phi !! Toroidal angles of each flux surface polygon set real(FP), dimension(:,:,:,:), allocatable, private :: rho_polygons !! Vertices of polygons which describe flux surfaces real(FP), dimension(:), allocatable, private :: rho_array !! Values of rho on each flux surface from file real(FP), dimension(:), allocatable, private :: axis_array_x !! Magnetic axis x-coordinate in the first field period (equidistant in phi) real(FP), dimension(:), allocatable, private :: axis_array_y !! Magnetic axis y-coordinate in the first field period (equidistant in phi) contains ! Getter routines procedure, public :: get_l_pol procedure, public :: get_m_tor_consecutive procedure, public :: get_num_field_periods 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, private :: check_fitting_coef procedure, private :: init_CD_CN procedure, private :: init_Imn procedure, private :: init_B_norm procedure, private :: read_bndry_polygons procedure, private :: read_rho_polygons procedure, private :: init_rho_array procedure, private :: init_mag_axis_loc end type dommaschk_t interface module subroutine read_bndry_polygons(self, filename, variable, & phi_array, polygon_vertices) !! Reads in boundary / exclusion polygon data from NetCDF class(dommaschk_t), intent(in) :: self !! Instance of the type character(len=*), intent(in) :: filename !! NetCDF filename character(len=*), intent(in) :: variable !! NetCDF variable name real(FP), dimension(:), allocatable :: phi_array !! Toroidal angles of each polygon real(FP), dimension(:,:,:), allocatable, intent(out) :: polygon_vertices !! Closed polygon vertices at each toroidal angle end subroutine module subroutine read_rho_polygons(self, filename) !! Reads in flux surface polygon data from NetCDF class(dommaschk_t), intent(inout) :: self !! Instance of the type character(len=*), intent(in) :: filename !! NetCDF filename end subroutine end interface contains pure integer function get_l_pol(self) class(dommaschk_t), intent(in) :: self get_l_pol = self%l_pol end function get_l_pol pure integer function get_m_tor_consecutive(self) class(dommaschk_t), intent(in) :: self get_m_tor_consecutive = self%m_tor_consecutive end function get_m_tor_consecutive pure integer function get_num_field_periods(self) class(dommaschk_t), intent(in) :: self get_num_field_periods = self%num_field_periods end function get_num_field_periods subroutine init(self, filename) class(dommaschk_t), intent(inout) :: self character(len=*), intent(in), optional :: filename if (present(filename)) then call read_params_dommaschk(filename) end if self%x0 = get_dommaschk_x0() self%y0 = get_dommaschk_y0() self%phi0 = get_dommaschk_phi0() self%xmin = get_dommaschk_xmin() self%xmax = get_dommaschk_xmax() self%ymin = get_dommaschk_ymin() self%ymax = get_dommaschk_ymax() self%rhomin = get_dommaschk_rhomin() self%rhomax = get_dommaschk_rhomax() self%m_tor_consecutive = get_dommaschk_m_tor_consecutive() self%l_pol = get_dommaschk_l_pol() self%num_field_periods = get_dommaschk_num_field_periods() self%bndry_from_file = get_dommaschk_bndry_from_file() self%exclusion_from_file = get_dommaschk_exclusion_from_file() self%rho_from_file = get_dommaschk_rho_from_file() self%fitting_coef => get_dommaschk_fitting_coef() call self%check_fitting_coef() call self%init_CD_CN() call self%init_Imn() call self%init_B_norm() if (self%bndry_from_file) then ! Outer boundary is required call self%read_bndry_polygons(get_dommaschk_outer_bndry_file(), & get_dommaschk_outer_bndry_var(), & self%outer_bndry_polygons_phi, & self%outer_bndry_polygons) ! Inner boundary is optional if (get_dommaschk_inner_bndry_file() /= "") then self%use_inner_bndry = .true. call self%read_bndry_polygons(get_dommaschk_inner_bndry_file(), & get_dommaschk_inner_bndry_var(), & self%inner_bndry_polygons_phi, & self%inner_bndry_polygons) else self%use_inner_bndry = .false. endif endif if (self%exclusion_from_file) then call self%read_bndry_polygons(get_dommaschk_exclusion_file(), & get_dommaschk_exclusion_var(), & self%exclusion_polygons_phi, & self%exclusion_polygons) endif if (self%rho_from_file) then call self%read_rho_polygons(get_dommaschk_rho_file()) call self%init_rho_array() endif call self%init_mag_axis_loc() self%initialized = .true. end subroutine init subroutine check_fitting_coef(self) !! Check fitting_coef for consistency. The first condition (from Eq. 12 !! in [*]) is strict, while the second and third (from Eq. 13a) only !! enforce stellarator symmetry, which can be violated. Hence, only a !! warning is provided there. class(dommaschk_t), intent(inout) :: self if (any(self%fitting_coef(3:4, :, 0) /= 0.0_FP)) then call handle_error("Dommaschk fitting coefficients are & &inconsistent!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__, & additional_info=& error_info_t("Violated c_{m,l=0} = & &d_{m,l=0} = 0 (Eq. 12)")) end if if (any(self%fitting_coef([1, 4], :, 0:self%l_pol:2) /= 0.0_FP)) then call handle_error("Dommaschk fitting coefficients violate & &stellarator symmetry!", & PARALLAX_WRN_GENERAL, __LINE__, __FILE__, & additional_info=& error_info_t("Violated a_{m,l=0} = & &d_{m,l=0} = 0 for even l & &(Eq. 13a)")) end if if (any(self%fitting_coef([2, 3], :, 1:self%l_pol:2) /= 0.0_FP)) then call handle_error("Dommaschk fitting coefficients violate & &stellarator symmetry!", & PARALLAX_WRN_GENERAL, __LINE__, __FILE__, & additional_info=& error_info_t("Violated b_{m,l=0} = & &c_{m,l=0} = 0 for odd l & &(Eq. 13a)")) end if end subroutine check_fitting_coef subroutine init_CD_CN(self) !! Compute and store C^D_{m,k}(R) and C^N_{m,k}(R) (Eq. 31 and 32) class(dommaschk_t), intent(inout) :: self integer :: mi, m, k, j integer :: ubound_kloop, ubound_r, ubound_lnr integer :: n_r, n_lnr ! Maximum upper bound of k sum in Eq. 3 ubound_kloop = int(self%l_pol / 2) ! Maximum number of R terms from Eqs. 31 and 32 ubound_r = 2 * (ubound_kloop + 1) - 1 ! Maximum number of log(R) terms from Eqs. 31 and 32 ubound_lnr = ubound_kloop ! Allocate coefficient and power arrays. The first dimension gives m, ! the second k, and the third contains the R and log(R) terms. allocate(self%CD_r_coef(0:self%m_tor_consecutive, & 0:ubound_kloop, & 0:ubound_r)) allocate(self%CN_r_coef, & self%C_r_power, & mold=self%CD_r_coef) self%CD_r_coef = 0.0_FP self%CN_r_coef = 0.0_FP self%C_r_power = 0.0_FP allocate(self%CD_lnr_coef(0:self%m_tor_consecutive, & 0:ubound_kloop, & 0:ubound_lnr)) allocate(self%CN_lnr_coef, & self%C_lnr_power, & mold=self%CD_lnr_coef) self%CD_lnr_coef = 0.0_FP self%CN_lnr_coef = 0.0_FP self%C_lnr_power = 0.0_FP ! Upon differentiation, a logarithm term produces a regular polynomial ! term and another logarithm term, due to the product rule. Thus the ! number of R terms in the derivative is the sum of the number of R ! terms and log(R) terms. allocate(self%dCD_dR_r_coef(0:self%m_tor_consecutive, & 0:ubound_kloop, & 0:ubound_r+ubound_lnr+1)) allocate(self%dCN_dR_r_coef, & self%dC_dR_r_power, & mold=self%dCD_dR_r_coef) self%dCD_dR_r_coef = 0.0_FP self%dCN_dR_r_coef = 0.0_FP self%dC_dR_r_power = 0.0_FP allocate(self%dCD_dR_lnr_coef(0:self%m_tor_consecutive, & 0:ubound_kloop, & 0:ubound_lnr)) allocate(self%dCN_dR_lnr_coef, & self%dC_dR_lnr_power, & mold=self%dCD_dR_lnr_coef) self%dCD_dR_lnr_coef = 0.0_FP self%dCN_dR_lnr_coef = 0.0_FP self%dC_dR_lnr_power = 0.0_FP do mi = 0, self%m_tor_consecutive m = mi * self%num_field_periods do k = 0, ubound_kloop do j = 0, k ! Eq. 31, first term self%CD_lnr_coef(mi, k, j) = -alpha_n(m, j) & * alpha_n_st(m, k - m - j) self%C_lnr_power(mi, k, j) = 2.0_FP * j + m self%CD_r_coef(mi, k, 2 * j) = & -( alpha_n(m, j) * ( gamma_n_st(m, k - m - j) & - alpha_n(m, k - m - j)) & - gamma_n(m, j) * alpha_n_st(m, k - m - j) & + alpha_n(m, j) * beta_n_st(m, k - j)) self%C_r_power(mi, k, 2 * j) = 2.0_FP * j + m ! Second term self%CD_r_coef(mi, k, 2 * j + 1) = beta_n(m, j) & * alpha_n_st(m, k - j) self%C_r_power(mi, k, 2 * j + 1) = 2.0_FP * j - m ! Eq. 32, first term. Note that the powers for CN and CD ! are the same self%CN_lnr_coef(mi, k, j) = alpha_n(m, j) & * alpha_n(m, k - m - j) self%CN_r_coef(mi, k, 2 * j) = & alpha_n(m, j) * gamma_n(m, k - m - j) & - gamma_n(m, j) * alpha_n(m, k - m - j) & + alpha_n(m, j) * beta_n(m, k - j) ! Second term self%CN_r_coef(mi, k, 2 * j + 1) = -beta_n(m, j) & * alpha_n(m, k - j) end do ! Differentiate polynomial terms n_lnr = k n_r = 2 * (n_lnr + 1) - 1 self%dCD_dR_r_coef(mi, k, :n_r) = self%CD_r_coef(mi, k, :n_r) & * self%C_r_power(mi, k, :n_r) self%dCN_dR_r_coef(mi, k, :n_r) = self%CN_r_coef(mi, k, :n_r) & * self%C_r_power(mi, k, :n_r) self%dC_dR_r_power(mi, k, :n_r) = self%C_r_power(mi, k, :n_r) & - 1.0_FP ! Polynomial terms from logarithm differentiation self%dCD_dR_r_coef(mi, k, n_r+1:n_r+n_lnr+1) = & self%CD_lnr_coef(mi, k, :n_lnr) self%dCN_dR_r_coef(mi, k, n_r+1:n_r+n_lnr+1) = & self%CN_lnr_coef(mi, k, :n_lnr) self%dC_dR_r_power(mi, k, n_r+1:n_r+n_lnr+1) = & self%C_lnr_power(mi, k, :n_lnr) - 1.0_FP ! Logarithm terms from logarithm differentiation self%dCD_dR_lnr_coef(mi, k, :) = self%CD_lnr_coef(mi, k, :) & * self%C_lnr_power(mi, k, :) self%dCN_dR_lnr_coef(mi, k, :) = self%CN_lnr_coef(mi, k, :) & * self%C_lnr_power(mi, k, :) self%dC_dR_lnr_power(mi, k, :) = self%C_lnr_power(mi, k, :) & - 1.0_FP end do end do contains pure real(FP) function alpha_n(m, n) !! Eq. 27 in [*] integer, intent(in) :: m, n real(FP) :: num, denom if (n < 0) then ! See notes under Eq. 31 and 32 alpha_n = 0.0_FP else num = (-1.0_FP)**n ! Gamma(n) = factorial(n-1) if n is a positive integer denom = factorial(m + n) * factorial(n) * 2.0_FP**(2 * n + m) alpha_n = num / denom end if end function alpha_n pure real(FP) function alpha_n_st(m, n) !! Eq. 27 in [*] integer, intent(in) :: m, n alpha_n_st = (2 * n + m) * alpha_n(m, n) end function alpha_n_st pure real(FP) function beta_n(m, n) !! Eq. 28 in [*] integer, intent(in) :: m, n integer :: pwr if ((n < 0) .or. (n >= m)) then ! See notes under Eq. 31 and 32 beta_n = 0.0_FP else pwr = 2 * n - m + 1 ! Gamma(n) = (n-1)! if n is a positive integer beta_n = factorial(m - n - 1) / (factorial(n) * 2.0_FP**pwr) end if end function beta_n pure real(FP) function beta_n_st(m, n) !! Eq. 28 in [*] integer, intent(in) :: m, n beta_n_st = (2 * n - m) * beta_n(m, n) end function beta_n_st pure real(FP) function gamma_n(m, n) result(res) !! Eq. 33 in [*] integer, intent(in) :: m, n integer :: i res = 0.0_FP do i = 1, n res = res + 1.0_FP / i + 1.0_FP / (m + i) end do res = 0.5_FP * alpha_n(m, n) * res end function gamma_n pure real(FP) function gamma_n_st(m, n) !! Eq. 33 in [*] integer, intent(in) :: m, n gamma_n_st = (2 * n + m) * gamma_n(m, n) end function gamma_n_st end subroutine init_CD_CN subroutine init_Imn(self) !! Compute and store D_{m,n} and N_{m,n} (as used in Eq. 12 in [*]) !! via Eq. 3 class(dommaschk_t), intent(inout) :: self integer :: mi, n, k integer :: ubound_kloop ! Maximum upper bound of k sum in Eq. 3 ubound_kloop = int(self%l_pol / 2) ! Allocate coefficient and power arrays. The first dimension gives m, ! the second n, and the third contains the Z terms. ! In Eq. 12, N_{m,n} is used until n = l_pol - 1, and the entire N_{m,n} ! term should be zero when l = 0. Thus the array is indexed from -1 in ! the n-dimension and left to zero for the -1 index. allocate(self%I_mn_coef( 0:self%m_tor_consecutive, & -1:self%l_pol, & 0:ubound_kloop)) allocate(self%I_mn_power, & self%dI_dZ_coef, & self%dI_dZ_power, & mold=self%I_mn_coef) self%I_mn_coef = 0.0_FP self%I_mn_power = 0.0_FP self%dI_dZ_coef = 0.0_FP self%dI_dZ_power = 0.0_FP do mi = 0, self%m_tor_consecutive do n = 0, self%l_pol do k = 0, int(n / 2) ! Eq. 3 self%I_mn_coef(mi, n, k) = 1.0_FP / factorial(n - 2 * k) self%I_mn_power(mi, n, k) = n - 2.0_FP * k self%dI_dZ_coef(mi, n, k) = self%I_mn_coef(mi, n, k) & * self%I_mn_power(mi, n, k) self%dI_dZ_power(mi, n, k) = self%I_mn_power(mi, n, k) & - 1.0_FP ! Correct for negative powers of Z when the leading ! coefficient is zero. Otherwise for the case Z = 0.0, NaN ! will be returned if (self%dI_dZ_power(mi, n, k) < 0.0_FP .and. & abs(self%dI_dZ_coef(mi, n, k)) < FP_EPS) then self%dI_dZ_power(mi, n, k) = 0.0_FP elseif (self%dI_dZ_power(mi, n, k) < 0.0_FP .and. & abs(self%dI_dZ_coef(mi, n, k)) > FP_EPS) then call handle_error("A negative power of Z exists!", & PARALLAX_ERR_EQUILIBRIUM, & __LINE__, __FILE__, & additional_info=& error_info_t("dI_dZ_power:", & [self%dI_dZ_power(mi, n, k)])) endif end do end do end do end subroutine init_Imn subroutine init_B_norm(self) !! Determine the magnetic field normalization with the un-normalized !! value of btor evaluated at the location of the magnetic axis. This is !! used to normalize all subsequent magnetic field calculations. class(dommaschk_t), intent(inout) :: self real(FP) :: btor_unnormalized ! Set the normalization factor to 1 to calculate the un-normalized ! value of btor self%B_norm_axis = 1.0_FP btor_unnormalized = self%btor(self%x0, self%y0, self%phi0) self%B_norm_axis = btor_unnormalized end subroutine init_B_norm subroutine init_rho_array(self) !! Calculates the value of rho for every surface in the given rho file. !! Used later in rho calculation, to interpolate between surfaces class(dommaschk_t), intent(inout) :: self integer :: n_vertices, n_surfaces, i type(closed_polygon2d_t) :: surf_polygon real(FP) :: area_last_surface n_vertices = size(self%rho_polygons, 1) n_surfaces = size(self%rho_polygons, 4) surf_polygon%n_pts = n_vertices allocate(surf_polygon%pts(n_vertices, 2)) allocate(self%rho_array(n_surfaces)) ! For the first plane (phi = 0.0), iterate over all surfaces from file ! and calculate rho as the ratio of the surface areas of the current ! surface and the last flux surface do i = n_surfaces, 1, -1 surf_polygon%pts = self%rho_polygons(:, :, 1, i) if (i == n_surfaces) then area_last_surface = abs(surf_polygon%signed_area()) self%rho_array(i) = 1.0_FP else self%rho_array(i) = abs(surf_polygon%signed_area()) & / area_last_surface endif enddo end subroutine init_rho_array subroutine init_mag_axis_loc(self) !! Initializes two arrays (equidistant in phi) of x- and !! y-coordinates of magnetic axis within the first field period !! calculated via field line tracing. These are used as data !! points for fast interpolation in 'mag_axis_loc' class(dommaschk_t), intent(inout) :: self !! Instance of class real(FP) :: phi_max, phi_spacing, dphi_tmp, phi_trace_start integer :: i, trace_status integer, parameter :: num_of_data_points = 128 ! Hardcoded number of data points for the interpolation phi_max = TWO_PI / self%num_field_periods phi_spacing = phi_max / num_of_data_points allocate(self%axis_array_x(num_of_data_points)) allocate(self%axis_array_y(num_of_data_points)) ! Make sure that tracing happens within the first field period phi_trace_start = modulo(self%phi0, phi_max) do i = 1, num_of_data_points dphi_tmp = (i - 1) * phi_spacing - phi_trace_start call trace(x_init=self%x0, & y_init=self%y0, & phi_init=phi_trace_start, & dphi=dphi_tmp, & equi=self, & x_end=self%axis_array_x(i), & y_end=self%axis_array_y(i), & istat=trace_status) if (trace_status /= PARALLAX_SUCCESS) then exit endif enddo ! NOTE: This step of checking the status is currently just necessary for ! the Dommaschk 'base_testcase' which is purely for checking ! Dommaschk coefficients. This equilibrium does not have nested ! flux surfaces. ! For any realistic equilibrium the tracing routine must not fail. if (trace_status /= PARALLAX_SUCCESS) then do i = 1, num_of_data_points self%axis_array_x(i) = self%x0 self%axis_array_y(i) = self%y0 enddo call handle_error("Tracing the magnetic axis failed!", & PARALLAX_WRN_GENERAL, __LINE__, __FILE__, & additional_info=error_info_t("Returning coords x0 and y0 & &instead of actual axis location.")) endif end subroutine init_mag_axis_loc subroutine display(self) class(dommaschk_t), intent(in) :: self write(get_stdout(), 99) & self%num_field_periods, self%m_tor_consecutive, self%l_pol 99 format( / & 1X,'DOMMASCHK equilibrium:'/, & 1X,'num_field_periods = ',I4 /, & 1X,'m_tor_consecutive = ',I4 /, & 1X,'l_pol = ',I4 / & ) end subroutine display subroutine debug(self) class(dommaschk_t), intent(in) :: self write(get_stdout(), 100) & self%x0, self%y0, self%phi0, self%rhomin, self%rhomax, & self%xmin, self%xmax, self%ymin, self%ymax, & self%num_field_periods, self%m_tor_consecutive, self%l_pol, & self%B_norm_axis, & self%bndry_from_file, self%exclusion_from_file, self%rho_from_file 100 format( / & 1X,'DOMMASCHK equilibrium:'/, & 1X,'x0 = ',ES14.6E3 /, & 1X,'y0 = ',ES14.6E3 /, & 1X,'phi0 = ',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,'num_field_periods = ',I4 /, & 1X,'m_tor_consecutive = ',I4 /, & 1X,'l_pol = ',I4 /, & 1X,'B_norm_axis = ',ES14.6E3 /, & 1X,'bndry_from_file = ',L /, & 1X,'exclusion_from_file = ',L /, & 1X,'rho_from_file = ',L / & ) end subroutine debug logical function is_axisymmetric(self) class(dommaschk_t), intent(in) :: self is_axisymmetric = .false. end function real(FP) function rho(self, x, y, phi) class(dommaschk_t), intent(in) :: self real(FP), intent(in) :: x, y, phi logical, save :: warning_printed = .false. real(FP) :: phi_mod type(closed_polygon2d_t), dimension(:), allocatable :: & polygons_interp_phi real(FP), dimension(:,:), allocatable :: tmp_vertices integer :: n_planes, n_vertices, n_surfaces, i real(FP) :: dphi, interp_frac integer :: plane_left, plane_right real(FP) :: dist_inner, dist_outer if(.not. self%rho_from_file) then if(.not. warning_printed) then call handle_error("No Dommaschk rho file provided!", & PARALLAX_WRN_GENERAL, __LINE__, __FILE__, & additional_info=error_info_t("Returning rho = 2.0")) warning_printed = .true. endif rho = 2.0_FP return endif phi_mod = modulo(phi, TWO_PI / self%num_field_periods) n_vertices = size(self%rho_polygons, 1) n_planes = size(self%rho_polygons, 3) n_surfaces = size(self%rho_polygons, 4) allocate(polygons_interp_phi(n_surfaces)) allocate(tmp_vertices(n_vertices, 2)) ! Interpolate flux surface polygons at given phi do i = n_surfaces, 1, -1 ! If only one polygon specified, no need to interpolate in phi if (size(self%rho_polygons_phi) == 1) then tmp_vertices = self%rho_polygons(:, :, 1, i) ! Interpolate polygons to get one located at the given phi else ! For polygon array, assume constant spacing between planes dphi = self%rho_polygons_phi(2) - self%rho_polygons_phi(1) ! Index of nearest plane to the left of the given phi, assuming ! that the planes are ordered and have even spacing plane_left = int((phi_mod - self%rho_polygons_phi(1)) / dphi) + 1 ! Index of nearest plane to the right, wrapping around to first plane_right = modulo(plane_left, n_planes) + 1 ! How close given phi is to left plane interp_frac = (phi_mod - self%rho_polygons_phi(plane_left)) / dphi ! Basic linear interpolation between right and left polygon tmp_vertices = self%rho_polygons(:, :, plane_left, i) & + ( self%rho_polygons(:, :, plane_right, i) & - self%rho_polygons(:, :, plane_left, i)) * interp_frac endif ! Need to call create_polygon2d here to generate internal ! arrays used in dist_to_polyedge call polygons_interp_phi(i)%create_polygon2d(n_vertices, & tmp_vertices(:, 1), tmp_vertices(:, 2)) ! Break loop at the first surface which does not enclose the given ! point. This is called the inner surface - the one just previous ! is called the outer surface if(.not. polygons_interp_phi(i)%pt_inside(x, y)) then exit endif enddo if(i == n_surfaces) then ! Point is outside of last flux surface given in file ! NOTE: Depending on the wall polygon, open field lines can have ! rho < 1! rho = 1.01_FP return endif if(i == 0) then ! Point is inside first flux surface (magnetic axis) rho = 0.0_FP return endif dist_inner = polygons_interp_phi(i )%dist_to_polyedge(x, y) dist_outer = polygons_interp_phi(i + 1)%dist_to_polyedge(x, y) interp_frac = dist_inner / (dist_inner + dist_outer) ! Interpolate rho between inner and outer surface rho = self%rho_array(i) & + (self%rho_array(i + 1) - self%rho_array(i)) * interp_frac end function rho real(FP) function bx(self, x, y, phi) !! Evaluates the radial component of the vacuum magnetic field B !! according to B_R = dV/dR class(dommaschk_t), intent(in) :: self real(FP), intent(in) :: x, y, phi integer :: mi, m, l, k, j real(FP) :: logR, cos_m_phi, sin_m_phi real(FP) :: dV_dR, dD_ml_dR, dN_ml_1_dR, dCD_mk_dR, dCN_mk_dR logR = log(x) bx = 0.0_FP do mi = 0, self%m_tor_consecutive m = mi * self%num_field_periods cos_m_phi = cos(m * phi) sin_m_phi = sin(m * phi) do l = 0, self%l_pol dD_ml_dR = 0.0_FP do k = 0, int(l / 2) dCD_mk_dR = 0.0_FP do j = 0, k dCD_mk_dR = dCD_mk_dR & + self%dCD_dR_r_coef(mi, k, 3 * j) & * x**self%dC_dR_r_power(mi, k, 3 * j) dCD_mk_dR = dCD_mk_dR & + self%dCD_dR_r_coef(mi, k, 3 * j + 1) & * x**self%dC_dR_r_power(mi, k, 3 * j + 1) dCD_mk_dR = dCD_mk_dR & + self%dCD_dR_r_coef(mi, k, 3 * j + 2) & * x**self%dC_dR_r_power(mi, k, 3 * j + 2) dCD_mk_dR = dCD_mk_dR & + self%dCD_dR_lnr_coef(mi, k, j) & * x**self%dC_dR_lnr_power(mi, k, j) * logR end do dD_ml_dR = dD_ml_dR + self%I_mn_coef(mi, l, k) & * y**self%I_mn_power(mi, l, k) & * dCD_mk_dR end do dN_ml_1_dR = 0.0_FP do k = 0, int((l - 1) / 2) dCN_mk_dR = 0.0_FP do j = 0, k dCN_mk_dR = dCN_mk_dR & + self%dCN_dR_r_coef(mi, k, 3 * j) & * x**self%dC_dR_r_power(mi, k, 3 * j) dCN_mk_dR = dCN_mk_dR & + self%dCN_dR_r_coef(mi, k, 3 * j + 1) & * x**self%dC_dR_r_power(mi, k, 3 * j + 1) dCN_mk_dR = dCN_mk_dR & + self%dCN_dR_r_coef(mi, k, 3 * j + 2) & * x**self%dC_dR_r_power(mi, k, 3 * j + 2) dCN_mk_dR = dCN_mk_dR & + self%dCN_dR_lnr_coef(mi, k, j) & * x**self%dC_dR_lnr_power(mi, k, j) * logR end do dN_ml_1_dR = dN_ml_1_dR + self%I_mn_coef(mi, l - 1, k) & * y**self%I_mn_power(mi, l - 1, k) & * dCN_mk_dR end do ! Eq. 12 (differentiated with respect to R) dV_dR = ( self%fitting_coef(1, mi, l) * cos_m_phi & + self%fitting_coef(2, mi, l) * sin_m_phi & ) * dD_ml_dR & + ( self%fitting_coef(3, mi, l) * cos_m_phi & + self%fitting_coef(4, mi, l) * sin_m_phi & ) * dN_ml_1_dR bx = bx + dV_dR end do end do ! Normalize field magnitude bx = bx / self%B_norm_axis end function bx real(FP) function by(self, x, y, phi) !! Evaluates the radial component of the vacuum magnetic field B !! according to B_Z = dV/dZ class(dommaschk_t), intent(in) :: self real(FP), intent(in) :: x, y, phi integer :: mi, m, l, k, j real(FP) :: logR, cos_m_phi, sin_m_phi real(FP) :: dV_dZ, dD_ml_dZ, dN_ml_1_dZ, CD_mk, CN_mk logR = log(x) by = 0.0_FP do mi = 0, self%m_tor_consecutive m = mi * self%num_field_periods cos_m_phi = cos(m * phi) sin_m_phi = sin(m * phi) do l = 0, self%l_pol dD_ml_dZ = 0.0_FP do k = 0, int(l / 2) CD_mk = 0.0_FP do j = 0, k CD_mk = CD_mk & + self%CD_r_coef(mi, k, 2 * j) & * x**self%C_r_power(mi, k, 2 * j) CD_mk = CD_mk & + self%CD_r_coef(mi, k, 2 * j + 1) & * x**self%C_r_power(mi, k, 2 * j + 1) CD_mk = CD_mk & + self%CD_lnr_coef(mi, k, j) & * x**self%C_lnr_power(mi, k, j) * logR end do dD_ml_dZ = dD_ml_dZ + self%dI_dZ_coef(mi, l, k) & * y**self%dI_dZ_power(mi, l, k) * CD_mk end do dN_ml_1_dZ = 0.0_FP do k = 0, int((l - 1) / 2) CN_mk = 0.0_FP do j = 0, k CN_mk = CN_mk & + self%CN_r_coef(mi, k, 2 * j) & * x**self%C_r_power(mi, k, 2 * j) CN_mk = CN_mk & + self%CN_r_coef(mi, k, 2 * j + 1) & * x**self%C_r_power(mi, k, 2 * j + 1) CN_mk = CN_mk & + self%CN_lnr_coef(mi, k, j) & * x**self%C_lnr_power(mi, k, j) * logR end do dN_ml_1_dZ = dN_ml_1_dZ + self%dI_dZ_coef(mi, l - 1, k) & * y**self%dI_dZ_power(mi, l - 1, k) & * CN_mk end do ! Eq. 12 (differentiated with respect to Z) dV_dZ = ( self%fitting_coef(1, mi, l) * cos_m_phi & + self%fitting_coef(2, mi, l) * sin_m_phi & ) * dD_ml_dZ & + ( self%fitting_coef(3, mi, l) * cos_m_phi & + self%fitting_coef(4, mi, l) * sin_m_phi & ) * dN_ml_1_dZ by = by + dV_dZ end do end do ! Normalize field magnitude by = by / self%B_norm_axis end function by real(FP) function btor(self, x, y, phi) !! Evaluates the poloidal component of the vacuum magnetic field B !! according to B_phi = 1/R dV/dphi class(dommaschk_t), intent(in) :: self real(FP), intent(in) :: x, y, phi integer :: mi, m, l, k, j real(FP) :: logR, cos_m_phi, sin_m_phi real(FP) :: dV_dphi, D_ml, N_ml_1, CD_mk, CN_mk logR = log(x) ! Derivative of phi term in Eq. 1 btor = 1.0_FP do mi = 0, self%m_tor_consecutive m = mi * self%num_field_periods cos_m_phi = cos(m * phi) sin_m_phi = sin(m * phi) do l = 0, self%l_pol D_ml = 0.0_FP do k = 0, int(l / 2) CD_mk = 0.0_FP do j = 0, k CD_mk = CD_mk & + self%CD_r_coef(mi, k, 2 * j) & * x**self%C_r_power(mi, k, 2 * j) CD_mk = CD_mk & + self%CD_r_coef(mi, k, 2 * j + 1) & * x**self%C_r_power(mi, k, 2 * j + 1) CD_mk = CD_mk & + self%CD_lnr_coef(mi, k, j) & * x**self%C_lnr_power(mi, k, j) * logR end do D_ml = D_ml + self%I_mn_coef(mi, l, k) & * y**self%I_mn_power(mi, l, k) * CD_mk end do N_ml_1 = 0.0_FP do k = 0, int((l - 1) / 2) CN_mk = 0.0_FP do j = 0, k CN_mk = CN_mk & + self%CN_r_coef(mi, k, 2 * j) & * x**self%C_r_power(mi, k, 2 * j) CN_mk = CN_mk & + self%CN_r_coef(mi, k, 2 * j + 1) & * x**self%C_r_power(mi, k, 2 * j + 1) CN_mk = CN_mk & + self%CN_lnr_coef(mi, k, j) & * x**self%C_lnr_power(mi, k, j) * logR end do N_ml_1 = N_ml_1 + self%I_mn_coef(mi, l - 1, k) & * y**self%I_mn_power(mi, l - 1, k) * CN_mk end do ! Eq. 12 (differentiated with respect to phi) dV_dphi = (-self%fitting_coef(1, mi, l) * m * sin_m_phi & + self%fitting_coef(2, mi, l) * m * cos_m_phi & ) * D_ml & + (-self%fitting_coef(3, mi, l) * m * sin_m_phi & + self%fitting_coef(4, mi, l) * m * cos_m_phi & ) * N_ml_1 btor = btor + dV_dphi end do end do btor = btor / x ! Normalize field magnitude btor = btor / self%B_norm_axis end function btor real(FP) function jacobian(self, x, y, phi) class(dommaschk_t), intent(in) :: self real(FP), intent(in) :: x, y, phi jacobian = x end function jacobian subroutine epol(self, x, y, phi, epolx, epoly) class(dommaschk_t), intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP), intent(out) :: epolx, epoly real(FP) :: bx, by, bpol bx = self%bx(x, y, phi) by = self%by(x, y, phi) bpol = sqrt(bx * bx + by * by) epolx = bx / bpol epoly = by / bpol end subroutine epol subroutine erad(self, x, y, phi, eradx, erady) class(dommaschk_t), intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP), intent(out) :: eradx, erady real(FP) :: bx, by, bpol bx = self%bx(x, y, phi) by = self%by(x, y, phi) bpol = sqrt(bx * bx + by * by) eradx = by / bpol erady = -bx / bpol end subroutine erad integer function district(self, x, y, phi) class(dommaschk_t), intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP) :: phi_mod phi_mod = modulo(phi, TWO_PI / self%num_field_periods) ! Check for domain exclusion if (self%exclusion_from_file) then if (.not. point_inside_polygons(x, y, phi_mod, & self%exclusion_polygons_phi, & self%exclusion_polygons) ) then district = DISTRICT_OUT return endif else ! Naive box limits if ( (x <= self%xmin) .or. (x >= self%xmax) .or. & (y <= self%ymin) .or. (y >= self%ymax) ) then district = DISTRICT_OUT return endif endif if (.not.self%bndry_from_file) then district = DISTRICT_CLOSED return endif if (self%use_inner_bndry) then if (point_inside_polygons(x, y, phi_mod, & self%inner_bndry_polygons_phi, & self%inner_bndry_polygons)) then district = DISTRICT_CORE return endif endif if (point_inside_polygons(x, y, phi_mod, & self%outer_bndry_polygons_phi, & self%outer_bndry_polygons)) then district = DISTRICT_CLOSED else district = DISTRICT_WALL endif end function district logical function in_vessel(self, x, y, phi) class(dommaschk_t), intent(in) :: self real(FP), intent(in) :: x, y, phi integer :: district district = self%district(x, y, phi) if ( (district == DISTRICT_CLOSED) .or. & (district == DISTRICT_CORE) ) then in_vessel = .true. else in_vessel = .false. endif end function in_vessel subroutine mag_axis_loc(self, phi, axis_x, axis_y) !! Returns the coordinates of magnetic axis class(dommaschk_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 real(FP) :: phi_max, phi_spacing, phi_mod integer :: num_of_samples integer, parameter :: intorder=3 num_of_samples = size(self%axis_array_x) phi_max = TWO_PI / self%num_field_periods phi_spacing = phi_max / num_of_samples phi_mod = modulo(phi, phi_max) ! Interpolate coordinates of the axis axis_x = interpol1d(intorder, num_of_samples, phi_spacing, & self%axis_array_x, phi_mod) axis_y = interpol1d(intorder, num_of_samples, phi_spacing, & self%axis_array_y, phi_mod) end subroutine mag_axis_loc logical function point_inside_polygons(x, y, phi_mod, & phi_array, polygon_array) !! Determine if a given point (x, y, phi_mod) (phi_mod lies within !! the first field period) lies within the given 2D polygons. If the !! polygon array size is larger than one, this is done using simple !! linear interpolation between adjacent planes. real(FP), intent(in) :: x !! x-coordinate real(FP), intent(in) :: y !! y-coordinate real(FP), intent(in) :: phi_mod !! Phi_coordinate mapped back to first field period real(FP), dimension(:), intent(in) :: phi_array !! Phi locations where polygons are defined real(FP), dimension(:,:,:), intent(in) :: polygon_array !! Polygon vertices type(closed_polygon2d_t) :: interp_polygon integer :: n_vertices real(FP) :: dphi, interp_frac integer :: plane_left, plane_right n_vertices = size(polygon_array, 1) ! If only one polygon specified, no need to interpolate in phi if (size(phi_array) == 1) then interp_polygon%N_pts = n_vertices allocate(interp_polygon%pts(n_vertices, 2)) interp_polygon%pts = polygon_array(:, :, 1) ! Interpolate polygons to get one located at the given phi else interp_polygon%N_pts = n_vertices allocate(interp_polygon%pts(n_vertices, 2)) ! For polygon array, assume constant spacing between planes dphi = phi_array(2) - phi_array(1) ! Index of nearest plane to the left of the given phi, assuming ! that the planes are ordered and have even spacing plane_left = int((phi_mod - phi_array(1)) / dphi) + 1 ! Index of nearest plane to the right, wrapping around to first plane_right = modulo(plane_left, size(phi_array)) + 1 ! How close given phi is to left plane interp_frac = (phi_mod - phi_array(plane_left)) / dphi ! Basic linear interpolation between right and left polygon interp_polygon%pts = polygon_array(:, :, plane_left) & + ( polygon_array(:, :, plane_right) & - polygon_array(:, :, plane_left)) * interp_frac endif point_inside_polygons = interp_polygon%pt_inside(x, y) end function end module dommaschk_equilibrium_m