module divertor_equilibrium_m use equilibrium_m, only: equilibrium_t use polygon_m, only: polygon2d_t, closed_polygon2d_t, limiting_polygon2d_t, & create_closed_polygon2d_t, create_limiting_polygon2d_t use precision_m, only: FP use descriptors_m use screen_io_m, only: get_stderr use error_handling_m, only: handle_error, error_info_t use status_codes_m, only: PARALLAX_ERR_EQUILIBRIUM implicit none integer :: max_polygon_N_pts = 100 !! How many points can be manually supplied as polygon points in the parameter file type, abstract, extends(equilibrium_t) :: divertor_equilibrium_t logical, public :: flip_Z = .false. !! Inverts the Z direction, which is equivalent to reversing the field direction !! Default false, use true to flip the equilibrium vertically (magnetic axis position unchanged). real(FP) :: poloidal_field_factor = 1.0_FP !! Poloidal field factor, which can be used to reverse the helicity type(closed_polygon2d_t) :: divertor_polygon !! Closed polygon which defines the divertor type(closed_polygon2d_t) :: exclusion_polygon !! Closed polygon which marks points as being off the grid ! logical, private :: set_divertor_polygon !! Use default first-wall + divertor polygon if false (default), or from parameter !! file if true ! logical, private :: set_exclusion_polygon !! Use default exclusion polygon if false (default), or from parameter !! file if true logical, public :: invert_divertor_polygon !! Flips the divertor polygon from being assumed concave to convex. logical, public :: invert_exclusion_polygon !! Flips the exclusion polygon from being assumed concave to convex. real(FP) :: axis_Btor !! Toroidal field at the the magnetic axis -- used to normalise magnetic field real(FP) :: O_point_psi !! Poloidal flux of the magnetic axis -- used to define the 0-point of rho real(FP) :: X_point_psi !! Poloidal flux of the (primary) seperatrix -- used to define the 1-point of rho character(len=:), allocatable, public :: district_definition !! Indicates, how districts are defined !! 'flux': Based on the limiting flux surface rho_max (default) !! 'wall': Based on divertor and exclusion polygons contains procedure, public :: is_axisymmetric procedure, public :: rho procedure, public :: jacobian procedure, public :: epol procedure, public :: erad procedure, public, nopass :: make_polygon procedure, public :: make_polygon_from_params procedure, public :: polygon_projection procedure, public :: in_vessel procedure, public :: mag_axis_loc procedure, public :: on_grid procedure, public :: district procedure, private :: district_flux_defined procedure, private :: district_wall_defined procedure(psi), public, deferred, pass(self) :: psi procedure(check_privflux_regions), public, deferred, pass(self) :: check_privflux_regions end type abstract interface real(FP) function psi(self, x, y, phi) !! Poloidal flux (in Weber) import :: FP, divertor_equilibrium_t class(divertor_equilibrium_t), intent(in) :: self real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) end function psi subroutine check_privflux_regions(self, x, y, phi, local_rhomin, local_rhomin_exists, local_rhomax, local_rhomax_exists) !! Returns the flux limits to be applied at (x, y, phi), taking into account region-specific limits import :: FP, divertor_equilibrium_t class(divertor_equilibrium_t), intent(in) :: self real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) real(FP), intent(inout) :: local_rhomax !! rhomax to be applied at (x, y, phi). Pass in global rhomax: the most restrictive limit is used real(FP), intent(inout) :: local_rhomin !! rhomin to be applied at (x, y, phi). Pass in global rhomin: the most restrictive limit is used logical, intent(out) :: local_rhomin_exists, local_rhomax_exists !! flags to mark whether a local limit is used, as opposed to the global limits passed in !! This denotes the "private flux" integer :: region_it end subroutine end interface contains logical function is_axisymmetric(self) class(divertor_equilibrium_t), intent(in) :: self is_axisymmetric = .true. end function real(FP) function rho(self, x, y, phi) !! Flux surface label normalised such that rho = 0 at the magnetic axis, and rho = 1 at the seperatrix class(divertor_equilibrium_t), intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP) :: psi_pt, rho_squared psi_pt = self%psi(x, y, phi) rho_squared = (psi_pt - self%O_point_psi)/(self%X_point_psi - self%O_point_psi) if (rho_squared > 0.0_FP) then rho = sqrt(rho_squared) else ! Probably comes from a wire (external to the device) -- set as zero ! N.b. this requires you to have been fairly careful in your selection of magnetic axis rho = 0.0_FP endif end function rho real(FP) function jacobian(self, x, y, phi) ! jacobian, i.e. J=R class(divertor_equilibrium_t), intent(in) :: self real(FP), intent(in) :: x, y, phi jacobian = x end function jacobian subroutine epol(self, x, y, phi, epolx, epoly) ! unit vector along poloidal (along flux surface) direction ! epolx: x-component of poloidal unit vector ! epoly: y-component of poloidal unit vector class(divertor_equilibrium_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 = self%bpol(x, y, phi) epolx = Bx / Bpol epoly = By / Bpol end subroutine epol subroutine erad(self, x, y, phi, eradx, erady) ! unit vector along radial (across flux surface) direction ! epolx: x-component of radial unit vector ! epoly: y-component of radial unit vector class(divertor_equilibrium_t), intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP), intent(out) :: eradx, erady real(FP) :: Bx, By, Bpol Bx = self%poloidal_field_factor * self%bx(x, y, phi) By = self%poloidal_field_factor * self%by(x, y, phi) Bpol = self%bpol(x, y, phi) eradx = By / Bpol erady = -Bx / Bpol end subroutine erad subroutine make_polygon(N_points, X_points, Y_points, polygon, & local_min, local_max, use_local_min, use_local_max) !! Initialises a polygon object "polygon" with points {X_points, Y_points} !! If local_min and local_max are given, a limiting_polygon2d_t is created, which stores these limits for later recall !! N.b. X_points and Y_points should be given normalised class(polygon2d_t), intent(out) :: polygon integer, intent(in) :: N_points real(FP), dimension(N_points), intent(in) :: X_points, Y_points real(FP), intent(in), optional :: local_min, local_max logical, intent(in), optional :: use_local_min, use_local_max select type (polygon) type is (closed_polygon2d_t) polygon = create_closed_polygon2d_t(N_points, X_points, Y_points) type is (limiting_polygon2d_t) if(present(local_min) .and. present(local_max) .and. present(use_local_max) .and. present(use_local_max)) then ! All must be defined, otherwise will get a segmentation fault polygon = create_limiting_polygon2d_t(N_points, X_points, Y_points, & local_min, local_max, use_local_min, use_local_max) else write(get_stderr(),*) & "error: limiting_polygon2d_t called without & &defined limits" ERROR STOP ERR_UNHANDLED endif class default write(get_stderr(),*) "error: polygon class not recognised" ERROR STOP ERR_UNHANDLED end select end subroutine make_polygon subroutine make_polygon_from_params(self, N_points, X_points, Y_points, polygon) !! Makes polygon from class(divertor_equilibrium_t), intent(in) :: self class(polygon2d_t), intent(out) :: polygon integer, intent(in) :: N_points real(FP), dimension(N_points), intent(in) :: X_points, Y_points if (self%flip_Z) then ! Reverse order to keep anti-clockwise, Z -> -Z call self%make_polygon(N_points = N_points, & X_points = X_points(N_points:1:-1), & Y_points = -1 * Y_points(N_points:1:-1), & polygon = polygon) else call self%make_polygon(N_points = N_points, & X_points = X_points(1:N_points), & Y_points = Y_points(1:N_points), & polygon = polygon) endif if (polygon%signed_area() <= 0.0_FP) then write(get_stderr(),*) & "Error: polygon must be defined anticlockwise. & &Signed area was ", polygon%signed_area() write(get_stderr(),*) "x", X_points, "y", Y_points ERROR STOP ERR_UNHANDLED endif end subroutine make_polygon_from_params subroutine polygon_projection(self, x, y, phi, projection_x, projection_y) !! Find the x and y components of a unit vector along the nearest polygon edge !! The polygon boundary dV is expected to have positive (counterclockwise) orientation !! Can check this with signed_area (in initialisation) class(divertor_equilibrium_t), intent(in) :: self real(FP), intent(in) :: x, y, phi real(FP), intent(out) :: projection_x, projection_y integer :: ipolyedge real(FP) :: dVx, dVy ! Find the closest polygon edge (by straight-line distance: can be improved -- i.e. using parallel connected points) ipolyedge = self%divertor_polygon%closest_polyedge(x,y) ! Find the normalised tangent vector of the nearest edge dVx = self%divertor_polygon%p2p(ipolyedge,1) / self%divertor_polygon%p2p_mag(ipolyedge) dVy = self%divertor_polygon%p2p(ipolyedge,2) / self%divertor_polygon%p2p_mag(ipolyedge) ! Normalise the components projection_x = dVx / abs(self%bx(x,y,phi) * dVy - self%by(x,y,phi) * dVx) projection_y = dVy / abs(self%bx(x,y,phi) * dVy - self%by(x,y,phi) * dVx) end subroutine polygon_projection logical function in_vessel(self, x, y, phi) !! Returns whether a point (x, y, phi) lies inside the divertor vessel !! Not especially performant, would be better to check if the point is outside of flux limits !! rather than checking the full district routine class(divertor_equilibrium_t), intent(in) :: self real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) integer :: district_pt district_pt = self%district(x, y, phi) if ((district_pt == DISTRICT_DOME) .or. (district_pt == DISTRICT_WALL)) then ! Outside vessel in_vessel = .false. else if (self%divertor_polygon%pt_inside(x, y) .neqv. self%invert_divertor_polygon) then ! Inside vessel in_vessel = .true. else ! In vessel 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(divertor_equilibrium_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 axis_x = self%x0 axis_y = self%y0 end subroutine mag_axis_loc logical function on_grid(self, x, y, phi) !! Returns whether a point (x, y, phi) lies on the computational grid class(divertor_equilibrium_t), intent(in) :: self real(FP), intent(in) :: x, y, phi !! 3D position (x and y normalized) if ( (x < self%xmin) .or. (x > self%xmax) .or. (y < self%ymin) .or. (y > self%ymax) ) then !Check if point is outside the grid extents on_grid = .false. elseif (.not.(self%exclusion_polygon%pt_inside(x, y)) .neqv. self%invert_exclusion_polygon) then ! Check if the point is inside the exclusion polygon specified on_grid = .false. else on_grid = .true. endif end function on_grid integer function district(self, x, y, phi) class(divertor_equilibrium_t), intent(in) :: self real(FP), intent(in) :: x real(FP), intent(in) :: y real(FP), intent(in) :: phi select case(self%district_definition) case('flux') district = self%district_flux_defined(x, y, phi) case('wall') district = self%district_wall_defined(x, y, phi) case default call handle_error('district_definition not valid', & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__, & additional_info=& error_info_t(self%district_definition)) end select end function integer function district_flux_defined(self, x, y, phi) !! Returns a descriptor for a given point, where outer boundary !! is based on limiting flux surface(s) ! The algorithm only applies local limits towards the separatrix. ! This is because in some cases (i.e. snowflake) it is ! possible to have both rho < 1.0 and rho > 1.0 private flux regions. class(divertor_equilibrium_t), intent(in) :: self !! Instance of type real(FP), intent(in) :: x !! x-coordinate real(FP), intent(in) :: y !! y-coordinate real(FP), intent(in) :: phi !! Toroidal angle real(FP) :: rho_pt, local_rhomax, local_rhomin logical :: local_rhomin_exists, local_rhomax_exists ! Check if point is outside the grid extents if (.not.(self%on_grid(x, y, phi))) then district_flux_defined = DISTRICT_OUT return endif ! If not a DISTRICT_OUT point, find the district of the point ! from the rho-label rho_pt = self%rho(x, y, phi) ! Default to using the global minimum or maximum as flux limits local_rhomin = self%rhomin local_rhomax = self%rhomax ! Check if there is a local rho limit at (x, y, phi) call self%check_privflux_regions(x, y, phi, & local_rhomin, local_rhomin_exists, & local_rhomax, local_rhomax_exists) ! If rhomax is less than 1 then only want to simulate closed field lines ! - don't consider special flux-limited regions if (self%rhomax <= 1.0_FP) then if ( (rho_pt > self%rhomax) .or. & local_rhomin_exists .or. & local_rhomax_exists ) then ! This deals with an issue where there may be regions ! of very low psi in the private flux district_flux_defined = DISTRICT_WALL return else if (rho_pt < self%rhomin) then district_flux_defined = DISTRICT_CORE return else district_flux_defined = DISTRICT_CLOSED return endif endif if (rho_pt < 1.0_FP) then if (local_rhomin_exists) then ! Check whether rho(x, y, phi) is within the local rho limits if ( (rho_pt > local_rhomin) .and. (rho_pt < self%rhomax) ) then district_flux_defined = DISTRICT_PRIVFLUX return else district_flux_defined = DISTRICT_DOME return endif else if (rho_pt < self%rhomin) then district_flux_defined = DISTRICT_CORE return else if (rho_pt > self%rhomax) then district_flux_defined = DISTRICT_WALL return else district_flux_defined = DISTRICT_CLOSED return endif end if else ! Open field-line region if (local_rhomax_exists) then ! Check whether rho(x, y, phi) is within the local rho limits if ((rho_pt > self%rhomin) .and. (rho_pt < local_rhomax)) then district_flux_defined = DISTRICT_PRIVFLUX return else district_flux_defined = DISTRICT_DOME return endif else if ((rho_pt < self%rhomin) .or. (rho_pt > self%rhomax)) then district_flux_defined = DISTRICT_WALL return else district_flux_defined = DISTRICT_SOL return end if endif endif ! Should never get to this point call handle_error('district_flux_defined not assignable', & PARALLAX_ERR_EQUILIBRIUM, __LINE__, __FILE__) end function integer function district_wall_defined(self, x, y, phi) !! Returns a descriptor for a given point, where outer boundary !! is based on divertor and exclusion polygons class(divertor_equilibrium_t), intent(in) :: self !! Instance of type real(FP), intent(in) :: x !! x-coordinate real(FP), intent(in) :: y !! y-coordinate real(FP), intent(in) :: phi !! Toroidal angle logical :: in_div, in_exc real(FP) :: rho, local_rhomin, local_rhomax logical :: local_rhomin_exists, local_rhomax_exists in_exc = self%exclusion_polygon%pt_inside(x, y) if (.not. in_exc) then district_wall_defined = DISTRICT_OUT return endif in_div = self%divertor_polygon%pt_inside(x, y) if (.not. in_div) then district_wall_defined = DISTRICT_WALL return endif ! For check, if point inside private flux region local_rhomin = self%rhomin local_rhomax = self%rhomax call self%check_privflux_regions(x, y, phi, & local_rhomin, local_rhomin_exists, & local_rhomax, local_rhomax_exists) rho = self%rho(x, y, phi) if (rho > 1.0_FP) then if (local_rhomax_exists) then district_wall_defined = DISTRICT_PRIVFLUX return else district_wall_defined = DISTRICT_SOL return endif endif if (local_rhomin_exists) then district_wall_defined = DISTRICT_PRIVFLUX else if (rho < self%rhomin) then district_wall_defined = DISTRICT_CORE else district_wall_defined = DISTRICT_CLOSED endif endif end function end module divertor_equilibrium_m