Commit 548748f5 authored by Juan Caneses Marin (nfc)'s avatar Juan Caneses Marin (nfc)
Browse files

Demonstrated ExtrapolateMomentsToMesh, we are now going to interpolate density... 

Demonstrated ExtrapolateMomentsToMesh, we are now going to interpolate density and flow to particle positions prior to the collision operator
parent 47a87d46

InputFiles/xp_TestDrag.in

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

! Simulation name:

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

params%fileDescriptor = '2021_03_10a',

params%fileDescriptor = '2021_03_11a',



! Magnetic field input data:

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

src/Modules.f90

+60 −11
Original line number Diff line number Diff line
@@ -284,6 +284,7 @@ TYPE outputTYP 
 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 :: U, Ppar, Pper, Tpar, Tper

 REAL(r8)   , DIMENSION(:)  , ALLOCATABLE :: zm



END TYPE outputTYP
@@ -292,6 +293,7 @@ END TYPE outputTYP 
TYPE meshTYP

 REAL(r8) , DIMENSION(:), ALLOCATABLE :: zm

 REAL(r8) , DIMENSION(:), ALLOCATABLE :: n, nU, nUE, P11, P22

 REAL(r8) , DIMENSION(:), ALLOCATABLE :: U, Ppar, Pper, Tpar, Tper

 REAL(r8) , DIMENSION(:), ALLOCATABLE :: E, B, dB, ddB

 REAL(r8) , DIMENSION(:), ALLOCATABLE :: unU

 REAL(r8) :: LZ, zmin, zmax, dzm
@@ -321,11 +323,16 @@ SUBROUTINE AllocateMesh(mesh,params) 
   ALLOCATE(mesh%unU(NZmesh + 4))

   ALLOCATE(mesh%P11(NZmesh + 4))

   ALLOCATE(mesh%P22(NZmesh + 4))

   ALLOCATE(mesh%nuE(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))

   ALLOCATE(mesh%U(NZmesh + 4))

   ALLOCATE(mesh%Ppar(NZmesh + 4))

   ALLOCATE(mesh%Pper(NZmesh + 4))

   ALLOCATE(mesh%Tpar(NZmesh + 4))

   ALLOCATE(mesh%Tper(NZmesh + 4))



END SUBROUTINE AllocateMesh
@@ -364,6 +371,11 @@ SUBROUTINE InitializeMesh(mesh,params) 
   mesh%E    = 0.

   mesh%dB   = 0.

   mesh%ddB  = 0.

   mesh%U    = 0.

   mesh%Ppar = 0.

   mesh%Pper = 0.

   mesh%Tpar = 0.

   mesh%Tper = 0.



END SUBROUTINE InitializeMesh
@@ -549,6 +561,11 @@ SUBROUTINE AllocateOutput(output,params) 
  ALLOCATE(output%E(NZmesh +4   ,jsize))

  ALLOCATE(output%dB(NZmesh + 4 ,jsize))

  ALLOCATE(output%ddB(NZmesh + 4,jsize))

  ALLOCATE(output%U(NZmesh + 4,jsize))

  ALLOCATE(output%Ppar(NZmesh + 4,jsize))

  ALLOCATE(output%Pper(NZmesh + 4,jsize))

  ALLOCATE(output%Tpar(NZmesh + 4,jsize))

  ALLOCATE(output%Tper(NZmesh + 4,jsize))



END SUBROUTINE AllocateOutput
@@ -638,6 +655,38 @@ SUBROUTINE SaveData(output,dir1) 
    OPEN(unit=8,file=fileName,form="unformatted",status="unknown")

    WRITE(8) output%n

    CLOSE(unit=8)

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

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

    WRITE(8) output%nU

    CLOSE(unit=8)

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

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

    WRITE(8) output%unU

    CLOSE(unit=8)

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

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

    WRITE(8) output%nUE

    CLOSE(unit=8)

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

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

    WRITE(8) output%Ppar

    CLOSE(unit=8)

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

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

    WRITE(8) output%Pper

    CLOSE(unit=8)

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

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

    WRITE(8) output%Tpar

    CLOSE(unit=8)

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

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

    WRITE(8) output%Tper

    CLOSE(unit=8)    

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

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

    WRITE(8) output%U

    CLOSE(unit=8)





    ! Saving pcount to file:

src/MoveParticlePack.f90

+33 −15
Original line number Diff line number Diff line
@@ -288,33 +288,37 @@ TYPE(paramsTYP), INTENT(IN) :: params 


! Define local variables:

INTEGER(i4) :: i, ix, frame, NZ

INTEGER(i4), DIMENSION(3) :: ixLeft, ixRight

REAL(r8) :: vpar, Ma, Ep, alpha, xip, vper, v, a

REAL(r8), DIMENSION(mesh%NZmesh + 4) :: n, nU, unU, P11, P22, nUE



! Initialize mesh quantities:

! Initialize local mesh quantities:

n = 0.

nU = 0.

unU = 0.

mesh%P11 = 0.

mesh%P22 = 0.

mesh%nUE = 0.

P11 = 0.

P22 = 0.

nUE = 0.



! Species "a" mass:

Ma = params%Ma

! Scaling factor:

alpha = plasma%alpha

! Size of mesh quantities:

NZ = mesh%NZmesh + 4

NZ = mesh%NZmesh

! Set the edge cell number:

ixLeft  = (/1,2,3/) 

ixRight = (/2,3,4/) + NZ



! Calculate moments and extrapolate to mesh points

!$OMP PARALLEL DO PRIVATE(ix,a,Ep,v,xip,vpar,vper) REDUCTION(+:n, nU, unU, P11, P22, nuE)

!$OMP PARALLEL DO PRIVATE(a,Ep,v,xip,vpar,vper,ix) REDUCTION(+:n, nU, unU, P11, P22, nUE)

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)

		v    = sqrt(2.*Ep/Ma)

		xip  = plasma%xip(i)

		vpar = v*xip

		vper = v*sqrt(1 - xip**2.)
@@ -356,9 +360,14 @@ DO i = 1,params%NC 
END DO

!$OMP END PARALLEL DO



! Correct edge cells of particle density:

n(ixLeft)  = n(4)

n(ixRight) = n(NZ+1)



! Apply magnetic compression:

mesh%n   =  n/params%Area0

mesh%nU  = nU/params%Area0

mesh%unU = unU

mesh%P11 = P11/params%Area0

mesh%P22 = P22/params%Area0

mesh%nUE = nUE/params%Area0
@@ -374,15 +383,24 @@ mesh%nUE = alpha*mesh%nUE/mesh%dzm 
! Apply smoothing:

frame = 9

!$OMP PARALLEL

IF (OMP_GET_THREAD_NUM() .EQ. 0) CALL MovingMean(mesh%n  ,NZ,frame)

IF (OMP_GET_THREAD_NUM() .EQ. 1) CALL MovingMean(mesh%nU ,NZ,frame)

IF (OMP_GET_THREAD_NUM() .EQ. 2) CALL MovingMean(mesh%unU,NZ,frame)

IF (OMP_GET_THREAD_NUM() .EQ. 3) CALL MovingMean(mesh%P11,NZ,frame)

IF (OMP_GET_THREAD_NUM() .EQ. 4) CALL MovingMean(mesh%P22,NZ,frame)

IF (OMP_GET_THREAD_NUM() .EQ. 5) CALL MovingMean(mesh%nUE,NZ,frame)

IF (OMP_GET_THREAD_NUM() .EQ. 0) CALL MovingMean(mesh%n(3:(NZ+2))  ,NZ,frame)

IF (OMP_GET_THREAD_NUM() .EQ. 1) CALL MovingMean(mesh%nU(3:(NZ+2)) ,NZ,frame)

IF (OMP_GET_THREAD_NUM() .EQ. 2) CALL MovingMean(mesh%unU(3:(NZ+2)),NZ,frame)

IF (OMP_GET_THREAD_NUM() .EQ. 3) CALL MovingMean(mesh%P11(3:(NZ+2)),NZ,frame)

IF (OMP_GET_THREAD_NUM() .EQ. 4) CALL MovingMean(mesh%P22(3:(NZ+2)),NZ,frame)

IF (OMP_GET_THREAD_NUM() .EQ. 5) CALL MovingMean(mesh%nUE(3:(NZ+2)),NZ,frame)

!$OMP END PARALLEL



! Calculate U, Tpar, Tper:

! Calculate U:

mesh%U    = mesh%nU/mesh%n



! Calculate Ppar and Pper:

mesh%Ppar = mesh%P11 - Ma*mesh%nU*mesh%U

mesh%Pper = mesh%P22



! Calculate Tpar and Tper:

mesh%Tpar = mesh%Ppar/(e_c*mesh%n)

mesh%Tper = mesh%Pper/(e_c*mesh%n)



END SUBROUTINE ExtrapolateMomentsToMesh

src/linFP.f90

+11 −3
Original line number Diff line number Diff line
@@ -226,9 +226,17 @@ NS_loop: DO j = 1,params%NS 
          END DO

         !$OMP END PARALLEL DO

         ! Record mesh defined quantities:

	 DO i = 1,(mesh%NZmesh + 4)

	     output%n(i,k) = mesh%n(i)

	 END DO

	 !DO i = 1,(mesh%NZmesh + 4)

	 output%n(:,k)    = mesh%n

	 output%nU(:,k)   = mesh%nU

	 output%unU(:,k)  = mesh%unU

	 output%nUE(:,k)  = mesh%nUE

	 output%Ppar(:,k) = mesh%Ppar

	 output%Pper(:,k) = mesh%Pper

	 output%Tpar(:,k) = mesh%Tpar

	 output%Tper(:,k) = mesh%Tper

	 output%U(:,k)    = mesh%U

	 !END DO

         ! Increment counter:

         k = k + 1

      END IF