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 moose_configparser use moose_utils, only: isdir, join, startswith use moose_error use FLARE_control, only : FLARE_init, FLARE_set_stdout => set_stdout use FLARE_model, only : FLARE_TYPE_EQUI3D_MGRID => TYPE_EQUI3D_MGRID, & FLARE_TYPE_EQUI3D_COILSET => TYPE_EQUI3D_COILSET, & FLARE_TYPE_EQUI3D_HINT => TYPE_EQUI3D_HINT, & FLARE_TYPE_EQUI2D_GEQDSK => TYPE_EQUI2D_GEQDSK, & FLARE_alloc_equi3d => alloc_equi3d, & FLARE_load_equi3d_mgrid => load_equi3d_mgrid, & FLARE_load_equi3d_coilset => load_equi3d_coilset, & FLARE_load_equi3d_hint => load_equi3d_hint, & FLARE_alloc_equi2d => alloc_equi2d, & FLARE_load_geqdsk => load_geqdsk, & FLARE_setup_model => setup_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_equi3d_path_bfield, Rmin, Rmax, Zmin, Zmax, & flare_equi3d_path_hint, group, bmax, & 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 call FLARE_set_stdout(get_stdout()) ! 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_MGRID) 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_EQUI3D_COILSET) call load_bfield_coilsets(flare_equi3d_path_bfield) self%axisymmetric = .false. case(FLARE_TYPE_EQUI3D_HINT) ! group is 0 for vacuum and 1 for plasma, ! in the python interface it is called plasma instead. call FLARE_alloc_equi3d(1, flare_equi3d_path_axis) call FLARE_load_equi3d_hint(1, flare_equi3d_path_hint, & group, bmax) 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 if (equi_flare_type == FLARE_TYPE_EQUI3D_COILSET) then FLARE_bfield%lb(1) = Rmin FLARE_bfield%ub(1) = Rmax FLARE_bfield%lb(2) = Zmin FLARE_bfield%ub(2) = Zmax endif ! 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() if (rho_flare_type == 'kisslinger') then do is = 1, self%nsurfs write(get_stdout(),*)'flux surface rho = ', self%rho_surfs(is) call self%kiss_rho(is)%display() enddo endif 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 write(get_stdout(),100) & self%axis_r_flare, & self%axis_z_flare, & self%axis_btor_flare, & self%rhomin, & self%rhomax, & self%xmin, & self%xmax, & self%ymin, & self%ymax 100 FORMAT( / & 1X,'FLARE Equilibrium:'/, & 1X,'axis_r_flare = ',ES14.6E3 /, & 1X,'axis_z_flare = ',ES14.6E3 /, & 1X,'axis_btor_flare = ',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 & ) 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 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 subroutine load_bfield_coilsets(bfield_path) !! Function to initialize bfield coilset equilibrium with multiple coils !! Parses the contents of the .bfield file !! Loads all coilset files listed in the file character(len=*), intent(in) :: bfield_path !! Path to .bfield config file type(configparser) :: config type(cp_section), pointer :: section character(len=:), allocatable :: config_path, filename, dtype character(len=:), allocatable :: units, coilfile, sval integer :: nequi3d, i, ios logical :: ex real(FP) :: amplitude logical :: have_amp if (isdir(bfield_path)) then config_path = bfield_path filename = join(config_path, ".bfield") else config_path = dirname(bfield_path) filename = bfield_path end if inquire(file=filename, exist=ex) if (.not.ex) call ERROR(& "configuration file "//trim(filename)//" does not exist") config = configparser() call config%read(filename) ! Count number of equi3d_coilset sections and allocate once nequi3d = 0 section => config%first_section() do if (.not.associated(section)) exit if (get_dtype(section%key) == "equi3d_coilset") & nequi3d = nequi3d + 1 section => config%next_section() end do if (nequi3d == 0) call ERROR("no equi3d_coilset sections found") call FLARE_alloc_equi3d(nequi3d, flare_equi3d_path_axis) i = 1 section => config%first_section() do if (.not.associated(section)) exit if (get_dtype(section%key) == "equi3d_coilset") then coilfile = trim(section%get("filename", & raw=.false., fallback="")) if (len_trim(coilfile) == 0) & & call ERROR( & "section '"//trim(section%key)//"' missing 'filename'") if (.not.is_absolute_path(coilfile)) coilfile = join(config_path, & coilfile) units = trim(section%get("units", raw=.false., fallback="m")) sval = trim(section%get("amplitude", raw=.false., fallback="")) have_amp = .false. if (len_trim(sval) > 0) then read(sval, *, iostat=ios) amplitude if (ios == 0) have_amp = .true. end if if (have_amp) then call FLARE_load_equi3d_coilset(i, coilfile, units=units, & amplitude=amplitude) else call FLARE_load_equi3d_coilset(i, coilfile, units=units, & amplitude=1.0_FP) end if i = i + 1 end if section => config%next_section() end do end subroutine load_bfield_coilsets pure logical function is_absolute_path(p) result(abs) !! Function to return true if p is an absolute path character(len=*), intent(in) :: p !! Filepath integer :: n n = len_trim(p) abs = .false. if (n >= 1) then if (p(1:1) == '/' .or. p(1:1) == '\') abs = .true. if (n >= 2) then if (p(2:2) == ':' .and. ((p(1:1) >= 'A' .and. p(1:1) <= 'Z') & .or. (p(1:1) >= 'a' .and. p(1:1) <= 'z'))) abs = .true. end if end if end function is_absolute_path pure function dirname(path) result(dir) !! Resolves the directory name for a given path character(len=*), intent(in) :: path !! Filepath or file name character(len=:), allocatable :: dir !! Path to directory integer :: i, n n = len_trim(path) do i = n, 1, -1 if (path(i:i) == '/' .or. path(i:i) == '\') then dir = path(1:i-1) return end if end do dir = "." end function dirname pure function get_dtype(key) result(dtype) !! Gets data value from the given key character(len=*), intent(in) :: key !! equi3d_coilset section key character(len=:), allocatable :: dtype !! Resulting data e.g. equi3d_coilset integer :: i, n n = len_trim(key) do i = n, 1, -1 if (key(i:i) == ":") then dtype = adjustl(key(i+1:n)) return end if end do dtype = adjustl(key(1:n)) end function get_dtype end module