Commit 5f6a289f authored by Juan Caneses Marin (nfc)'s avatar Juan Caneses Marin (nfc)
Browse files

Demonstrated Electromagnetic field Interpolation for B, dB, ddB

parent 40f270e3

InputFiles/xp_TestDrag.in

+6 −6
Original line number Diff line number Diff line 
&params_nml

! Simulation name:

! ===============

params%fileDescriptor = '2021_03_11h',

params%fileDescriptor = '2021_03_15b',



! Magnetic field input data:

! =========================

params%repoDir       ="/home/nfc/myRepos/LinearFokkerPlanck_Axisymmetric",

params%BFieldFileDir = "/BfieldData",

params%BFieldFile    = "/Bfield_a.txt",

params%BFieldFile    = "/Bfield_b.txt",

params%nz            = 501,                                                     ! Number of points in BFieldFile



! Simulation conditions:

! =====================

params%NC = 90000,                                                              ! Total number of particles

params%NS = 50000,                                                              ! Total number of time steps

params%NC = 200000,                                                              ! Total number of particles

params%NS = 20000,                                                              ! Total number of time steps

params%dt = 0.5E-7,                                                             ! Time step in [s]

params%G  = 2.5E+20,                                                            ! Real particle injection rate 



! Time steps to record:

! ====================

params%jstart = 1,                                                              ! start frame

params%jend   = 50000,                                                          ! end frame

params%jend   = 20000,                                                          ! end frame

params%jincr  = 500,                                                            ! increment



! Domain geometry:
@@ -57,7 +57,7 @@ params%CollOperType = 2, 
params%BC_Type    = 2,                                                          ! 1: Isotropic plasma source, 2: NBI, 3: Periodic

params%BC_zp_mean = 0.0,                                                        ! Mean particle injection

params%BC_zp_std  = 0.3,                                                        ! STD of particle injection

params%BC_Ep      = 250.,                                                      ! Drift kinetic energy

params%BC_Ep      = 2000.,                                                      ! Drift kinetic energy

params%BC_Tp      = 10,                                                         ! Thermal kinetic energy

params%BC_xip     = 0.707,                                                      ! Mean pitch angle where xip = cos(theta) = vpar/v

src/Modules.f90

+74 −19
Original line number Diff line number Diff line
@@ -275,7 +275,7 @@ END TYPE plasmaTYP 


! -----------------------------------------------------------------------------

TYPE fieldSplineTYP

 TYPE(splTYP) :: B, dB, ddB, V, dV

 TYPE(splTYP) :: B, dB, ddB, E

END TYPE fieldSplineTYP



! -----------------------------------------------------------------------------
@@ -284,6 +284,7 @@ TYPE outputTYP 
 REAL(r8)   , DIMENSION(:,:), ALLOCATABLE :: zp, kep, xip, a

 INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: m

 REAL(r8)   , DIMENSION(:,:), ALLOCATABLE :: np, Up, Tparp, Tperp

 REAL(r8)   , DIMENSION(:,:), ALLOCATABLE :: Ep, Bp, dBp, ddBp

 REAL(r8)   , DIMENSION(:)  , ALLOCATABLE :: tp, jrng

! Global quantities:

 REAL(r8)   , DIMENSION(:)  , ALLOCATABLE :: NR, NSP, ER
@@ -355,6 +356,8 @@ SUBROUTINE InitializeMesh(mesh,params,fieldspline) 
   ! Declare local variables:

   INTEGER(i4), DIMENSION(params%NZmesh) :: m

   INTEGER(i4) :: i

   INTEGER(i4), DIMENSION(params%NZmesh+4) :: mg

   REAL(r8), DIMENSION(params%NZmesh+4) :: zmg



   ! Populate fields:

   mesh%NZmesh   = params%NZmesh
@@ -376,16 +379,34 @@ SUBROUTINE InitializeMesh(mesh,params,fieldspline) 
   mesh%nUE  = 0.

   mesh%P11  = 0.

   mesh%P22  = 0.

   mesh%B    = 0. ! Needs to be interpolated from fieldspline%B 

   mesh%B    = 0. 

   mesh%E    = 0.

   mesh%dB   = 0. ! Needs to be interpolated from fieldspline%dB

   mesh%ddB  = 0. ! Needs to be interpolated from fieldspline%ddB

   mesh%dB   = 0. 

   mesh%ddB  = 0. 

   mesh%U    = 0.

   mesh%Ppar = 0.

   mesh%Pper = 0.

   mesh%Tpar = 0.

   mesh%Tper = 0.



   ! Populate E, B, dB and ddB with fieldspline data:

   mg = (/ (i, i=1,mesh%NZmesh+4, 1) /)

   zmg = (mg-1)*mesh%dzm + 0.5*mesh%dzm + mesh%zmin - 2.*mesh%dzm

   DO i = 1,SIZE(zmg)

	CALL Interp1(zmg(i),mesh%B(i)  ,fieldspline%B  )

	CALL Interp1(zmg(i),mesh%dB(i) ,fieldspline%dB )

	CALL Interp1(zmg(i),mesh%ddB(i),fieldspline%ddB)

   END DO

   

   ! Output data to test:

   IF (.TRUE.) THEN

	   OPEN(unit=8,file="meshB.txt",form="formatted",status="unknown")

	   DO i = 1,SIZE(zmg)

	 	 WRITE(8,*) zmg(i), mesh%B(i), mesh%dB(i), mesh%ddB(i)

	   END DO

	   CLOSE(unit=8) 

   END IF



END SUBROUTINE InitializeMesh



! ----------------------------------------------------------------------------
@@ -415,21 +436,26 @@ SUBROUTINE AllocatePlasma(plasma,params) 
END SUBROUTINE AllocatePlasma



! --------------------------------------------------------------------------

SUBROUTINE InitializePlasma(plasma,params)

SUBROUTINE InitializePlasma(plasma,mesh,params)

   USE PhysicalConstants

   USE OMP_LIB



   IMPLICIT NONE



   ! Declare interface variables:

   TYPE(plasmaTYP), INTENT(INOUT) :: plasma

   TYPE(meshTYP)  , INTENT(IN)    :: mesh

   TYPE(paramsTYP), INTENT(INOUT) :: params

   

   ! Declare local variables:

   INTEGER(i4) :: i



   ! Derived parameters: Reference cross sectional area

   ! 1- Initialize Parametes:

   ! ===================================================================

   ! Reference cross sectional area

   params%Area0 = 0.5*params%dtheta*( params%r2**2. - params%r1**2.)



   ! Derived parameters: Select test species

   ! Select test species

   IF (params%species_a .EQ. 1) THEN

       params%qa = -e_c

       params%Ma = m_e
@@ -438,13 +464,13 @@ SUBROUTINE InitializePlasma(plasma,params) 
       params%Ma = params%Aion*m_p

   END IF

   

   ! Initialize plasma: scalar quantities:

   ! 2- Initialize plasma:

   ! ===================================================================

   ! 2.1- Global quantities:

   plasma%tp      = 0.

   plasma%NR      = params%ne0*params%Area0*(params%zmax - params%zmin)

   plasma%NSP     = params%NC

   plasma%alpha   = plasma%NR/plasma%NSP

   plasma%Eplus   = 0.

   plasma%Eminus  = 0.

   plasma%Ndot1   = 0.

   plasma%Ndot2   = 0.

   plasma%Ndot3   = 0.
@@ -455,27 +481,40 @@ SUBROUTINE InitializePlasma(plasma,params) 
   plasma%Edot4   = 0.

   plasma%uEdot3  = 0.



   ! 2.2- Initialize plasma: zp, kep, xip and a:

   !$OMP PARALLEL DO

   ! Initialize plasma: zp, kep, xip and a:

   DO i = 1,params%NC

     plasma%a(i)    = 1.

     CALL loadParticles(i,plasma,params)

   END DO

   !$OMP END PARALLEL DO



   ! 2.3- Assign particles to cells:

   CALL AssignCell(plasma,mesh,params)



   ! 2.4- Initialize moments:

   ! CALL InterpolateMoments(plasma,mesh,params)

   !$OMP PARALLEL DO

   DO i = 1,params%NC

     plasma%np(i)    = 0.

     plasma%Up(i)    = 0.

     plasma%Tparp(i) = 0.

     plasma%Tperp(i) = 0.

     plasma%Bp(i)    = 0. ! Needs to be PIC interpolated from mesh%B

     plasma%Ep(i)    = 0.

     plasma%dBp(i)   = 0. ! Needs to be PIC interpolated from mesh%dB

     plasma%ddBp(i)  = 0. ! Needs to be PIC interpolated from mesh%ddB

   END DO

   !$OMP END PARALLEL DO



   ! 2.5- Initialize electromagnetic field via interpolation: 

   CALL InterpolateElectromagneticFields(plasma,mesh,params)



   ! Test interpolate EM fields:

   IF (.TRUE.) THEN

         OPEN(unit=8,file="InterpB.txt",form="formatted",status="unknown")

         DO i = 1,params%NC

                 WRITE(8,*) plasma%zp(i), plasma%Bp(i), plasma%dBp(i), plasma%ddBp(i)

         END DO

         CLOSE(unit=8)

   END IF



END SUBROUTINE InitializePlasma



! --------------------------------------------------------------------------
@@ -520,8 +559,8 @@ SUBROUTINE AllocateFieldSpline(fieldspline,params) 
   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%E  ,NZ,0._8,0._8)

!   CALL AllocateSpline(fieldspline%dV ,NZ,0._8,0._8)



END SUBROUTINE AllocateFieldSpline
@@ -536,8 +575,8 @@ SUBROUTINE ComputeFieldSpline(fieldspline) 
   CALL ComputeSpline(fieldspline%B)

   CALL ComputeSpline(fieldspline%dB)

   CALL ComputeSpline(fieldspline%ddB)

   CALL ComputeSpline(fieldspline%V)

   CALL ComputeSpline(fieldspline%dV)

   CALL ComputeSpline(fieldspline%E)

!   CALL ComputeSpline(fieldspline%dV)



END SUBROUTINE ComputeFieldSpline
@@ -565,6 +604,10 @@ SUBROUTINE AllocateOutput(output,params) 
   ALLOCATE(output%Up(params%NC  ,jsize))

   ALLOCATE(output%Tparp(params%NC  ,jsize))

   ALLOCATE(output%Tperp(params%NC  ,jsize))

   ALLOCATE(output%Ep(params%NC  ,jsize))

   ALLOCATE(output%Bp(params%NC  ,jsize))

   ALLOCATE(output%dBp(params%NC  ,jsize))

   ALLOCATE(output%ddBp(params%NC  ,jsize))

   ALLOCATE(output%tp(jsize))



   ! Create array with the indices of the time steps to save:
@@ -688,6 +731,18 @@ SUBROUTINE SaveData(output,dir1) 
    OPEN(unit=8,file=fileName,form="unformatted",status="unknown")

    WRITE(8) output%Tperp

    CLOSE(unit=8)

    fileName = trim(trim(dir1)//'/'//'Bp.out')

    OPEN(unit=8,file=fileName,form="unformatted",status="unknown")

    WRITE(8) output%Bp

    CLOSE(unit=8)

    fileName = trim(trim(dir1)//'/'//'dBp.out')

    OPEN(unit=8,file=fileName,form="unformatted",status="unknown")

    WRITE(8) output%dBp

    CLOSE(unit=8)

    fileName = trim(trim(dir1)//'/'//'ddBp.out')

    OPEN(unit=8,file=fileName,form="unformatted",status="unknown")

    WRITE(8) output%ddBp

    CLOSE(unit=8)



    ! Saving mesh-defined quantities:

    ! -------------------------------------------------------------------------

src/MoveParticlePack.f90

+90 −12
Original line number Diff line number Diff line 
! =======================================================================================================

SUBROUTINE AdvanceParticles(plasma,fieldspline,params)

SUBROUTINE AdvanceParticles(plasma,mesh,fieldspline,params)

! =======================================================================================================

USE local

USE dataTYP
@@ -9,6 +9,7 @@ IMPLICIT NONE 


! Declare interface variables:

TYPE(plasmaTYP)     , INTENT(INOUT) :: plasma

TYPE(meshTYP)       , INTENT(IN)    :: mesh

TYPE(paramsTYP)     , INTENT(IN)    :: params

TYPE(fieldSplineTYP), INTENT(IN)    :: fieldspline
@@ -21,9 +22,9 @@ dresNum = 0. 
resNum0 = 0.

resNum1 = 0.



! Advance particle positions and velocities:

!$OMP PARALLEL DO FIRSTPRIVATE(dresNum, resNum0, resNum1) SCHEDULE(STATIC) 

DO i = 1,params%NC



	! 2.1 - Reset flags:

	CALL ResetFlags(i,plasma)
@@ -61,6 +62,84 @@ END DO 
RETURN

END SUBROUTINE AdvanceParticles



! ======================================================================================================

SUBROUTINE InterpolateElectromagneticFields(plasma,mesh,params)

! ======================================================================================================

USE LOCAL

USE dataTYP



IMPLICIT NONE



! Declare interface variables:

TYPE(plasmaTYP), INTENT(INOUT) :: plasma

TYPE(meshTYP)  , INTENT(IN)    :: mesh

TYPE(paramsTYP), INTENT(IN)    :: params



! Declare local variables:

INTEGER(i4) :: i, ix

REAL(r8), DIMENSION(3) ::  w, E, B, dB, ddB



!$OMP PARALLEL DO PRIVATE(ix, w, E, B, dB, ddB)

DO i = 1,params%NC

	IF (plasma%f1(i) .EQ. 0 .AND. plasma%f2(i) .EQ. 0) THEN

                ! Get nearest grid point:

                ix = plasma%m(i) + 2



                ! Assignment function:

                w(1) = plasma%wL(i)

                w(2) = plasma%wC(i)

                w(3) = plasma%wR(i)



                ! Electric field:

                E(1) = mesh%E(ix - 1)

                E(2) = mesh%E(ix)

                E(3) = mesh%E(ix + 1)

                plasma%Ep(i) = w(1)*E(1) + w(2)*E(2) + w(3)*E(3)



                ! Magnetic field:

                B(1) = mesh%B(ix - 1)

                B(2) = mesh%B(ix)

                B(3) = mesh%B(ix + 1)

                plasma%Bp(i) = w(1)*B(1) + w(2)*B(2) + w(3)*B(3)



                ! Magnetic field: 1st derivative

                dB(1) = mesh%dB(ix - 1)

                dB(2) = mesh%dB(ix)

                dB(3) = mesh%dB(ix + 1)

                plasma%dBp(i) = w(1)*dB(1) + w(2)*dB(2) + w(3)*dB(3)



                ! Magnetic field: 2nd derivative

		!IF (params%iHeat) THEN

		IF (.TRUE.) THEN

	                ddB(1) = mesh%ddB(ix - 1)

        	        ddB(2) = mesh%ddB(ix)

	        	ddB(3) = mesh%ddB(ix + 1)

                	plasma%ddBp(i) = w(1)*ddB(1) + w(2)*ddB(2) + w(3)*ddB(3)

		END IF

	END IF

END DO

!$OMP END PARALLEL DO



END SUBROUTINE InterpolateElectromagneticFields



! =======================================================================================================

SUBROUTINE InterpolateMoments(plasma,mesh,params)

! =======================================================================================================

USE LOCAL

USE dataTYP



IMPLICIT NONE



! Declare interface variables:

TYPE(plasmaTYP), INTENT(INOUT) :: plasma

TYPE(meshTYP)  , INTENT(INOUT) :: mesh

TYPE(paramsTYP), INTENT(IN)    :: params



! Declare local variables:

INTEGER(i4) :: i



END SUBROUTINE InterpolateMoments



! =======================================================================================================

SUBROUTINE CheckBoundary(i,plasma,params)

! =======================================================================================================
@@ -198,7 +277,7 @@ TYPE(paramsTYP) , INTENT(IN) :: params 
TYPE(fieldSplineTYP), INTENT(IN) :: fieldspline



! Local variables:

REAL(r8) :: dB, dV

REAL(r8) :: dB, E

REAL(r8) :: Ma, qa



! Test particle mass:
@@ -210,13 +289,13 @@ qa = params%qa 
! Calculate the magnetic field gradient at zp0:

CALL Interp1(zp0,dB,fieldspline%dB)



! Calculate electric potential gradient at zp0:

CALL Interp1(zp0,dV,fieldspline%dV)

! Calculate electric field at zp0:

CALL Interp1(zp0,E,fieldspline%E)



! Assign values to output variables:

K = upar0

L = -(1/Ma)*(mu0*dB + qa*dV)

M = 0

L = (-mu0*dB + qa*E)/Ma

M = 0.



RETURN

END SUBROUTINE RightHandSide
@@ -293,7 +372,6 @@ unU = 0. 
P11 = 0.

P22 = 0.

nUE = 0.



! Species "a" mass:

Ma = params%Ma

! Scaling factor:
@@ -587,7 +665,7 @@ TYPE(fieldSplineTYP) , INTENT(IN) :: fieldspline 
REAL(r8) :: zp0, kep0, xip0

REAL(r8) :: u0, upar0, uper0

REAL(r8) :: kep_par0, kep_per0

REAL(r8) :: dB, ddB, dV

REAL(r8) :: dB, ddB, E

REAL(r8) :: Bf, Omega, dOmega, ddOmega

REAL(r8) :: Omega_dot, Omega_ddot, tau_rf

REAL(r8) :: rl, flr, besselterm
@@ -615,7 +693,7 @@ kep_per0 = kep0*(1. - xip0**2.) 
! Gradients:

CALL Interp1(zp0,dB ,fieldspline%dB )

CALL Interp1(zp0,ddB,fieldspline%ddB)

CALL Interp1(zp0,dV ,fieldspline%dV )

CALL Interp1(zp0,E  ,fieldspline%E  )



! Spatial derivatives of the magnetic field:

CALL Interp1(zp0,Bf,fieldspline%B)
@@ -625,7 +703,7 @@ ddOmega = params%n_harmonic*e_c*ddB/Ma 


! Calculate the first and second time derivative of Omega:

Omega_dot = upar0*dOmega

Omega_ddot = (upar0**2.)*ddOmega  - (uper0**2.)*dOmega*dOmega/(2.*Omega) - qa*dV*dOmega/Ma

Omega_ddot = (upar0**2.)*ddOmega  - (uper0**2.)*dOmega*dOmega/(2.*Omega) + qa*E*dOmega/Ma



! Calculate the interaction time (tau_RF):

IF ( (Omega_ddot**2.) .GT. 4.8175*ABS(Omega_dot**3.) )  then

src/linFP.f90

+17 −10
Original line number Diff line number Diff line
@@ -79,22 +79,20 @@ CALL ReadSpline(fieldspline%B,fileName) 
CALL diffSpline(fieldspline%B ,fieldspline%dB )

CALL diffSpline(fieldspline%dB,fieldspline%ddB)



! 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)

CALL diffSpline(fieldspline%V,fieldspline%dV)

! Electric field E:

fieldspline%E%x = fieldspline%B%x

fieldspline%E%y = 0

CALL ComputeSpline(fieldspline%E)



! Complete setting up the spline data :

CALL ComputeFieldSpline(fieldspline)



! Initialize simulation variables:

! ==============================================================================

WRITE(*,*) "Initialize Mesh"

CALL InitializeMesh(mesh,params,fieldspline)

CALL InitializePlasma(plasma,params)

WRITE(*,*) "Initialize Plasma"

CALL InitializePlasma(plasma,mesh,params)



!$OMP PARALLEL

IF (OMP_GET_THREAD_NUM() .EQ. 0) params%threads_given = OMP_GET_NUM_THREADS()
@@ -130,7 +128,7 @@ NS_loop: DO j = 1,params%NS 
    dresNum = 0.; resNum0 = 0.; resNum1 = 0.



    ! Advance particles:

    CALL AdvanceParticles(plasma,fieldspline,params)

    CALL AdvanceParticles(plasma,mesh,fieldspline,params)



    ! Assign charge to cells:

    CALL AssignCell(plasma,mesh,params)
@@ -186,6 +184,12 @@ NS_loop: DO j = 1,params%NS 
    END DO NC_loop4

    !$OMP END PARALLEL DO



    ! Field solve:

    ! AdvanceEfield(plasma,mesh,params)



    ! Interpolate electromagnetic fields to particle positions:

    CALL InterpolateElectromagneticFields(plasma,mesh,params)



    ! Calculate new number of real particles NR and super-particles NSP:

    ! ==============================================================================

    dNR  = a_new*(uN1 + uN2)*dt - (N1 + N2)*dt
@@ -230,6 +234,9 @@ NS_loop: DO j = 1,params%NS 
             output%Up(i,k)    = plasma%Up(i)

             output%Tparp(i,k) = plasma%Tparp(i)

             output%Tperp(i,k) = plasma%Tperp(i)

             output%Bp(i,k)    = plasma%Bp(i)

             output%dBp(i,k)   = plasma%dBp(i)

             output%ddBp(i,k)  = plasma%ddBp(i)

          END DO

         !$OMP END PARALLEL DO

src/make_linFP.f

+2 −5
Original line number Diff line number Diff line
@@ -2,7 +2,7 @@ 
COMPILER = gfortran

OPTFLAGS = -O3

DBGFLAGS = -g

OBJ_1 = Modules.o linFP.o PotentialProfile.o

OBJ_1 = Modules.o linFP.o 

OBJ_2 = fitpack.o MoveParticlePack.o CoulombCollisions.o



All: $(OBJ_1) $(OBJ_2)
@@ -10,14 +10,11 @@ All: $(OBJ_1) $(OBJ_2) 
	rm *.o *.mod



Modules.o: Modules.f90

	$(COMPILER) $(OPTFLAGS) -c Modules.f90

	$(COMPILER) $(OPTFLAGS) -fopenmp -c Modules.f90



linFP.o: linFP.f90

	$(COMPILER) $(OPTFLAGS) -fopenmp -c linFP.f90



PotentialProfile.o: PotentialProfile.f90

	$(COMPILER) $(OPTFLAGS) -c PotentialProfile.f90



MoveParticlePack.o: MoveParticlePack.f90

	$(COMPILER) $(OPTFLAGS) -c -fopenmp MoveParticlePack.f90