Skip to content

Commit

Permalink
added collimation with round or rectangular aperture option by C-H Hu…
Browse files Browse the repository at this point in the history 
…ang.
  • Loading branch information
@qianglbl
qianglbl committed Aug 17, 2024
1 parent 4978d39 commit 3b228d9
Show file tree
Hide file tree
Showing 3 changed files with 222 additions and 29 deletions.
130 changes: 124 additions & 6 deletions src/Appl/BeamBunch.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1123,8 +1123,9 @@ subroutine kick2t_BeamBunch(innp,innx,inny,innz,rays,exg,&

end subroutine kick2t_BeamBunch

!check the particles outside the computational domain for each bunch/bin
subroutine lost_BeamBunch(this,xrad,yrad,zleng,zcent,nplc,nptot)
!------------------------------------------------------------------
!revised by C-H. Huang
subroutine lostREC_BeamBunch(this,xrad,yrad,zleng,zcent,nplc,nptot)
implicit none
include 'mpif.h'
type (BeamBunch), intent(inout) :: this
Expand Down Expand Up @@ -1154,9 +1155,9 @@ subroutine lost_BeamBunch(this,xrad,yrad,zleng,zcent,nplc,nptot)
this%Pts1(5,i) = this%Pts1(5,i0)
this%Pts1(6,i) = this%Pts1(6,i0)
rr = this%Pts1(1,i0)**2+this%Pts1(3,i0)**2
if(rr.ge.rrap) then
ilost = ilost + 1
else if(this%Pts1(1,i0).le.ptrange(1)) then
!if(rr.ge.rrap) then
! ilost = ilost + 1
if(this%Pts1(1,i0).le.ptrange(1)) then
ilost = ilost + 1
else if(this%Pts1(1,i0).ge.ptrange(2)) then
ilost = ilost + 1
Expand Down Expand Up @@ -1189,7 +1190,76 @@ subroutine lost_BeamBunch(this,xrad,yrad,zleng,zcent,nplc,nptot)
this%Current = this%Current*nptot/this%Npt
this%Npt = nptot

end subroutine lost_BeamBunch
end subroutine lostREC_BeamBunch


!check the particles outside the computational domain for each bunch/bin
subroutine lostDFO_BeamBunch(this,xrad,yrad,zleng,zcent,nplc,nptot)
implicit none
include 'mpif.h'
type (BeamBunch), intent(inout) :: this
double precision, intent(in) :: xrad,yrad,zleng,zcent
integer, intent(out) :: nplc,nptot
integer :: i
integer :: ilost,i0,ierr
double precision, dimension(6) :: ptrange
double precision :: rr,rrap

ptrange(1) = -xrad/Scxlt
ptrange(2) = xrad/Scxlt
ptrange(3) = -yrad/Scxlt
ptrange(4) = yrad/Scxlt
!ptrange(5) = zcent-0.5d0*zleng/Scxlt
!ptrange(6) = zcent+0.5d0*zleng/Scxlt
ptrange(5) = 0.0
ptrange(6) = zleng/Scxlt
rrap = ptrange(1)**2
ilost = 0
do i0 = 1, this%Nptlocal
i = i0 - ilost
this%Pts1(1,i) = this%Pts1(1,i0)
this%Pts1(2,i) = this%Pts1(2,i0)
this%Pts1(3,i) = this%Pts1(3,i0)
this%Pts1(4,i) = this%Pts1(4,i0)
this%Pts1(5,i) = this%Pts1(5,i0)
this%Pts1(6,i) = this%Pts1(6,i0)
rr = this%Pts1(1,i0)**2+this%Pts1(3,i0)**2
if(rr.ge.rrap) then
ilost = ilost + 1
!if(this%Pts1(1,i0).le.ptrange(1)) then
! ilost = ilost + 1
!else if(this%Pts1(1,i0).ge.ptrange(2)) then
! ilost = ilost + 1
!else if(this%Pts1(3,i0).le.ptrange(3)) then
! ilost = ilost + 1
!else if(this%Pts1(3,i0).ge.ptrange(4)) then
! ilost = ilost + 1
else if(this%Pts1(5,i0).le.ptrange(5) .and. &
this%Pts1(6,i0).lt.0.0) then
ilost = ilost + 1
else if(this%Pts1(5,i0).ge.ptrange(6)) then
ilost = ilost + 1
! else if(this%Pts1(6,i0).lt.0.0) then !this does not allow particles move in negative direction
! ilost = ilost + 1
else
endif
enddo
do i = this%Nptlocal - ilost + 1, this%Nptlocal
this%Pts1(1,i) = 0.0
this%Pts1(2,i) = 0.0
this%Pts1(3,i) = 0.0
this%Pts1(4,i) = 0.0
this%Pts1(5,i) = -1.0e5
this%Pts1(6,i) = 0.0
enddo
this%Nptlocal = this%Nptlocal - ilost
nplc = this%Nptlocal
call MPI_ALLREDUCE(nplc,nptot,1,MPI_INTEGER,&
MPI_SUM,MPI_COMM_WORLD,ierr)
this%Current = this%Current*nptot/this%Npt
this%Npt = nptot

end subroutine lostDFO_BeamBunch


!check the particles outside the computational domain for each bunch/bin
Expand Down Expand Up @@ -1249,6 +1319,54 @@ subroutine lostXY_BeamBunch(this,xradmin,xradmax,yradmin,yradmax,nplc,nptot)
this%Npt = nptot

end subroutine lostXY_BeamBunch


subroutine lostROUND_BeamBunch(this,xradmin,xradmax,yradmin,yradmax,nplc,nptot)
implicit none
include 'mpif.h'
type (BeamBunch), intent(inout) :: this
double precision, intent(in) :: xradmin,xradmax,yradmin,yradmax
integer, intent(out) :: nplc,nptot
integer :: i
integer :: ilost,i0,ierr
double precision, dimension(6) :: ptrange
double precision :: rr,rrap
ptrange(1) = xradmin/Scxlt
ptrange(2) = xradmax/Scxlt
ptrange(3) = yradmin/Scxlt
ptrange(4) = yradmax/Scxlt
rrap = ptrange(1)**2
ilost = 0
do i0 = 1, this%Nptlocal
i = i0 - ilost
this%Pts1(1,i) = this%Pts1(1,i0)
this%Pts1(2,i) = this%Pts1(2,i0)
this%Pts1(3,i) = this%Pts1(3,i0)
this%Pts1(4,i) = this%Pts1(4,i0)
this%Pts1(5,i) = this%Pts1(5,i0)
this%Pts1(6,i) = this%Pts1(6,i0)
if((this%Pts1(1,i0)**2 + this%Pts1(3,i0)**2).ge.rrap) then
ilost = ilost + 1
else
endif
enddo
do i = this%Nptlocal - ilost + 1, this%Nptlocal
this%Pts1(1,i) = 0.0
this%Pts1(2,i) = 0.0
this%Pts1(3,i) = 0.0
this%Pts1(4,i) = 0.0
this%Pts1(5,i) = -1.0e15
this%Pts1(6,i) = 0.0
enddo
this%Nptlocal = this%Nptlocal - ilost
nplc = this%Nptlocal
call MPI_ALLREDUCE(nplc,nptot,1,MPI_INTEGER,&
MPI_SUM,MPI_COMM_WORLD,ierr)
this%Current = this%Current*nptot/this%Npt
this%Npt = nptot

end subroutine lostROUND_BeamBunch
!------------------------------------------------------------------

!//Calculate the space-charge E and B forces from point-to-point
!//summation including relativistic effects, and update the particle
Expand Down
3 changes: 2 additions & 1 deletion src/Appl/BeamLineElem.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -579,7 +579,8 @@ subroutine getradius_BeamLineElem(this,piperadius)
elseif(associated(this%psc)) then
call getparam_SC(this%psc,6,piperadius)
elseif(associated(this%pbpm)) then
call getparam_BPM(this%pbpm,2,piperadius)
piperadius = 1e8
!call getparam_BPM(this%pbpm,2,piperadius)
elseif(associated(this%pcf)) then
call getparam_ConstFoc(this%pcf,5,piperadius)
elseif(associated(this%pslrf)) then
Expand Down
118 changes: 96 additions & 22 deletions src/Contrl/AccSimulator.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -644,7 +644,7 @@ subroutine run_AccSimulator()
integer :: nsubstep,integerSamplePeriod
double precision :: zcent,distance,blnLength,dzz
integer, allocatable, dimension(:,:) :: idrfile
integer :: ibend,ibstart,isw,ibinit,ibendold,iifile,ii,ibinitold,idbeamln
integer :: ibend,ibstart,isw,ibinit,ibendold,iifile,ii,ibinitold,idbeamln,ibel
double precision :: zmax,t,dtless,zshift,gammazavg,curr
integer :: tmpflag,ib,ibb,ibunch,inib,nplctmp,nptmp,nptottmp
double precision, allocatable, dimension(:) :: gammaz
Expand All @@ -656,6 +656,8 @@ subroutine run_AccSimulator()
double precision :: zz,zorgin,zorgin2,vref,gamin,gam
double precision, dimension(6) :: ptref
integer :: idbd,idbend,flagbctmp
double precision :: radmin,piperadius
integer :: iflagdmshp
integer :: npttmplc,npttmp
double precision :: ztmp1,ztmp2,deltaz
integer :: ipt,iptnew
Expand Down Expand Up @@ -717,8 +719,8 @@ subroutine run_AccSimulator()
integer :: totnpts
real*8 :: qchg
!for collimator
integer :: icol
real*8, dimension(101) :: tcol,xradmin,xradmax,yradmin,yradmax
integer :: icol,iflagcol
real*8, dimension(101) :: tcol,xradmin,xradmax,yradmin,yradmax,flagcol
!for instant applying linear transfer matrix to the beam
real*8, dimension(101) :: tmap
integer :: imap
Expand Down Expand Up @@ -888,6 +890,7 @@ subroutine run_AccSimulator()
isteer = 0
icol = 0
irstart = 0
iflagdmshp = 0

!idrfile is used to store the <element type>, <external data file name>,
!and <id> for the internal data storage of each beamline element
Expand Down Expand Up @@ -1039,6 +1042,7 @@ subroutine run_AccSimulator()
call getparam_BeamLineElem(Blnelem(i),4,xradmax(icol))
call getparam_BeamLineElem(Blnelem(i),5,yradmin(icol))
call getparam_BeamLineElem(Blnelem(i),6,yradmax(icol))
call getparam_BeamLineElem(Blnelem(i),7,flagcol(icol))
endif

!transfer matrix information
Expand Down Expand Up @@ -1464,18 +1468,35 @@ subroutine run_AccSimulator()
endif
exit
endif

!check the particles outside the computational domain
do ib = 1, ibunch
call lost_BeamBunch(Ebunch(ib),xrad,yrad,Perdlen,zcent,&
nplctmp,nptmp)
Nplocal(ib) = nplctmp
Np(ib) = nptmp
!**
!check the particles outside the computational domain
!do ib = 1, ibunch
! call lostREC_BeamBunch(Ebunch(ib),xrad,yrad,Perdlen,zcent,&
! nplctmp,nptmp)
! Nplocal(ib) = nplctmp
! Np(ib) = nptmp
!! print*,"npt: ",ib,myid,nplctmp,nptmp
enddo
!enddo

! if(iflagdmshp.eq.0) then
! do ib = 1, ibunch
! call lostREC_BeamBunch(Ebunch(ib),xrad,yrad,Perdlen,zcent,&
! nplctmp,nptmp)
! Nplocal(ib) = nplctmp
! Np(ib) = nptmp
! enddo
! else if(iflagdmshp.eq.1) then
! do ib = 1, ibunch
! call lostDFO_BeamBunch(Ebunch(ib),xrad,yrad,Perdlen,zcent,&
! nplctmp,nptmp)
! Nplocal(ib) = nplctmp
! Np(ib) = nptmp
! enddo
! endif

totnpts = sum(Np)
if(totnpts.le.0) exit
!**

!find the beginning and end beam line elements that the effective
!bunch particles occupy
Expand Down Expand Up @@ -1604,6 +1625,69 @@ subroutine run_AccSimulator()
ibinitold = ibinit
ibendold = ibend

!-------------------------------------------------
!modified by C-H. Huang
radmin = xrad
!print*, 'ibinit',ibinit
!print*, 'ibend',ibend
do ibel = ibinit,ibend
call getradius_BeamLineElem(Blnelem(ibel),piperadius)
if(piperadius.le.1e-8) then
piperadius = xrad
endif
if(abs(piperadius).lt.radmin) then
radmin = abs(piperadius)
! write(97,*) radmin
endif
end do

!collimation
if(distance.le.tcol(icol+1) .and. (distance+dzz).ge.tcol(icol+1)) then
icol = icol + 1
iflagcol = flagcol(icol)+0.0001
if (iflagcol.le.10) then
do ib = 1, ibunch
call lostXY_BeamBunch(Ebunch(ib),xradmin(icol),xradmax(icol),&
yradmin(icol),yradmax(icol),nplctmp,nptmp)
Nplocal(ib) = nplctmp
Np(ib) = nptmp
iflagdmshp = 0
enddo
!else if (iflagcol.gt.10) then
else
do ib = 1, ibunch
call lostROUND_BeamBunch(Ebunch(ib),xradmin(icol),xradmax(icol),&
yradmin(icol),yradmax(icol),nplctmp,nptmp)
Nplocal(ib) = nplctmp
Np(ib) = nptmp
iflagdmshp = 1
enddo
!else
endif
endif

!check the particles outside the computational domain
if(iflagdmshp.eq.0) then
do ib = 1, ibunch
call lostREC_BeamBunch(Ebunch(ib),radmin,radmin,Perdlen,zcent,&
nplctmp,nptmp)
Nplocal(ib) = nplctmp
Np(ib) = nptmp
!print*, 'Nplocal', nplctmp
!print*, 'Np', nptmp
enddo
else if(iflagdmshp.eq.1) then
do ib = 1, ibunch
call lostDFO_BeamBunch(Ebunch(ib),radmin,radmin,Perdlen,zcent,&
nplctmp,nptmp)
Nplocal(ib) = nplctmp
Np(ib) = nptmp
!print*, 'Nplocal', nplctmp
!print*, 'Np', nptmp
enddo
endif
!-------------------------------------------------

if(idbend.ne.1) then !not bending magnet

if(zmax.le.0.0) goto 1000 !if no particles emittted pass field calcuation
Expand Down Expand Up @@ -2475,17 +2559,7 @@ subroutine run_AccSimulator()
enddo
endif

!collimation
if(distance.le.tcol(icol+1) .and. (distance+dzz).ge.tcol(icol+1)) then
icol = icol + 1
do ib = 1, ibunch
call lostXY_BeamBunch(Ebunch(ib),xradmin(icol),xradmax(icol),&
yradmin(icol),yradmax(icol),nplctmp,nptmp)
Nplocal(ib) = nplctmp
Np(ib) = nptmp
enddo
endif


!meager multiple bins into one bin
!if(t.le.tmger .and. (t+dtless*Dt).ge.tmger) then
if(distance.le.tmger .and. (distance+dzz).ge.tmger) then
Expand Down

0 comments on commit 3b228d9

Please sign in to comment.