module connection_length_m !! Connection length and shortest distance to the target routines use precision_m, only : FP, FP_LARGEST use equilibrium_m, only: equilibrium_t use fieldline_tracer_m, only : trace use descriptors_m, only : DISTRICT_OUT use error_handling_m, only: handle_error use status_codes_m, only: PARALLAX_WRN_GENERAL use screen_io_m, only: get_stdout implicit none public :: connection_length public :: shortest_dist_to_bound public :: distance_to_boundary contains subroutine distance_to_boundary(equi, x, y, phi, dphi_max, arclength_fwd, & arclength_bwd, target_reached_fwd, target_reached_bwd, maxstepsize, & dphi_fwd, dphi_bwd) !! Computes the distance to the target in each direction !! If the (x,y,phi) point is located outside of the vessel, the !! distances are defined negative class(equilibrium_t), target, intent(inout) :: equi !! Equilibrium defining the field line real(FP), intent(in) :: x !! Spatial coordinate x at which the field line is located real(FP), intent(in) :: y !! Spatial coordinate y at which the field line is located real(FP), intent(in) :: phi !! Spatial coordinate phi at which the field line is located real(FP), intent(in) :: dphi_max !! Maximum toroidal angle to be traced real(FP), intent(out) :: arclength_fwd !! For points that are inside vessel (based on equi%in_vessel): !! Arclength in the direction of B until the target is reached !! For points that are not inside vessel: !! Negative arclength in the direction of -B until target is reached real(FP), intent(out) :: arclength_bwd !! For points that are inside vessel (based on equi%in_vessel): !! Arclength in the direction of -B until the target is reached !! For points that are not inside vessel: !! Negative arclength in the direction of B until target is reached logical, intent(inout) :: target_reached_fwd !! True, if the target is reached during the forward (as defined above) !! tracing logical, intent(inout) :: target_reached_bwd !! True, if the target is reached during the backward (as defined above) !! tracing real(FP), intent(in) :: maxstepsize !! Maximum stepsize for the DOP853 integrator real(FP), intent(out), optional :: dphi_fwd !! Toroidal angle traced in direction of B until target is reached !! It follows the same sign convention as arclength_fwd real(FP), intent(out), optional :: dphi_bwd !! Toroidal angle traced in direction of -B until target is reached !! It follows the same sign convention as arclength_bwd real(FP) :: dphi_out ! Toroidal angle to end of domain (DISTRICT_OUT) real(FP), parameter :: trace_large = sqrt(FP_LARGEST) ! Large value returned if "OUT" region is reached real(FP) :: x_target_fwd, y_target_fwd, phi_target_fwd real(FP) :: x_target_bwd, y_target_bwd, phi_target_bwd ! Compute arclengths up to the plasma-wall boundaries ----------------- if (.not. condition_in_target(x, y, phi)) then ! Trace until forward / backward target (or maximum distance) call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = dphi_max, & equi = equi, & x_end = x_target_fwd, & y_end = y_target_fwd, & phi_end = phi_target_fwd, & arclen = arclength_fwd, & maxstep = maxstepsize, & stop_int = .true., & tracing_terminated = target_reached_fwd, & condition = condition_in_target) if (present(dphi_fwd)) then dphi_fwd = abs(phi_target_fwd - phi) endif call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = -dphi_max, & equi = equi, & x_end = x_target_bwd, & y_end = y_target_bwd, & phi_end = phi_target_bwd, & arclen = arclength_bwd, & maxstep = maxstepsize, & stop_int = .true., & tracing_terminated = target_reached_bwd, & condition = condition_in_target) if (present(dphi_bwd)) then dphi_bwd = abs(phi_target_bwd - phi) endif if((.not. target_reached_fwd) .eqv. target_reached_bwd) then call handle_error(& "Wall reached only in one direction in the SOL", & PARALLAX_WRN_GENERAL, __LINE__, __FILE__) endif else ! For points in the wall, firstly one has to trace until dphi_max, ! but with the condition that the state vector is still in the ! domain. The endpoint of that trace is dphi_out, which can be ! used for the new maximum angle to be traced in the direction ! in question. ! For the forward arclength one traces backward, and for the backward ! target one traces forward, since these distances are defined ! negative in the wall region. ! Trace toroidal distance until field line leaves domain call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = -dphi_max, & equi = equi, & phi_end = dphi_out, & maxstep = maxstepsize, & stop_int = .true., & condition = condition_not_in_domain) dphi_out = dphi_out - phi ! Trace until in reverse direction until target is hit (up to at ! most phi_out) call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = dphi_out, & equi = equi, & phi_end = phi_target_fwd, & arclen = arclength_fwd, & maxstep = maxstepsize, & stop_int = .true., & tracing_terminated = target_reached_fwd, & condition = condition_not_in_target) ! Set large values if out region is reached if (.not.target_reached_fwd) then arclength_fwd = trace_large phi_target_fwd = trace_large endif if (present(dphi_fwd)) then dphi_fwd = -abs(phi_target_fwd - phi) endif ! Same for the other direction call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = dphi_max, & equi = equi, & phi_end = dphi_out, & maxstep = maxstepsize, & stop_int = .true., & condition = condition_not_in_domain) dphi_out = dphi_out - phi call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = dphi_out, & equi = equi, & phi_end = phi_target_bwd, & arclen = arclength_bwd, & maxstep = maxstepsize, & stop_int = .true., & tracing_terminated = target_reached_bwd, & condition = condition_not_in_target) ! Set large values if out region is reached if (.not.target_reached_bwd) then arclength_bwd = trace_large phi_target_bwd = -trace_large endif if (present(dphi_bwd)) then dphi_bwd = -abs(phi_target_bwd - phi) endif arclength_fwd = - arclength_fwd arclength_bwd = - arclength_bwd endif contains logical function condition_in_target(x, y, phi) real(FP), intent(in) :: x real(FP), intent(in) :: y real(FP), intent(in) :: phi condition_in_target = (.not.equi%in_vessel(x, y, phi)) end function logical function condition_not_in_target(x, y, phi) real(FP), intent(in) :: x real(FP), intent(in) :: y real(FP), intent(in) :: phi condition_not_in_target = equi%in_vessel(x, y, phi) end function logical function condition_not_in_domain(x, y, phi) real(FP), intent(in) :: x real(FP), intent(in) :: y real(FP), intent(in) :: phi integer :: district district = equi%district(x, y, phi) if (district == DISTRICT_OUT) then condition_not_in_domain = .TRUE. else condition_not_in_domain = .FALSE. endif end function end subroutine real(FP) function connection_length(equi, x, y, phi, dphi_max, maxstep) !! It computes the connection length for the field line defined !! by x, y, phi (the distance from boundary to boundary). !! If (x, y, phi) is located inside the wall, the output is negative. !! If (x, y, phi) is located in the closed flux surface region, the !! output is the largest allowable number for PARALLAX precision. !! If the boundary is reached along the field line, but only in one !! direction, then the routine will return with the arclength to the !! found boundary point plus the arclength up to dphi_max in the other !! direction. class(equilibrium_t), target, intent(inout) :: equi !! Equilibrium defining the field line real(FP), intent(in) :: x !! Spatial coordinate x at which the field line is located real(FP), intent(in) :: y !! Spatial coordinate y at which the field line is located real(FP), intent(in) :: phi !! Spatial coordinate phi at which the field line is located real(FP), intent(in) :: dphi_max !! Maximum toroidal angle to be traced real(FP), intent(in), optional :: maxstep !! Maximum stepsize for the DOP853 integrator real(FP) :: arclength_fwd !! Arclength in the positive direction until the target is reached real(FP) :: arclength_bwd !! Arclength in the negative direction until the target is reached logical:: target_reached_fwd !! True, if the target is reached during the forward tracing logical:: target_reached_bwd !! True, if the target is reached during the backward tracing if(present(maxstep)) then call distance_to_boundary(equi, x, y, phi, dphi_max, arclength_fwd, & arclength_bwd, target_reached_fwd, target_reached_bwd, maxstep) else call distance_to_boundary(equi, x, y, phi, dphi_max, arclength_fwd, & arclength_bwd, target_reached_fwd, target_reached_bwd, dphi_max) endif if((.not. target_reached_fwd) .and. (.not. target_reached_bwd)) then connection_length = FP_LARGEST return endif connection_length = arclength_fwd + arclength_bwd end function subroutine shortest_dist_to_bound(equi, x, y, phi, dphi_max, & shortest_dist, maxstep, dirind) !! It computes the shortest distance to the boundary along B. !! If (x, y, phi) is located inside the wall, the output is negative. !! If (x, y, phi) is located in the closed flux surface region, the !! output is the largest allowable number for PARALLAX precision. class(equilibrium_t), target, intent(inout) :: equi !! Equilibrium defining the field line real(FP), intent(in) :: x !! Spatial coordinate x at which the field line is located real(FP), intent(in) :: y !! Spatial coordinate y at which the field line is located real(FP), intent(in) :: phi !! Spatial coordinate phi at which the field line is located real(FP), intent(in) :: dphi_max !! Maximum toroidal angle to be traced real(FP), intent(out) :: shortest_dist !! Shortest distance to boundary along the field line from (x,y,phi) real(FP), intent(in), optional :: maxstep !! Maximum stepsize for the DOP853 integrator integer, intent(out), optional :: dirind !! Refers to the direction to the closest boundary along B: !! = 1, if the closest transition is in the +phi direction !! = -1, if the closest transition is in the -phi direction !! = 0, if the it could not found any boundary along the field line. real(FP) :: arclength_fwd !! Arclength in the positive direction until the target is reached real(FP) :: arclength_bwd !! Arclength in the negative direction until the target is reached logical:: target_reached_fwd !! True, if the target is reached during the forward tracing logical:: target_reached_bwd !! True, if the target is reached during the backward tracing if(present(maxstep)) then call distance_to_boundary(equi, x, y, phi, dphi_max, arclength_fwd, & arclength_bwd, target_reached_fwd, target_reached_bwd, maxstep) else call distance_to_boundary(equi, x, y, phi, dphi_max, arclength_fwd, & arclength_bwd, target_reached_fwd, target_reached_bwd, dphi_max) endif if((.not. target_reached_fwd) .and. (.not. target_reached_bwd)) then shortest_dist = FP_LARGEST if(present(dirind)) dirind = 0 return endif if(abs(arclength_fwd) < abs(arclength_bwd)) then shortest_dist = arclength_fwd if(present(dirind)) dirind = 1 else shortest_dist = arclength_bwd if(present(dirind)) dirind = -1 endif end subroutine end module