From 3a10b455953b6a192a17ae8cbbf1d54eed821cea Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Tue, 30 Jul 2024 10:57:07 +0200 Subject: [PATCH] src/limiter.f90: remove 1. Use the `contour` type for limiter and plasma boundary (rlim, zlim, rbnd, zbnd) 2. Replace `inside` with `contour%contains` 3. Replace `range2rect` with a `contour` interface 4. Remove the limiter module which just re-exports the limiter as a global; instead just pass the contour object around --- src/equilibrium.f90 | 63 ++++++++++++++--------------- src/gray_core.f90 | 20 ++++----- src/gray_jetto1beam.f90 | 16 +++----- src/gray_params.f90 | 8 ++-- src/gray_tables.f90 | 8 ++-- src/limiter.f90 | 42 ------------------- src/magsurf_data.f90 | 27 ++++++------- src/main.f90 | 20 +++------ src/multipass.f90 | 27 +++++++------ src/reflections.f90 | 79 +++++++++++++++++++----------------- src/types.f90 | 89 +++++++++++++++++++++++++++++++++++++++++ 11 files changed, 215 insertions(+), 184 deletions(-) delete mode 100644 src/limiter.f90 diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index b3cf91c..4e386b9 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -7,6 +7,7 @@ ! the data is interpolated using splines. module equilibrium use const_and_precisions, only : wp_ + use types, only : contour use splines, only : spline_simple, spline_1d, spline_2d, linear_1d use gray_params, only : EQ_VACUUM, EQ_ANALYTICAL, & EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL @@ -24,12 +25,6 @@ module equilibrium real(wp_) :: B0 ! Magnetic field at the magnetic axis (T) end type - ! A 2D contour in the (R,z) plane - type contour - real(wp_), allocatable :: R(:) - real(wp_), allocatable :: z(:) - end type - ! Order of the splines integer, parameter :: kspl=3, ksplp=kspl + 1 @@ -105,6 +100,7 @@ contains character(len=48) :: string real(wp_) :: dr, dz, dps, rleft, zmid, zleft, psiedge, psiaxis real(wp_) :: xdum ! dummy variable, used to discard data + real(wp_), allocatable :: R(:), z(:) u = get_free_unit(unit) @@ -161,30 +157,30 @@ contains end if ! Get size of boundary and limiter arrays and allocate them - read (u,*) nbnd, nlim - if (allocated(data%rbnd)) deallocate(data%rbnd) - if (allocated(data%zbnd)) deallocate(data%zbnd) - if (allocated(data%rlim)) deallocate(data%rlim) - if (allocated(data%zlim)) deallocate(data%zlim) + read (u, *) nbnd, nlim ! Load plasma boundary data if(nbnd > 0) then - allocate(data%rbnd(nbnd), data%zbnd(nbnd)) + allocate(R(nbnd), z(nbnd)) if (params%ifreefmt) then - read(u, *) (data%rbnd(i), data%zbnd(i), i=1,nbnd) + read(u, *) (R(i), z(i), i=1,nbnd) else - read(u, '(5e16.9)') (data%rbnd(i), data%zbnd(i), i=1,nbnd) + read(u, '(5e16.9)') (R(i), z(i), i=1,nbnd) end if + data%boundary = contour(R, z) + deallocate(R, z) end if ! Load limiter data if(nlim > 0) then - allocate(data%rlim(nlim), data%zlim(nlim)) + allocate(R(nlim), z(nlim)) if (params%ifreefmt) then - read(u, *) (data%rlim(i), data%zlim(i), i=1,nlim) + read(u, *) (R(i), z(i), i=1,nlim) else - read(u, '(5e16.9)') (data%rlim(i), data%zlim(i), i=1,nlim) + read(u, '(5e16.9)') (R(i), z(i), i=1,nlim) end if + data%limiter = contour(R, z) + deallocate(R, z) end if ! End of G-EQDSK file @@ -229,6 +225,7 @@ contains ! local variables integer :: i, u, nlim + real(wp_), allocatable :: R(:), z(:) u = get_free_unit(unit) @@ -247,13 +244,12 @@ contains if(ipass >= 2) then ! get size of boundary and limiter arrays and allocate them read (u,*) nlim - if (allocated(data%rlim)) deallocate(data%rlim) - if (allocated(data%zlim)) deallocate(data%zlim) ! store boundary and limiter data - if(nlim > 0) then - allocate(data%rlim(nlim), data%zlim(nlim)) - read(u,*) (data%rlim(i), data%zlim(i), i = 1, nlim) + if (nlim > 0) then + allocate(R(nlim), z(nlim)) + read(u,*) (R(i), z(i), i = 1, nlim) + data%limiter = contour(R, z) end if end if close(u) @@ -442,7 +438,7 @@ contains case (EQ_EQDSK_PARTIAL) ! Data valid only inside boundary (data%psin=0 outside), ! presence of boundary anticipated here to filter invalid data - nbnd = min(size(data%rbnd), size(data%zbnd)) + nbnd = size(data%boundary%R) ! allocate knots and spline coefficients arrays if (allocated(psi_spline%knots_x)) deallocate(psi_spline%knots_x) @@ -456,7 +452,7 @@ contains do j=1,nz do i=1,nr if (nbnd.gt.0) then - if(.not.inside(data%rbnd,data%zbnd,data%rv(i),data%zv(j))) cycle + if (.not. data%boundary%contains(data%rv(i), data%zv(j))) cycle else if(data%psin(i,j).le.0.0d0) cycle end if @@ -470,7 +466,7 @@ contains do j=1,nz do i=1,nr if (nbnd.gt.0) then - if(.not.inside(data%rbnd,data%zbnd,data%rv(i),data%zv(j))) cycle + if (.not. data%boundary%contains(data%rv(i), data%zv(j))) cycle else if(data%psin(i,j).le.0.0d0) cycle end if @@ -574,12 +570,12 @@ contains ! Use provided boundary to set an initial guess ! for the search of O/X points - nbnd=min(size(data%rbnd), size(data%zbnd)) - if (nbnd>0) then - call vmaxmini(data%zbnd,nbnd,zbinf,zbsup,ibinf,ibsup) - rbinf=data%rbnd(ibinf) - rbsup=data%rbnd(ibsup) - call vmaxmin(data%rbnd,nbnd,rbmin,rbmax) + nbnd = size(data%boundary%R) + if (nbnd > 0) then + call vmaxmini(data%boundary%z, nbnd, zbinf, zbsup, ibinf, ibsup) + rbinf = data%boundary%R(ibinf) + rbsup = data%boundary%R(ibsup) + call vmaxmin(data%boundary%R, nbnd, rbmin, rbmax) else zbinf=data%zv(2) zbsup=data%zv(nz-1) @@ -673,8 +669,7 @@ contains call set_rho_spline(sqrt(data%psinr), rhotn) ! Compute the domain of the ψ mapping - psi_domain%R = data%rbnd - psi_domain%z = data%zbnd + psi_domain = data%boundary call rescale_boundary(psi_domain, O=[rmaxis, zmaxis], t0=1.5_wp_) end subroutine set_equil_spline @@ -1015,7 +1010,7 @@ contains case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) ! Numerical data - if (inside(psi_domain%R, psi_domain%z, R, z)) then + if (psi_domain%contains(R, z)) then ! Within the interpolation range if (present(psi_n)) psi_n = psi_spline%eval(R, z) if (present(dpsidr)) dpsidr = psi_spline%deriv(R, z, 1, 0) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 78020ca..d671bc9 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -53,7 +53,7 @@ contains integer :: iox,nharm,nhf,iokhawa,ierrn,ierrhcd, index_rt, parent_index_rt integer :: ip,ib,iopmin,ipar,child_index_rt - integer :: igrad_b,istop_pass,nbeam_pass,nlim + integer :: igrad_b,istop_pass,nbeam_pass logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff ! i: integration step, jk: global ray index @@ -103,9 +103,6 @@ contains real(wp_), dimension(params%output%nrho) :: currins_beam, dpdv_beam, jcd_beam ! ======== set environment BEGIN ======== - ! Number of limiter contourn points - nlim = size(data%equilibrium%zlim) - ! Compute X=ω/ω_ce and Y=(ω/ω_pe)² (with B=1) call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn) @@ -121,7 +118,7 @@ contains if (params%equilibrium%iequil /= EQ_VACUUM) then ! Initialise the magsurf_data module - call compute_flux_averages(params%equilibrium, results%tables) + call compute_flux_averages(params, results%tables) ! Initialise the output profiles call pec_init(params%output, rhout) @@ -145,7 +142,7 @@ contains ! print ψ surface for q=3/2 and q=2/1 call find_flux_surfaces( & - qvals=[1.5_wp_, 2.0_wp_], params=params%equilibrium, & + qvals=[1.5_wp_, 2.0_wp_], params=params, & tbl=results%tables%flux_surfaces) ! print initial position @@ -231,8 +228,7 @@ contains zzm = yw(3,jk)*0.01_wp_ rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_ - if(inside(data%equilibrium%rlim, data%equilibrium%zlim, & - rrm, zzm)) then ! * start propagation in/outside vessel? + if (data%equilibrium%limiter%contains(rrm, zzm)) then ! * start propagation in/outside vessel? iow(jk) = 1 ! + inside else iow(jk) = 0 ! + outside @@ -295,8 +291,7 @@ contains rrm = sqrt(xv(1)*xv(1)+xv(2)*xv(2))*0.01_wp_ ins_pl = (psinv>=zero .and. psinvin else if(ext_wl) then ! ray exit vessel - call wall_out(jk, ins_pl, xv, anv, params%raytracing%dst, & - bres, sox, psipol, chipol, iow, iop, ext, eyt) + call wall_out(jk, data%equilibrium%limiter, ins_pl, xv, anv, & + params%raytracing%dst, bres, sox, psipol, chipol, & + iow, iop, ext, eyt) yw(:,jk) = [xv, anv] ! * updated coordinates (reflected) igrad_b = 0 ! * switch to ray-tracing diff --git a/src/gray_jetto1beam.f90 b/src/gray_jetto1beam.f90 index f8c1f2b..e96abc2 100644 --- a/src/gray_jetto1beam.f90 +++ b/src/gray_jetto1beam.f90 @@ -48,21 +48,17 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & ! Set a simple limiter following the boundary of the data grid simple_limiter: block - use utils, only : range2rect - use limiter, only : limiter_set_globals=>set_globals use const_and_precisions, only : cm + use types, only : contour ! 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), & - x=data%equilibrium%rlim, y=data%equilibrium%zlim) - - ! Set the global variables of `limiter` - call limiter_set_globals(data%equilibrium) + ! Convert to a closed contour + data%equilibrium%limiter = contour( & + Rmin=params%misc%rwall, Rmax=R_max, zmin=z(1), zmax=z(mz)) end block simple_limiter end if first_call @@ -70,6 +66,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & ! Set MHD equilibrium data init_equilibrium: block use equilibrium, only : set_equil_spline + use types, only : contour ! Copy argument arrays data%equilibrium%rv = r @@ -79,11 +76,10 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & data%equilibrium%zax = zax data%equilibrium%psinr = psrad data%equilibrium%fpol = fpol - data%equilibrium%rbnd = rbnd - data%equilibrium%zbnd = zbnd data%equilibrium%psia = psia data%equilibrium%psin = psin data%equilibrium%qpsi = qpsi + data%equilibrium%boundary = contour(rbnd, zbnd) ! Compute splines call set_equil_spline(params%equilibrium, data%equilibrium, err) diff --git a/src/gray_params.f90 b/src/gray_params.f90 index 108df58..97210d7 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -1,7 +1,7 @@ module gray_params use const_and_precisions, only : wp_ - use types, only : table + use types, only : table, contour implicit none @@ -171,10 +171,8 @@ module gray_params type equilibrium_data real(wp_), allocatable :: rv(:) ! R of the uniform grid real(wp_), allocatable :: zv(:) ! Z of the uniform grid - real(wp_), allocatable :: rlim(:) ! R of the limiter contour (wall) - real(wp_), allocatable :: zlim(:) ! Z of the limiter contour - real(wp_), allocatable :: rbnd(:) ! R of the boundary contour (plasma) - real(wp_), allocatable :: zbnd(:) ! Z of the boundary contour + type(contour) :: limiter ! limiter contour (wall) + type(contour) :: boundary ! boundary contour (plasma) real(wp_), allocatable :: fpol(:) ! Poloidal current function real(wp_), allocatable :: qpsi(:) ! Safety factor on the flux grid real(wp_), allocatable :: psin(:,:) ! Poloidal flux on a uniform grid diff --git a/src/gray_tables.f90 b/src/gray_tables.f90 index fa5848c..7c9fc75 100644 --- a/src/gray_tables.f90 +++ b/src/gray_tables.f90 @@ -288,7 +288,7 @@ contains ! Finds the ψ for a set of values of q and stores the ! associated surfaces to the flux surface table - use gray_params, only : equilibrium_parameters + use gray_params, only : gray_parameters use equilibrium, only : fq, frhotor, rmaxis, zmaxis, zbsup, zbinf use magsurf_data, only : contours_psi, npoints use logger, only : log_info @@ -296,9 +296,9 @@ contains use types, only : table ! subroutine arguments - real(wp_), intent(in) :: qvals(:) - type(equilibrium_parameters), intent(in) :: params - type(table), intent(inout) :: tbl + real(wp_), intent(in) :: qvals(:) + type(gray_parameters), intent(in) :: params + type(table), intent(inout) :: tbl ! local variables integer :: i diff --git a/src/limiter.f90 b/src/limiter.f90 deleted file mode 100644 index b0b9ce7..0000000 --- a/src/limiter.f90 +++ /dev/null @@ -1,42 +0,0 @@ -module limiter - use const_and_precisions, only : wp_ - - implicit none - - ! Inner wall radius - real(wp_), save :: rwallm - - ! Limiter contourn - integer, public, save :: nlim - real(wp_), dimension(:), allocatable, save :: rlim, zlim - -contains - - subroutine set_globals(data) - ! Set global variables exposed by this module. - use gray_params, only : equilibrium_data - - ! subroutine arguments - type(equilibrium_data), intent(in) :: data - - if (allocated(rlim)) deallocate(rlim) - if (allocated(zlim)) deallocate(zlim) - nlim = size(data%rlim) - allocate(rlim(nlim), zlim(nlim)) - rlim = data%rlim - zlim = data%zlim - rwallm = minval(rlim) - end subroutine set_globals - - - subroutine unset_globals - ! Unset global variables exposed by this module. - use const_and_precisions, only : zero - - if(allocated(rlim)) deallocate(rlim) - if(allocated(zlim)) deallocate(zlim) - nlim = 0 - rwallm = zero - end subroutine unset_globals - -end module limiter diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index a35f33f..027e54e 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -101,11 +101,11 @@ contains bfield,frhotor,fq,tor_curr,psia,pol_flux use dierckx, only : regrid,coeff_parder use types, only : table, wrap - use gray_params, only : equilibrium_parameters, gray_tables + use gray_params, only : gray_parameters, gray_tables ! subroutine arguments - type(equilibrium_parameters), intent(in) :: params - type(gray_tables), intent(inout), optional :: tables + type(gray_parameters), intent(in) :: params + type(gray_tables), intent(inout), optional :: tables ! local constants integer, parameter :: nnintp=101,ncnt=100,nlam=101,ksp=3, & @@ -210,7 +210,7 @@ contains if (rhotjp /= rhotjp) print *, 'Nan!!!!!!!!!!' psicon(jp)=height - call contours_psi(params,psinjp,rctemp,zctemp,rup,zup,rlw,zlw) + call contours_psi(params, psinjp, rctemp, zctemp, rup, zup, rlw, zlw) rcon(:,jp) = rctemp zcon(:,jp) = zctemp @@ -446,23 +446,22 @@ contains - subroutine contours_psi(params, h,rcn,zcn,rup,zup,rlw,zlw) - use const_and_precisions, only : wp_,pi - use gray_params, only : equilibrium_parameters + subroutine contours_psi(params, h, rcn, zcn, rup, zup, rlw, zlw) + use const_and_precisions, only : wp_, pi + use gray_params, only : gray_parameters use logger, only : log_warning use dierckx, only : profil, sproota use equilibrium, only : model, frhotor, psi_spline, & kspl, points_tgo - use limiter, only : rwallm ! local constants integer, parameter :: mest=4 ! subroutine arguments - type(equilibrium_parameters), intent(in) :: params - real(wp_), intent(in) :: h - real(wp_), intent(out) :: rcn(:), zcn(:) - real(wp_), intent(inout) :: rup, zup, rlw, zlw + type(gray_parameters), intent(in) :: params + real(wp_), intent(in) :: h + real(wp_), intent(out) :: rcn(:), zcn(:) + real(wp_), intent(inout) :: rup, zup, rlw, zlw ! local variables integer :: npoints,np,info,ic,ier,iopt,m @@ -474,7 +473,7 @@ contains np=(npoints-1)/2 th=pi/dble(np) - if (params%iequil<2) then + if (params%equilibrium%iequil<2) then block real(wp_) :: r_p ! poloidal radius r_p = sqrt(h) * model%a @@ -515,7 +514,7 @@ contains end if call sproota(h, psi_spline%knots_x, psi_spline%nknots_x, & czc, zeroc, mest, m, ier) - if (zeroc(1).gt.rwallm) then + if (zeroc(1) > params%misc%rwall) then rcn(ic)=zeroc(1) rcn(npoints+1-ic)=zeroc(2) else diff --git a/src/main.f90 b/src/main.f90 index 0f70415..98e0b96 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -138,11 +138,9 @@ program main free_mem: block use coreprofiles, only : unset_profiles_spline use equilibrium, only : unset_equil_spline - use limiter, only : limiter_unset_globals=>unset_globals call unset_profiles_spline call unset_equil_spline - call limiter_unset_globals end block free_mem end block no_save @@ -296,8 +294,7 @@ contains ! Performs miscellanous initial tasks, before the gray_main subroutine. use gray_params, only : EQ_VACUUM, EQ_ANALYTICAL, & EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL - use utils, only : range2rect - use limiter, only : limiter_set_globals=>set_globals + use types, only : contour use const_and_precisions, only : cm, comp_huge ! subroutine arguments @@ -305,14 +302,12 @@ contains type(gray_data), intent(inout) :: data ! Build a basic limiter if one is not provided by equilibrium.txt - if (.not. allocated(data%equilibrium%rlim) & + if (.not. allocated(data%equilibrium%limiter%R) & .or. params%raytracing%ipass < 0) then ! 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 @@ -353,14 +348,11 @@ contains ! ║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, & - x=data%equilibrium%rlim, y=data%equilibrium%zlim) + data%equilibrium%limiter = contour( & + Rmin=params%misc%rwall, Rmax=R_max, & + zmin=z_launch - delta_z, zmax=z_launch + delta_z) end block end if - - ! Set the global variables of the `limiter` module - call limiter_set_globals(data%equilibrium) end subroutine init_misc @@ -393,7 +385,7 @@ contains integer, allocatable :: beams(:) ! Initialise the magsurf_data module - call compute_flux_averages(params%equilibrium) + call compute_flux_averages(params) ! Initialise the output profiles call pec_init(params%output) diff --git a/src/multipass.f90 b/src/multipass.f90 index 2771def..ef1e123 100644 --- a/src/multipass.f90 +++ b/src/multipass.f90 @@ -99,21 +99,23 @@ contains end subroutine plasma_out - subroutine wall_out(i, ins, xv, anv, dst, bres, sox, & + subroutine wall_out(i, wall, ins, xv, anv, dst, bres, sox, & psipol1, chipol1, iow, iop, ext, eyt) ! Ray exits vessel + use types, only : contour ! subroutine arguments - integer, intent(in) :: i ! ray index - logical, intent(in) :: ins ! inside plasma? (ins=1 plasma/wall overlap) - real(wp_), intent(inout) :: xv(3), anv(3) - real(wp_), intent(in) :: dst ! step size - real(wp_), intent(in) :: bres - integer, intent(in) :: sox - real(wp_), intent(out) :: psipol1, chipol1 - integer, intent(inout) :: iow(:) ! in/out vessel and plasma flags - integer, intent(inout) :: iop(:) ! in/out vessel and plasma flags - complex(wp_), intent(inout) :: ext(:), eyt(:) + integer, intent(in) :: i ! ray index + type(contour), intent(in) :: wall ! wall contour + logical, intent(in) :: ins ! inside plasma? (ins=1 plasma/wall overlap) + real(wp_), intent(inout) :: xv(3), anv(3) + real(wp_), intent(in) :: dst ! step size + real(wp_), intent(in) :: bres + integer, intent(in) :: sox + real(wp_), intent(out) :: psipol1, chipol1 + integer, intent(inout) :: iow(:) ! in/out vessel and plasma flags + integer, intent(inout) :: iop(:) ! in/out vessel and plasma flags + complex(wp_), intent(inout) :: ext(:), eyt(:) ! local variables integer :: irfl @@ -122,7 +124,8 @@ contains iow(i)=iow(i)+1 ! out->in if(ins) call plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) ! plasma-wall overlapping - call wall_refl(xv-dst*anv,anv,ext(i),eyt(i),xvrfl,anvrfl,ext1,eyt1,walln,irfl) ! ray reflects at wall + call wall_refl(wall, xv-dst*anv, anv, ext(i), eyt(i), & + xvrfl, anvrfl, ext1, eyt1, walln, irfl) ! ray reflects at wall ext(i) = ext1 ! save parameters at wall reflection eyt(i) = eyt1 xv = xvrfl diff --git a/src/reflections.f90 b/src/reflections.f90 index 9f2cf8f..2c63f8f 100644 --- a/src/reflections.f90 +++ b/src/reflections.f90 @@ -26,22 +26,27 @@ end subroutine reflect -subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw) - real(wp_), intent(in), dimension(3) :: xv,kv - integer, intent(in) :: nw - real(wp_), dimension(nw), intent(in) :: rw,zw - real(wp_), intent(out) :: sint - real(wp_), dimension(3), intent(out) :: normw - integer :: i,j,ni,iint,first - real(wp_), dimension(2) :: si,ti - real(wp_) :: drw,dzw,xint,yint,rint,l,kxy +subroutine inters_linewall(xv, kv, wall, sint, normw) + use types, only : contour + + ! subroutine arguments + real(wp_), intent(in) :: xv(3), kv(3) + type(contour), intent(in) :: wall + real(wp_), intent(out) :: sint + real(wp_), intent(out) :: normw(3) + + ! local variables + integer :: i, j, ni, iint, first + real(wp_) :: si(2), ti(2) + real(wp_) :: drw, dzw, xint, yint, rint, l, kxy + sint=comp_huge iint=0 normw=zero - do i=1,nw-1 - !search intersections with i-th wall segment - call linecone_coord(xv,kv,rw(i:i+1),zw(i:i+1),si,ti,ni) - !discard solutions with s<=0 + do i=1, size(wall%R)-1 + ! search intersections with i-th wall segment + call linecone_coord(xv, kv, wall%R(i:i+1), wall%z(i:i+1), si, ti, ni) + ! discard solutions with s<=0 first=ni+1 do j=1,ni if (si(j)>zero) then @@ -51,16 +56,16 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw) end do do j=first,ni if ((si(j)=zero .and. ti(j)<=one) then - !check intersection is in r,z range and keep the closest + ! check intersection is in r,z range and keep the closest sint = si(j) iint = i end if end do end do if (iint==0) return - !calculate wall normal at intersection point - drw = rw(iint+1)-rw(iint) - dzw = zw(iint+1)-zw(iint) + ! calculate wall normal at intersection point + drw = wall%R(iint+1) - wall%R(iint) + dzw = wall%z(iint+1) - wall%z(iint) xint = xv(1) + sint*kv(1) yint = xv(2) + sint*kv(2) rint = hypot(xint, yint) @@ -74,7 +79,7 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw) normw(1) = kv(1)/kxy*dzw/l normw(2) = kv(2)/kxy*dzw/l end if - !reverse normal if k.n>0 + ! reverse normal if k.n>0 if (dot_product(normw,kv)>zero) normw=-normw end subroutine inters_linewall @@ -187,30 +192,31 @@ function interssegm(xa,ya,xb,yb) end function interssegm -subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln,irfl) - use limiter, only : rlim,zlim,nlim - use utils, only : inside -! arguments +subroutine wall_refl(wall, xv, anv, ext, eyt, xvrfl, & + anvrfl, extr, eytr, walln, irfl) + use types, only : contour + + ! subroutine arguments + type(contour), intent(in) :: wall integer :: irfl real(wp_), dimension(3) :: xv,anv,xvrfl,anvrfl,walln complex(wp_) :: ext,eyt,extr,eytr -! local variables + ! local variables real(wp_) :: smax,rrm,zzm real(wp_), dimension(3) :: anv0,vv1,vv2,vv3 complex(wp_) :: eztr complex(wp_), dimension(3) :: evin,evrfl -! + anv0=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2) rrm=1.0e-2_wp_*sqrt(xv(1)**2+xv(2)**2) zzm=1.0e-2_wp_*xv(3) -! -! computation of reflection coordinates and normal to the wall - call inters_linewall(xv/1.0e2_wp_,anv0,rlim(1:nlim),zlim(1:nlim), & - nlim,smax,walln) + + ! computation of reflection coordinates and normal to the wall + call inters_linewall(xv/1.0e2_wp_, anv0, wall, smax, walln) smax=smax*1.0e2_wp_ xvrfl=xv+smax*anv0 irfl=1 - if (.not.inside(rlim,zlim,rrm,zzm)) then + if (.not. wall%contains(rrm,zzm)) then ! first wall interface is outside-inside if (dot_product(walln,walln) null() ! pointer to the next item end type + type queue ! A queue is a list of items with O(1) insertion and extraction. ! The first item inserted (`put` operation) is the first to be @@ -24,6 +25,7 @@ module types procedure :: empty => queue_empty ! checks whether the queue is empty end type + type table ! A type for storing tabular data before serialisation integer :: id ! the table unique ID @@ -39,6 +41,7 @@ module types procedure :: save => table_save ! processes and writes the table to file end type + type wrap ! A wrapper type for storing (references of) heterogeneous ! values in an array or other homogeneuous container type @@ -50,6 +53,21 @@ module types procedure :: wrap_init end interface + + type contour + ! A closed contour in the (R,z) plane + real(wp_), allocatable :: R(:) + real(wp_), allocatable :: z(:) + contains + procedure :: contains => contour_contains ! test if contour contains a point + end type + + ! Interface for custom type constructor + interface contour + procedure :: contour_init + procedure :: contour_init_rect + end interface + contains subroutine queue_put(self, val) @@ -242,6 +260,77 @@ contains end subroutine table_save + function contour_init(R, z) result(self) + ! Creates a contour + + ! functions arguments + real(wp_), intent(in) :: R(:), z(:) + type(contour) :: self + + ! local variables + integer :: n + + ! Ensure the first and last point are the same + n = size(R) + self%R = [R, R(1)] + self%z = [z, z(1)] + end function contour_init + + + function contour_init_rect(Rmin, Rmax, zmin, zmax) result(self) + ! Given two ranges [Rmin, Rmax], [zmin, zmax] creates a + ! rectangular contour as follows: + ! + ! (Rmin, zmax)╔═════╗(Rmax, zmax) + ! ║4 3║ + ! ║ ║ + ! ║1 2║ + ! (Rmin, zmin)╚═════╝(Rmax, zmax) + ! + + ! subroutine arguments + real(wp_), intent(in) :: Rmin, Rmax, zmin, zmax + type(contour) :: self + + self = contour_init([Rmin, Rmax, Rmax, Rmin], [zmin, zmin, zmax, zmax]) + end function contour_init_rect + + + function contour_contains(self, R0, z0) result(inside) + ! Tests whether the point (`R`, `z`) lies inside the 2D contour + + use utils, only : intlinf, locate_unord + + ! subroutine arguments + class(contour), intent(in) :: self + real(wp_), intent(in) :: R0, z0 + logical :: inside + + ! local variables + integer :: seg(size(self%R)), i, nsegs + real(wp_) :: R + + inside = .false. + + ! Find the `nsegs` segments that intersect the horizontal + ! line z=z0, i.e. `self%z(seg(i)) < z0 < self%z(seg(i)+1)` + call locate_unord(self%z, z0, seg, size(self%R), nsegs) + + ! No intersections, it must be outside (above or below) + if (nsegs == 0) return + + ! Count the number of intersections that lie to the left + ! (equivalently, to the right) of the point. An even number + ! means that the point is outside the polygon. + do i = 1, nsegs + ! coordinate of the intersection between segment and z=z0 + R = intlinf(self%z(seg(i)), self%R(seg(i)), & + self%z(seg(i)+1), self%R(seg(i)+1), z0) + if (R < R0) inside = .not. inside + end do + end function contour_contains + + subroutine test_queue() integer, target :: x=1, y=2 real, target :: z = 1.231