Made output variables into a structure (4e191980) · Commits · Kumar, Atul / LinearFokkerPlanck Axisymmetric

InputFiles/xp_TestDrag.in

+1 −1

Original line number	Diff line number	Diff line
		&params_nml
		! Simulation name:
		! ===============
		params%fileDescriptor = '2021_02_17c',
		params%fileDescriptor = '2021_02_19a',

		! Magnetic field input data:
		! =========================

src/Modules.f90

+102 −23

Original line number	Diff line number	Diff line
		@@ -40,7 +40,7 @@ END TYPE splTYP

		CONTAINS
		! ----------------------------------------------------------------------------
		SUBROUTINE InitSpline(spline0,n,slp1,slpn)
		SUBROUTINE AllocateSpline(spline0,n,slp1,slpn)
		IMPLICIT NONE
		TYPE(splTYP) :: spline0
		INTEGER(i4) :: n
		@@ -51,7 +51,7 @@ CONTAINS
		spline0%n = n
		spline0%slp1 = slp1
		spline0%slpn = slpn
		END SUBROUTINE InitSpline
		END SUBROUTINE AllocateSpline

		! ----------------------------------------------------------------------------
		SUBROUTINE ReadSpline(spline0,fileName)
		@@ -208,13 +208,16 @@ END MODULE spline_fits

		! MODULE dataTYP
		! =============================================================================
		! Module containing definition of an object to contain all data used in
		! computation
		! The following types are prototyped:
		! paramsTYP: Contain simulation basic and derived parameters
		! plasmaTYP: Contain arrays of simulation data
		! fieldSplineTYP: Contain splines for fields
		MODULE dataTYP
		USE local
		USE spline_fits

		IMPLICIT NONE
		! -----------------------------------------------------------------------------
		TYPE paramsTYP
		! Simulation name:
		CHARACTER*150 :: fileDescriptor
		@@ -251,6 +254,7 @@ TYPE paramsTYP

		END TYPE paramsTYP

		! -----------------------------------------------------------------------------
		TYPE plasmaTYP
		REAL(r8) , DIMENSION(:), ALLOCATABLE :: zp, kep, xip, a
		INTEGER(i4), DIMENSION(:), ALLOCATABLE :: f1 , f2 , f3 , f4
		@@ -262,13 +266,21 @@ TYPE plasmaTYP
		REAL(r8) , DIMENSION(:), ALLOCATABLE :: Edot1, Edot2, Edot3, Edot4, Edot3_hat
		END TYPE plasmaTYP

		! -----------------------------------------------------------------------------
		TYPE fieldSplineTYP
		TYPE(splTYP) :: B, dB, ddB, V, dV, U
		END TYPE fieldSplineTYP

		! -----------------------------------------------------------------------------
		TYPE outputTYP
		REAL(r8) , DIMENSION(:,:), ALLOCATABLE :: zp, kep, xip, a
		REAL(r8) , DIMENSION(:) , ALLOCATABLE :: tp, jrng
		INTEGER(i4) :: jsize
		END TYPE outputTYP

		CONTAINS
		! ----------------------------------------------------------------------------
		SUBROUTINE InitPlasma(plasma,params)
		SUBROUTINE AllocatePlasma(plasma,params)
		IMPLICIT NONE
		! Declare interface variables:
		TYPE(plasmaTYP), INTENT(INOUT) :: plasma
		@@ -277,31 +289,74 @@ SUBROUTINE InitPlasma(plasma,params)
		! Declare local variables:
		INTEGER(i4) :: NC, NS

		! NC: Number of computational particles
		NC = params%NC
		! NS: Number of time steps
		NS = params%NS

		! Allocate memory:
		! Allocate memory: For all computational particles
		ALLOCATE(plasma%zp(NC) ,plasma%kep(NC) ,plasma%xip(NC) ,plasma%a(NC))
		ALLOCATE(plasma%f1(NC) ,plasma%f2(NC) ,plasma%f3(NC) ,plasma%f4(NC))
		ALLOCATE(plasma%dE1(NC) ,plasma%dE2(NC) ,plasma%dE3(NC) ,plasma%dE4(NC))
		ALLOCATE(plasma%dErf_hat(NC),plasma%doppler(NC))

		! Allocate memory: For all time steps
		ALLOCATE(plasma%NR(NS) ,plasma%NSP(NS))
		ALLOCATE(plasma%Eplus(NS),plasma%Eminus(NS))
		ALLOCATE(plasma%Ndot1(NS),plasma%Ndot2(NS) ,plasma%Ndot3(NS),plasma%Ndot4(NS))
		ALLOCATE(plasma%Edot1(NS),plasma%Edot2(NS) ,plasma%Edot3(NS),plasma%Edot4(NS))

		! Allocate memory: RF heating operator terms
		ALLOCATE(plasma%dE3_hat(NC),plasma%Edot3_hat(NS))
		ALLOCATE(plasma%dErf_hat(NC),plasma%doppler(NC))

		END SUBROUTINE AllocatePlasma

		! --------------------------------------------------------------------------
		SUBROUTINE InitializePlasma(plasma,params)
		USE PhysicalConstants
		USE OMP_LIB
		IMPLICIT NONE
		! Declare interface variables:
		TYPE(plasmaTYP), INTENT(INOUT) :: plasma
		TYPE(paramsTYP), INTENT(INOUT) :: params
		INTEGER(i4) :: i

		! Derived parameters: Reference cross sectional area
		params%Area0 = 0.5params%dtheta( params%r22. - params%r12.)

		! Derived parameters: Select test species
		if (params%species_a .EQ. 1) then
		params%qa = -e_c
		params%Ma = m_e
		else
		params%qa = +params%Zion*e_c
		params%Ma = params%Aion*m_p
		end if

		! Initialize plasma: NR and NC
		! When source contrained fueling is used, these quantities
		! can be time depedent
		plasma%NR = params%ne0params%Area0(params%zmax - params%zmin)
		plasma%NSP = params%NC

		! Initialize plasma: particle weights
		plasma%a = 1.

		! Initialize variables:
		!$OMP PARALLEL DO
		DO i = 1,params%NC
		CALL loadParticles(i,plasma,params)
		END DO
		!$OMP END PARALLEL DO

		! Initialize plasma: flags and Energy increments
		plasma%f1 = 0.; plasma%f2 = 0.; plasma%f3 = 0.; plasma%f4 = 0.
		plasma%dE1 = 0.; plasma%dE2 = 0.; plasma%dE3 = 0.; plasma%dE4 = 0.
		plasma%dErf_hat = 0.; plasma%doppler = 0.;

		END SUBROUTINE InitPlasma
		END SUBROUTINE InitializePlasma

		! --------------------------------------------------------------------------
		SUBROUTINE InitFieldSpline(fieldspline,params)
		SUBROUTINE AllocateFieldSpline(fieldspline,params)
		USE spline_fits

		IMPLICIT NONE
		@@ -316,14 +371,14 @@ SUBROUTINE InitFieldSpline(fieldspline,params)
		NZ = params%nz

		! Allocate memory:
		CALL InitSpline(fieldspline%B ,NZ,0._8,0._8)
		CALL InitSpline(fieldspline%dB ,NZ,0._8,0._8)
		CALL InitSpline(fieldspline%ddB,NZ,0._8,0._8)
		CALL InitSpline(fieldspline%V ,NZ,0._8,0._8)
		CALL InitSpline(fieldspline%dV ,NZ,0._8,0._8)
		CALL InitSpline(fieldspline%U ,NZ,0._8,0._8)
		CALL AllocateSpline(fieldspline%B ,NZ,0._8,0._8)
		CALL AllocateSpline(fieldspline%dB ,NZ,0._8,0._8)
		CALL AllocateSpline(fieldspline%ddB,NZ,0._8,0._8)
		CALL AllocateSpline(fieldspline%V ,NZ,0._8,0._8)
		CALL AllocateSpline(fieldspline%dV ,NZ,0._8,0._8)
		CALL AllocateSpline(fieldspline%U ,NZ,0._8,0._8)

		END SUBROUTINE InitFieldSpline
		END SUBROUTINE AllocateFieldSpline

		SUBROUTINE ComputeFieldSpline(fieldspline)
		USE spline_fits
		@@ -342,4 +397,28 @@ SUBROUTINE ComputeFieldSpline(fieldspline)

		END SUBROUTINE ComputeFieldSpline

		!------------------------------------------------------------------------
		SUBROUTINE AllocateOutput(output,params)
		IMPLICIT none
		! Declare interface variables:
		TYPE(outputTYP), INTENT(INOUT) :: output
		TYPE(paramsTYP), INTENT(IN) :: params
		INTEGER(i4) :: j

		! Determine size of temporal snapshots to record:
		output%jsize = (params%jend-params%jstart+1)/params%jincr

		! Allocate memory:
		ALLOCATE(output%jrng(output%jsize))
		ALLOCATE(output%zp(params%NC ,output%jsize))
		ALLOCATE(output%kep(params%NC,output%jsize))
		ALLOCATE(output%xip(params%NC,output%jsize))
		ALLOCATE(output%a(params%NC ,output%jsize))
		ALLOCATE(output%tp(output%jsize))

		! Create array with the indices of the time steps to save:
		output%jrng = (/ (j, j=params%jstart, params%jend, params%jincr) /)

		END SUBROUTINE AllocateOutput

		END MODULE dataTYP

src/linFP.f90

+50 −98

Original line number	Diff line number	Diff line
		@@ -21,12 +21,9 @@ IMPLICIT NONE
		TYPE(paramsTYP) :: params
		TYPE(plasmaTYP) :: plasma
		TYPE(fieldSplineTYP) :: fieldspline
		TYPE(outputTYP) :: output
		! DO loop indices:
		INTEGER(i4) :: i,j,k
		! Size of time interval:
		INTEGER(i4) :: jsize
		! Indices of time steps to save:
		INTEGER(i4), DIMENSION(:), ALLOCATABLE :: jrng
		! Thread ID:
		INTEGER(i4) :: id
		! OPENMP computational time:
		@@ -36,10 +33,6 @@ REAL(r8) :: dresNum, resNum0, resNum1
		! Main simulation variables:
		! simulation time:
		REAl(r8) :: tp
		! subset of zp, kep and xip to save:
		REAL(r8), DIMENSION(:,:), ALLOCATABLE :: zp_hist, kep_hist, xip_hist, a_hist
		! Subset of tp:
		REAl(r8), DIMENSION(:) , ALLOCATABLE :: t_hist
		! Simulation diagnotics:
		! Count the number of particles incident on (1) dump, (2) target, (3) cyclotron resonance, (4) slowing down
		REAL(r8), DIMENSION(:) , ALLOCATABLE :: pcount1, pcount2, pcount3, pcount4
		@@ -52,40 +45,21 @@ REAL(r8) :: pcnt1, pcnt2, pcnt3, pcnt4
		CHARACTER*300 :: command, mpwd
		INTEGER(i4) :: n_mpwd, STATUS
		CHARACTER*300 :: inputFileDir, inputFile, fileName, xpSelector, repoDir, dir0, dir1

		! Create input namelist from the user-defined structures:
		! ==============================================================================
		! Namelists:
		namelist/params_nml/params

		! Get root directory:
		! Environment:
		! ==============================================================================
		CALL GET_ENVIRONMENT_VARIABLE('REPO_DIR',repoDir)

		! Get input file name and directory:
		! =============================================================================
		CALL GET_ENVIRONMENT_VARIABLE('INPUT_FILE',inputFile)
		CALL GET_ENVIRONMENT_VARIABLE('INPUT_FILE_DIR',inputFileDir)

		! Read input data into "params" structure:
		! Read input parameters into "params" structure:
		! ==============================================================================
		OPEN(unit=4,file=inputFileDir,status='old',form='formatted')
		read(4,params_nml)
		CLOSE(unit=4)

		! Calculate cross sectional area of plasma at reference location:
		! ==============================================================================
		params%Area0 = 0.5params%dtheta( params%r22. - params%r1*1.)

		! Select the test species:
		! ==============================================================================
		if (params%species_a .eq. 1) then
		params%qa = -e_c
		params%Ma = m_e
		else
		params%qa = +params%Zion*e_c
		params%Ma = params%Aion*m_p
		end if

		! Print to the terminal:
		! ==============================================================================
		WRITE(,) ''
		@@ -116,27 +90,17 @@ WRITE(,) ''

		! Allocate memory to splines:
		! ===========================================================================
		CALL InitFieldSpline(fieldspline,params)
		CALL AllocateFieldSpline(fieldspline,params)

		! Allocate memory to main simulation variables:
		! ==============================================================================
		CALL InitPlasma(plasma,params)
		CALL AllocatePlasma(plasma,params)
		ALLOCATE(pcount1(params%NS) ,pcount2(params%NS) ,pcount3(params%NS) ,pcount4(params%NS))
		ALLOCATE(ecount1(params%NS) ,ecount2(params%NS) ,ecount3(params%NS) ,ecount4(params%NS))

		! Allocate memory to output variables:
		! ==============================================================================
		! Determine size of temporal snapshots to record:
		jsize = (params%jend-params%jstart+1)/params%jincr
		ALLOCATE(jrng(jsize) )
		ALLOCATE(zp_hist(params%NC,jsize) )
		ALLOCATE(kep_hist(params%NC,jsize))
		ALLOCATE(xip_hist(params%NC,jsize))
		ALLOCATE(a_hist(params%NC,jsize) )
		ALLOCATE(t_hist(jsize) )

		! Create array with the indices of the time steps to save:
		jrng = (/ (j, j=params%jstart, params%jend, params%jincr) /)
		CALL AllocateOutput(output,params)

		! Create splines of electromagnetic fields:
		! ===========================================================================
		@@ -145,56 +109,39 @@ fileName = trim(adjustl(params%BFieldFile))
		fileName = trim(adjustl(params%BFieldFileDir))//fileName
		fileName = trim(adjustl(params%repoDir))//fileName

		! Populate spline profile data using external file:
		! Magnetic field B:
		! Magnetic field B, dB and ddB:
		CALL ReadSpline(fieldspline%B,fileName)
		! dB and ddB:
		CALL diffSpline(fieldspline%B ,fieldspline%dB )
		CALL diffSpline(fieldspline%dB,fieldspline%ddB)
		! Electric potential V:

		! Electric potential V, dV:
		fieldspline%V%x = fieldspline%B%x
		fieldspline%V%y = 0
		if (params%iPotential) then
		CALL PotentialProfile(fieldspline,params)
		end if
		CALL ComputeSpline(fieldspline%V)
		! dV:
		CALL diffSpline(fieldspline%V,fieldspline%dV)

		! Plasma flow U:
		fieldspline%U%x = fieldspline%B%x
		fieldspline%U%y = 0
		fieldspline%U%y = 0.

		! Complete setting up the spline data :
		CALL ComputeFieldSpline(fieldspline)

		! Inititalize zp, kep, xip
		! Initialize simulation variables:
		! ==============================================================================
		WRITE(,) "Initializing PDF..."
		CALL InitializePlasma(plasma,params)

		!$OMP PARALLEL
		if (OMP_GET_THREAD_NUM() .EQ. 0) then
		params%threads_given = OMP_GET_NUM_THREADS()
		if (OMP_GET_THREAD_NUM() .EQ. 0) params%threads_given = OMP_GET_NUM_THREADS()
		!$OMP END PARALLEL
		WRITE(,) ''
		WRITE(,) '*********************************************************************'
		WRITE(,) "Number of threads given: ", params%threads_given
		WRITE(,) '*********************************************************************'
		WRITE(,) ''
		end if
		!$OMP DO
		DO i = 1,params%NC
		CALL loadParticles(i,plasma,params)
		! Need to initialize "a"
		END DO
		!$OMP END DO
		!$OMP END PARALLEL
		WRITE(,) "Initialization complete"

		! Test initial distribution:
		! ==========================================================================
		fileName = "LoadParticles.dat"
		OPEN(unit=8,file=fileName,form="formatted",status="unknown")
		do i = 1,params%NC
		WRITE(8,*) plasma%zp(i), plasma%kep(i), plasma%xip(i)
		end do
		CLOSE(unit=8)

		! Initialize simulation diagnostics:
		! ==============================================================================
		@@ -206,9 +153,6 @@ tp = 0
		pcount1 = 0; pcount2 = 0; pcount3 = 0; pcount4 = 0
		! Energy leak diagnotics:
		ecount1 = 0; ecount2 = 0; ecount3 = 0; ecount4 = 0
		! NR and NSP:
		plasma%NR(1) = params%ne0paramsArea0*(params%zmax - params%zmin)
		plasma%NC(1) = params%NC

		! Record start time:
		! ==============================================================================
		@@ -230,36 +174,44 @@ AllTime: do j = 1,params%NS
		! Loop over particles:
		! ==============================================================================
		!$OMP DO SCHEDULE(STATIC)
		AllParticles: do i = 1,params%NC
		AllParticles: DO i = 1,params%NC

		! Calculate Cyclotron resonance number:
		! ------------------------------------------------------------------------
		if (params%iHeat) CALL CyclotronResonanceNumber(i,plasma,resNum0,fieldspline,params)
		IF (params%iHeat) CALL CyclotronResonanceNumber(i,plasma,resNum0,fieldspline,params)

		! Push particles adiabatically:
		! ------------------------------------------------------------------------
		if (params%iPush) CALL MoveParticle(i,plasma,fieldspline,params)
		IF (params%iPush) CALL MoveParticle(i,plasma,fieldspline,params)

		! Re-inject particles:
		! ------------------------------------------------------------------------
		if (plasma%zp(i) .GE. params%zmax) CALL ReinjectParticles(i,plasma,params,ecnt2,pcnt2)
		if (plasma%zp(i) .LE. params%zmin) CALL ReinjectParticles(i,plasma,params,ecnt1,pcnt1)
		IF (plasma%zp(i) .GE. params%zmax) THEN
		plasma%f2(i) = 1
		plasma%dE2(i) = plasma%kep(i)
		CALL ReinjectParticles(i,plasma,params,ecnt2,pcnt2)
		ELSE IF (plasma%zp(i) .LE. params%zmin) THEN
		plasma%f1(i) = 1
		plasma%dE1(i) = plasma%kep(i)
		CALL ReinjectParticles(i,plasma,params,ecnt1,pcnt1)
		END IF

		! Apply Coulomb collision operator:
		! ------------------------------------------------------------------------
		if (params%iColl) CALL collisionOperator(i,plasma,ecnt4,pcnt4,params)
		IF (params%iColl) CALL collisionOperator(i,plasma,ecnt4,pcnt4,params)

		! Apply RF heating operator:
		! ------------------------------------------------------------------------
		if (params%iHeat) then
		IF (params%iHeat) THEN
		CALL CyclotronResonanceNumber(i,plasma,resNum1,fieldspline,params)
		dresNum = dsign(1.d0,resNum0*resNum1)
		if (dresNum .LT. 0 .AND. plasma%zp(i) .GT. params%zRes1 .AND. plasma%zp(i) .LT. params%zRes2) then
		IF (dresNum .LT. 0 .AND. plasma%zp(i) .GT. params%zRes1 .AND. plasma%zp(i) .LT. params%zRes2) THEN
		plasma%f3(i) = 1
		CALL RFHeatingOperator(i,plasma,ecnt3,pcnt3,fieldspline,params)
		end if
		end if
		END IF
		END IF

		end do AllParticles
		END DO AllParticles
		!$OMP END DO

		!$OMP CRITICAL
		@@ -279,14 +231,14 @@ AllTime: do j = 1,params%NS
		! =====================================================================
		! Check if data is to be saved
		if (params%iSave) then
		if (j .EQ. jrng(k)) then
		t_hist(k) = tp
		if (j .EQ. output%jrng(k)) then
		output%tp(k) = tp
		!$OMP PARALLEL DO
		do i = 1,params%NC
		! Record "ith" particle at "kth" time
		zp_hist(i,k) = plasma%zp(i)
		kep_hist(i,k) = plasma%kep(i)
		xip_hist(i,k) = plasma%xip(i)
		output%zp(i,k) = plasma%zp(i)
		output%kep(i,k) = plasma%kep(i)
		output%xip(i,k) = plasma%xip(i)
		end do
		!$OMP END PARALLEL DO
		k = k + 1
		@@ -348,28 +300,28 @@ if (params%iSave) then
		! --------------------------------------------------------------------------
		fileName = trim(trim(dir1)//'/'//'zp.out')
		OPEN(unit=8,file=fileName,form="unformatted",status="unknown")
		WRITE(8) zp_hist
		WRITE(8) output%zp
		CLOSE(unit=8)

		! Saving kep_hist to file:
		! --------------------------------------------------------------------------
		fileName = trim(trim(dir1)//'/'//'kep.out')
		OPEN(unit=8,file=fileName,form="unformatted",status="unknown")
		WRITE(8) kep_hist
		WRITE(8) output%kep
		CLOSE(unit=8)

		! Saving xip_hist to file:
		! --------------------------------------------------------------------------
		fileName = trim(trim(dir1)//'/'//'xip.out')
		OPEN(unit=8,file=fileName,form="unformatted",status="unknown")
		WRITE(8) xip_hist
		WRITE(8) output%xip
		CLOSE(unit=8)

		! Saving t_hist to file:
		! --------------------------------------------------------------------------
		fileName = trim(trim(dir1)//'/'//'tp.out')
		OPEN(unit=8,file=fileName,form="unformatted",status="unknown")
		WRITE(8) t_hist
		WRITE(8) output%tp
		CLOSE(unit=8)

		! Saving pcount to file: