Loading InputFiles/xp_TestDrag.in +6 −6 Original line number Diff line number Diff line ¶ms_nml ! Simulation name: ! =============== params%fileDescriptor = '2021_03_08a', params%fileDescriptor = '2021_03_09b', ! Magnetic field input data: ! ========================= Loading @@ -12,16 +12,16 @@ params%nz = 501, ! Simulation conditions: ! ===================== params%NC = 50000, ! Total number of particles params%NS = 90000, ! Total number of time steps params%NC = 90000, ! Total number of particles params%NS = 50000, ! 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 = 90000, ! end frame params%jincr = 600, ! increment params%jend = 50000, ! end frame params%jincr = 500, ! increment ! Domain geometry: ! ==================== Loading Loading @@ -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 = 5000., ! 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 Loading src/Modules.f90 +194 −45 Original line number Diff line number Diff line Loading @@ -259,14 +259,14 @@ END TYPE paramsTYP ! ----------------------------------------------------------------------------- TYPE plasmaTYP REAL(r8) , DIMENSION(:), ALLOCATABLE :: zp, kep,xip, a INTEGER(i4), DIMENSION(:), ALLOCATABLE :: f1, f2, f3, f4 REAL(r8) , DIMENSION(:), ALLOCATABLE :: wL, wC, wR INTEGER(i4), DIMENSION(:), ALLOCATABLE :: m, f1, f2, f3, f4 REAL(r8) , DIMENSION(:), ALLOCATABLE :: dE1, dE2, dE3, dE4, dE5 REAL(r8) , DIMENSION(:), ALLOCATABLE :: udE3, udErf, doppler REAL(r8) :: NR , NSP, alpha REAL(r8) :: NR , NSP, alpha, tp REAL(r8) :: Eplus, Eminus REAL(r8) :: Ndot1, Ndot2, Ndot3, Ndot4 REAL(r8) :: Edot1, Edot2, Edot3, Edot4, uEdot3 REAL(r8) , DIMENSION(:), ALLOCATABLE :: n, nU, unU, U END TYPE plasmaTYP ! ----------------------------------------------------------------------------- Loading @@ -277,42 +277,29 @@ END TYPE fieldSplineTYP ! ----------------------------------------------------------------------------- TYPE outputTYP REAL(r8) , DIMENSION(:,:), ALLOCATABLE :: zp, kep, xip, a INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: m REAL(r8) , DIMENSION(:) , ALLOCATABLE :: tp, jrng REAL(r8) , DIMENSION(:) , ALLOCATABLE :: NR, NSP, ER REAL(r8) , DIMENSION(:) , ALLOCATABLE :: Eplus, Eminus REAL(r8) , DIMENSION(:) , ALLOCATABLE :: Ndot1, Ndot2, Ndot3, Ndot4, Ndot5 REAL(r8) , DIMENSION(:) , ALLOCATABLE :: Edot1, Edot2, Edot3, Edot4, Edot5 REAL(r8) , DIMENSION(:,:), ALLOCATABLE :: n, nU, unU, nUE, P11, P22, E, B, dB, ddB REAL(r8) , DIMENSION(:) , ALLOCATABLE :: zm END TYPE outputTYP ! ----------------------------------------------------------------------------- TYPE meshTYP REAL(r8) , DIMENSION(:), ALLOCATABLE :: zm REAL(r8) :: LZ, zmin, zmax, dz INTEGER(i4) :: NZ REAL(r8) , DIMENSION(:), ALLOCATABLE :: n, nU, nUE, P11, P22 REAL(r8) , DIMENSION(:), ALLOCATABLE :: E, B, dB, ddB REAL(r8) , DIMENSION(:), ALLOCATABLE :: unU REAL(r8) :: LZ, zmin, zmax, dzm INTEGER(i4) :: NZmesh END TYPE meshTYP TYPE fieldsTYP REAL(r8) , DIMENSION(:), ALLOCATABLE :: E, B REAL(r8) , DIMENSION(:), ALLOCATABLE :: V, dB, ddB END TYPE fieldsTYP CONTAINS ! ---------------------------------------------------------------------------- SUBROUTINE AllocateFields(fields,params) IMPLICIT NONE ! Declare interface variables: TYPE(fieldsTYP), INTENT(INOUT) :: fields TYPE(paramsTYP), INTENT(IN) :: params ! Declare local variables: INTEGER(i4) :: NZ ! Allocate memory to mesh: NZ = params%NZmesh ALLOCATE(fields%E(NZ),fields%B(NZ),fields%dB(NZ),fields%ddB(NZ)) END SUBROUTINE AllocateFields ! ---------------------------------------------------------------------------- SUBROUTINE AllocateMesh(mesh,params) IMPLICIT NONE Loading @@ -321,11 +308,25 @@ SUBROUTINE AllocateMesh(mesh,params) TYPE(paramsTYP), INTENT(IN) :: params ! Declare local variables: INTEGER(i4) :: NZ INTEGER(i4) :: NZmesh ! Allocate memory to mesh: NZ = params%NZmesh ALLOCATE(mesh%zm(NZ)) NZmesh = params%NZmesh ALLOCATE(mesh%zm(NZmesh)) ! Allocate memory for mesh defined quantities: ! Ghost cells are used: ALLOCATE(mesh%n(NZmesh + 4)) ALLOCATE(mesh%nU(NZmesh + 4)) ALLOCATE(mesh%unU(NZmesh + 4)) ALLOCATE(mesh%P11(NZmesh + 4)) ALLOCATE(mesh%P22(NZmesh + 4)) ALLOCATE(mesh%nuE(NZmesh + 4)) ALLOCATE(mesh%B(NZmesh + 4)) ALLOCATE(mesh%dB(NZmesh + 4)) ALLOCATE(mesh%ddB(NZmesh + 4)) ALLOCATE(mesh%E(NZmesh + 4)) END SUBROUTINE AllocateMesh ! ---------------------------------------------------------------------------- Loading @@ -336,21 +337,33 @@ SUBROUTINE InitializeMesh(mesh,params) TYPE(paramsTYP), INTENT(IN) :: params ! Declare local variables: INTEGER(i4), DIMENSION(params%NZ) :: m INTEGER(i4), DIMENSION(params%NZmesh) :: m INTEGER(i4) :: i ! Populate fields: mesh%NZ = params%NZmesh mesh%NZmesh = params%NZmesh mesh%zmin = params%zmin mesh%zmax = params%zmax ! Derived quantities: mesh%LZ = params%zmax - params%zmin mesh%dz = mesh%LZ/mesh%NZ m = (/ (i, i=1,mesh%NZ, 1) /) mesh%dzm = mesh%LZ/mesh%NZmesh m = (/ (i, i=1,mesh%NZmesh, 1) /) ! Create mesh: mesh%zm = (m-1)*mesh%dz + 0.5*mesh%dz + mesh%zmin mesh%zm = (m-1)*mesh%dzm + 0.5*mesh%dzm + mesh%zmin ! Initialize all mesh-defined quantities: mesh%n = 0. mesh%nU = 0. mesh%unU = 0. mesh%nUE = 0. mesh%P11 = 0. mesh%P22 = 0. mesh%B = 0. mesh%E = 0. mesh%dB = 0. mesh%ddB = 0. END SUBROUTINE InitializeMesh Loading @@ -367,8 +380,9 @@ SUBROUTINE AllocatePlasma(plasma,params) ! NC: Number of computational particles NC = params%NC ! Allocate memory: For all computational particles ! Allocate memory: for all computational particles ALLOCATE(plasma%zp(NC) ,plasma%kep(NC) ,plasma%xip(NC) ,plasma%a(NC)) ALLOCATE(plasma%m(NC) ,plasma%wL(NC) ,plasma%wC(NC) ,plasma%wR(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), plasma%dE5(NC)) ALLOCATE(plasma%udErf(NC)) Loading Loading @@ -402,6 +416,7 @@ SUBROUTINE InitializePlasma(plasma,params) END IF ! Initialize plasma: scalar quantities: plasma%tp = 0. plasma%NR = params%ne0*params%Area0*(params%zmax - params%zmin) plasma%NSP = params%NC plasma%alpha = plasma%NR/plasma%NSP Loading Loading @@ -498,7 +513,7 @@ SUBROUTINE AllocateOutput(output,params) TYPE(paramsTYP), INTENT(IN) :: params ! Declare local variables: INTEGER(i4) :: j, jsize, NS INTEGER(i4) :: j, jsize, NS, NZmesh ! Determine size of temporal snapshots to record: jsize = (params%jend-params%jstart+1)/params%jincr Loading @@ -509,6 +524,7 @@ SUBROUTINE AllocateOutput(output,params) ALLOCATE(output%kep(params%NC,jsize)) ALLOCATE(output%xip(params%NC,jsize)) ALLOCATE(output%a(params%NC ,jsize)) ALLOCATE(output%m(params%NC ,jsize)) ALLOCATE(output%tp(jsize)) ! Create array with the indices of the time steps to save: Loading @@ -520,9 +536,23 @@ SUBROUTINE AllocateOutput(output,params) ALLOCATE(output%Ndot1(NS),output%Ndot2(NS),output%Ndot3(NS),output%Ndot4(NS) ,output%Ndot5(NS)) ALLOCATE(output%Edot1(NS),output%Edot2(NS),output%Edot3(NS),output%Edot4(NS) ,output%Edot5(NS)) ! Allocate memory for mesh-defined quantities: NZmesh = params%NZmesh ALLOCATE(output%zm(NZmesh)) ALLOCATE(output%n(NZmesh + 4 ,jsize)) ALLOCATE(output%nU(NZmesh + 4 ,jsize)) ALLOCATE(output%unU(NZmesh + 4,jsize)) ALLOCATE(output%nUE(NZmesh + 4,jsize)) ALLOCATE(output%P11(NZmesh + 4,jsize)) ALLOCATE(output%P22(NZmesh + 4,jsize)) ALLOCATE(output%B(NZmesh + 4 ,jsize)) ALLOCATE(output%E(NZmesh +4 ,jsize)) ALLOCATE(output%dB(NZmesh + 4 ,jsize)) ALLOCATE(output%ddB(NZmesh + 4,jsize)) END SUBROUTINE AllocateOutput ! --------------------------------------------------------------------------- ! ================================================================================= SUBROUTINE PrintParamsToTerminal(params,inputFile) IMPLICIT NONE ! Declare interface variables: Loading Loading @@ -558,4 +588,123 @@ SUBROUTINE PrintParamsToTerminal(params,inputFile) END SUBROUTINE PrintParamsToTerminal ! ================================================================================= SUBROUTINE SaveData(output,dir1) IMPLICIT NONE !Declare interface variables: TYPE(outputTYP), INTENT(IN) :: output CHARACTER*300 , INTENT(IN) :: dir1 ! Declare local variables: CHARACTER*300 :: fileName WRITE(*,*) "Saving data ..." ! Save plasma quantities: ! -------------------------------------------------------------------------- fileName = TRIM(TRIM(dir1)//'/'//'zp.out') OPEN(UNIT=8,FILE=fileName,FORM="unformatted",STATUS="unknown") WRITE(8) output%zp CLOSE(UNIT=8) fileName = trim(trim(dir1)//'/'//'kep.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%kep CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'xip.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%xip CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'a.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%a CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'tp.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%tp CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'m.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%m CLOSE(unit=8) ! Saving mesh-defined quantities: ! ------------------------------------------------------------------------- fileName = trim(trim(dir1)//'/'//'z_mesh.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%zm CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'n_mesh.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%n CLOSE(unit=8) ! Saving pcount to file: ! -------------------------------------------------------------------------- fileName = trim(trim(dir1)//'/'//'pcount1.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Ndot1 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'pcount2.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Ndot2 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'pcount3.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Ndot3 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'pcount4.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Ndot4 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'pcount5.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Ndot5 CLOSE(unit=8) ! Saving ecount to file: ! -------------------------------------------------------------------------- fileName = trim(trim(dir1)//'/'//'ecount1.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Edot1 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'ecount2.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Edot2 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'ecount3.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Edot3 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'ecount4.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Edot4 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'ecount5.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Edot5 CLOSE(unit=8) ! Save ER, NR and NSP: ! ------------------------------------------------------------------------- fileName = trim(trim(dir1)//'/'//'ER.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%ER CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'NR.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%NR CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'NSP.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%NSP CLOSE(unit=8) WRITE(*,*) "Data saving complete" END SUBROUTINE SaveData END MODULE dataTYP src/MoveParticlePack.f90 +210 −0 Original line number Diff line number Diff line ! ======================================================================================================= SUBROUTINE AdvanceParticles(plasma,fieldspline,params) ! ======================================================================================================= USE local USE dataTYP USE OMP_LIB IMPLICIT NONE ! Declare interface variables: TYPE(plasmaTYP) , INTENT(INOUT) :: plasma TYPE(paramsTYP) , INTENT(IN) :: params TYPE(fieldSplineTYP), INTENT(IN) :: fieldspline ! Declare local variables: REAL(r8) :: dresNum, resNum0, resNum1 INTEGER(i4) :: i ! Initialize resonance number: dresNum = 0. resNum0 = 0. resNum1 = 0. !$OMP PARALLEL DO FIRSTPRIVATE(dresNum, resNum0, resNum1) SCHEDULE(STATIC) DO i = 1,params%NC ! 2.1 - Reset flags: CALL ResetFlags(i,plasma) ! 2.2- Compute resonance number: IF (params%iHeat) THEN CALL CyclotronResonanceNumber(i,plasma,fieldspline,params,resNum0) END IF ! 2.3- Advance zp, kep, xip: IF (params%iPush) THEN CALL MoveParticle(i,plasma,fieldspline,params) END IF ! 2.4- Check boundaries: set f1 and f2 flags IF (params%iPush) THEN CALL CheckBoundary(i,plasma,params) END IF ! Apply Coulomb collision operator: ! ------------------------------------------------------------------------ IF (params%iColl) THEN CALL collisionOperator(i,plasma,params) END IF ! 2.5- Compute resonance number: IF (params%iHeat) THEN CALL CyclotronResonanceNumber(i,plasma,fieldspline,params,resNum0) END IF ! 2.6- Check resonance: set f3 flag IF (params%iHeat) THEN 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 plasma%f3(i) = 1 CALL RFoperatorTerms(i,plasma,fieldspline,params) END IF END IF END DO !$OMP END PARALLEL DO RETURN END SUBROUTINE AdvanceParticles ! ======================================================================================================= SUBROUTINE CheckBoundary(i,plasma,params) ! ======================================================================================================= USE LOCAL USE dataTYP IMPLICIT NONE ! Declare interface variables: INTEGER(i4), INTENT(IN) :: i TYPE(plasmaTYP), INTENT(INOUT) :: plasma TYPE(paramsTYP), INTENT(IN) :: params IF (plasma%zp(i) .GE. params%zmax) THEN plasma%f2(i) = 1 plasma%dE2(i) = plasma%kep(i) ELSE IF (plasma%zp(i) .LE. params%zmin) THEN plasma%f1(i) = 1 plasma%dE1(i) = plasma%kep(i) END IF RETURN END SUBROUTINE CheckBoundary ! ======================================================================================================= SUBROUTINE MoveParticle(i,plasma,fieldspline,params) ! ======================================================================================================= Loading Loading @@ -133,6 +227,122 @@ M = 0 RETURN END SUBROUTINE RightHandSide ! ======================================================================================================= SUBROUTINE AssignCell(plasma,mesh,params) ! ======================================================================================================= USE LOCAL USE dataTYP USE OMP_LIB IMPLICIT NONE ! Define interface variables: TYPE(plasmaTYP), INTENT(INOUT) :: plasma TYPE(meshTYP) , INTENT(IN) :: mesh TYPE(paramsTYP), INTENT(IN) :: params ! Declare local variables: INTEGER(i4) :: i REAL(r8) :: Z, dz, zi, zoffset ! Mesh grid size: dz = mesh%dzm ! Mesh offset: zoffset = mesh%zmin !$OMP PARALLEL DO PRIVATE(zi,Z) DO i = 1,params%NC IF (plasma%f1(i) .EQ. 0 .AND. plasma%f2(i) .EQ. 0) THEN ! Using zm and zp, find nearest grid point (NGP): zi = plasma%zp(i) plasma%m(i) = NINT(0.5 + (zi - zoffset)/dz) ! Compute wL(i), wC(i), wR(i) Z = mesh%zm(plasma%m(i)) - zi plasma%wC(i) = 0.75 - (Z/dz)**2. plasma%wL(i) = 0.5*(1.5 + ((Z - dz)/dz) )**2. plasma%wR(i) = 0.5*(1.5 - ((Z + dz)/dz) )**2. END IF END DO !$OMP END PARALLEL DO RETURN END SUBROUTINE AssignCell ! ======================================================================================================= SUBROUTINE ExtrapolateMomentsToMesh(plasma,mesh,params) ! ======================================================================================================= USE LOCAL USE PhysicalConstants USE dataTYP USE OMP_LIB IMPLICIT NONE ! Define interface variables: TYPE(plasmaTYP), INTENT(IN) :: plasma TYPE(meshTYP) , INTENT(INOUT) :: mesh TYPE(paramsTYP), INTENT(IN) :: params ! Define local variables: INTEGER(i4) :: i, ix, frame REAL(r8) :: vpar, Ma, Ep, alpha, xip, vper, v, a ! Initialize mesh quantities: mesh%n = 0. mesh%nU = 0. mesh%unU = 0. mesh%P11 = 0. mesh%P22 = 0. mesh%nuE = 0. ! Species "a" mass: Ma = params%Ma ! Scaling factor: alpha = plasma%alpha ! Calculate moments and extrapolate to mesh points !$OMP PARALLEL DO private(ix,a,Ep,v,xip,vpar,vper) DO i = 1,params%NC IF (plasma%f1(i) .EQ. 0 .AND. plasma%f2(i) .EQ. 0) THEN ! Derived quantities: a = plasma%a(i) Ep = e_c*plasma%kep(i) v = sqrt(2*Ep/Ma) xip = plasma%xip(i) vpar = v*xip vper = v*sqrt(1 - xip**2.) ! Mesh point with ghost cells: ix = plasma%m(i) + 2 ! Plasma density: [NSP m^-3] mesh%n(ix-1) = mesh%n(ix-1) + plasma%wL(i)*a; mesh%n(ix) = mesh%n(ix) + plasma%wC(i)*a; mesh%n(ix+1) = mesh%n(ix+1) + plasma%wR(i)*a; END IF END DO !$OMP END PARALLEL DO ! Apply magnetic compression: mesh%n = mesh%n/params%Area0 ! Apply scaling factor: mesh%n = alpha*mesh%n/mesh%dzm ! Apply smoothing: frame = 9 ! CALL MovingMean(mesh%n ,frame) ! Calculate U, Tpar, Tper: END SUBROUTINE ExtrapolateMomentsToMesh ! ======================================================================================================= SUBROUTINE ReinjectParticles(i,plasma,params) ! ======================================================================================================= Loading src/linFP.f90 +27 −167 File changed.Preview size limit exceeded, changes collapsed. Show changes src/make_linFP.f +2 −2 Original line number Diff line number Diff line Loading @@ -6,7 +6,7 @@ OBJ_1 = Modules.o linFP.o PotentialProfile.o OBJ_2 = fitpack.o MoveParticlePack.o CoulombCollisions.o All: $(OBJ_1) $(OBJ_2) gfortran $(OPTFLAGS) -fopenmp $(OBJ_1) $(OBJ_2) -o linFP $(COMPILER) $(OPTFLAGS) -fopenmp $(OBJ_1) $(OBJ_2) -o linFP rm *.o *.mod Modules.o: Modules.f90 Loading @@ -22,7 +22,7 @@ fitpack.o: fitpack.f $(COMPILER) $(OPTFLAGS) -c -w fitpack.f MoveParticlePack.o: MoveParticlePack.f90 $(COMPILER) $(OPTFLAGS) -c MoveParticlePack.f90 $(COMPILER) $(OPTFLAGS) -c -fopenmp MoveParticlePack.f90 CoulombCollisions.o: CoulombCollisions.f90 $(COMPILER) $(OPTFLAGS) -fopenmp -c CoulombCollisions.f90 Loading Loading
InputFiles/xp_TestDrag.in +6 −6 Original line number Diff line number Diff line ¶ms_nml ! Simulation name: ! =============== params%fileDescriptor = '2021_03_08a', params%fileDescriptor = '2021_03_09b', ! Magnetic field input data: ! ========================= Loading @@ -12,16 +12,16 @@ params%nz = 501, ! Simulation conditions: ! ===================== params%NC = 50000, ! Total number of particles params%NS = 90000, ! Total number of time steps params%NC = 90000, ! Total number of particles params%NS = 50000, ! 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 = 90000, ! end frame params%jincr = 600, ! increment params%jend = 50000, ! end frame params%jincr = 500, ! increment ! Domain geometry: ! ==================== Loading Loading @@ -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 = 5000., ! 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 Loading
src/Modules.f90 +194 −45 Original line number Diff line number Diff line Loading @@ -259,14 +259,14 @@ END TYPE paramsTYP ! ----------------------------------------------------------------------------- TYPE plasmaTYP REAL(r8) , DIMENSION(:), ALLOCATABLE :: zp, kep,xip, a INTEGER(i4), DIMENSION(:), ALLOCATABLE :: f1, f2, f3, f4 REAL(r8) , DIMENSION(:), ALLOCATABLE :: wL, wC, wR INTEGER(i4), DIMENSION(:), ALLOCATABLE :: m, f1, f2, f3, f4 REAL(r8) , DIMENSION(:), ALLOCATABLE :: dE1, dE2, dE3, dE4, dE5 REAL(r8) , DIMENSION(:), ALLOCATABLE :: udE3, udErf, doppler REAL(r8) :: NR , NSP, alpha REAL(r8) :: NR , NSP, alpha, tp REAL(r8) :: Eplus, Eminus REAL(r8) :: Ndot1, Ndot2, Ndot3, Ndot4 REAL(r8) :: Edot1, Edot2, Edot3, Edot4, uEdot3 REAL(r8) , DIMENSION(:), ALLOCATABLE :: n, nU, unU, U END TYPE plasmaTYP ! ----------------------------------------------------------------------------- Loading @@ -277,42 +277,29 @@ END TYPE fieldSplineTYP ! ----------------------------------------------------------------------------- TYPE outputTYP REAL(r8) , DIMENSION(:,:), ALLOCATABLE :: zp, kep, xip, a INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: m REAL(r8) , DIMENSION(:) , ALLOCATABLE :: tp, jrng REAL(r8) , DIMENSION(:) , ALLOCATABLE :: NR, NSP, ER REAL(r8) , DIMENSION(:) , ALLOCATABLE :: Eplus, Eminus REAL(r8) , DIMENSION(:) , ALLOCATABLE :: Ndot1, Ndot2, Ndot3, Ndot4, Ndot5 REAL(r8) , DIMENSION(:) , ALLOCATABLE :: Edot1, Edot2, Edot3, Edot4, Edot5 REAL(r8) , DIMENSION(:,:), ALLOCATABLE :: n, nU, unU, nUE, P11, P22, E, B, dB, ddB REAL(r8) , DIMENSION(:) , ALLOCATABLE :: zm END TYPE outputTYP ! ----------------------------------------------------------------------------- TYPE meshTYP REAL(r8) , DIMENSION(:), ALLOCATABLE :: zm REAL(r8) :: LZ, zmin, zmax, dz INTEGER(i4) :: NZ REAL(r8) , DIMENSION(:), ALLOCATABLE :: n, nU, nUE, P11, P22 REAL(r8) , DIMENSION(:), ALLOCATABLE :: E, B, dB, ddB REAL(r8) , DIMENSION(:), ALLOCATABLE :: unU REAL(r8) :: LZ, zmin, zmax, dzm INTEGER(i4) :: NZmesh END TYPE meshTYP TYPE fieldsTYP REAL(r8) , DIMENSION(:), ALLOCATABLE :: E, B REAL(r8) , DIMENSION(:), ALLOCATABLE :: V, dB, ddB END TYPE fieldsTYP CONTAINS ! ---------------------------------------------------------------------------- SUBROUTINE AllocateFields(fields,params) IMPLICIT NONE ! Declare interface variables: TYPE(fieldsTYP), INTENT(INOUT) :: fields TYPE(paramsTYP), INTENT(IN) :: params ! Declare local variables: INTEGER(i4) :: NZ ! Allocate memory to mesh: NZ = params%NZmesh ALLOCATE(fields%E(NZ),fields%B(NZ),fields%dB(NZ),fields%ddB(NZ)) END SUBROUTINE AllocateFields ! ---------------------------------------------------------------------------- SUBROUTINE AllocateMesh(mesh,params) IMPLICIT NONE Loading @@ -321,11 +308,25 @@ SUBROUTINE AllocateMesh(mesh,params) TYPE(paramsTYP), INTENT(IN) :: params ! Declare local variables: INTEGER(i4) :: NZ INTEGER(i4) :: NZmesh ! Allocate memory to mesh: NZ = params%NZmesh ALLOCATE(mesh%zm(NZ)) NZmesh = params%NZmesh ALLOCATE(mesh%zm(NZmesh)) ! Allocate memory for mesh defined quantities: ! Ghost cells are used: ALLOCATE(mesh%n(NZmesh + 4)) ALLOCATE(mesh%nU(NZmesh + 4)) ALLOCATE(mesh%unU(NZmesh + 4)) ALLOCATE(mesh%P11(NZmesh + 4)) ALLOCATE(mesh%P22(NZmesh + 4)) ALLOCATE(mesh%nuE(NZmesh + 4)) ALLOCATE(mesh%B(NZmesh + 4)) ALLOCATE(mesh%dB(NZmesh + 4)) ALLOCATE(mesh%ddB(NZmesh + 4)) ALLOCATE(mesh%E(NZmesh + 4)) END SUBROUTINE AllocateMesh ! ---------------------------------------------------------------------------- Loading @@ -336,21 +337,33 @@ SUBROUTINE InitializeMesh(mesh,params) TYPE(paramsTYP), INTENT(IN) :: params ! Declare local variables: INTEGER(i4), DIMENSION(params%NZ) :: m INTEGER(i4), DIMENSION(params%NZmesh) :: m INTEGER(i4) :: i ! Populate fields: mesh%NZ = params%NZmesh mesh%NZmesh = params%NZmesh mesh%zmin = params%zmin mesh%zmax = params%zmax ! Derived quantities: mesh%LZ = params%zmax - params%zmin mesh%dz = mesh%LZ/mesh%NZ m = (/ (i, i=1,mesh%NZ, 1) /) mesh%dzm = mesh%LZ/mesh%NZmesh m = (/ (i, i=1,mesh%NZmesh, 1) /) ! Create mesh: mesh%zm = (m-1)*mesh%dz + 0.5*mesh%dz + mesh%zmin mesh%zm = (m-1)*mesh%dzm + 0.5*mesh%dzm + mesh%zmin ! Initialize all mesh-defined quantities: mesh%n = 0. mesh%nU = 0. mesh%unU = 0. mesh%nUE = 0. mesh%P11 = 0. mesh%P22 = 0. mesh%B = 0. mesh%E = 0. mesh%dB = 0. mesh%ddB = 0. END SUBROUTINE InitializeMesh Loading @@ -367,8 +380,9 @@ SUBROUTINE AllocatePlasma(plasma,params) ! NC: Number of computational particles NC = params%NC ! Allocate memory: For all computational particles ! Allocate memory: for all computational particles ALLOCATE(plasma%zp(NC) ,plasma%kep(NC) ,plasma%xip(NC) ,plasma%a(NC)) ALLOCATE(plasma%m(NC) ,plasma%wL(NC) ,plasma%wC(NC) ,plasma%wR(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), plasma%dE5(NC)) ALLOCATE(plasma%udErf(NC)) Loading Loading @@ -402,6 +416,7 @@ SUBROUTINE InitializePlasma(plasma,params) END IF ! Initialize plasma: scalar quantities: plasma%tp = 0. plasma%NR = params%ne0*params%Area0*(params%zmax - params%zmin) plasma%NSP = params%NC plasma%alpha = plasma%NR/plasma%NSP Loading Loading @@ -498,7 +513,7 @@ SUBROUTINE AllocateOutput(output,params) TYPE(paramsTYP), INTENT(IN) :: params ! Declare local variables: INTEGER(i4) :: j, jsize, NS INTEGER(i4) :: j, jsize, NS, NZmesh ! Determine size of temporal snapshots to record: jsize = (params%jend-params%jstart+1)/params%jincr Loading @@ -509,6 +524,7 @@ SUBROUTINE AllocateOutput(output,params) ALLOCATE(output%kep(params%NC,jsize)) ALLOCATE(output%xip(params%NC,jsize)) ALLOCATE(output%a(params%NC ,jsize)) ALLOCATE(output%m(params%NC ,jsize)) ALLOCATE(output%tp(jsize)) ! Create array with the indices of the time steps to save: Loading @@ -520,9 +536,23 @@ SUBROUTINE AllocateOutput(output,params) ALLOCATE(output%Ndot1(NS),output%Ndot2(NS),output%Ndot3(NS),output%Ndot4(NS) ,output%Ndot5(NS)) ALLOCATE(output%Edot1(NS),output%Edot2(NS),output%Edot3(NS),output%Edot4(NS) ,output%Edot5(NS)) ! Allocate memory for mesh-defined quantities: NZmesh = params%NZmesh ALLOCATE(output%zm(NZmesh)) ALLOCATE(output%n(NZmesh + 4 ,jsize)) ALLOCATE(output%nU(NZmesh + 4 ,jsize)) ALLOCATE(output%unU(NZmesh + 4,jsize)) ALLOCATE(output%nUE(NZmesh + 4,jsize)) ALLOCATE(output%P11(NZmesh + 4,jsize)) ALLOCATE(output%P22(NZmesh + 4,jsize)) ALLOCATE(output%B(NZmesh + 4 ,jsize)) ALLOCATE(output%E(NZmesh +4 ,jsize)) ALLOCATE(output%dB(NZmesh + 4 ,jsize)) ALLOCATE(output%ddB(NZmesh + 4,jsize)) END SUBROUTINE AllocateOutput ! --------------------------------------------------------------------------- ! ================================================================================= SUBROUTINE PrintParamsToTerminal(params,inputFile) IMPLICIT NONE ! Declare interface variables: Loading Loading @@ -558,4 +588,123 @@ SUBROUTINE PrintParamsToTerminal(params,inputFile) END SUBROUTINE PrintParamsToTerminal ! ================================================================================= SUBROUTINE SaveData(output,dir1) IMPLICIT NONE !Declare interface variables: TYPE(outputTYP), INTENT(IN) :: output CHARACTER*300 , INTENT(IN) :: dir1 ! Declare local variables: CHARACTER*300 :: fileName WRITE(*,*) "Saving data ..." ! Save plasma quantities: ! -------------------------------------------------------------------------- fileName = TRIM(TRIM(dir1)//'/'//'zp.out') OPEN(UNIT=8,FILE=fileName,FORM="unformatted",STATUS="unknown") WRITE(8) output%zp CLOSE(UNIT=8) fileName = trim(trim(dir1)//'/'//'kep.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%kep CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'xip.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%xip CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'a.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%a CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'tp.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%tp CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'m.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%m CLOSE(unit=8) ! Saving mesh-defined quantities: ! ------------------------------------------------------------------------- fileName = trim(trim(dir1)//'/'//'z_mesh.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%zm CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'n_mesh.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%n CLOSE(unit=8) ! Saving pcount to file: ! -------------------------------------------------------------------------- fileName = trim(trim(dir1)//'/'//'pcount1.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Ndot1 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'pcount2.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Ndot2 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'pcount3.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Ndot3 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'pcount4.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Ndot4 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'pcount5.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Ndot5 CLOSE(unit=8) ! Saving ecount to file: ! -------------------------------------------------------------------------- fileName = trim(trim(dir1)//'/'//'ecount1.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Edot1 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'ecount2.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Edot2 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'ecount3.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Edot3 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'ecount4.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Edot4 CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'ecount5.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%Edot5 CLOSE(unit=8) ! Save ER, NR and NSP: ! ------------------------------------------------------------------------- fileName = trim(trim(dir1)//'/'//'ER.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%ER CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'NR.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%NR CLOSE(unit=8) fileName = trim(trim(dir1)//'/'//'NSP.out') OPEN(unit=8,file=fileName,form="unformatted",status="unknown") WRITE(8) output%NSP CLOSE(unit=8) WRITE(*,*) "Data saving complete" END SUBROUTINE SaveData END MODULE dataTYP
src/MoveParticlePack.f90 +210 −0 Original line number Diff line number Diff line ! ======================================================================================================= SUBROUTINE AdvanceParticles(plasma,fieldspline,params) ! ======================================================================================================= USE local USE dataTYP USE OMP_LIB IMPLICIT NONE ! Declare interface variables: TYPE(plasmaTYP) , INTENT(INOUT) :: plasma TYPE(paramsTYP) , INTENT(IN) :: params TYPE(fieldSplineTYP), INTENT(IN) :: fieldspline ! Declare local variables: REAL(r8) :: dresNum, resNum0, resNum1 INTEGER(i4) :: i ! Initialize resonance number: dresNum = 0. resNum0 = 0. resNum1 = 0. !$OMP PARALLEL DO FIRSTPRIVATE(dresNum, resNum0, resNum1) SCHEDULE(STATIC) DO i = 1,params%NC ! 2.1 - Reset flags: CALL ResetFlags(i,plasma) ! 2.2- Compute resonance number: IF (params%iHeat) THEN CALL CyclotronResonanceNumber(i,plasma,fieldspline,params,resNum0) END IF ! 2.3- Advance zp, kep, xip: IF (params%iPush) THEN CALL MoveParticle(i,plasma,fieldspline,params) END IF ! 2.4- Check boundaries: set f1 and f2 flags IF (params%iPush) THEN CALL CheckBoundary(i,plasma,params) END IF ! Apply Coulomb collision operator: ! ------------------------------------------------------------------------ IF (params%iColl) THEN CALL collisionOperator(i,plasma,params) END IF ! 2.5- Compute resonance number: IF (params%iHeat) THEN CALL CyclotronResonanceNumber(i,plasma,fieldspline,params,resNum0) END IF ! 2.6- Check resonance: set f3 flag IF (params%iHeat) THEN 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 plasma%f3(i) = 1 CALL RFoperatorTerms(i,plasma,fieldspline,params) END IF END IF END DO !$OMP END PARALLEL DO RETURN END SUBROUTINE AdvanceParticles ! ======================================================================================================= SUBROUTINE CheckBoundary(i,plasma,params) ! ======================================================================================================= USE LOCAL USE dataTYP IMPLICIT NONE ! Declare interface variables: INTEGER(i4), INTENT(IN) :: i TYPE(plasmaTYP), INTENT(INOUT) :: plasma TYPE(paramsTYP), INTENT(IN) :: params IF (plasma%zp(i) .GE. params%zmax) THEN plasma%f2(i) = 1 plasma%dE2(i) = plasma%kep(i) ELSE IF (plasma%zp(i) .LE. params%zmin) THEN plasma%f1(i) = 1 plasma%dE1(i) = plasma%kep(i) END IF RETURN END SUBROUTINE CheckBoundary ! ======================================================================================================= SUBROUTINE MoveParticle(i,plasma,fieldspline,params) ! ======================================================================================================= Loading Loading @@ -133,6 +227,122 @@ M = 0 RETURN END SUBROUTINE RightHandSide ! ======================================================================================================= SUBROUTINE AssignCell(plasma,mesh,params) ! ======================================================================================================= USE LOCAL USE dataTYP USE OMP_LIB IMPLICIT NONE ! Define interface variables: TYPE(plasmaTYP), INTENT(INOUT) :: plasma TYPE(meshTYP) , INTENT(IN) :: mesh TYPE(paramsTYP), INTENT(IN) :: params ! Declare local variables: INTEGER(i4) :: i REAL(r8) :: Z, dz, zi, zoffset ! Mesh grid size: dz = mesh%dzm ! Mesh offset: zoffset = mesh%zmin !$OMP PARALLEL DO PRIVATE(zi,Z) DO i = 1,params%NC IF (plasma%f1(i) .EQ. 0 .AND. plasma%f2(i) .EQ. 0) THEN ! Using zm and zp, find nearest grid point (NGP): zi = plasma%zp(i) plasma%m(i) = NINT(0.5 + (zi - zoffset)/dz) ! Compute wL(i), wC(i), wR(i) Z = mesh%zm(plasma%m(i)) - zi plasma%wC(i) = 0.75 - (Z/dz)**2. plasma%wL(i) = 0.5*(1.5 + ((Z - dz)/dz) )**2. plasma%wR(i) = 0.5*(1.5 - ((Z + dz)/dz) )**2. END IF END DO !$OMP END PARALLEL DO RETURN END SUBROUTINE AssignCell ! ======================================================================================================= SUBROUTINE ExtrapolateMomentsToMesh(plasma,mesh,params) ! ======================================================================================================= USE LOCAL USE PhysicalConstants USE dataTYP USE OMP_LIB IMPLICIT NONE ! Define interface variables: TYPE(plasmaTYP), INTENT(IN) :: plasma TYPE(meshTYP) , INTENT(INOUT) :: mesh TYPE(paramsTYP), INTENT(IN) :: params ! Define local variables: INTEGER(i4) :: i, ix, frame REAL(r8) :: vpar, Ma, Ep, alpha, xip, vper, v, a ! Initialize mesh quantities: mesh%n = 0. mesh%nU = 0. mesh%unU = 0. mesh%P11 = 0. mesh%P22 = 0. mesh%nuE = 0. ! Species "a" mass: Ma = params%Ma ! Scaling factor: alpha = plasma%alpha ! Calculate moments and extrapolate to mesh points !$OMP PARALLEL DO private(ix,a,Ep,v,xip,vpar,vper) DO i = 1,params%NC IF (plasma%f1(i) .EQ. 0 .AND. plasma%f2(i) .EQ. 0) THEN ! Derived quantities: a = plasma%a(i) Ep = e_c*plasma%kep(i) v = sqrt(2*Ep/Ma) xip = plasma%xip(i) vpar = v*xip vper = v*sqrt(1 - xip**2.) ! Mesh point with ghost cells: ix = plasma%m(i) + 2 ! Plasma density: [NSP m^-3] mesh%n(ix-1) = mesh%n(ix-1) + plasma%wL(i)*a; mesh%n(ix) = mesh%n(ix) + plasma%wC(i)*a; mesh%n(ix+1) = mesh%n(ix+1) + plasma%wR(i)*a; END IF END DO !$OMP END PARALLEL DO ! Apply magnetic compression: mesh%n = mesh%n/params%Area0 ! Apply scaling factor: mesh%n = alpha*mesh%n/mesh%dzm ! Apply smoothing: frame = 9 ! CALL MovingMean(mesh%n ,frame) ! Calculate U, Tpar, Tper: END SUBROUTINE ExtrapolateMomentsToMesh ! ======================================================================================================= SUBROUTINE ReinjectParticles(i,plasma,params) ! ======================================================================================================= Loading
src/make_linFP.f +2 −2 Original line number Diff line number Diff line Loading @@ -6,7 +6,7 @@ OBJ_1 = Modules.o linFP.o PotentialProfile.o OBJ_2 = fitpack.o MoveParticlePack.o CoulombCollisions.o All: $(OBJ_1) $(OBJ_2) gfortran $(OPTFLAGS) -fopenmp $(OBJ_1) $(OBJ_2) -o linFP $(COMPILER) $(OPTFLAGS) -fopenmp $(OBJ_1) $(OBJ_2) -o linFP rm *.o *.mod Modules.o: Modules.f90 Loading @@ -22,7 +22,7 @@ fitpack.o: fitpack.f $(COMPILER) $(OPTFLAGS) -c -w fitpack.f MoveParticlePack.o: MoveParticlePack.f90 $(COMPILER) $(OPTFLAGS) -c MoveParticlePack.f90 $(COMPILER) $(OPTFLAGS) -c -fopenmp MoveParticlePack.f90 CoulombCollisions.o: CoulombCollisions.f90 $(COMPILER) $(OPTFLAGS) -fopenmp -c CoulombCollisions.f90 Loading