Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow tracking of subset of particles filtered by energy range #16

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Allow tracking of subset of particles filtered by energy range #16

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
aafe992
Reinstate Doxygen-style comments
mattachu Oct 14, 2020
a572779
Allow calculation of distance and dzz using filtered energy
mattachu Mar 26, 2021
4781d70
Add code to read in filter parameters from the input file
mattachu Mar 26, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
118 changes: 116 additions & 2 deletions src/Appl/Ranger.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,9 @@
!----------------------------------------------------------------
module Rangerclass
use Pgrid2dclass
interface singlerange
module procedure singlerange_standard, singlerange_filtered
end interface
interface globalrange
module procedure globalrange1,globalrange2,globalrange3,globalrange4
end interface
Expand All @@ -21,7 +24,7 @@ module Rangerclass
!> find the range (minimum and maxium of a single bunch) and the centroid
!> of a single beam bunch.
!--------------------------------------------------------------------------------------
subroutine singlerange(partcls,nplc,nptot,range,center)
subroutine singlerange_standard(partcls,nplc,nptot,range,center)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The singlerange subroutine has not been modified apart from being renamed.
The new functionality is all in the new singlerange_filtered subroutine.

implicit none
include 'mpif.h'
integer, intent(in) :: nplc, nptot
Expand Down Expand Up @@ -79,7 +82,118 @@ subroutine singlerange(partcls,nplc,nptot,range,center)
enddo
endif

end subroutine singlerange
end subroutine singlerange_standard

!--------------------------------------------------------------------------------------
!> @brief
!> find the range (minimum and maxium of a single bunch) and the centroid
!> of a single beam bunch, filtering for a certain momentum range.
!--------------------------------------------------------------------------------------
subroutine singlerange_filtered(partcls,nplc,nptot,range,center,&
fmin, fmax, frange, fcenter, fnptot)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The new subroutine contains all the same code as the original subroutine, but with additional ranges and particle counts calculated for the particles filtered by energy range.

The results for singlerange_filtered(partcls,nplc,nptot,range,center,…) should match exactly the results for singlerange_standard(partcls,nplc,nptot,range,center), but also return the additional values fmin, fmax, frange, fcenter, fnptot.

implicit none
include 'mpif.h'
integer, intent(in) :: nplc, nptot
double precision, dimension(:,:) :: partcls
double precision, dimension(6), intent(out) :: range,center
double precision, intent(in) :: fmin, fmax
double precision, dimension(6), intent(out) :: frange, fcenter
integer, intent(out) :: fnptot
double precision, dimension(3) :: localrange1,localrange2,temp1,temp2
double precision, dimension(6) :: centerlc
double precision, dimension(3) :: frange_local_min, frange_local_max
double precision, dimension(6) :: fcenter_local
integer :: fnptot_local
integer :: i, j, ierr

do i = 1, 3
localrange1(i) = 1000000.0
localrange2(i) = -1000000.0
frange_local_min(i) = 1000000.0
frange_local_max(i) = -1000000.0
centerlc(2*i-1) = 0.0
centerlc(2*i) = 0.0
fcenter_local(2*i-1) = 0.0
fcenter_local(2*i) = 0.0
enddo
fnptot_local = 0

! print*,"nplc, partcls: ",nplc,partcls
do i = 1, nplc
do j = 1, 3
if(localrange1(j) > partcls(2*j-1,i)) then
localrange1(j) = partcls(2*j-1,i)
endif
if(localrange2(j) < partcls(2*j-1,i)) then
localrange2(j) = partcls(2*j-1,i)
endif
centerlc(2*j-1) = centerlc(2*j-1) + partcls(2*j-1,i)
centerlc(2*j) = centerlc(2*j) + partcls(2*j,i)
enddo
if(partcls(6,i) >= fmin .and. partcls(6,i) <= fmax) then
fnptot_local = fnptot_local + 1
do j = 1, 3
if(partcls(2*j-1,i) < frange_local_min(j)) then
frange_local_min(j) = partcls(2*j-1,i)
endif
if(partcls(2*j-1,i) > frange_local_max(j)) then
frange_local_max(j) = partcls(2*j-1,i)
endif
fcenter_local(2*j-1) = fcenter_local(2*j-1) + partcls(2*j-1,i)
fcenter_local(2*j) = fcenter_local(2*j) + partcls(2*j,i)
enddo
endif
enddo

! print*,"localrange1: ",localrange1
! print*,"localrange2: ",localrange1
! print*,"centerlc: ",centerlc

call MPI_ALLREDUCE(localrange1,temp1,3,MPI_DOUBLE_PRECISION,MPI_MIN,MPI_COMM_WORLD,&
ierr);
call MPI_ALLREDUCE(localrange2,temp2,3,MPI_DOUBLE_PRECISION,MPI_MAX,MPI_COMM_WORLD,&
ierr);
call MPI_ALLREDUCE(centerlc,center,6,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,&
ierr);

if(nptot.gt.0) then
do i = 1, 3
range(2*i-1) = temp1(i)
range(2*i) = temp2(i)
center(2*i-1) = center(2*i-1)/nptot
center(2*i) = center(2*i)/nptot
enddo
else
do i = 1, 3
range(2*i-1) = temp1(i)
range(2*i) = temp2(i)
center(2*i-1) = center(2*i-1)
center(2*i) = center(2*i)
enddo
endif

call MPI_ALLREDUCE(frange_local_min, temp1, 3, MPI_DOUBLE_PRECISION, &
MPI_MIN, MPI_COMM_WORLD, ierr);
call MPI_ALLREDUCE(frange_local_max, temp2, 3, MPI_DOUBLE_PRECISION, &
MPI_MAX, MPI_COMM_WORLD, ierr);
call MPI_ALLREDUCE(fcenter_local, fcenter, 6, MPI_DOUBLE_PRECISION, &
MPI_SUM, MPI_COMM_WORLD, ierr);
call MPI_ALLREDUCE(fnptot_local, fnptot, 1, MPI_INTEGER, &
MPI_SUM, MPI_COMM_WORLD, ierr);

if(fnptot == 0) then
frange = 0.0d0
fcenter = 0.0d0
else
do i = 1, 3
frange(2*i-1) = temp1(i)
frange(2*i) = temp2(i)
fcenter(2*i-1) = fcenter(2*i-1)/fnptot
fcenter(2*i) = fcenter(2*i)/fnptot
enddo
end if

end subroutine singlerange_filtered

!--------------------------------------------------------------------------------------
!> @brief
Expand Down
45 changes: 41 additions & 4 deletions src/Contrl/AccSimulator.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,9 +63,12 @@ module AccSimulatorclass
!! \# of beam elems, type of integrator.
!! FlagImage: switch flag for image space-charge force calculation: "1" for yes,
!! otherwise for no.
!! FlagEnergyRange: calculate distances for a specific subset of particles
!! based on energy range rather than for all particles
!> @{
integer :: Nx,Ny,Nz,Nxlocal,Nylocal,Nzlocal,Flagbc,&
Nblem,Flagmap,Flagdiag,FlagImage
integer :: FlagEnergyRange = 0
!> @}

!> @name
Expand All @@ -85,6 +88,12 @@ module AccSimulatorclass
Perdlen,dt,xrad,yrad
!> @}

!> @name
!! Energy range minimum and maximum for filtering large energy spread
!> @{
double precision :: filter_min, filter_max
!> @}

!> @name
!! conts. in init. dist.
!> @{
Expand Down Expand Up @@ -197,6 +206,11 @@ subroutine init_AccSimulator(time)
Flagsubstep,phsini,dt,ntstep,Nbunch,FlagImage,Nemission,&
temission,zimage)

! Hard-code energy range for testing - should be included in input file
FlagEnergyRange = 1
filter_min = 95.0e6 ! eV
filter_max = 105.0e6 ! eV

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The reading of values from input file is in the next commit.

! print*,"Np: ",Np,dt,ntstep,Nx,Ny,Nz
! print*,"Bcurr: ",Bcurr,Bkenergy,Bmass,Bcharge,Bfreq
!-------------------------------------------------------------------
Expand Down Expand Up @@ -679,6 +693,11 @@ subroutine run_AccSimulator()
integer :: tmpflag,ib,ibb,ibunch,inib,nplctmp,nptmp,nptottmp
double precision, allocatable, dimension(:) :: gammaz
double precision, allocatable, dimension(:,:) :: brange
double precision :: fmin, fmax
double precision, allocatable, dimension(:,:) :: frange
double precision, dimension(6) :: fgrange
double precision :: fgammaz, fzcent
integer, allocatable, dimension(:) :: fnp
double precision :: dGspread
integer, dimension(Maxoverlap) :: tmpfile
double precision :: tmpcur,totchrg,r0
Expand Down Expand Up @@ -1149,6 +1168,8 @@ subroutine run_AccSimulator()

allocate(gammaz(Nbunch))
allocate(brange(12,Nbunch))
allocate(frange(12,Nbunch))
allocate(fnp(Nbunch))
!count the total current and # of particles and local # of particles for each
!bunch or bin
curr = 0.0d0
Expand Down Expand Up @@ -1431,8 +1452,16 @@ subroutine run_AccSimulator()
!only the bunch with zmax>0 is counted as an effective bunch
do ib = 1, ibunch
!//find the range and center of each bunch/bin
call singlerange(Ebunch(ib)%Pts1,Nplocal(ib),Np(ib),&
ptrange,sgcenter)
if(FlagEnergyRange == 1) then
fmin = sqrt((1.0d0 + filter_min/Ebunch(ib)%mass)**2 - 1.0d0)
fmax = sqrt((1.0d0 + filter_max/Ebunch(ib)%mass)**2 - 1.0d0)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This converts the energy range to a momentum range.

call singlerange(Ebunch(ib)%Pts1, Nplocal(ib), Np(ib), &
ptrange, sgcenter, fmin, fmax, &
frange(1:6,ib), frange(7:12,ib), fnp(ib))
else
call singlerange(Ebunch(ib)%Pts1,Nplocal(ib),Np(ib),&
ptrange,sgcenter)
endif
!Ebunch(ib)%refptcl(5) = sgcenter(5) + zorgin
Ebunch(ib)%refptcl(5) = sgcenter(5)
gammaz(ib) = sqrt(1.0+sgcenter(6)**2)!//gammaz from <gamma_i betaz_i>
Expand Down Expand Up @@ -1468,8 +1497,14 @@ subroutine run_AccSimulator()

! print*,"gamz: ",gammaz(1)
!get the distance of the center of all effective bunches/bins
distance = zcent*Scxlt
dzz = sqrt(1.0d0-1.0d0/gammazavg**2)*Clight*dtless*Dt
if(FlagEnergyRange == 1) then
call globalrange(frange, fgrange, fgammaz, fzcent, fnp, ibunch)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Main values grange, gammazavg, zcent are not modified, here we use globalrange() to calculate new filtered values just for the purpose of setting distance and zzz

distance = fzcent*Scxlt
dzz = sqrt(1.0d0-1.0d0/fgammaz**2)*Clight*dtless*Dt
else
distance = zcent*Scxlt
dzz = sqrt(1.0d0-1.0d0/gammazavg**2)*Clight*dtless*Dt
endif
nptottmp = sum(Np)
!exit if the beam is outside the beamline
if(distance.gt.blnLength .or. distance.gt.tstop .or. nptottmp.lt.1) then
Expand Down Expand Up @@ -2560,6 +2595,8 @@ subroutine run_AccSimulator()
deallocate(idrfile)
deallocate(gammaz)
deallocate(brange)
deallocate(frange)
deallocate(fnp)
if(Flagbc.eq.3) then
deallocate(besscoef)
deallocate(bessnorm)
Expand Down