Subrotuine that checks if particles in the simulation are within the spatial domain where interpolants and fields are known.
External fields and interpolants can have different spatial domains where they are defined. Therefore, it is necessary to check if a given particle has left these spatial domains to stop following it, otherwise this will cause an error in the simulation.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(FIELDS), | intent(in) | :: | F | |||
| real(kind=rp), | intent(in), | DIMENSION(:,:), ALLOCATABLE | :: | Y | Particles' position in cylindrical coordinates, Y(1,:) = , Y(2,:) = , and Y(3,:) = . |
|
| integer(kind=is), | intent(inout), | DIMENSION(:), ALLOCATABLE | :: | flag | Flag that determines whether particles are followed in the simulation (flag=1), or not (flag=0). |
subroutine check_if_in_fields_domain(F,Y,flag)
!! @note Subrotuine that checks if particles in the simulation are within
!! the spatial domain where interpolants and fields are known. @endnote
!! External fields and interpolants can have different spatial domains where
!! they are defined. Therefore, it is necessary to
!! check if a given particle has left these spatial domains to stop
!! following it, otherwise this will cause an error in the simulation.
TYPE(FIELDS), INTENT(IN) :: F
REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: Y
!! Particles' position in cylindrical coordinates,
!! Y(1,:) = \(R\), Y(2,:) = \(\phi\), and Y(3,:) = \(Z\).
INTEGER(is), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: flag
!! Flag that determines whether particles are followed in the
!! simulation (flag=1), or not (flag=0).
INTEGER :: IR
!! Variable used to localize the grid cell in the \((R,\phi,Z)\)
!! or \((R,Z)\) grid containing the fields data that corresponds
!! to the radial position of the particles.
INTEGER :: IPHI
!! Variable used to localize the grid cell in the \((R,\phi,Z)\)
!! or \((R,Z)\) grid containing the fields data that corresponds
!! to the azimuthal position of the particles.
INTEGER :: IZ
!! Variable used to localize the grid cell in the \((R,\phi,Z)\)
!! or \((R,Z)\) grid containing the fields data that corresponds
!! to the vertical position of the particles.
INTEGER(ip) :: pp
!! Particle iterator.
INTEGER(ip) :: ss
!! Species iterator.
REAL(rp) :: Rwall
if (size(Y,1).eq.1) then
ss = size(Y,1)
else
if (Y(2,1).eq.0) then
ss=1_idef
else
ss = size(Y,1)
end if
endif
! write(output_unit_write,'("R: ",E15.10)') Y(1,1)
! write(output_unit_write,'("PHI: ",E15.10)') Y(2,1)
! write(output_unit_write,'("Z: ",E15.10)') Y(1,3)
! write(output_unit_write,*) 'Flag',flag(1)
if (ALLOCATED(fields_domain%FLAG3D)) then
if (F%Dim2x1t) then
!$OMP PARALLEL DO FIRSTPRIVATE(ss) PRIVATE(pp,IR,IZ) &
!$OMP& SHARED(Y,flag,fields_domain,bfield_2X1T)
do pp=1_idef,ss
IR = INT(FLOOR((Y(pp,1) - fields_domain%Ro + 0.5_rp* &
fields_domain%DR)/fields_domain%DR) + 1.0_rp,idef)
IZ = INT(FLOOR((Y(pp,3) + ABS(fields_domain%Zo) + 0.5_rp* &
fields_domain%DZ)/fields_domain%DZ) + 1.0_rp,idef)
if ((fields_domain%FLAG3D(IR,1,IZ).NE.1_is).OR. &
((IR.GT.bfield_2X1T%NR).OR.(IZ.GT.bfield_2X1T%NZ))) then
if (F%Analytic_IWL.eq.'NONE') then
flag(pp) = 0_is
else if (F%Analytic_IWL.eq.'D3D') then
if ((IR.lt.floor(bfield_2d%NR/6._rp)).and. &
(IZ.gt.floor(bfield_2d%NZ/5._rp)).and. &
(IZ.lt.floor(4._rp*bfield_2d%NZ/5._rp))) then
Rwall=F%circumradius*cos(C_PI/F%ntiles)/ &
(cos((mod(Y(pp,2),2*C_PI/F%ntiles))-C_PI/F%ntiles))
if (Y(pp,1).lt.Rwall) flag(pp) = 0_is
else
flag(pp) = 0_is
endif
endif
!write(output_unit_write,'("YR:",E17.10)') Y(1,1)
!write(output_unit_write,'("YZ:",E17.10)') Y(1,3)
!write(output_unit_write,'("IR: ",I16)') IR
!write(output_unit_write,'("IZ: ",I16)') IZ
end if
end do
!$OMP END PARALLEL DO
else
!$OMP PARALLEL DO FIRSTPRIVATE(ss) PRIVATE(pp,IR,IPHI,IZ) &
!$OMP& SHARED(Y,flag,fields_domain,bfield_3d)
do pp=1_idef,ss
IR = INT(FLOOR((Y(pp,1) - fields_domain%Ro + 0.5_rp* &
fields_domain%DR)/fields_domain%DR) + 1.0_rp,idef)
IPHI = INT(FLOOR((Y(pp,2) + 0.5_rp*fields_domain%DPHI)/ &
fields_domain%DPHI) + 1.0_rp,idef)
IZ = INT(FLOOR((Y(pp,3) + ABS(fields_domain%Zo) + 0.5_rp* &
fields_domain%DZ)/fields_domain%DZ) + 1.0_rp,idef)
if ((fields_domain%FLAG3D(IR,IPHI,IZ).NE.1_is).OR. &
((IR.GT.bfield_3d%NR).OR.(IZ.GT.bfield_3d%NZ))) then
flag(pp) = 0_is
!write(output_unit_write,'("YR:",E17.10)') Y(1,1)
!write(output_unit_write,'("YZ:",E17.10)') Y(1,3)
!write(output_unit_write,'("IR: ",I16)') IR
!write(output_unit_write,'("IZ: ",I16)') IZ
end if
end do
!$OMP END PARALLEL DO
end if
else if (ALLOCATED(fields_domain%FLAG2D)) then
!$OMP PARALLEL DO FIRSTPRIVATE(ss) PRIVATE(pp,IR,IZ) &
!$OMP& SHARED(Y,flag,fields_domain,bfield_2d)
do pp=1_idef,ss
IR = INT(FLOOR((Y(pp,1) - fields_domain%Ro + 0.5_rp* &
fields_domain%DR)/fields_domain%DR) + 1.0_rp,idef)
IZ = INT(FLOOR((Y(pp,3) + ABS(fields_domain%Zo) + 0.5_rp* &
fields_domain%DZ)/fields_domain%DZ) + 1.0_rp,idef)
#if DBG_CHECK
if ((IR.lt.0).or.(IZ.lt.0).or.(IR.GT. &
bfield_2d%NR).OR.(IZ.GT.bfield_2d%NZ)) then
write(6,'("YR:",E17.10)') Y(1,1)
write(6,'("YZ:",E17.10)') Y(1,3)
write(6,'("IR: ",I16)') IR
write(6,'("IZ: ",I16)') IZ
!call KORC_ABORT(23)
end if
#endif
!write(output_unit_write,'("IR: ",I16)') IR
!write(output_unit_write,'("IZ: ",I16)') IZ
if (((IR.lt.0).or.(IZ.lt.0).or. &
(IR.GT.bfield_2d%NR).OR.(IZ.GT.bfield_2d%NZ)).or. &
(fields_domain%FLAG2D(IR,IZ).NE.1_is)) then
!write(output_unit_write,*) 'here'
if (F%Analytic_IWL.eq.'NONE') then
flag(pp) = 0_is
else if (F%Analytic_IWL.eq.'D3D') then
if ((IR.lt.floor(bfield_2d%NR/6._rp)).and. &
(IZ.gt.floor(bfield_2d%NZ/5._rp)).and. &
(IZ.lt.floor(4._rp*bfield_2d%NZ/5._rp))) then
Rwall=F%circumradius*cos(C_PI/F%ntiles)/ &
(cos((mod(Y(pp,2),2*C_PI/F%ntiles))-C_PI/F%ntiles))
if (Y(pp,1).lt.Rwall) flag(pp) = 0_is
else
flag(pp) = 0_is
endif
endif
end if
end do
!$OMP END PARALLEL DO
end if
end subroutine check_if_in_fields_domain