module numerical_equilibrium_m !****************************************************************************** ! Defines the `numerical` equilibrium subclass ! N.b. this is axisymmetric; phi, when required, is not used. ! ! The object structure is quite large -- it is designed to be fully self-descriptive, such ! that multiple equilibria may be developed independently, without resulting in conflicts between ! the equilibria objects. As a result, most variables are stored within the object (i.e. as ! object%variable) rather than shared across the module. ! ! The `numerical` equilibrium object develops an equilibrium object from a NetCDF, specified ! in the params.in file, as "${path_to_netcdf}/${equilibrium_case}.nc". This file should define ! ! * Global attributes "description", "history", "author" ! * A group "Magnetic_geometry" with ! * dimensions "R", "Z" ! * variables "R", "Z" (corresponding dimension values) and "psi" (R,Z grid). The ! "btor" (R,Z grid) variable is optional (unused) if use_spline_for_btor == .false. ! (in params.in), required if use_spline_for_btor == .true. ! * attributes "magnetic_axis_R", "magnetic_axis_Z", "x_point_R", "x_point_Z", "axis_Btor" ! * A group "Psi_limits" with ! * attributes "global_min" and "global_min" ! * subgroups (name not used) corresponding to regions (polygons) where a different psi limit ! should be applied. Each subgroup should have ! * attributes "local_min" and "local_max", which may be NaN (IEEE conformant) if undefined ! * a nc_polygon (see below) ! * A group "Box_limits" with a corresponding nc_polygon ! * A group "FW_and_divertor_polygon" with a corresponding nc_polygon ! * A group "Shadow_polygon" with a corresponding nc_polygon ! ! Groups which contain a nc_polygon should have the following information ! * dimension "N_points" which gives the number of polygon points ! * variables "R_points" and "Z_points" which gives the (R, Z) points ! Note that the R_points and Z_points will be automatically normalised to the magnetic_axis_R ! ! N.b. this module relies on bspline_fortran module from ! https://github.com/jacobwilliams/bspline-fortran ! !****************************************************************************** use, intrinsic :: iso_fortran_env use netcdf use bspline_module use polygon_m, only: closed_polygon2d_t, limiting_polygon2d_t use precision_m, only : FP, FP_EPS use screen_io_m, only : get_stdout, get_stderr, nf90_handle_err use comm_handling_m, only: is_master use constants_m, only : pi, two_pi use descriptors_m use divertor_equilibrium_m, only : divertor_equilibrium_t use elementary_functions_m, only : smoothstep => step_hermite implicit none private integer, parameter :: string_buffer_length = 200 ! Maximum length of variable-length strings integer, parameter :: flux_regions_buffer = 10 ! Maximum number of flux-limit polygon regions real(FP), parameter :: spline_edge_buffer = 20.0_FP ! Number of points to go to zero poloidal field over, as the trace approaches the spline edge type, public, extends(divertor_equilibrium_t) :: numerical_t ! Private class variables private character(len = string_buffer_length) :: equilibrium_case character(len = string_buffer_length) :: path_to_netcdf character(len = string_buffer_length) :: params_filename real(FP), public :: R0 ! radial position of the magnetic axis from the NetCDF, used to normalise spatial dimension integer :: spline_order ! What order spline to use for interpolation -- needs to be high enough to allow for derivatives. Set via params.in type(bspline_2d) :: psi_spl2d ! NOTE: Units must be Wb=V*s=T*m2 logical :: psi_spl2d_extrap = .false. ! Spline interpolator (2D) for poloidal flux. Dimensions are normalised to self%R0 real(FP) :: spl_limits_R(2), spl_limits_Z(2) ! Extent of interpolator grid (extrapolation is not allowed) as (min, max) array integer :: spl_size_R, spl_size_Z ! Number of points on the interpolator grid in each dimension real(FP) :: rhomin_privflux_params = 0.0_FP ! Allows the user to specify a region-specific for any specially handled regions. Only recommended for cases ! with a single X-point. Default to minimum possible value (rho is always normalised to positive values). type(limiting_polygon2d_t), dimension(:), allocatable :: flux_regions ! Array of polygons, for holding polygons defining regions which require specific flux limiting integer :: num_flux_regions ! Actual number of flux-limit polygon regions logical :: set_divertor_polygon ! Use divertor polygon from NetCDF if false (default), or from parameter ! file if true logical :: set_exclusion_polygon ! Use exclusion polygon from NetCDF if false (default), or from parameter ! file if true type(closed_polygon2d_t) :: spline_edge_polygon ! Closed polygon which identifies the edge of the spline grid character(len=string_buffer_length) :: ncdf_description, ncdf_history, ncdf_author ! Strings which contains meta-information about the netcdf logical :: set_rhomin, set_rhomax, set_rhomin_privflux ! Controls for whether rho limits should be taken from the parameter file logical :: raise_exception_off_grid ! Default false. If true, if a point off-grid is queried then an exception will be raised contains procedure, public :: init !! Initialises equilibrium, i.e.~reads any required parameters from file. procedure, public :: display !! Print to console information about the equilibrium procedure, public :: debug !! Print to console extended information about the equilibrium procedure, public :: psi !! Non-normalised poloidal flux function (in Wb=V*s=T*m2) procedure, public :: bx !! Radial field component normalized to on axis field strength. procedure, public :: by !! Vertical field component normalized to on axis field strength. procedure, public :: btor !! Toroidal field component normalized to on axis field strength. procedure, public :: check_privflux_regions ! procedure, public :: rholimit_privflux !! Returns the equivalent to 'rholimit_privflux' for a given (x,y) point ! procedure, public, nopass :: convert_psi_to_rho !! Convert psi value to normalised value (rho) ! procedure, public, nopass :: convert_rho_to_psi !! Convert normalised psi value (rho) to non-normalised psi end type numerical_t ! Define initialisation routines in initialise_numerical_equilibrium submodule interface module subroutine read_params_in(self, filename) class(numerical_t), intent(inout) :: self character(len = *), intent(in), optional :: filename end subroutine read_params_in module subroutine read_from_ncdf(self) class(numerical_t), intent(inout) :: self end subroutine read_from_ncdf end interface contains subroutine init(self, filename) ! Initialises equilibrium from a NetCDF file ! ! Reads/sets equilibrium parameters from params.in ! Reads NetCDF specified and builds the equilibrium from this class(numerical_t), intent(inout) :: self character(len = *), intent(in), optional :: filename if (present(filename)) then call read_params_in(self, filename) else call read_params_in(self) endif call read_from_ncdf(self) if (self%exclusion_polygon%signed_area() <= 0.0_FP) then write(get_stderr(), *) & "Error: exclusion polygon is defined clockwise, or is empty", & self%exclusion_polygon%signed_area() ERROR STOP ERR_UNHANDLED endif if (self%divertor_polygon%signed_area() <= 0.0_FP) then write(get_stderr(), *) & "Error: divertor polygon is defined clockwise, or is empty", & self%divertor_polygon%signed_area() ERROR STOP ERR_UNHANDLED endif ! Mark the equilibrium as initialized and ready to use self%initialized = .true. end subroutine init real(FP) function psi(self, x, y, phi) ! Evaluate the poloidal flux, actual value (don't take derivatives so supply dx = 0, dy = 0) ! If evaluating off-grid, use a nearest-neighbour approximation to psi class(numerical_t), intent(in) :: self real(FP), intent(in) :: x, y, phi ! 3D position (x and y normalized) integer :: status CHARACTER(len=:), ALLOCATABLE :: msg real(FP) :: x_in_grid, y_in_grid real(FP), parameter :: eps = 1E-10 !Small number which ensures that the end result is actually on the grid (rather than a border point) if (self%raise_exception_off_grid) then ! Don't use the nearest-neighbour back to the spline grid. Instead, if a point off grid is queried, raise an error x_in_grid = x y_in_grid = y else ! Use nearest-neighbour extrapolation to find the nearest on-grid point ! self%spl_limits_R(1) = minval(R), self%spl_limits_R(2) = maxval(R) ! self%spl_limits_Z(1) = minval(Z), self%spl_limits_Z(2) = maxval(Z) ! Find the nearest (x, y) point in grid x_in_grid = min(max(x, self%spl_limits_R(1) + eps), self%spl_limits_R(2) - eps) y_in_grid = min(max(y, self%spl_limits_Z(1) + eps), self%spl_limits_Z(2) - eps) endif call self%psi_spl2d%evaluate(x_in_grid, y_in_grid, 0, 0, psi, status) if (status /= 0) then call get_status_message(status, msg) write(get_stderr(), *) & 'Error evaluating 2D spline (numerical_equilibrium_t::psi): '& // msg // "at ", x, y ERROR STOP ERR_UNHANDLED endif end function psi subroutine check_privflux_regions(self, x, y, phi, local_rhomin, local_rhomin_exists, local_rhomax, local_rhomax_exists) class(numerical_t), intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP), intent(inout) :: local_rhomax, local_rhomin logical, intent(out) :: local_rhomin_exists, local_rhomax_exists integer :: region_it local_rhomin_exists = .false. local_rhomax_exists = .false. ! Iterate over all the regions, check whether the point is inside the polygon ! N.b. this can be expensive if the polygon is complicated! do region_it = 1, self%num_flux_regions if (self%flux_regions(region_it)%pt_inside(x, y)) then if (self%flux_regions(region_it)%use_local_min) then ! Use the maximum of the global and local maximum (i.e. more restrictive limit) local_rhomin = max(self%flux_regions(region_it)%local_min, local_rhomin) local_rhomin_exists = .true. endif if (self%flux_regions(region_it)%use_local_max) then ! Use the minimum of the global and local maximum (i.e. more restrictive limit) local_rhomax = min(self%flux_regions(region_it)%local_max, local_rhomax) local_rhomax_exists = .true. endif endif end do end subroutine check_privflux_regions real(FP) function bx(self, x, y, phi) ! Evaluate the poloidal flux, taking the 0th x derivative and the 1st y derivative ! Returns zero if the point is off the spline grid class(numerical_t), intent(in) :: self real(FP), intent(in) :: x, y, phi integer :: status real(FP) :: psi_dz CHARACTER(len=:), ALLOCATABLE :: msg call self%psi_spl2d%evaluate(x, y, 0, 1, psi_dz, status) if (((status == 601) .or. (status == 602)) .and. (.not.(self%raise_exception_off_grid))) then ! Evaluating off the spline grid bx = 0 elseif (status /= 0) then call get_status_message(status, msg) write(get_stderr(), *) & 'Error evaluating 2D spline (numerical_equilibrium_t::bx): '& // msg // "at ", x, y, status ERROR STOP ERR_UNHANDLED else ! N.b. One factor of self%R0 comes from evaluating the derivative on the ! (x', y') normalised coordinate and the second comes from using ! normalised x', y' for the axis bx = -self%poloidal_field_factor * psi_dz / (two_pi * x * self%R0 * self%R0 * self%axis_Btor) ! Reduce to toroidal-only field as trace approaches spline edge bx = bx * (1.0_FP - & smoothstep(spline_edge_buffer/2.0_FP, self%spline_edge_polygon%dist_to_polyedge(x, y), spline_edge_buffer)) endif end function bx real(FP) function by(self, x, y, phi) ! Evaluate the poloidal flux, taking the 1st x derivative and the 0th y derivative ! Returns zero if the point is off the spline grid class(numerical_t), intent(in) :: self real(FP), intent(in) :: x, y, phi integer :: status real(FP) :: psi_dr CHARACTER(len=:), ALLOCATABLE :: msg call self%psi_spl2d%evaluate(x, y, 1, 0, psi_dr, status) if (((status == 601) .or. (status == 602)) .and. (.not.(self%raise_exception_off_grid))) then ! Evaluating off the spline grid by = 0 elseif (status /= 0) then call get_status_message(status, msg) write(get_stderr(), *) & 'Error evaluating 2D spline (numerical_equilibrium_t::by): '& // msg // "at ", x, y, status ERROR STOP ERR_UNHANDLED else ! N.b. One factor of self%R0 comes from evaluating the derivative on the ! (x', y') normalised coordinate and the second comes from using ! normalised x', y' for the axis by = self%poloidal_field_factor * psi_dr / (two_pi * x * self%R0 * self%R0 * self%axis_Btor) ! Reduce to toroidal-only field as trace approaches spline edge by = by * (1.0_FP - & smoothstep(spline_edge_buffer/2.0_FP, self%spline_edge_polygon%dist_to_polyedge(x, y), spline_edge_buffer)) endif end function by real(FP) function btor(self, x, y, phi) ! TODO: The option due to use a spline interpolator has been removed ! Even when self%use_spline_for_btor = .false., this function caused a severe ! performance hit. Possible that spline was being evaluated via branch prediction ! resulting in poor thread synchronisation. Should investigate. class(numerical_t), intent(in) :: self real(FP), intent(in) :: x, y, phi ! integer :: status btor = 1.0_FP / x ! call self%btor_spl2d%evaluate(x, y, 0, 0, btor_spl, status) ! if (status /= 0) then ! write(stdout, *)'Error evaluating 2D spline (numerical_t_equilibrium::btor_spl): '& ! //get_status_message(status) // "at ", x, y ! ERROR STOP ERR_UNHANDLED ! endif ! btor = btor_spl / self%axis_Btor end function btor subroutine display(self) class(numerical_t), intent(in) :: self character(len=500) :: flux_regions_print character(len = string_buffer_length) :: region_it_print integer :: region_it flux_regions_print = "" do region_it = 1, self%num_flux_regions write(region_it_print, '(I4)') region_it write(flux_regions_print, '(A)') & trim(flux_regions_print) // trim(region_it_print) // ": " // & & trim(self%flux_regions(region_it)%polygon_string()) // ", " end do write(get_stdout(), 99) & trim(self%equilibrium_case) ,& self%ncdf_history ,& self%params_filename ,& self%rhomax, self%rhomin ,& trim(flux_regions_print) ,& trim(self%divertor_polygon%polygon_string()) ,& trim(self%exclusion_polygon%polygon_string()) 99 FORMAT(& 1X, 'NUMERICAL equilibrium ' /, & 1X, 'equilibrium case: ', 3X, A /, & 1X, 'equi netcdf history: ', 3X, A /, & 1X, 'params read from file: ', 3X, A /, & 1X, 'rho limits -- max:', ES10.3E2, ', min:', ES10.3E2 /, & 1X, '<<Psi limit regions>>' /, & 1X, 'Regions: ', 3X, A /, & 1X, '<<First wall + divertor polygon>>' /, & 1X, A /, & 1X, '<<Exclusion polygon>>' /, & 1X, A / & /) end subroutine display subroutine debug(self) class(numerical_t), intent(in) :: self character(len=500) :: flux_regions_print character(len = string_buffer_length) :: region_it_print integer :: region_it, spline_size_bits call self%psi_spl2d%size_of(spline_size_bits) flux_regions_print = "" do region_it = 1, self%num_flux_regions write(region_it_print, '(I4)') region_it write(flux_regions_print, '(A)') trim(flux_regions_print) // trim(region_it_print) // ": " // & & trim(self%flux_regions(region_it)%polygon_string()) // ", " end do write(get_stdout(), 101) & "NUMERICAL_t" ,& self%params_filename ,& trim(self%equilibrium_case) ,& trim(self%path_to_netcdf) ,& self%ncdf_description ,& self%ncdf_history ,& self%ncdf_author ,& self%flip_Z ,& self%poloidal_field_factor ,& self%spline_order ,& self%x0, self%y0 ,& self%R0 ,& self%O_point_psi ,& self%X_point_psi ,& self%axis_Btor ,& spline_size_bits/8.0 ,& self%spl_limits_R(1), self%spl_limits_R(2), self%spl_size_R ,& self%spl_limits_Z(1), self%spl_limits_Z(2), self%spl_size_Z ,& self%xmin, self%xmax ,& self%ymin, self%ymax ,& self%rhomax, self%rhomin ,& trim(flux_regions_print) ,& trim(self%divertor_polygon%polygon_string()) ,& trim(self%exclusion_polygon%polygon_string()) 101 FORMAT(& /& 1X, '<<Equilibrium>> ' /, & 1X, 'equilibrium type: ', 3X, A /, & 1X, 'params read from file: ', 3X, A /, & 1X, 'equilibrium case: ', 3X, A /, & 1X, 'path to equilibrium data: ', 3X, A /, & 1X, '<<NetCDF attributes>> ' /, & 1X, 'description: ', 3X, A /, & 1X, 'history: ', 3X, A /, & 1X, 'author: ', 3X, A /, & 1X, '<<From params.in>> ' /, & 1X, 'flip_Z (used for toroidal field reversal) = ', L10 /, & 1X, 'poloidal_field_factor (for helicity) = ', ES10.3E2 /, & 1X, 'spline_order = ', I10 /, & 1X, '<<Magnetic axis (normalised)>> ' /, & 1X, 'self%x0: ', ES10.3E2, '(m), self%y0: ', ES10.3E2, '(m)' /, & 1X, '<<Normalisation values>>' /, & 1X, 'self%R0: = ', ES10.3E2 /, & 1X, 'O_point_psi: = ', ES10.3E2 /, & 1X, 'X_point_psi: = ', ES10.3E2 /, & 1X, 'axis_Btor = ', ES10.3E2 /, & 1X, '<<Interpolators>>' /, & 1X, 'Size of pol. flx. interpolator (bytes): ', ES10.3E2 /, & 1X, 'R/self%R0 spl. dimension -- min:', ES10.3E2, ', max:', ES10.3E2, ', pts:', I10 /, & 1X, 'Z/self%R0 spl. dimension -- min:', ES10.3E2, ', max:', ES10.3E2, ', pts:', I10 /, & 1X, '<<Grid limits (normalised)>>' /, & 1X, 'R/self%R0 limits -- min:', ES10.3E2, ', max:', ES10.3E2 /, & 1X, 'Z/self%R0 limits -- min:', ES10.3E2, ', max:', ES10.3E2 /, & 1X, '<<Poloidal-flux global limits>>' /, & 1X, 'Poloidal flux limits (normalised, n.b. rhomin/max correspond to rhomin/min)'/, & 1X, 'rho limits -- max:', ES10.3E2, ', min:', ES10.3E2 /, & 1X, '<<Psi limit regions>>' /, & 1X, 'Regions: ', 3X, A /, & 1X, '<<First wall + divertor polygon>>' /, & 1X, A /, & 1X, '<<Exclusion polygon>>' /, & 1X, A /, & 1X, '<<End equilibrium>>' / & /) end subroutine debug end module numerical_equilibrium_m