Verified the collision operator; however, collision operator does not scale well with openMP (e4ef9d6c) · Commits · Kumar, Atul / LinearFokkerPlanck Axisymmetric

InputFiles/xp_IonSlowingDown.in

+3 −3

Original line number	Diff line number	Diff line
		&in_nml
		! Simulation name:
		! ===============
		in%fileDescriptor = '2021_02_07b',
		in%fileDescriptor = '2021_02_07d',

		! Magnetic field input data:
		! =========================
		@@ -19,7 +19,7 @@ in%zmax = +5.0,
		in%zmin = -5.0,

		in%iSave = .true.,
		in%iPush = .true., ! Enable RK4 particle pusher
		in%iPush = .false., ! Enable RK4 particle pusher
		in%iColl = .true., ! Collisions or no collisions
		in%iHeat = .false., ! Enable RF heating operator
		in%iPotential = .false., ! Enable electric potential
		@@ -53,7 +53,7 @@ in%zp_InitType = 2, !
		in%zp_init = 0.0, ! Mean particle injection
		in%zp_init_std = 0.3, ! STD of particle injection
		in%Ep_init = 1000., ! Drift kinetic energy
		in%Tp_init = 1.; ! Thermal kinetic energy
		in%Tp_init = 0.1, ! Thermal kinetic energy
		in%xip_init = 0.9240, ! Mean pitch angle where xip = cos(theta) = vpar/v

		! RF heating operator conditions:

src/linFP.f90

+86 −98

Original line number	Diff line number	Diff line
		@@ -34,8 +34,8 @@ INTEGER(i4), DIMENSION(:), ALLOCATABLE :: jrng
		INTEGER(i4) :: id
		! OPENMP computational time:
		DOUBLE PRECISION :: ostart, oend, oend_estimate
		! Cyclotron resonance number change:
		REAL(r8) :: df
		! Cyclotron resonance numbers:
		REAL(r8) :: df, f0, f1
		! Main simulation variables:
		! simulation time:
		REAl(r8) :: tp
		@@ -50,8 +50,6 @@ REAl(r8), DIMENSION(:) , ALLOCATABLE :: t_hist
		REAL(r8), DIMENSION(:) , ALLOCATABLE :: pcount1, pcount2, pcount3, pcount4
		! Record the total energy of particle incident on (1) dump, (2) target, (3) cyclotron resonance, (4) slowing down
		REAL(r8), DIMENSION(:) , ALLOCATABLE :: ecount1, ecount2, ecount3, ecount4
		! Cyclotron resonance number:
		REAL(r8), DIMENSION(:) , ALLOCATABLE :: fcurr, fnew
		! Local simulation diagnostic counters:
		REAL(r8) :: ecnt1, ecnt2, ecnt3, ecnt4
		REAL(r8) :: pcnt1, pcnt2, pcnt3, pcnt4
		@@ -144,11 +142,6 @@ ALLOCATE(zp_hist(in%Nparts,jsize),kep_hist(in%Nparts,jsize),xip_hist(in%Nparts,j
		! Create array with the indices of the time steps to save:
		jrng = (/ (j, j=in%jstart, in%jend, in%jincr) /)

		! Allocate memory to local variables:
		! ==============================================================================
		! Cyclotron resonance number:
		ALLOCATE(fcurr(in%Nparts),fnew(in%Nparts))

		! Create splines of electromagnetic fields:
		! ===========================================================================
		! Magnetic field data:
		@@ -213,6 +206,8 @@ CLOSE(unit=8)

		! Initialize simulation diagnostics:
		! ==============================================================================
		! Record time index:
		k = 1
		! Time array:
		tp = 0
		! Leak current diagnostics:
		@@ -228,7 +223,7 @@ ostart = OMP_GET_WTIME()
		! ==============================================================================
		AllTime: do j = 1,in%Nsteps

		!$OMP PARALLEL PRIVATE(pcnt1,pcnt2,pcnt3,pcnt4,ecnt1,ecnt2,ecnt3,ecnt4,df)
		!$OMP PARALLEL PRIVATE(pcnt1,pcnt2,pcnt3,pcnt4,ecnt1,ecnt2,ecnt3,ecnt4,df,f0,f1)

		! Initialize private particle counters:
		pcnt1 = 0; pcnt2 = 0; pcnt3 = 0; pcnt4 = 0
		@@ -244,8 +239,7 @@ AllTime: do j = 1,in%Nsteps

		! Calculate Cyclotron resonance number:
		! ------------------------------------------------------------------------
		if (in%iHeat) CALL CyclotronResonanceNumber(zp(i),kep(i),xip(i),fcurr(i),in,spline_B)
		! fcurr and fnew could be declared private
		if (in%iHeat) CALL CyclotronResonanceNumber(zp(i),kep(i),xip(i),f0,in,spline_B)

		! Push particles adiabatically:
		! ------------------------------------------------------------------------
		@@ -270,8 +264,8 @@ AllTime: do j = 1,in%Nsteps
		! Apply RF heating operator:
		! ------------------------------------------------------------------------
		if (in%iHeat) then
		CALL CyclotronResonanceNumber(zp(i),kep(i),xip(i),fnew(i),in,spline_B)
		df = dsign(1.d0,fcurr(i)*fnew(i))
		CALL CyclotronResonanceNumber(zp(i),kep(i),xip(i),f1,in,spline_B)
		df = dsign(1.d0,f0*f1)
		if (df .LT. 0 .AND. zp(i) .GT. in%zRes1 .AND. zp(i) .LT. in%zRes2) then
		CALL RFHeatingOperator(zp(i),kep(i),xip(i),ecnt3,pcnt3,in,spline_B,spline_dB,spline_ddB,spline_dV)
		end if
		@@ -293,40 +287,34 @@ AllTime: do j = 1,in%Nsteps

		!$OMP END PARALLEL

		! Update time array:
		! =========================================================================
		tp = tp + in%dt

		! Select data to save:
		! =====================================================================
		! Check if data is to be saved
		if (in%iSave) then
		do k = 1,jsize
		if (j .EQ. jrng(k)) then
		t_hist(k) = tp
		!$OMP PARALLEL DO PRIVATE(i)
		!$OMP PARALLEL DO
		do i = 1,in%Nparts
		! Record "ith" particle position at "kth" time
		! Record "ith" particle at "kth" time
		zp_hist(i,k) = zp(i)
		! Record "ith" particle KE at "kth" time
		kep_hist(i,k) = kep(i)
		! Record "ith" particle pitch angle at "kth" time
		xip_hist(i,k) = xip(i)
		end do
		!$OMP END PARALLEL DO
		k = k + 1
		end if
		end do
		end if

		! Update time array:
		! =========================================================================
		tp = tp + in%dt

		! Estimate computational time:
		! =====================================================================
		id = OMP_GET_THREAD_NUM()
		if (id .EQ. 0) then
		if (j .EQ. 150) then
		oend_estimate = OMP_GET_WTIME()
		WRITE(,) 'Estimated compute time: ', in%Nsteps*(oend_estimate-ostart)/j,' [s]'
		end if
		end if

		end do AllTime

src/run_linFP.sh

+1 −1

Original line number	Diff line number	Diff line
		@@ -7,7 +7,7 @@ INPUT_FILE="xp_IonSlowingDown.in"

		# Set number of threads:
		# ===================================================
		export OMP_NUM_THREADS=48
		export OMP_NUM_THREADS=1

		# Set processor binding for openMP:
		# ===================================================