Commit 2e134be9 authored by Juan Caneses Marin (nfc)'s avatar Juan Caneses Marin (nfc)
Browse files

Extracted the collision operator from AdvanceParticles and sucessfully made it... 

Extracted the collision operator from AdvanceParticles and sucessfully made it a separate step. No mesh interpolated quantities used yet
parent e201619c

InputFiles/xp_TestDrag.in

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

! Simulation name:

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

params%fileDescriptor = '2021_03_11a',

params%fileDescriptor = '2021_03_11c',



! Magnetic field input data:

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

src/CoulombCollisions.f90

+48 −4
Original line number Diff line number Diff line
@@ -15,10 +15,55 @@ TYPE(paramsTYP), INTENT(IN) :: params 


! Define local variables:

INTEGER(i4) :: i, ix

REAL(r8), DIMENSION (3) :: w, n, U, Tpar, Tper



! Interpolate n ,U ,Tpar and Tper to particle positions:

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

!$OMP PARALLEL DO PRIVATE(ix, w, n, U, Tpar, Tper)

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)

		

		! Plasma density:

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

		n(2) = mesh%n(ix)

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

		plasma%np(i) = w(1)*n(1) + w(2)*n(2) + w(3)*n(3)  



		! Parallel drift velocity:

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

		U(2) = mesh%U(ix)

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

		plasma%Up(i) = w(1)*U(1) + w(2)*U(2) + w(3)*U(3)  



		! Parallel Temperature:

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

		Tpar(2) = mesh%Tpar(ix)

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

		plasma%Tparp(i) = w(1)*Tpar(1) + w(2)*Tpar(2) + w(3)*Tpar(3) 



 		! Perpendicular Temperature:

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

		Tper(2) = mesh%Tper(ix)

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

		plasma%Tperp(i) = w(1)*Tper(1) + w(2)*Tper(2) + w(3)*Tper(3) 

	END IF

END DO

!$OMP END PARALLEL DO



! Interpolate n and U to particle positions:





! Apply collision operator:

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

!$OMP PARALLEL DO

 DO i = 1,params%NC

	 CALL CollisionOperator(i,plasma,params)

 END DO

!$OMP END PARALLEL DO



END SUBROUTINE ApplyCollisionOperator
@@ -36,7 +81,6 @@ IMPLICIT NONE 
INTEGER(i4)    , INTENT(IN)    :: i

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

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

!REAL(r8)       , INTENT(INOUT) :: ecnt, pcnt



! Define local variables:

REAL(r8) :: zp0, xip0, kep0

src/Modules.f90

+51 −24
Original line number Diff line number Diff line
@@ -258,12 +258,15 @@ END TYPE paramsTYP 


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

TYPE plasmaTYP

 ! Particle quantities:

 REAL(r8)   , DIMENSION(:), ALLOCATABLE :: zp, kep,xip, a

 REAL(r8)   , DIMENSION(:), ALLOCATABLE :: wL, wC, wR

 REAL(r8)   , DIMENSION(:), ALLOCATABLE :: n, U, Tpar, Tper, B, E, dB, ddB

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

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

 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

 ! Global quantities (scalars):

 REAL(r8) :: NR , NSP, alpha, tp

 REAL(r8) :: Eplus, Eminus

 REAL(r8) :: Ndot1, Ndot2, Ndot3, Ndot4
@@ -277,13 +280,17 @@ END TYPE fieldSplineTYP 


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

TYPE outputTYP

 ! Particle quantities:

 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 :: tp, jrng

! Global quantities:

 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

! Mesh quantities:

 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
@@ -363,16 +370,16 @@ SUBROUTINE InitializeMesh(mesh,params,fieldspline) 
   mesh%zm = (m-1)*mesh%dzm + 0.5*mesh%dzm + mesh%zmin



   ! Initialize all mesh-defined quantities:

   mesh%n    = 0. ! interpolation?

   mesh%n    = 0.

   mesh%nU   = 0.

   mesh%unU  = 0.

   mesh%nUE  = 0.

   mesh%P11  = 0.

   mesh%P22  = 0.

   mesh%B    = 0. ! Interpolation?

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

   mesh%E    = 0.

   mesh%dB   = 0. ! Interpolation?

   mesh%ddB  = 0. ! Interpolation?

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

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

   mesh%U    = 0.

   mesh%Ppar = 0.

   mesh%Pper = 0.
@@ -394,11 +401,11 @@ SUBROUTINE AllocatePlasma(plasma,params) 
   ! NC: Number of computational particles

   NC = params%NC

 

   ! Allocate memory: for all computational particles

   ! Allocate memory: Particle quantities: 

   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%n(NC)   ,plasma%U(NC)   ,plasma%Tpar(NC),plasma%Tper(NC))

   ALLOCATE(plasma%B(NC)   ,plasma%E(NC)   ,plasma%dB(NC)  ,plasma%ddB(NC))

   ALLOCATE(plasma%np(NC)  ,plasma%Up(NC)  ,plasma%Tparp(NC),plasma%Tperp(NC))

   ALLOCATE(plasma%Bp(NC)  ,plasma%Ep(NC)  ,plasma%dBp(NC)  ,plasma%ddBp(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))
@@ -458,14 +465,14 @@ SUBROUTINE InitializePlasma(plasma,params) 


   !$OMP PARALLEL DO

   DO i = 1,params%NC

     plasma%n(i)    = params%ne0

     plasma%U(i)    = 0.

     plasma%Tpar(i) = params%Ti0

     plasma%Tper(i) = params%Ti0

     plasma%B(i)    = 0.

     plasma%E(i)    = 0.

     plasma%dB(i)   = 0.

     plasma%ddB(i)  = 0.

     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
@@ -547,25 +554,29 @@ SUBROUTINE AllocateOutput(output,params) 
   ! Determine size of temporal snapshots to record:

   jsize = (params%jend-params%jstart+1)/params%jincr

   

   ! Allocate memory:

   ! Allocate memory: Particle quantities

   ALLOCATE(output%jrng(jsize))

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

   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%np(params%NC  ,jsize))

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

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

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

   ALLOCATE(output%tp(jsize))



   ! Create array with the indices of the time steps to save:

   output%jrng = (/ (j, j=params%jstart, params%jend, params%jincr) /)



  ! Allocate memory: For all time steps

  ! Allocate memory: Global quantities

   NS = params%NS

   ALLOCATE(output%NR(NS)   ,output%NSP(NS)  ,output%Eplus(NS),output%Eminus(NS),output%ER(NS))

   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:

  ! Allocate memory: mesh-defined quantities:

  NZmesh = params%NZmesh

  ALLOCATE(output%zm(NZmesh))

  ALLOCATE(output%n(NZmesh + 4  ,jsize))
@@ -635,7 +646,7 @@ SUBROUTINE SaveData(output,dir1) 
 

 WRITE(*,*) "Saving data ..."



    ! Save plasma quantities:

    ! Save particle quantities:

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

    fileName = TRIM(TRIM(dir1)//'/'//'zp.out')

    OPEN(UNIT=8,FILE=fileName,FORM="unformatted",STATUS="unknown")
@@ -661,6 +672,22 @@ SUBROUTINE SaveData(output,dir1) 
    OPEN(unit=8,file=fileName,form="unformatted",status="unknown")

    WRITE(8) output%m

    CLOSE(unit=8)

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

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

    WRITE(8) output%np

    CLOSE(unit=8)

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

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

    WRITE(8) output%Up

    CLOSE(unit=8)

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

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

    WRITE(8) output%Tparp

    CLOSE(unit=8)

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

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

    WRITE(8) output%Tperp

    CLOSE(unit=8)



    ! Saving mesh-defined quantities:

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

src/MoveParticlePack.f90

+0 −6
Original line number Diff line number Diff line
@@ -42,12 +42,6 @@ DO i = 1,params%NC 
		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)

src/linFP.f90

+13 −9
Original line number Diff line number Diff line
@@ -217,19 +217,23 @@ NS_loop: DO j = 1,params%NS 
    ! Check if data is to be saved

    IF (params%iSave) THEN

       IF (j .EQ. output%jrng(k)) THEN

	  ! Record particle quantities:

          output%tp(k) = plasma%tp

          !$OMP PARALLEL DO

          DO i = 1,params%NC

             ! Record "ith" particle at "kth" time

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

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

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

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

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

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

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

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

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

          END DO

         !$OMP END PARALLEL DO

         ! Record mesh defined quantities:

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



	 ! Record mesh quantities:

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

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

	 output%unU(:,k)  = mesh%unU
@@ -239,7 +243,7 @@ NS_loop: DO j = 1,params%NS 
	 output%Tpar(:,k) = mesh%Tpar

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

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

	 !END DO



         ! Increment counter:

         k = k + 1

      END IF