From 3c925eb64a6608ce1cc74e729e897a658d746858 Mon Sep 17 00:00:00 2001 From: Eda Huang Date: Tue, 13 Aug 2024 14:37:03 -0700 Subject: [PATCH] check the radius of elements and collimate the transport beam profile into round or rectangular shape --- src/Appl/BeamBunch.f90 | 127 ++++++++++++++++++++++++++++++++++-- src/Appl/BeamLineElem.f90 | 3 +- src/Contrl/AccSimulator.f90 | 115 +++++++++++++++++++++++++------- 3 files changed, 216 insertions(+), 29 deletions(-) diff --git a/src/Appl/BeamBunch.f90 b/src/Appl/BeamBunch.f90 index 6b224be..b090e35 100644 --- a/src/Appl/BeamBunch.f90 +++ b/src/Appl/BeamBunch.f90 @@ -1123,8 +1123,7 @@ 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) + subroutine lostREC_BeamBunch(this,xrad,yrad,zleng,zcent,nplc,nptot) implicit none include 'mpif.h' type (BeamBunch), intent(inout) :: this @@ -1154,9 +1153,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 @@ -1189,7 +1188,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 @@ -1249,6 +1317,53 @@ 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 diff --git a/src/Appl/BeamLineElem.f90 b/src/Appl/BeamLineElem.f90 index e93e08e..3cef44c 100644 --- a/src/Appl/BeamLineElem.f90 +++ b/src/Appl/BeamLineElem.f90 @@ -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 diff --git a/src/Contrl/AccSimulator.f90 b/src/Contrl/AccSimulator.f90 index c64215c..85ae4ac 100644 --- a/src/Contrl/AccSimulator.f90 +++ b/src/Contrl/AccSimulator.f90 @@ -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 @@ -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 @@ -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 @@ -888,6 +890,7 @@ subroutine run_AccSimulator() isteer = 0 icol = 0 irstart = 0 + iflagdmshp = 0 !idrfile is used to store the , , !and for the internal data storage of each beamline element @@ -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 @@ -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 @@ -1604,6 +1625,66 @@ subroutine run_AccSimulator() ibinitold = ibinit ibendold = ibend + !* + 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) + 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 + 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 @@ -2475,17 +2556,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