From dcc832ce65c23901ffcf98c631ec7454f84dca7e Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Wed, 20 Sep 2023 15:12:17 +0200 Subject: [PATCH] fix simple limiter in analytic equilibrium In the case of analytic equilibrium without a limiter contour, the simple limiter was built incorrectly due to an unnecessary conversion from cm (the equilibrium data are already in metre). --- src/const_and_precisions.f90 | 1 + src/coreprofiles.f90 | 4 +-- src/gray_jetto1beam.f90 | 20 ++++++++----- src/main.f90 | 57 +++++++++++++++++++++++++----------- 4 files changed, 56 insertions(+), 26 deletions(-) diff --git a/src/const_and_precisions.f90 b/src/const_and_precisions.f90 index 644a8f6..bcfa581 100644 --- a/src/const_and_precisions.f90 +++ b/src/const_and_precisions.f90 @@ -57,5 +57,6 @@ module const_and_precisions real(wp_), parameter :: mc2_ = me_*c_**2/kev_ ! Electron rest energy, keV real(wp_), parameter :: mu0inv = 1._wp_/mu0_ ! Inverse magnetic permeability of vacuum, m/H real(wp_), parameter :: wce1_ = e_/me_ ! ECR (angular) frequency / magnetic field, rad/s/T + real(wp_), parameter :: cm = 0.01_wp_ ! Converts cm to m end module const_and_precisions diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90 index 014b7d6..b7926f6 100644 --- a/src/coreprofiles.f90 +++ b/src/coreprofiles.f90 @@ -451,9 +451,9 @@ contains ! invalid as they are given assuming a vacuum (N=1) if (present(launch_pos)) then block - use equilibrium, only : equinum_psi + use equilibrium, only : equinum_psi + use const_and_precisions, only : cm real(wp_) :: R, Z, psi - real(wp_), parameter :: cm = 1.0e-2_wp_ ! Convert from cartesian to cylindrical coordinates R = hypot(launch_pos(1), launch_pos(2)) * cm ! R = √(x²+y²) diff --git a/src/gray_jetto1beam.f90 b/src/gray_jetto1beam.f90 index d1a66bc..2251cb4 100644 --- a/src/gray_jetto1beam.f90 +++ b/src/gray_jetto1beam.f90 @@ -60,15 +60,21 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & call set_globals(params) end block init_params - ! Set a simple limiter + ! Set a simple limiter following the boundary of the data grid simple_limiter: block - use reflections, only : range2rect - use limiter, only : limiter_set_globals=>set_globals + use reflections, only : range2rect + use limiter, only : limiter_set_globals=>set_globals + use const_and_precisions, only : cm + + ! Avoid clipping out the launcher + real(wp_) :: R_launch, R_max + R_launch = norm2(params%antenna%pos(1:2)) * cm + R_max = max(R_launch, R(mr)) + + ! Convert to a list of R,z + call range2rect(xmin=params%misc%rwall, xmax=R_max, ymin=z(1), ymax=z(mz), & + xv=data%equilibrium%rlim, yv=data%equilibrium%zlim) - real(wp_) :: r0m - r0m = norm2(params%antenna%pos(1:2)) * 0.01_wp_ - call range2rect(params%misc%rwall, max(r0m, r(mr)), z(1), z(mz), & - data%equilibrium%rlim, data%equilibrium%zlim) ! Set the global variables of `limiter` call limiter_set_globals(data%equilibrium) end block simple_limiter diff --git a/src/main.f90 b/src/main.f90 index 2575b33..af9e3dc 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -384,8 +384,9 @@ contains subroutine init_misc(params, data) ! Performs miscellanous initial tasks, before the gray_main subroutine. - use reflections, only : range2rect - use limiter, only : limiter_set_globals=>set_globals + use reflections, only : range2rect + use limiter, only : limiter_set_globals=>set_globals + use const_and_precisions, only : cm implicit none @@ -393,27 +394,49 @@ contains type(gray_parameters), intent(inout) :: params type(gray_data), intent(inout) :: data - ! Build a basic limiter when one is not provided by the EQDSK + ! Build a basic limiter if one is not provided by equilibrium.txt if (.not. allocated(data%equilibrium%rlim) & .or. params%raytracing%ipass < 0) then - block - real(wp_) :: rmxm, r0m, z0m, dzmx - r0m = norm2(params%antenna%pos(1:2)) * 0.01_wp_ - dzmx = abs(params%raytracing%ipass) * & - params%raytracing%dst * params%raytracing%nstep * 0.01_wp_ - z0m = params%antenna%pos(3) * 0.01_wp_ - allocate(data%equilibrium%rlim(5)) - allocate(data%equilibrium%zlim(5)) - params%raytracing%ipass = abs(params%raytracing%ipass) + ! Restore sign of ipass + params%raytracing%ipass = abs(params%raytracing%ipass) + + allocate(data%equilibrium%rlim(5)) + allocate(data%equilibrium%zlim(5)) + block + real(wp_) :: R_launch, z_launch, R_max, delta_z + + ! Launcher coordinates + R_launch = norm2(params%antenna%pos(1:2)) * cm + z_launch = params%antenna%pos(3) * cm + + ! Max vertical distance a ray can travel + delta_z = abs(params%raytracing%ipass) * & + params%raytracing%dst * params%raytracing%nstep * cm + + ! Max radius, either due to the plasma extent or equilibrium grid if (params%equilibrium%iequil < 2) then - rmxm = (data%equilibrium%rv(1) + data%equilibrium%rv(2)) * 0.01_wp_ + ! analytic equilibrium, use R₀+a + R_max = data%equilibrium%rv(1) + data%equilibrium%rv(2) else - rmxm = data%equilibrium%rv(size(data%equilibrium%rv)) + ! numeric equilibrium, use max R of the grid + R_max = data%equilibrium%rv(size(data%equilibrium%rv)) end if - call range2rect(params%misc%rwall, max(r0m, rmxm), & - z0m - dzmx, z0m + dzmx, & - data%equilibrium%rlim, data%equilibrium%zlim) + + ! Avoid clipping out the launcher + R_max = max(R_launch, R_max) + + ! Convert to a list of R,z: + ! + ! (R_wall, z_launch+Δz)╔═════╗(R_max, z_launch+Δz) + ! ║4 3║ + ! ║ ║ + ! ║1 2║ + ! (R_wall, z_launch-Δz)╚═════╝(R_max, z_launch-Δz) + ! + call range2rect(xmin=params%misc%rwall, xmax=R_max, & + ymin=z_launch - delta_z, ymax=z_launch + delta_z, & + xv=data%equilibrium%rlim, yv=data%equilibrium%zlim) end block end if