module kisslinger_m !! Contains boundary data for 3D magnetic configurations provided in !! Kisslinger's format, see !! [Y Feng and W7-X-team 2022 Plasma Phys. Control. Fusion 64 125012, !! DOI 10.1088/1361-6587/ac9ed9], Figure A3 !! !! The boundary is prescribed as a set of polygons on phi cross sections use precision_m, only : FP, FP_EPS, FP_NAN use constants_m, only : TWO_PI, PI use screen_io_m, only : get_stdout use error_handling_m, only: handle_error, error_info_t use status_codes_m, only: PARALLAX_ERR_EQUILIBRIUM use list_operations_m, only : rank_list implicit none private integer, public, parameter :: TILT_SYMMETRY = 1 type, public :: kisslinger_t character(len=256), private :: header_description !! Description of header in File integer, private :: nphi !! Number of toroidal cross sections integer, private :: np !! Number of Polygon points per cross section integer, private :: mtor !! Toroidal symmetry real(FP), private, dimension(:), allocatable :: phi !! Toroidal angles of cross sections real(FP), private, dimension(:,:), allocatable :: xp !! x-coordinate of polygons real(FP), private, dimension(:,:), allocatable :: yp !! y-coordinate of polygons integer, private :: symm !! Type of symmetry that has been applied to data logical, private :: closed_toroidal !! = true if polygons are closed toroidally logical, private, dimension(:), allocatable :: closed_poloidal !! = true if polygon is closed poloidally contains procedure, public :: init => init_kisslinger procedure, public :: polygon_at_tor_xsection procedure, private :: read_from_file procedure, private :: apply_tilt_symmetry procedure, private :: sort_ascending_phi procedure, public :: display final :: destructor end type contains subroutine init_kisslinger(self, filename, axis_r, axis_z, symm) !! Initialises boundary data class(kisslinger_t), intent(inout) :: self !! Instance of type character(len=*), intent(in) :: filename !! Filename where data is read from real(FP), intent(in) :: axis_r !! R-location of magnetic axis used for normalisation real(FP), intent(in) :: axis_z !! Z-location of magnetic axis used for normalisation integer, intent(in) :: symm !! Type of symmetry to be applied to data !! - NO_SYMMETRY (default) !! - TILT_SYMMETRY (x -> x, phi -> -phi, y -> -y) integer :: iphi, ip real(FP) :: phi_first, phi_last, dst ! Read raw data from file call self%read_from_file(filename) ! Normalise data to magnetic axis do iphi = 1, self%nphi self%phi(iphi) = self%phi(iphi) * PI / 180.0_FP do ip = 1, self%np self%xp(ip, iphi) = self%xp(ip,iphi) / axis_r self%yp(ip, iphi) = (self%yp(ip,iphi) - axis_z) / axis_r enddo enddo ! Apply symmetry self%symm = symm select case(self%symm) case(TILT_SYMMETRY) call self%apply_tilt_symmetry() case default end select ! Sort data ascending in phi call self%sort_ascending_phi() ! Determine closedness phi_first = modulo(self%phi(1), TWO_PI/self%mtor) phi_last = modulo(self%phi(self%nphi), TWO_PI/self%mtor) if (abs(phi_first-phi_last) < FP_EPS) then self%closed_toroidal = .true. else self%closed_toroidal = .false. endif allocate(self%closed_poloidal(self%nphi)) do iphi = 1, self%nphi dst = sqrt((self%xp(1,iphi) - self%xp(self%np,iphi))**2 + & (self%yp(1,iphi) - self%yp(self%np,iphi))**2 ) if (abs(dst) < FP_EPS) then self%closed_poloidal(iphi) = .true. else self%closed_poloidal(iphi) = .false. endif enddo end subroutine subroutine read_from_file(self, filename) !! Reads raw data from file in Kisslinger format class(kisslinger_t), intent(inout) :: self !! Instance of type character(len=*), intent(in) :: filename !! Filename where data is read from integer :: funit, io_error integer :: iphi, ip real(FP) :: roffset, zoffset ! Read raw data open(newunit=funit, file=filename, status='old', action='read', & iostat=io_error) if (io_error /= 0) then call handle_error("Opening Kisslinger file "//filename//" failed!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__, & additional_info= error_info_t("Iostat:", [io_error])) end if read(funit,'(A)',iostat=io_error)self%header_description if (io_error /= 0) then call handle_error("Reading Kisslinger file "//filename//" failed!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__) end if read(funit,*,iostat=io_error)self%nphi, self%np, self%mtor, roffset, zoffset if (io_error /= 0) then call handle_error("Reading Kisslinger file "//filename//" failed!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__) end if allocate(self%phi(self%nphi)) allocate(self%xp(self%np, self%nphi)) allocate(self%yp(self%np, self%nphi)) do iphi = 1, self%nphi read(funit,*,iostat=io_error)self%phi(iphi) if (io_error /= 0) then call handle_error("Reading Kisslinger file "//filename//" failed!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__) end if do ip = 1, self%np read(funit,*,iostat=io_error)self%xp(ip, iphi), self%yp(ip, iphi) self%xp(ip, iphi) = self%xp(ip, iphi) + roffset self%yp(ip, iphi) = self%yp(ip, iphi) + zoffset enddo enddo close(funit) end subroutine subroutine apply_tilt_symmetry(self) !! Applies tilt symmetry to data !! x -> x, phi -> -phi, y -> -y class(kisslinger_t), intent(inout) :: self !! Instance of type integer :: iphi_zero, nphi_sym, iphi, ir, ip, ip_rev real(FP), allocatable, dimension(:) :: wrk_phi real(FP), allocatable, dimension(:,:) :: wrk_xp, wrk_yp ! Find entry phi == 0, to which symmetry operation shall not be applied iphi_zero = 0 do iphi = 1, self%nphi if (abs(self%phi(iphi)) < FP_EPS) then iphi_zero = iphi exit endif enddo if (iphi_zero > 0) then nphi_sym = 2 * self%nphi - 1 else nphi_sym = 2 * self%nphi endif ! Copy data to work array and reallocate data with new size allocate(wrk_phi, source=self%phi) allocate(wrk_xp, source=self%xp) allocate(wrk_yp, source=self%yp) deallocate(self%yp) deallocate(self%xp) deallocate(self%phi) allocate(self%phi(nphi_sym)) allocate(self%xp(self%np, nphi_sym)) allocate(self%yp(self%np, nphi_sym)) ! Copy first nphi data do iphi = 1, self%nphi self%phi(iphi) = wrk_phi(iphi) self%xp(:,iphi) = wrk_xp(:,iphi) self%yp(:,iphi) = wrk_yp(:,iphi) enddo ! Apply tilt symmetry, exclude phi=0 ir = 0 do iphi = 1, self%nphi if (iphi == iphi_zero) then cycle endif ir = ir + 1 self%phi(self%nphi+ir) = -wrk_phi(iphi) do ip = 1, self%np ! Run through wrk in reversed direction to ! conserve directionality of polygon ip_rev = self%np - (ip-1) self%xp(ip,self%nphi+ir) = wrk_xp(ip_rev,iphi) self%yp(ip,self%nphi+ir) = -wrk_yp(ip_rev,iphi) enddo enddo self%nphi = nphi_sym deallocate(wrk_yp) deallocate(wrk_xp) deallocate(wrk_phi) end subroutine subroutine sort_ascending_phi(self) !! Sorts data ascending in phi class(kisslinger_t), intent(inout) :: self !! Instance of type integer :: iphi integer, dimension(self%nphi) :: perm real(FP), allocatable, dimension(:) :: wrk_phi real(FP), allocatable, dimension(:,:) :: wrk_xp, wrk_yp allocate(wrk_phi, source=self%phi) allocate(wrk_xp, source=self%xp) allocate(wrk_yp, source=self%yp) call rank_list(wrk_phi, perm) do iphi = 1, self%nphi self%phi(iphi) = wrk_phi(perm(iphi)) self%xp(:,iphi) = wrk_xp(:,perm(iphi)) self%yp(:,iphi) = wrk_yp(:,perm(iphi)) enddo deallocate(wrk_yp) deallocate(wrk_xp) deallocate(wrk_phi) end subroutine subroutine polygon_at_tor_xsection(self, phi_xsec, ni, xi, yi) !! Returns polygon at toroidal cross section, !! via interpolation along phi class(kisslinger_t), intent(in) :: self !! Instance of type real(FP), intent(in) :: phi_xsec !! Toroidal angle where polygon is requested integer, intent(out) :: ni !! Number of polygon points at requested cross section real(FP), allocatable, dimension(:), intent(out) :: xi !! x-coordinates of polygon interpolated to cross section at phi_xsec real(FP), allocatable, dimension(:), intent(out) :: yi !! y-coordinates of polygon interpolated to cross section at phi_xsec real(FP) :: phi, phi_mb1, phi_mb2, dst_adj(2) integer :: iphi, ip, iadj(2) if (allocated(xi)) deallocate(xi) if (allocated(yi)) deallocate(yi) ! Find phi mapped to toroidal segment, where Kisslinger is defined select case(self%symm) case(TILT_SYMMETRY) ! Kisslinger defined for phi in [-pi/mtor, pi/mtor[ phi_mb1 = modulo(phi_xsec, TWO_PI/self%mtor) phi_mb2 = modulo(phi_xsec, PI/self%mtor) if (phi_mb1 == phi_mb2) then phi = phi_mb1 else phi = -PI/self%mtor + phi_mb2 endif case default ! Kisslinger defined for phi in [0, 2pi/mtor[ phi = modulo(phi_xsec, TWO_PI/self%mtor) end select ! Return empty polygon if phi not within toroidal range of Kisslinger if (.not.self%closed_toroidal) then if ( (phi > self%phi(self%nphi)) .or. & (phi < self%phi(1)) ) then ni = 0 allocate(xi(ni)) allocate(yi(ni)) return endif if (self%symm == TILT_SYMMETRY) then if (abs(phi) < minval(abs(self%phi))) then ni = 0 allocate(xi(ni)) allocate(yi(ni)) return endif endif endif ! Find phi indices adjacent to phi iadj = 0 do iphi = 1, self%nphi if ((phi - self%phi(iphi)) >= 0.0_FP) then iadj(1) = iphi iadj(2) = iphi + 1 endif enddo if (iadj(2) > self%nphi) then if (self%closed_toroidal) then ! Ensure periodicity (Last phi point is the same as first) iadj(1) = self%nphi-1 iadj(2) = self%nphi else iadj(2) = 0 endif endif ! Distance to adjacent cross sections dst_adj = FP_NAN if (iadj(1) > 0) then dst_adj(1) = phi - self%phi(iadj(1)) endif if (iadj(2) > 0) then dst_adj(2) = phi - self%phi(iadj(2)) endif ! Perform interpolation along phi ni = self%np allocate(xi(ni)) allocate(yi(ni)) if (iadj(1) == 0) then call handle_error("Interpolation of Kisslinger failed!", & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__) endif if (iadj(2) == 0) then ! Nearest neighbour interpolation xi(:) = self%xp(:,iadj(1)) yi(:) = self%yp(:,iadj(1)) return endif ! Linear interpolation do ip = 1, ni xi(ip) = (self%xp(ip,iadj(2)) * dst_adj(1) - & self%xp(ip,iadj(1)) * dst_adj(2)) / & (dst_adj(1) - dst_adj(2)) yi(ip) = (self%yp(ip,iadj(2)) * dst_adj(1) - & self%yp(ip,iadj(1)) * dst_adj(2)) / & (dst_adj(1) - dst_adj(2)) enddo end subroutine subroutine display(self) !! Writes basic information to stdout class(kisslinger_t), intent(in) :: self !! Instance of type write(get_stdout(), 99) self%header_description, self%nphi, self%np, & self%mtor, self%closed_toroidal, self%symm 99 FORMAT(1X,'Kisslinger boundary polygon segments:'/, & 1X,'description = ',A64 /, & 1X,'nphi = ',I4 /, & 1X,'np = ',I4 /, & 1X,'mtor = ',I4 /, & 1X,'closed_toroidal = ',L1 /, & 1X,'tilt_symmetry = ',I3) end subroutine subroutine destructor(self) !! Destructor type(kisslinger_t), intent(inout) :: self !! Instance of type if (allocated(self%phi)) deallocate(self%phi) if (allocated(self%xp)) deallocate(self%xp) if (allocated(self%yp)) deallocate(self%yp) if (allocated(self%closed_poloidal)) deallocate(self%closed_poloidal) end subroutine end module