module flare_equilibrium_m !! Equilibrium based on FLARE: !! (field line analysis and reconstruction for 3D boundary plasma modeling) !! see H. Frerichs 2024 Nucl. Fusion 64 106034 DOI 10.1088/1741-4326/ad7303 !! !! The source code of FLARE needs to be compiled separately !! https://gitlab.com/hfrerichs/flare !! (stable with #4351aff722f9d2b7495211f09f774bc406064ac1) !! !! FLARE relies on MOOSE, which needs to be compile seperately, too. !! https://gitlab.com/hfrerichs/moose.git !! (stable with #4c4a046b615ba19f68e1c1003c975941db7c6cf7) !! !! All routines or variables that are provided externally from FLARE are !! imported with the prefix "FLARE_" use, intrinsic :: iso_fortran_env, only: IOSTAT_END, IOSTAT_EOR use error_handling_m, only : handle_error, error_info_t use status_codes_m, only : PARALLAX_WRN_NOT_IMPLEMENTED, & PARALLAX_ERR_PARAMETERS, & PARALLAX_WRN_GENERAL, & PARALLAX_ERR_EQUILIBRIUM use precision_m, only : FP use screen_io_m, only : get_stdout use comm_handling_m, only : is_master use descriptors_m, only : DISTRICT_CORE, DISTRICT_CLOSED, DISTRICT_WALL, & DISTRICT_SOL use equilibrium_m, only : equilibrium_t use FLARE_control, only : FLARE_init use FLARE_model, only : FLARE_TYPE_EQUI3D => TYPE_EQUI3D, & FLARE_TYPE_EQUI2D_GEQDSK => TYPE_EQUI2D_GEQDSK, & FLARE_alloc_equi3d => alloc_equi3d, & FLARE_load_equi3d_mgrid => load_equi3d_mgrid, & FLARE_alloc_equi2d => alloc_equi2d, & FLARE_load_geqdsk => load_geqdsk, & FLARE_setup_model => setup_model, & FLARE_broadcast_model => broadcast_model, & FLARE_assert_model => assert_model, & FLARE_bfield => bfield use kisslinger_m, only : kisslinger_t use params_equi_flare_m, only : & read_equi_params_flare_type, write_equi_params_flare_type, & read_flare_equi_params, write_flare_equi_params, & read_flare_bnd_params, write_flare_bnd_params, & read_flare_rho_params, write_flare_rho_params, & equi_flare_type, boundary_flare_type, rho_flare_type, & amplitudes, spline_order, bspline_data, scale_bt, scale_ip, & rhomin, rhomax, & flare_equi3d_path_axis, flare_equi3d_path_mgrid, flare_geqdsk_path, & flare_rho_dirpath, flare_rho_prefix, rho_symms, & flare_units_system, & nbnd_seg, bnd_paths, bnd_symms use euclidean_geo_m, only : intersection_lines use polygon_m, only : closed_polygon2d_t, create_closed_polygon2d_t implicit none private type, public, extends(equilibrium_t) :: flare_t !! Equilibrium type based on FLARE magnetic reconstruction real(FP), private :: axis_r_flare !! R-coordinate of the magnetic axis in FLARE normalisation real(FP), private :: axis_z_flare !! Z-coordinate of the magnetic axis in FLARE normalisation real(FP), private :: axis_btor_flare !! Toroidal magnetic field strength on axis in FLARE normalisation logical, private :: axisymmetric type(kisslinger_t), private, allocatable, dimension(:) :: kiss_bnd !! Description of boundary in Kisslinger's format integer, private :: nsurfs !! Number of flux surface polygons real(FP), private, allocatable, dimension(:) :: rho_surfs !! rho values on flux surfaces type(kisslinger_t), private, allocatable, dimension(:) :: kiss_rho !! Description of flux surfaces in Kisslinger's format real(FP), private :: lenconvfac_flare_kiss !! Conversion factor for length scales from flare units (SI or CGS) !! to Kisslinger units (cm) 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, private :: read_flux_surfaces end type contains subroutine init(self, filename) !! Initialises a FLARE equilibrium class(flare_t), intent(inout) :: self !! Instance of class character(len=*), intent(in), optional :: filename integer :: ibnd real(FP), dimension(2) :: axis_flare real(FP), dimension(3) :: baxis_flare ! Read parameters if (present(filename)) then call read_equi_params_flare_type(filename) call read_flare_equi_params(filename) endif call FLARE_init(standalone=.false., greeting=.true.) ! Initialise FLARE equilibrium select case(equi_flare_type) case(FLARE_TYPE_EQUI3D) call FLARE_alloc_equi3d(1, flare_equi3d_path_axis) call FLARE_load_equi3d_mgrid(1, flare_equi3d_path_mgrid, & amplitudes, bspline_data, spline_order) self%axisymmetric = .false. case(FLARE_TYPE_EQUI2D_GEQDSK) call FLARE_alloc_equi2d() call FLARE_load_geqdsk(flare_geqdsk_path, scale_bt, scale_ip, & spline_order) self%axisymmetric = .true. case default call handle_error(& 'equi_type_flare='//equi_flare_type//' not valid', & PARALLAX_ERR_PARAMETERS, __LINE__, __FILE__) end select call FLARE_setup_model() call FLARE_assert_model() ! Get magnetic axis self%phi0 = 0.0_FP axis_flare = FLARE_bfield%equi%magnetic_axis(self%phi0) self%axis_r_flare = axis_flare(1) self%axis_z_flare = axis_flare(2) baxis_flare = FLARE_bfield%eval(& [self%axis_r_flare, self%axis_z_flare, self%phi0]) self%axis_btor_flare = baxis_flare(3) self%x0 = 1.0_FP self%y0 = 0.0_FP ! Set bounding box self%xmin = FLARE_bfield%lb(1) / self%axis_r_flare self%ymin = (FLARE_bfield%lb(2) - self%axis_z_flare) / self%axis_r_flare self%xmax = FLARE_bfield%ub(1) / self%axis_r_flare self%ymax = (FLARE_bfield%ub(2) - self%axis_z_flare) / self%axis_r_flare ! Unit conversion factor for length scales between ! flare and kisslinger (cm) if (flare_units_system == 'CGS') then self%lenconvfac_flare_kiss = 1.0_FP elseif (flare_units_system == 'SI') then self%lenconvfac_flare_kiss = 100.0_FP else call handle_error(& 'flare_units_system='//flare_units_system//' not valid', & PARALLAX_ERR_PARAMETERS, __LINE__, __FILE__) endif ! Initialise description of boundaries if (present(filename)) then call read_flare_bnd_params(filename) endif select case(boundary_flare_type) case('none') ! Do nothing case('kisslinger') allocate(self%kiss_bnd(nbnd_seg)) do ibnd = 1, nbnd_seg call self%kiss_bnd(ibnd)%init(bnd_paths(ibnd), & self%axis_r_flare*self%lenconvfac_flare_kiss, & self%axis_z_flare*self%lenconvfac_flare_kiss, & bnd_symms(ibnd)) enddo case default call handle_error(& 'boundary_flare_type='//boundary_flare_type//' not valid', & PARALLAX_ERR_PARAMETERS, __LINE__, __FILE__) end select ! Initialise description of flux surfaces (rho) if (present(filename)) then call read_flare_rho_params(filename) endif select case(rho_flare_type) case('none') case('kisslinger') call self%read_flux_surfaces() case default call handle_error(& 'rho_flare_type='//rho_flare_type//' not valid', & PARALLAX_ERR_PARAMETERS, __LINE__, __FILE__) end select self%rhomin = rhomin self%rhomax = rhomax self%initialized = .true. end subroutine subroutine display(self) !! Prints to console information about the equilibrium class(flare_t), intent(in) :: self !! Instance of class integer :: ibnd, is if (is_master()) then call write_equi_params_flare_type() call write_flare_equi_params() call write_flare_bnd_params() do ibnd = 1, nbnd_seg call self%kiss_bnd(ibnd)%display() enddo call write_flare_rho_params() do is = 1, self%nsurfs write(get_stdout(),*)'flux surface rho = ', self%rho_surfs(is) call self%kiss_rho(is)%display() enddo endif write(*,*)'axis_r_flare', self%axis_r_flare write(*,*)'axis_z_flare', self%axis_z_flare write(*,*)'axis_btor_flare', self%axis_btor_flare write(*,*)'xmin', self%xmin write(*,*)'ymin', self%ymin write(*,*)'xmax', self%xmax write(*,*)'ymax', self%ymax end subroutine subroutine debug(self) !! Prints to console extended debug information about the equilibrium class(flare_t), intent(in) :: self !! Instance of class call handle_error("FLARE: debug not yet implemented", & PARALLAX_WRN_NOT_IMPLEMENTED, __LINE__, __FILE__) end subroutine logical function is_axisymmetric(self) !! Returns if the equilibrium is axisymmetric or not class(flare_t), intent(in) :: self !! Instance of class is_axisymmetric = self%axisymmetric end function real(FP) function rho(self, x, y, phi) !! Returns the radial coordinate, i.e. distance to torus axis class(flare_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) integer :: is, npoly, ip, ipvt real(FP), allocatable, dimension(:) :: xpoly, ypoly type(closed_polygon2d_t), dimension(self%nsurfs) :: polys real(FP) :: dist_inner, dist_outer, interp_frac, axis_x, axis_y if (rho_flare_type == 'kisslinger') then ! Create closed_polygon2d_t polys from kisslinger's polygon ! Functionality for closed_polygon2d_t can then be used below do is = 1, self%nsurfs call self%kiss_rho(is)%polygon_at_tor_xsection( & phi, npoly, xpoly, ypoly) polys(is) = create_closed_polygon2d_t(npoly, xpoly, ypoly) enddo ! Capture case, where point inside first flux surface: ! Interpolate between magnetic axis and first flux surface. if (polys(1)%pt_inside(x, y)) then call self%mag_axis_loc(phi, axis_x, axis_y) dist_inner = sqrt((x - axis_x)**2 + (y - axis_y)**2) dist_outer = polys(1)%dist_to_polyedge(x, y) interp_frac = dist_inner / (dist_inner + dist_outer) rho = interp_frac * self%rho_surfs(1) return endif ! Get indices of two adjacent flux surfaces do is = self%nsurfs, 1, -1 if(.not. polys(is)%pt_inside(x, y)) then if(is == self%nsurfs) then ipvt = is - 1 else ipvt = is endif exit endif enddo ! Linear Inter/extrapolation of flux surface label to point dist_inner = polys(ipvt)%dist_to_polyedge(x, y) dist_outer = polys(ipvt+1)%dist_to_polyedge(x, y) interp_frac = dist_inner / (dist_inner + dist_outer) rho = self%rho_surfs(ipvt) & + (self%rho_surfs(ipvt+1) - self%rho_surfs(ipvt)) * interp_frac return endif ! If no description for rho is provided, just return rho = 0.0 rho = 0.0_FP end function real(FP) function bx(self, x, y, phi) !! Returns the x-component of the magnetic field class(flare_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) real(FP), dimension(3) :: coord_flare, b_flare !$omp critical ! Conversion to FLARE coordinates coord_flare(1) = x * self%axis_r_flare coord_flare(2) = y * self%axis_r_flare + self%axis_z_flare coord_flare(3) = phi b_flare = FLARE_bfield%eval(coord_flare) bx = b_flare(1) / self%axis_btor_flare !$omp end critical end function real(FP) function by(self, x, y, phi) !! Returns the y-component of the magnetic field class(flare_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) real(FP), dimension(3) :: coord_flare, b_flare !$omp critical ! Conversion to FLARE coordinates coord_flare(1) = x * self%axis_r_flare coord_flare(2) = y * self%axis_r_flare + self%axis_z_flare coord_flare(3) = phi b_flare = FLARE_bfield%eval(coord_flare) by = b_flare(2) / self%axis_btor_flare !$omp end critical end function real(FP) function btor(self, x, y, phi) !! Returns the toroidal component btor of the magnetic field class(flare_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) real(FP), dimension(3) :: coord_flare, b_flare !$omp critical ! Conversion to FLARE coordinates coord_flare(1) = x * self%axis_r_flare coord_flare(2) = y * self%axis_r_flare + self%axis_z_flare coord_flare(3) = phi b_flare = FLARE_bfield%eval(coord_flare) btor = b_flare(3) / self%axis_btor_flare !$omp end critical end function real(FP) function jacobian(self, x, y, phi) !! Returns the Jacobian of the geometry class(flare_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(flare_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 call handle_error("FLARE: epol not yet implemented", & PARALLAX_WRN_NOT_IMPLEMENTED, __LINE__, __FILE__) 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(flare_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 call handle_error("FLARE: erad not yet implemented", & PARALLAX_WRN_NOT_IMPLEMENTED, __LINE__, __FILE__) end subroutine integer function district(self, x, y, phi) !! Returns in which district point (x, y) is (see module descriptors_m) class(flare_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) real(FP) :: rho_p rho_p = self%rho(x, y, phi) if (rho_p < self%rhomin) then district = DISTRICT_CORE return endif if (rho_p > self%rhomax) then district = DISTRICT_WALL return endif if (rho_p > 1.0_FP) then district = DISTRICT_SOL return endif district = DISTRICT_CLOSED end function logical function in_vessel(self, x, y, phi) !! Returns whether a given point is within the vessel or not class(flare_t), intent(in) :: self !! Instance of class real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) real(FP) :: axis_x, axis_y 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. call self%mag_axis_loc(phi, axis_x, axis_y) 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) = axis_x pa(2) = axis_y 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(in_vessel .eqv. .false.) then exit endif else if (info == 1 .or. info == -2) then call handle_error(& "FLARE: the 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(flare_t), intent(in) :: self !! Instance of class real(kind=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), dimension(2) :: ax_flare !$omp critical ax_flare = FLARE_bfield%equi%magnetic_axis(phi) axis_x = ax_flare(1) / self%axis_r_flare axis_y = (ax_flare(2) - self%axis_z_flare) / self%axis_r_flare !$omp end critical end subroutine subroutine read_flux_surfaces(self) !! Reads flux surfaces for computing rho class(flare_t), intent(inout) :: self !! Instance of class character(len=:), allocatable :: flux_label_file, flux_surface_file character(len=256) :: io_errmsg integer :: funit, io_error, is character(len=3) :: is_c real(FP) :: mean_flux_area ! Read number of flux surfaces and (unnormalised) flux surface labels flux_label_file = trim(adjustl(flare_rho_dirpath)) // '/' //& trim(adjustl(flare_rho_prefix)) // & '_flux_labels.txt' open(newunit=funit, file=flux_label_file, & status='old', action='read', iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, PARALLAX_ERR_EQUILIBRIUM, & __LINE__, __FILE__) endif read(funit,*, iostat=io_error, iomsg=io_errmsg)self%nsurfs if (io_error /= 0) then call handle_error(io_errmsg, PARALLAX_ERR_EQUILIBRIUM, & __LINE__, __FILE__) endif allocate(self%rho_surfs(self%nsurfs)) do is = 1, self%nsurfs read(funit,*,iostat=io_error, iomsg=io_errmsg) & mean_flux_area, self%rho_surfs(is) if (io_error /= 0) then call handle_error(io_errmsg, PARALLAX_ERR_EQUILIBRIUM, & __LINE__, __FILE__) endif enddo close(funit) ! Read flux surfaces allocate(self%kiss_rho(self%nsurfs)) do is = 1, self%nsurfs write(is_c,'(I3.3)')is flux_surface_file = trim(adjustl(flare_rho_dirpath)) // '/' //& trim(adjustl(flare_rho_prefix)) // & '_isurf_' // is_c // '.txt' call self%kiss_rho(is)%init(flux_surface_file, & self%axis_r_flare*self%lenconvfac_flare_kiss, & self%axis_z_flare*self%lenconvfac_flare_kiss, & rho_symms) enddo end subroutine end module