Homec4science
Diffusion FENNECS: Finite Element Non-Neutral Electron Cloud Simulator (geom_testpardiso) Browse

fennecs/src/beam_mod.f90
f331bae49c3ageom_testpardiso

Last Change
Raw File

beam_mod.f90
View Options

MODULE beam
!
USE constants
USE mpi
USE mpihelper
USE basic, ONLY: mpirank, mpisize
IMPLICIT NONE
TYPE particles
INTEGER :: Nptot
INTEGER, DIMENSION(:), ALLOCATABLE :: Rindex, Zindex ! Indices of the particles on the grid
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: &
& R,Z, &
& WZ, WR, &
& pot, Er, Ez
DOUBLE PRECISION, DIMENSION(:),POINTER:: &
& UR, URold, &
& UTHET, UTHETold, &
& UZ, UZold, &
& Gamma, Gammaold
END TYPE particles
!
TYPE(particles) :: parts ! Storage for all the particles
SAVE :: parts
TYPE(particle), DIMENSION(:), ALLOCATABLE:: rrecvpartbuff, lrecvpartbuff, rsendpartbuff, lsendpartbuff ! buffer to send and receive particle from left and right processes
! Diagnostics (scalars)
DOUBLE PRECISION :: ekin ! Total kinetic energz (J)
DOUBLE PRECISION :: epot ! Total potential energy (J)
DOUBLE PRECISION :: etot, etot0 ! Current and Initial total energy (J)
INTEGER, DIMENSION(3), SAVE :: ireducerequest=MPI_REQUEST_NULL
!
INTEGER, DIMENSION(:), ALLOCATABLE :: Nplocs_all !F.B. get local numbers of particles
LOGICAL:: collected=.false. !Stores if the particles data have been collected to root process this timestep
!
CONTAINS
SUBROUTINE creat_parts(p, nparts)
TYPE(particles) :: p
INTEGER, INTENT(in) :: nparts
IF (.NOT. ALLOCATED(p%Z) ) THEN
p%Nptot = nparts
ALLOCATE(p%Z(nparts))
ALLOCATE(p%R(nparts))
ALLOCATE(p%WZ(nparts))
ALLOCATE(p%WR(nparts))
ALLOCATE(p%UR(nparts))
ALLOCATE(p%UZ(nparts))
ALLOCATE(p%UTHET(nparts))
ALLOCATE(p%URold(nparts))
ALLOCATE(p%UZold(nparts))
ALLOCATE(p%UTHETold(nparts))
ALLOCATE(p%Gamma(nparts))
ALLOCATE(p%Rindex(nparts))
ALLOCATE(p%Zindex(nparts))
ALLOCATE(p%pot(nparts))
ALLOCATE(p%Er(nparts))
ALLOCATE(p%Ez(nparts))
ALLOCATE(p%GAMMAold(nparts))
parts%URold=0
parts%UZold=0
parts%UTHETold=0
parts%rindex=0
parts%zindex=0
parts%WR=0
parts%WZ=0
parts%UR=0
parts%UZ=0
parts%UTHET=0
parts%Z=0
parts%R=0
parts%Gamma=1
parts%Er=0
parts%Ez=0
parts%pot=0
parts%gammaold=1
END IF
END SUBROUTINE creat_parts
!______________________________________________________________________________
SUBROUTINE load_parts
! Loads the particles at the beginning of the simulation and create the parts variable if necessary
USE basic, ONLY: nplasma, mpirank, ierr, distribtype, mpisize, nlclassical
USE mpi
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET
ALLOCATE(VZ(nplasma), VR(nplasma), VTHET(nplasma))
! Select case to define the type of distribution
SELECT CASE(distribtype)
CASE(1) ! Gaussian distribution in V and uniform in R
CALL loaduniformRZ(VR, VZ, VTHET)
CASE(2) !Stable distribution from Davidson 4.95 p.119
CALL loadDavidson(VR, VZ, VTHET, lodunir)
CASE(3) !Stable distribution from Davidson 4.95 p.119 but with distribution in R as 1/R**2
CALL loadDavidson(VR, VZ, VTHET, lodinvr)
CASE DEFAULT
IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of distribution:", distribtype
CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
END SELECT
IF(nlclassical) THEN
parts%Gamma=1
ELSE
parts%Gamma=sqrt(1/(1-VR**2+VZ**2+VTHET**2))
END IF
parts%UR=parts%Gamma*VR
parts%UZ=parts%Gamma*VZ
parts%UTHET=parts%Gamma*VTHET
DEALLOCATE(VZ, VR, VTHET)
CALL localisation
IF(mpisize .gt. 1) THEN
CALL distribpartsonprocs
END IF
END SUBROUTINE load_parts
!________________________________________________________________________________
SUBROUTINE bound
USE basic, ONLY: zgrid, nz, Zbounds, cstep, mpirank
INTEGER :: i, rsendnbparts=0, lsendnbparts=0
INTEGER, DIMENSION(parts%Nptot) :: sendhole
sendhole=0
lsendnbparts=0
rsendnbparts=0
IF (parts%Nptot .NE. 0) THEN
! Boundary condition at z direction
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
DO i=1,parts%Nptot
IF(cstep .ne. 0) THEN
IF (parts%Z(i) .ge. Zbounds(mpirank+1)) THEN
!$OMP CRITICAL (nbparts)
rsendnbparts=rsendnbparts+1
sendhole(lsendnbparts+rsendnbparts)=i
!$OMP END CRITICAL (nbparts)
ELSE IF (parts%Z(i) .lt. Zbounds(mpirank)) THEN
!$OMP CRITICAL (nbparts)
lsendnbparts=lsendnbparts+1
sendhole(lsendnbparts+rsendnbparts)=-i
!$OMP END CRITICAL (nbparts)
END IF
END IF
DO WHILE (parts%Z(i) .GT. zgrid(nz))
parts%Z(i) = parts%Z(i) - zgrid(nz) + zgrid(0)
END DO
DO WHILE (parts%Z(i) .LT. zgrid(0))
parts%Z(i) = parts%Z(i) + zgrid(nz) - zgrid(0)
END DO
END DO
!$OMP END PARALLEL DO
CALL particlescommunication(lsendnbparts, rsendnbparts, sendhole)
END IF
END subroutine bound
!________________________________________________________________________________
SUBROUTINE localisation
USE basic, ONLY: zgrid, rgrid, dz, nr, nnr, dr, rnorm
INTEGER :: j,i, nblostparts=0
INTEGER, DIMENSION(parts%Nptot) :: losthole
i=1
losthole=0
IF (parts%Nptot .NE. 0) THEN
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
DO i=1,parts%Nptot
parts%zindex(i)=(parts%Z(i)-zgrid(0))/dz
IF (parts%R(i) .GT. rgrid(0) .AND. parts%R(i) .LT. rgrid(nnr(1))) THEN
parts%rindex(i)=(parts%R(i)-rgrid(0))/dr(1)
ELSE IF(parts%R(i) .GT. rgrid(nnr(1)) .AND. parts%R(i) .LT. rgrid(nnr(1)+nnr(2))) THEN
parts%rindex(i)=(parts%R(i)-rgrid(nnr(1)))/dr(2)+nnr(1)
ELSE IF(parts%R(i) .GT. rgrid(nnr(1)+nnr(2)) .AND. parts%R(i) .LT. rgrid(nr)) THEN
parts%rindex(i)=(parts%R(i)-rgrid(nnr(1)+nnr(2)))/dr(3)+nnr(1)+nnr(2)
ELSE
WRITE(*,*) "Particle:",i , " of process:", mpirank, "is out of bound in r:", parts%R(i)*rnorm, rgrid(0)*rnorm, rgrid(nr)*rnorm
!nlend=.TRUE.
WRITE(*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Particle removed !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
!$OMP CRITICAL (lostparts)
nblostparts=nblostparts+1
losthole(nblostparts)=i
!$OMP END CRITICAL (lostparts)
END IF
END DO
!$OMP END PARALLEL DO
IF(nblostparts.gt.0) THEN
DO j=1,nblostparts
CALL move_part(parts%Nptot, losthole(j))
parts%Nptot=parts%Nptot-1
END DO
END IF
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
DO i=1,parts%Nptot
parts%WZ(i)=(parts%Z(i)-zgrid(parts%zindex(i)))/(zgrid(parts%zindex(i)+1)-zgrid(parts%zindex(i)));
parts%WR(i)=(parts%R(i)-rgrid(parts%rindex(i)))/(rgrid(parts%rindex(i)+1)-rgrid(parts%rindex(i)));
END DO
!$OMP END PARALLEL DO
END IF
END SUBROUTINE localisation
!________________________________________________________________________________________________________
SUBROUTINE comp_velocity
!
! Computes the new velocity of the particles due to Lorentz force
!
USE basic, ONLY : nlclassical
! Store old Velocities
CALL swappointer(parts%UZold, parts%UZ)
CALL swappointer(parts%URold, parts%UR)
CALL swappointer(parts%UTHETold, parts%UTHET)
CALL swappointer(parts%Gammaold, parts%Gamma)
IF (nlclassical) THEN
CALL comp_velocity_fun(gamma_classical)
ELSE
CALL comp_velocity_fun(gamma_relativistic)
END IF
END SUBROUTINE comp_velocity
!____________________________________________
SUBROUTINE comp_velocity_fun(gamma)
!
! Computes the new velocity of the particles due to Lorentz force
!
USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nz
interface
subroutine gamma(gam, UZ, UR, UTHET)
DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET
DOUBLE PRECISION, INTENT(OUT):: gam
end subroutine
end interface
DOUBLE PRECISION :: tau
DOUBLE PRECISION:: BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2
INTEGER:: J1, J2, J3, J4
INTEGER:: i
! Normalized time increment
tau=omegac/2/omegap*dt/tnorm
IF (parts%Nptot .NE. 0) THEN
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2)
DO i=1,parts%Nptot
J1=(parts%rindex(i))*(nz+1)+parts%zindex(i)+1
J2=J1+1
J3=(parts%rindex(i)+1)*(nz+1)+parts%zindex(i)+1
J4=J3+1
! Interpolation for magnetic field
BRZ=(1-parts%WZ(i))*(1-parts%WR(i))*Bz(J4) &
& +parts%WZ(i)*(1-parts%WR(i))*Bz(J3) &
& +(1-parts%WZ(i))*parts%WR(i)*Bz(J2) &
& +parts%WZ(i)*parts%WR(i)*Bz(J1)
BRR=(1-parts%WZ(i))*(1-parts%WR(i))*Br(J4) &
& +parts%WZ(i)*(1-parts%WR(i))*Br(J3) &
& +(1-parts%WZ(i))*parts%WR(i)*Br(J2) &
& +parts%WZ(i)*parts%WR(i)*Br(J1)
! First half of electric pulse
parts%UZ(i)=parts%UZold(i)+parts%Ez(i)*tau
parts%UR(i)=parts%URold(i)+parts%ER(i)*tau
CALL gamma(parts%Gamma(i), parts%UZ(i), parts%UR(i), parts%UTHETold(i))
! Rotation along magnetic field
ZBZ=tau*BRZ/parts%Gamma(i)
ZBR=tau*BRR/parts%Gamma(i)
ZPZ=parts%UZ(i)-ZBR*parts%UTHETold(i) !u'_{z}
ZPR=parts%UR(i)+ZBZ*parts%UTHETold(i) !u'_{r}
ZPTHET=parts%UTHETold(i)+(ZBR*parts%UZ(i)-ZBZ*parts%UR(i)) !u'_{theta}
SQR=1+ZBZ*ZBZ+ZBR*ZBR
ZBZ2=2*ZBZ/SQR
ZBR2=2*ZBR/SQR
parts%UZ(i)=parts%UZ(i)-ZBR2*ZPTHET !u+_{z}
parts%UR(i)=parts%UR(i)+ZBZ2*ZPTHET !u+_{r}
parts%UTHET(i)=parts%UTHETold(i)+(ZBR2*ZPZ-ZBZ2*ZPR) !u+_{theta}
! Second half of acceleration
parts%UZ(i)=parts%UZ(i)+parts%EZ(i)*tau
parts%UR(i)=parts%UR(i)+parts%ER(i)*tau
CALL gamma(parts%Gamma(i), parts%UZ(i), parts%UR(i), parts%UTHETold(i))
END DO
!$OMP END PARALLEL DO
END IF
collected=.false.
END SUBROUTINE comp_velocity_fun
!______________________________________________________________________________
SUBROUTINE gamma_classical(gamma, UZ, UR, UTHET)
DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET
DOUBLE PRECISION, INTENT(OUT):: gamma
gamma=1.0
END SUBROUTINE gamma_classical
!______________________________________________________________________________
SUBROUTINE gamma_relativistic(gamma, UZ, UR, UTHET)
DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET
DOUBLE PRECISION, INTENT(OUT):: gamma
gamma=sqrt(1+UZ**2+UR**2+UTHET**2)
END SUBROUTINE
!________________________________________________________________________________
SUBROUTINE push
Use basic, ONLY: dt, tnorm
DOUBLE PRECISION:: XP, YP, COSA, SINA, U1, U2
INTEGER :: i
IF (parts%Nptot .NE. 0) THEN
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i, XP, YP, COSA, SINA, U1, U2)
DO i=1,parts%Nptot
! Local Cartesian coordinates
XP=parts%R(i)+dt/tnorm*parts%UR(i)/parts%Gamma(i)
YP=dt/tnorm*parts%UTHET(i)/parts%Gamma(i)
! Conversion to cylindrical coordiantes
parts%Z(i)=parts%Z(i)+dt/tnorm*parts%UZ(i)/parts%Gamma(i)
parts%R(i)=sqrt(XP**2+YP**2)
! Velocity in rotated reference frame
IF (parts%R(i) .EQ. 0) THEN
COSA=1
SINA=0
ELSE
COSA=XP/parts%R(i)
SINA=YP/parts%R(i)
END IF
U1=COSA*parts%UR(i)+SINA*parts%UTHET(i)
U2=-SINA*parts%UR(i)+COSA*parts%UTHET(i)
parts%UR(i)=U1
parts%UTHET(i)=U2
END DO
!$OMP END PARALLEL DO
END IF
collected=.false.
END SUBROUTINE push
!_______________________________________________________________________________
!_______________________________________________________________________________
SUBROUTINE diagnostics
!
! Compute energies
!
!!$ CALL energies
USE constants, ONLY: vlight
USE basic, ONLY: qsim, phinorm, msim, cstep, nlclassical, nz, partdensity, ierr, step, it2d, nlend, nr, itparts, nplasma
DOUBLE PRECISION :: gammamid
INTEGER :: i
INTEGER, DIMENSION(:) :: stat(3*MPI_STATUS_SIZE)
CALL MPI_WAIT(ireducerequest(3), stat(2*MPI_STATUS_SIZE+1:3*MPI_STATUS_SIZE), ierr)
partdensity=0
ekin=0
epot=0
etot=0
!$OMP PARALLEL DO REDUCTION(+:epot, ekin, partdensity) DEFAULT(SHARED) PRIVATE(i,gammamid)
DO i=1,parts%Nptot
!
epot=epot+qsim*parts%pot(i)*phinorm
IF(.not. nlclassical) THEN
gammamid=1/sqrt(1-(abs(parts%URold(i)*parts%UR(i))+abs(parts%UZold(i)*parts%UZ(i))+abs(parts%UTHETold(i)*parts%UTHET(i)))/(parts%Gammaold(i)*parts%Gamma(i)))
ekin=ekin+msim*vlight**2*(gammamid-1)
ELSE
ekin=ekin+0.5*msim*vlight**2*(abs(parts%URold(i)*parts%UR(i))+abs(parts%UZold(i)*parts%UZ(i))+abs(parts%UTHETold(i)*parts%UTHET(i)))
END IF
IF(modulo(step,it2d).eq. 0 .or. nlend) THEN
partdensity(parts%zindex(i)+1+parts%rindex(i)*nz)=partdensity(parts%zindex(i)+1+parts%rindex(i)*nz)+1
END IF
END DO
!$OMP END PARALLEL DO
IF(mpirank .eq.0 ) THEN
CALL MPI_IREDUCE(MPI_IN_PLACE, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(1), ierr)
CALL MPI_IREDUCE(MPI_IN_PLACE, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(2), ierr)
ELSE
CALL MPI_IREDUCE(epot, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(1), ierr)
CALL MPI_IREDUCE(ekin, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(2), ierr)
END IF
! Stores the particle density on the grid
IF(modulo(step,it2d).eq. 0 .or. nlend) THEN
IF(mpirank .eq.0 ) THEN
CALL MPI_IREDUCE(MPI_IN_PLACE, partdensity, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(3), ierr)
ELSE
CALL MPI_IREDUCE(partdensity, partdensity, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(3), ierr)
END IF
END IF
IF(mpisize .gt. 1) THEN
Nplocs_all(mpirank)=parts%Nptot
IF(mpirank .eq.0 ) THEN
CALL MPI_gather(MPI_IN_PLACE, 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,&
& 0, MPI_COMM_WORLD, ierr)
ELSE
CALL MPI_gather(Nplocs_all(mpirank), 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,&
& 0, MPI_COMM_WORLD, ierr)
END IF
IF(mpirank .eq.0 .and. INT(SUM(Nplocs_all)).ne. nplasma) THEN
WRITE(*,*) "Error particle lost on some processus"
!CALL MPI_abort(MPI_COMM_WORLD,-1,ierr)
END IF
END IF
CALL MPI_WAIT(ireducerequest(1), stat(1:MPI_STATUS_SIZE), ierr)
CALL MPI_WAIT(ireducerequest(2), stat(MPI_STATUS_SIZE+1:2*MPI_STATUS_SIZE), ierr)
etot=epot+ekin
!
! Shift to Etot at cstep=1 (not valable yet at cstep=0!)
IF(cstep.EQ.1) THEN
etot0 = etot
END IF
IF(mpisize .gt. 1 .and. (modulo(step,itparts) .eq. 0 .or. nlend)) THEN
CALL collectparts
END IF
end subroutine diagnostics
SUBROUTINE collectparts
USE basic, ONLY: mpirank, mpisize, ierr
INTEGER, DIMENSION(:), ALLOCATABLE :: displs
INTEGER:: i
IF(collected) RETURN ! exit subroutine if particles have already been collected this time step
IF(mpirank .ne. 0) THEN
ALLOCATE(displs(0:mpisize-1))
displs=0
! Send Particles informations to root process
CALL MPI_Gatherv(parts%Z, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%Z, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(parts%R, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%R, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(parts%UR, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UR, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(parts%UZ, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UZ, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(parts%UTHET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UTHET, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(parts%pot, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%pot, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(parts%Rindex, Nplocs_all(mpirank), MPI_INTEGER, parts%Rindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(parts%Zindex, Nplocs_all(mpirank), MPI_INTEGER, parts%Zindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
ELSE
ALLOCATE(displs(0:mpisize-1))
displs(0)=0
Do i=1,mpisize-1
displs(i)=displs(i-1)+Nplocs_all(i-1)
END DO
! Receive particle information from all processes
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%Z, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%R, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UR, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UZ, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UTHET, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%pot, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%Rindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%Zindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
END IF
collected=.TRUE.
END SUBROUTINE collectparts
!_______________________________________________________________________________
SUBROUTINE adapt_vinit
!
! Computes the velocity at time -dt/2 from velocities computed at time 0
!
USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nlclassical, phinorm, nz, distribtype, H0, vnorm
DOUBLE PRECISION :: tau, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, &
& SQR, ZBZ2, ZBR2, Vperp, v2
INTEGER :: J1, J2, J3, J4, i
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET
IF(distribtype .EQ. 2 .OR. distribtype .EQ. 3) THEN
ALLOCATE(VR(parts%Nptot),VZ(parts%Nptot),VTHET(parts%Nptot))
CALL loduni(7,VZ)
VZ=VZ*2*pi
VTHET=parts%UTHET/parts%Gamma*vnorm
DO i=1,parts%Nptot
Vperp=sqrt(MAX(2*H0/me+2*elchar/me*parts%pot(i)*phinorm-VTHET(i)**2,0.0_db))
VR(i)=Vperp*sin(VZ(i))
VZ(i)=Vperp*cos(VZ(i))
IF(nlclassical) THEN
parts%Gamma(i)=1
ELSE
v2=VR(i)**2+VZ(i)**2+VTHET(i)**2
parts%Gamma(i)=sqrt(1/(1-v2))
END IF
parts%UR(i)=parts%Gamma(i)*VR(i)/vnorm
parts%UZ(i)=parts%Gamma(i)*VZ(i)/vnorm
parts%UTHET(i)=parts%Gamma(i)*VTHET(i)/vnorm
END DO
DEALLOCATE(VR,VZ,VTHET)
END IF
! Normalised time increment
tau=-omegac/2/omegap*dt/tnorm
! Store old Velocities
CALL swappointer(parts%UZold, parts%UZ)
CALL swappointer(parts%URold, parts%UR)
CALL swappointer(parts%UTHETold, parts%UTHET)
CALL swappointer(parts%Gammaold, parts%Gamma)
IF (parts%Nptot .NE. 0) THEN
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2)
DO i=1,parts%Nptot
J1=(parts%rindex(i))*(nz+1)+parts%zindex(i)+1
J2=J1+1
J3=(parts%rindex(i)+1)*(nz+1)+parts%zindex(i)+1
J4=J3+1
! Interpolation for magnetic field
BRZ=(1-parts%WZ(i))*(1-parts%WR(i))*Bz(J4) &
& +parts%WZ(i)*(1-parts%WR(i))*Bz(J3) &
& +(1-parts%WZ(i))*parts%WR(i)*Bz(J2) &
& +parts%WZ(i)*parts%WR(i)*Bz(J1)
BRR=(1-parts%WZ(i))*(1-parts%WR(i))*Br(J4) &
& +parts%WZ(i)*(1-parts%WR(i))*Br(J3) &
& +(1-parts%WZ(i))*parts%WR(i)*Br(J2) &
& +parts%WZ(i)*parts%WR(i)*Br(J1)
! Half inverse Rotation along magnetic field
ZBZ=tau*BRZ/parts%Gammaold(i)
ZBR=tau*BRR/parts%Gammaold(i)
SQR=1+ZBZ*ZBZ+ZBR*ZBR
ZPZ=(parts%UZold(i)-ZBR*parts%UTHETold(i))/SQR !u-_{z}
ZPR=(parts%URold(i)+ZBZ*parts%UTHETold(i))/SQR !u-_{r}
ZPTHET=parts%UTHETold(i)+(ZBR*parts%UZold(i)-ZBZ*parts%URold(i))/SQR !u-_{theta}
parts%UZ(i)=ZPZ
parts%UR(i)=ZPR
parts%UTHET(i)=ZPTHET
! half of decceleration
parts%UZ(i)=parts%UZ(i)+parts%Ez(i)*tau
parts%UR(i)=parts%UR(i)+parts%Er(i)*tau
IF(.not. nlclassical) THEN
parts%Gamma(i)=sqrt(1+parts%UZ(i)**2+parts%UR(i)**2+parts%UTHET(i)**2)
END IF
END DO
!$OMP END PARALLEL DO
END IF
END SUBROUTINE adapt_vinit
!_______________________________________________________________________________
SUBROUTINE distribpartsonprocs
! Computes the start and end indices for the Z boundaries on local processus
! Computes the particle indices from initial particle loading vector, that stay in current process
USE basic, ONLY: nz, nplasma, ierr, Zbounds
INTEGER:: k
DOUBLE PRECISION:: idealnbpartsperproc
INTEGER:: totparts ! Total number of particle from zgrid(0) to zgrid(k)
INTEGER:: partindexstart ! Starting index of local particle in the parts variable used later to move the local particles at the begining of the vector
INTEGER:: partindexlast ! Ending index of local particle in the parts variable used later to move the local particles at the begining of the vector
INTEGER, DIMENSION(0:nz):: partspercol ! Vector containing the number of particles between zgrid(n) and zgrid(n+1)
INTEGER:: Zmin, Zmax ! Minimum and maximum indices of particles in Z direction
INTEGER:: Zperproc
partspercol=0
idealnbpartsperproc = FLOOR(REAL(nplasma)/REAL(mpisize))
Zmin=MINVAL(parts%Zindex)
Zmax=MAXVAL(parts%Zindex)
Zperproc=(Zmax-Zmin)/mpisize
partindexstart=1
totparts=0
partindexlast=parts%Nptot
DO k=0,mpisize-2
Nplocs_all(k)=idealnbpartsperproc
END DO
NPlocs_all(mpisize-1)=idealnbpartsperproc+MODULO(nplasma,mpisize)
DO k=1, mpisize-1
Zbounds(k)=parts%Z(INT(k*idealnbpartsperproc))
END DO
! Store local end of particle vector
partindexlast=SUM(Nplocs_all(0:mpirank))
! Store local start of particle vector
IF(mpirank .ne. 0) partindexstart=SUM(Nplocs_all(0:mpirank-1))+1
IF(mpirank .ne. 0) THEN
DO k= partindexstart, partindexlast
CALL move_part(k,k-partindexstart+1)
END DO
END IF
parts%Nptot=Nplocs_all(mpirank)
WRITE(*,*) mpirank, " Zbounds: ", Zbounds(mpirank), Zbounds(mpirank+1), " indices ", partindexstart, partindexlast, " nptot", parts%Nptot
IF(INT(SUM(Nplocs_all)) .ne. nplasma) THEN
WRITE(*,*) "Error in particle counting on proc:", mpirank, " local sum:", INT(SUM(Nplocs_all)), " nplasma:", nplasma
CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
END IF
END SUBROUTINE distribpartsonprocs
!_______________________________________________________________________________
SUBROUTINE particlescommunication(lsendnbparts, rsendnbparts, sendholes)
USE mpihelper, ONLY: particle_type
USE basic, ONLY: rightproc, leftproc, ierr
INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts
INTEGER, INTENT(in) :: sendholes(parts%Nptot)
INTEGER, ASYNCHRONOUS :: rrecvnbparts=0, lrecvnbparts=0
INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendrequest, recvrequest
INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendstatus(2*MPI_STATUS_SIZE), recvstatus(2*MPI_STATUS_SIZE)
INTEGER :: lsentnbparts, rsentnbparts
INTEGER :: lreceivednbparts, rreceivednbparts
lsentnbparts=lsendnbparts
rsentnbparts=rsendnbparts
sendrequest=MPI_REQUEST_NULL
recvrequest=MPI_REQUEST_NULL
CALL MPI_IRECV(lrecvnbparts, 1, MPI_INTEGER, leftproc, 0, MPI_COMM_WORLD, recvrequest(1), ierr)
CALL MPI_IRECV(rrecvnbparts, 1, MPI_INTEGER, rightproc, 0, MPI_COMM_WORLD, recvrequest(2), ierr)
CALL MPI_ISEND(lsentnbparts, 1, MPI_INTEGER, leftproc, 0, MPI_COMM_WORLD, sendrequest(1), ierr)
CALL MPI_ISEND(rsentnbparts, 1, MPI_INTEGER, rightproc, 0, MPI_COMM_WORLD, sendrequest(2), ierr)
CALL MPI_Wait(recvrequest(1), recvstatus(1), ierr)
CALL MPI_Wait(recvrequest(2), recvstatus(2), ierr)
recvrequest=MPI_REQUEST_NULL
lreceivednbparts=lrecvnbparts
rreceivednbparts=rrecvnbparts
IF( ALLOCATED(rrecvpartbuff) .AND. (size(rrecvpartbuff) .lt. rrecvnbparts) ) DEALLOCATE(rrecvpartbuff)
IF(.not. ALLOCATED(rrecvpartbuff) .and. rrecvnbparts .gt. 0) ALLOCATE(rrecvpartbuff(INT(rreceivednbparts*1.3)))
IF( ALLOCATED(lrecvpartbuff) .AND. (size(lrecvpartbuff) .lt. lrecvnbparts) ) DEALLOCATE(lrecvpartbuff)
IF((.not. ALLOCATED(lrecvpartbuff) ).and. lrecvnbparts .gt. 0) ALLOCATE(lrecvpartbuff(INT(lreceivednbparts*1.3)))
! Receive particles from left and right processes
IF ( lrecvnbparts .gt. 0) THEN
CALL MPI_IRECV(lrecvpartbuff, lreceivednbparts, particle_type, leftproc, 1, MPI_COMM_WORLD, recvrequest(1), ierr)
END IF
IF( rrecvnbparts .gt. 0) THEN
CALL MPI_IRECV(rrecvpartbuff, rreceivednbparts, particle_type, rightproc, 1, MPI_COMM_WORLD, recvrequest(2), ierr)
END IF
IF ( (lsendnbparts + rsendnbparts) .gt. 0) THEN
CALL AddPartSendBuffers(lsendnbparts, rsendnbparts, sendholes)
END IF
CALL MPI_Wait(sendrequest(1), sendstatus(1), ierr)
CALL MPI_Wait(sendrequest(2), sendstatus(MPI_STATUS_SIZE+1), ierr)
IF( lsendnbparts .gt. 0) THEN
CALL MPI_ISEND(lsendpartbuff, lsendnbparts, particle_type, leftproc, 1, MPI_COMM_WORLD, sendrequest(1), ierr)
END IF
IF( rsendnbparts .gt. 0) THEN
CALL MPI_ISEND(rsendpartbuff, rsendnbparts, particle_type, rightproc, 1, MPI_COMM_WORLD, sendrequest(2), ierr)
END IF
IF ( lreceivednbparts .gt. 0) THEN
CALL MPI_Wait(recvrequest(1), recvstatus(1), ierr)
IF(ierr .ne. MPI_SUCCESS) THEN
WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(1)
CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
END IF
END IF
IF ( rreceivednbparts .gt. 0) THEN
CALL MPI_Wait(recvrequest(2), recvstatus(MPI_STATUS_SIZE+1), ierr)
IF(ierr .ne. MPI_SUCCESS) THEN
WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(2)
CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
END IF
END IF
CALL Addincomingparts(rreceivednbparts, lreceivednbparts, lsendnbparts+rsendnbparts, sendholes)
IF( lsendnbparts .gt. 0) THEN
CALL MPI_Wait(sendrequest(1), sendstatus(1), ierr)
END IF
IF( rsendnbparts .gt. 0) THEN
CALL MPI_Wait(sendrequest(2), sendstatus(2), ierr)
END IF
!
!
END SUBROUTINE particlescommunication
!_______________________________________________________________________________
SUBROUTINE Addincomingparts(rrecvnbparts, lrecvnbparts, sendnbparts, sendholes)
!
USE mpihelper
INTEGER, INTENT(in) :: rrecvnbparts, lrecvnbparts, sendnbparts
INTEGER, INTENT(in) :: sendholes(parts%Nptot)
INTEGER k,partpos, partdiff
! computes if we received less particles than we sent
partdiff=sendnbparts-rrecvnbparts-lrecvnbparts
IF(rrecvnbparts .gt. 0) THEN
Do k=1,rrecvnbparts
IF(k .le. sendnbparts) THEN
partpos=abs(sendholes(k))
ELSE
parts%Nptot=parts%Nptot+1
partpos=parts%Nptot
END IF
CALL Insertincomingpart(rrecvpartbuff, k, partpos)
END DO
END IF
IF(lrecvnbparts .gt. 0) THEN
Do k=1,lrecvnbparts
IF(k+rrecvnbparts .le. sendnbparts) THEN
partpos=abs(sendholes(k+rrecvnbparts))
ELSE
parts%Nptot=parts%Nptot+1
partpos=parts%Nptot
END IF
CALL Insertincomingpart(lrecvpartbuff, k, partpos)
END DO
END IF
IF(partdiff .gt. 0) THEN
DO k=rrecvnbparts+lrecvnbparts+1, sendnbparts
CALL move_part(parts%Nptot,abs(sendholes(k)))
parts%Nptot = parts%Nptot-1
END DO
END IF
!
END SUBROUTINE Addincomingparts
!_______________________________________________________________________________
SUBROUTINE Insertincomingpart(buffer, bufferindex, partsindex)
USE mpihelper
INTEGER, INTENT(in) :: bufferindex, partsindex
TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(in) :: buffer
parts%Rindex(partsindex) = buffer(bufferindex)%Rindex
parts%Zindex(partsindex) = buffer(bufferindex)%Zindex
parts%R(partsindex) = buffer(bufferindex)%R
parts%Z(partsindex) = buffer(bufferindex)%Z
parts%UZ(partsindex) = buffer(bufferindex)%UZ
parts%UR(partsindex) = buffer(bufferindex)%UR
parts%UTHET(partsindex) = buffer(bufferindex)%UTHET
parts%Gamma(partsindex) = buffer(bufferindex)%Gamma
!
!
END SUBROUTINE Insertincomingpart
!_______________________________________________________________________________
SUBROUTINE Insertsentpart(buffer, bufferindex, partsindex)
USE mpihelper
INTEGER, INTENT(in) :: bufferindex, partsindex
TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(inout) :: buffer
buffer(bufferindex)%Rindex = parts%Rindex(partsindex)
buffer(bufferindex)%Zindex = parts%Zindex(partsindex)
buffer(bufferindex)%R = parts%R(partsindex)
buffer(bufferindex)%Z = parts%Z(partsindex)
buffer(bufferindex)%UZ = parts%UZ(partsindex)
buffer(bufferindex)%UR = parts%UR(partsindex)
buffer(bufferindex)%UTHET = parts%UTHET(partsindex)
buffer(bufferindex)%Gamma = parts%Gamma(partsindex)
!
!
END SUBROUTINE Insertsentpart
!_______________________________________________________________________________
SUBROUTINE AddPartSendBuffers(lsendnbparts, rsendnbparts, sendholes)
!
USE mpihelper
INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts
INTEGER, INTENT(in) :: sendholes(parts%Nptot)
INTEGER:: partpos, k
INTEGER:: lsendpos, rsendpos
lsendpos=0
rsendpos=0
IF( ALLOCATED(rsendpartbuff) .AND. size(rsendpartbuff) .lt. rsendnbparts ) DEALLOCATE(rsendpartbuff)
IF( ALLOCATED(lsendpartbuff) .AND. size(lsendpartbuff) .lt. lsendnbparts ) DEALLOCATE(lsendpartbuff)
IF(.not. ALLOCATED(rsendpartbuff) .and. rsendnbparts .gt. 0) ALLOCATE(rsendpartbuff(INT(rsendnbparts*1.3)))
IF(.not. ALLOCATED(lsendpartbuff) .and. lsendnbparts .gt. 0) ALLOCATE(lsendpartbuff(INT(lsendnbparts*1.3)))
Do k=lsendnbparts+rsendnbparts,1,-1
partpos=abs(sendholes(k))
IF(sendholes(k) .GT. 0) THEN
rsendpos=rsendpos+1
CALL Insertsentpart(rsendpartbuff, rsendpos, partpos)
ELSE IF(sendholes(k) .LT. 0) THEN
lsendpos=lsendpos+1
CALL Insertsentpart(lsendpartbuff, lsendpos, partpos)
END IF
END DO
!
!
END SUBROUTINE AddPartSendBuffers
!_______________________________________________________________________________
SUBROUTINE exchange_parts(index1, index2)
INTEGER, INTENT(IN) :: index1, index2
DOUBLE PRECISION:: R,Z,UR, UZ, UTHET, Gamma
INTEGER :: Rindex, Zindex
!! Exchange particle at index1 with particle at index2
! Store part at index1 in temporary value
Gamma = parts%Gamma(index1)
R = parts%R(index1)
Z = parts%Z(index1)
UR = parts%UR(index1)
UTHET = parts%UTHET(index1)
UZ = parts%UZ(index1)
Rindex = parts%Rindex(index1)
Zindex = parts%Zindex(index1)
! Move part at index2 in part at index 1
parts%Gamma(index1) = parts%Gamma(index2)
parts%R(index1) = parts%R(index2)
parts%Z(index1) = parts%Z(index2)
parts%UR(index1) = parts%UR(index2)
parts%UTHET(index1) = parts%UTHET(index2)
parts%UZ(index1) = parts%UZ(index2)
parts%Rindex(index1) = parts%Rindex(index2)
parts%Zindex(index1) = parts%Zindex(index2)
! Move temporary values from part(index1) to part(index2)
parts%Gamma(index2) = Gamma
parts%R(index2) = R
parts%Z(index2) = Z
parts%UR(index2) = UR
parts%UTHET(index2) = UTHET
parts%UZ(index2) = UZ
parts%Rindex(index2) = Rindex
parts%Zindex(index2) = Zindex
END SUBROUTINE exchange_parts
!_______________________________________________________________________________
SUBROUTINE move_part(sourceindex, destindex)
!! Move particle with index sourceindex to particle with index destindex
!! This will destroy particle at destindex
INTEGER, INTENT(IN) :: destindex, sourceindex
! Move part at sourceindex in part at destindex
parts%Gamma(destindex) = parts%Gamma(sourceindex)
parts%R(destindex) = parts%R(sourceindex)
parts%Z(destindex) = parts%Z(sourceindex)
parts%UR(destindex) = parts%UR(sourceindex)
parts%UTHET(destindex) = parts%UTHET(sourceindex)
parts%UZ(destindex) = parts%UZ(sourceindex)
parts%Rindex(destindex) = parts%Rindex(sourceindex)
parts%Zindex(destindex) = parts%Zindex(sourceindex)
END SUBROUTINE move_part
!_______________________________________________________________________________
!________________________________________________________________________________
!
SUBROUTINE loduni(nbase, y)
!
! Load a uniform distribution using the Hammersley's sequence.
! (nbase = 0 => Random sampling)
!
INTEGER, INTENT(in) :: nbase
DOUBLE PRECISION, INTENT(out) :: y(:)
INTEGER :: n, i
!
n = SIZE(y)
!
SELECT CASE (nbase)
CASE(0)
CALL RANDOM_NUMBER(y)
CASE(1)
DO i=1,n
y(i) = (i-0.5_db)/n
END DO
CASE(2:)
DO i=1,n
y(i) = rev(nbase,i)
END DO
CASE default
WRITE(*,'(a,i5)') 'Invalid value of NBASE =', nbase
END SELECT
!
END SUBROUTINE loduni
!________________________________________________________________________________
SUBROUTINE lodgaus(nbase, y)
!
! Sample y from the Gauss distributrion
!
DOUBLE PRECISION, INTENT(out) :: y(:)
INTEGER, INTENT(in) :: nbase
!
INTEGER :: n, i, iflag
DOUBLE PRECISION :: r(SIZE(y))
DOUBLE PRECISION :: t
!
DOUBLE PRECISION :: c0, c1, c2
DOUBLE PRECISION :: d1, d2, d3
DATA c0, c1, c2/2.515517_db, 0.802853_db, 0.010328_db/
DATA d1, d2, d3/1.432788_db, 0.189269_db, 0.001308_db/
!
n = SIZE(y)
CALL loduni(nbase, r)
r = 1.0E-5_db + 0.99998_db*r
!
DO i=1,n
iflag = 1
IF (r(i) .GT. 0.5_db) THEN
r(i) = 1.0_db - r(i)
iflag = -1
END IF
t = SQRT(LOG(1.0_db/(r(i)*r(i))))
y(i) = t - (c0+c1*t+c2*t**2) / (1.0_db+d1*t+d2*t**2+d3*t**3)
y(i) = y(i) * iflag
END DO
y = y - SUM(y)/REAL(n,db)
END SUBROUTINE lodgaus
!________________________________________________________________________________
SUBROUTINE lodinvr(nbase, y, ra, rb)
!
! Sample y from the distribution f=1/(r)H(r-ra)H(rb-r) for where H is the heavyside function
!
DOUBLE PRECISION, INTENT(out) :: y(:)
INTEGER, INTENT(in) :: nbase
DOUBLE PRECISION, INTENT(in) :: rb, ra
!
DOUBLE PRECISION :: r(SIZE(y))
!
CALL loduni(nbase, r)
r=r*log(rb/ra)
y=ra*exp(r)
END SUBROUTINE lodinvr
!________________________________________________________________________________
SUBROUTINE lodunir(nbase, y, ra, rb)
DOUBLE PRECISION, INTENT(out) :: y(:)
INTEGER, INTENT(in) :: nbase
DOUBLE PRECISION, INTENT(in) :: rb, ra
!
DOUBLE PRECISION :: r(SIZE(y))
!
CALL loduni(nbase, r)
y=ra+(rb-ra)*r
END SUBROUTINE lodunir
!________________________________________________________________________________
REAL(db) FUNCTION rev(nbase,i)
!
! Return an element of the Hammersley's sequence
!
INTEGER, INTENT(IN) :: nbase ! Base of the Hammersley's sequence (IN)
INTEGER, INTENT(IN) :: i ! Index of the sequence (IN)
!
! Local vars
INTEGER :: j1, j2
REAL(db) :: xs, xsi
!-----------------------------------------------------------------------
xs = 0._db
xsi = 1.0_db
j2 = i
DO
xsi = xsi/nbase
j1 = j2/nbase
xs = xs + (j2-nbase*j1)*xsi
j2 = j1
IF( j2.LE.0 ) EXIT
END DO
rev = xs
END FUNCTION rev
!________________________________________________________________________________
!Modified Bessel functions of the first kind of the first order
FUNCTION bessi1(x)
DOUBLE PRECISION :: bessi1,x
DOUBLE PRECISION :: ax
DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/
DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, -0.3988024d-1, -0.362018d-2, 0.163801d-2,-0.1031555d-1, &
& 0.2282967d-1, -0.2895312d-1, 0.1787654d-1,-0.420059d-2/
if (abs(x).lt.3.75) then
y=(x/3.75)**2
bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
else
ax=abs(x)
y=3.75/ax
bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))))
if(x.lt.0.)bessi1=-bessi1
endif
return
END FUNCTION bessi1
!________________________________________________________________________________
SUBROUTINE destroy_parts(p)
TYPE(particles) :: p
p%nptot=0
IF(ALLOCATED(p%Z)) DEALLOCATE(p%Z)
IF(ALLOCATED(p%R)) DEALLOCATE(p%R)
IF(ALLOCATED(p%WZ)) DEALLOCATE(p%WZ)
IF(ALLOCATED(p%WR)) DEALLOCATE(p%WR)
IF(ASSOCIATED(p%UR)) DEALLOCATE(p%UR)
IF(Associated(p%URold)) DEALLOCATE(p%URold)
IF(Associated(p%UZ)) DEALLOCATE(p%UZ)
IF(Associated(p%UZold)) DEALLOCATE(p%UZold)
IF(Associated(p%UTHET)) DEALLOCATE(p%UTHET)
IF(Associated(p%UTHETold)) DEALLOCATE(p%UTHETold)
IF(Associated(p%Gamma)) DEALLOCATE(p%Gamma)
IF(Associated(p%Gammaold)) DEALLOCATE(p%Gammaold)
IF(ALLOCATED(p%Rindex)) DEALLOCATE(p%Rindex)
IF(ALLOCATED(p%Zindex)) DEALLOCATE(p%Zindex)
END SUBROUTINE
!________________________________________________________________________________
SUBROUTINE clean_beam
!
CALL destroy_parts(parts)
!
!
END SUBROUTINE clean_beam
!________________________________________________________________________________
RECURSIVE SUBROUTINE quicksortparts(leftlimit, rightlimit)
! Sorts the particle according to their Z position using quicksort algorithm
INTEGER,INTENT(IN):: leftlimit, rightlimit
DOUBLE PRECISION:: pivot
INTEGER::i, cnt, mid
IF(leftlimit .ge. rightlimit) RETURN
mid=(leftlimit+rightlimit)/2
IF(parts%Z(mid).lt.parts%Z(leftlimit)) CALL exchange_parts(leftlimit,mid)
IF(parts%Z(rightlimit).lt.parts%Z(leftlimit)) CALL exchange_parts(leftlimit,rightlimit)
IF(parts%Z(mid).lt.parts%Z(rightlimit)) CALL exchange_parts(rightlimit,mid)
! Store the pivot point for comparison
pivot=parts%Z(rightlimit)
cnt=leftlimit
! Move all parts with Z smaller than pivot to the left of pivot
DO i=leftlimit, rightlimit
IF(parts%Z(i) .le. pivot) THEN
CALL exchange_parts(i,cnt)
cnt=cnt+1
END IF
END DO
CALL quicksortparts(leftlimit,cnt-2)
CALL quicksortparts(cnt,rightlimit)
END SUBROUTINE quicksortparts
!________________________________________________________________________________
SUBROUTINE swappointer( pointer1, pointer2)
DOUBLE PRECISION, DIMENSION(:), POINTER, INTENT(inout):: pointer1, pointer2
DOUBLE PRECISION, DIMENSION(:), POINTER:: temppointer
temppointer=>pointer1
pointer1=>pointer2
pointer2=>temppointer
END SUBROUTINE swappointer
!_______________________________________________________________________________
SUBROUTINE loaduniformRZ(VR,VZ,VTHET)
USE basic, ONLY: plasmadim, rnorm, temp, vnorm
USE constants, ONLY: me, kb, elchar
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VZ, VR, VTHET
DOUBLE PRECISION:: vth
! Initial distribution in z with normalisation
CALL loduni(1,parts%Z(:))
parts%Z(:)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*parts%Z(:))/rnorm
! Initial distribution in r with normalisation
CALL loduni(2,parts%R(:))
parts%R=(plasmadim(3)+parts%R(:)*(plasmadim(4)-plasmadim(3)))/rnorm
! Initial velocities distribution
vth=sqrt(3*kb*temp/me)/vnorm !thermal velocity
CALL lodgaus(3,VZ(:))
CALL lodgaus(5,VR(:))
CALL lodgaus(7,VTHET(:))
VZ=VZ*vth
VR=VR*vth
VTHET=VTHET*vth
END SUBROUTINE loaduniformRZ
!_______________________________________________________________________________
SUBROUTINE loadDavidson(VR,VZ,VTHET,lodr)
USE constants, ONLY: me, kb, elchar
USE basic, ONLY: nplasma, rnorm, plasmadim, distribtype, H0, P0, Rcurv, B0, width, qsim, msim, vnorm, &
& omegac, zgrid, nz, rnorm, V, n0, nblock
interface
subroutine lodr(nbase,y,ra,rb)
DOUBLE PRECISION, INTENT(out) :: y(:)
INTEGER, INTENT(in) :: nbase
DOUBLE PRECISION, INTENT(in) :: rb, ra
end subroutine
end interface
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VZ, VR, VTHET
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE::ra, rb
DOUBLE PRECISION :: r0, athetpos, deltar2, rg, zg, halfLz, Mirrorratio, Le, Pcomp, Acomp, gamma, VOL
INTEGER :: i, j, n, blockstart, blockend, addedpart, remainparts
INTEGER, DIMENSION(:), ALLOCATABLE :: blocksize
Allocate(ra(nblock),rb(nblock))
!r0=(plasmadim(4)+plasmadim(3))/2
r0=sqrt(2*H0/(me*omegac**2))
halfLz=(zgrid(nz)+zgrid(0))/2
MirrorRatio=(Rcurv-1)/(Rcurv+1)
DO n=1,nblock
! Compute limits in radius and load radii for each part
Le=(plasmadim(2)-plasmadim(1))/nblock*(n-0.5)-halfLz*rnorm+plasmadim(1)
deltar2=1-MirrorRatio*cos(2*pi*Le/width)
rb(n)=r0*(1+sqrt(1-abs(omegac)*P0*deltar2/H0))/deltar2
ra(n)=r0*(1-sqrt(1-abs(omegac)*P0*deltar2/H0))/deltar2
WRITE(*,'(a,7(1pe12.4))') "ra,rb: ", omegac, H0, P0, r0, deltar2, ra(n), rb(n)
END DO
VOL=SUM(2*pi*ra*(rb-ra)*(plasmadim(2)-plasmadim(1))/nblock)
qsim=VOL*n0*elchar/nplasma
msim=abs(qsim)/elchar*me
V=VOL/rnorm**3
!H0=me*omegac**2*r0**2*0.5
gamma=1/sqrt(1-omegac**2*r0**2/vlight**2)
!P0=H0/omegac
blockstart=1
blockend=0
ALLOCATE(blocksize(nblock))
DO n=1,nblock
blocksize(n)=nplasma/VOL*2*pi*ra(n)*(rb(n)-ra(n))*(plasmadim(2)-plasmadim(1))/nblock
END DO
remainparts=parts%Nptot-SUM(blocksize)
addedpart=1
n=nblock/2
j=1
DO WHILE(remainparts .GT. 0)
blocksize(n)=blocksize(n)+addedpart
remainparts=remainparts-addedpart
n=n+j
j=-1*(j+SIGN(1,j))
END DO
DO n=1,nblock
blockstart=blockend+1
blockend=MIN(blockstart+blocksize(n)-1,parts%Nptot)
! Initial distribution in z with normalisation between magnetic mirrors
CALL loduni(1,parts%Z(blockstart:blockend))
parts%Z(blockstart:blockend)=(plasmadim(2)-plasmadim(1))/rnorm/nblock*((n-1)+parts%Z(blockstart:blockend))+plasmadim(1)/rnorm
CALL lodr(3,parts%R(blockstart:blockend),ra(n), rb(n))
parts%R(blockstart:blockend)=parts%R(blockstart:blockend)/rnorm
END DO
! Load velocities theta velocity
! Loading of r and z velocity is done in adapt_vinit to have
! access to parts%pot
DO i=1,parts%Nptot
! Interpolation for Magnetic potential
rg=parts%R(i)*rnorm
zg=(parts%Z(i)-halfLz)*rnorm
Athetpos=0.5*B0*(rg - width/pi*MirrorRatio*bessi1(2*pi*rg/width)*COS(2*pi*zg/width))
Pcomp=abs(P0/rg/me)
Acomp=-SIGN(elchar/me*Athetpos,qsim)
VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/me)),Pcomp+Acomp)
END DO
VTHET=VTHET/vnorm
VZ=0._db
VR=0._db
END SUBROUTINE loadDavidson
END MODULE beam
c4science · Help