Commit a486eec7 authored by Caneses Marin, Juan Francisco's avatar Caneses Marin, Juan Francisco
Browse files

Using flags to accumulate all events and has been shown to be identical to the... 

Using flags to accumulate all events and has been shown to be identical to the initial method used with pcount and ecount
parent 4e191980

src/CoulombCollisions.f90

+2 −0
Original line number Diff line number Diff line
@@ -149,9 +149,11 @@ if (kep_pf_1 .le. 0.) kep_pf_1 = kep_pf_0 
if (kep_pf_1 .GT. params%elevel*Tb) then

	! Record slowing down energy during time step dt

	ecnt = ecnt + dE_pf

        plasma%E4(i) = plasma%E4(i) +  dE_pf

	! Count how many particles are involved in the slowing down power calculation

	if(species_a .EQ. species_b) then

		pcnt = pcnt + 1

                plasma%f4(i) = 1

	end if

end if

src/Modules.f90

+53 −26
Original line number Diff line number Diff line
@@ -258,8 +258,8 @@ 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 :: dE1, dE2, dE3, dE4, dE3_hat

 REAL(r8)   , DIMENSION(:), ALLOCATABLE :: dErf_hat, doppler

 REAL(r8)   , DIMENSION(:), ALLOCATABLE :: E1, E2, E3, E4, E3_hat

 REAL(r8)   , DIMENSION(:), ALLOCATABLE :: Erf_hat, doppler

 REAL(r8)   , DIMENSION(:), ALLOCATABLE :: NR , NSP

 REAL(r8)   , DIMENSION(:), ALLOCATABLE :: Eplus, Eminus

 REAL(r8)   , DIMENSION(:), ALLOCATABLE :: Ndot1, Ndot2, Ndot3, Ndot4
@@ -297,17 +297,17 @@ SUBROUTINE AllocatePlasma(plasma,params) 
   ! 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%E1(NC) ,plasma%E2(NC)  ,plasma%E3(NC)  ,plasma%E4(NC))

   ALLOCATE(plasma%Erf_hat(NC))

   ALLOCATE(plasma%doppler(NC))

   ALLOCATE(plasma%E3_hat(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))

   ALLOCATE(plasma%Edot3_hat(NS))



END SUBROUTINE AllocatePlasma
@@ -319,7 +319,7 @@ SUBROUTINE InitializePlasma(plasma,params) 
   ! Declare interface variables:

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

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

   INTEGER(i4) :: i

   INTEGER(i4) :: i, j



   ! Derived parameters: Reference cross sectional area

    params%Area0 = 0.5*params%dtheta*( params%r2**2. - params%r1**2.)
@@ -333,28 +333,55 @@ SUBROUTINE InitializePlasma(plasma,params) 
       params%Ma = params%Aion*m_p

   end if

    

   !$OMP PARALLEL

   !$OMP DO

   ! Initialize plasma: time dependent quantities:

   DO j = 1,params%NS

     ! Initialize plasma: NR and NC

     ! When source contrained fueling is used, these quantities

     ! can be time depedent

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

   plasma%NSP = params%NC



   ! Initialize plasma: particle weights

   plasma%a  = 1.

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

     plasma%NSP(j)        = params%NC

     plasma%Eplus(j)      = 0. 

     plasma%Eminus(j)     = 0.

     plasma%Ndot1(j)      = 0.

     plasma%Ndot2(j)      = 0.

     plasma%Ndot3(j)      = 0.

     plasma%Ndot4(j)      = 0.

     plasma%Edot1(j)      = 0.

     plasma%Edot2(j)      = 0.

     plasma%Edot3(j)      = 0.

     plasma%Edot4(j)      = 0.

     plasma%Edot3_hat(j)  = 0.

   END DO

   !$OMP DO



   !$OMP PARALLEL DO

   !$OMP DO

   ! Initialize plasma: For all computational particles

   DO i = 1,params%NC  

     plasma%a(i)  = 1.

     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.;

   !$OMP END DO

   !$OMP END PARALLEL 



END SUBROUTINE InitializePlasma



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

SUBROUTINE ResetFlags(i,plasma)

   USE LOCAL

   IMPLICIT NONE

   ! Declare interface variables:

   INTEGER(i4)    , INTENT(IN)    :: i

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

   

   plasma%f1(i) = 0. ; plasma%f2(i) = 0. ; plasma%f3(i) = 0. ; plasma%f4(i) = 0.

   plasma%E1(i) = 0. ; plasma%E2(i) = 0. ; plasma%E3(i) = 0. ; plasma%E4(i) = 0.

   plasma%Erf_hat(i) = 0.

   plasma%doppler(i) = 0.;



END SUBROUTINE ResetFlags



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

SUBROUTINE AllocateFieldSpline(fieldspline,params)

   USE spline_fits

src/MoveParticlePack.f90

+2 −0
Original line number Diff line number Diff line
@@ -394,8 +394,10 @@ plasma%xip(i) = upar1/u1 


! Record resonance event:

pcnt = pcnt + 1

plasma%f3(i) = 1

! Record energy kick:

ecnt = ecnt + dkep

plasma%E3(i) = dkep



RETURN

END SUBROUTINE RFHeatingOperator

src/linFP.f90

+58 −26
Original line number Diff line number Diff line
@@ -22,6 +22,9 @@ TYPE(paramsTYP) :: params 
TYPE(plasmaTYP)      :: plasma

TYPE(fieldSplineTYP) :: fieldspline

TYPE(outputTYP)      :: output

! Accumulators:

REAL(i4) :: p1,p2,p3,p4

REAL(i4) :: q1,q2,q3,q4

! DO loop indices:

INTEGER(i4) :: i,j,k

! Thread ID:
@@ -161,6 +164,9 @@ ostart = OMP_GET_WTIME() 
! Loop over time:

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

AllTime: do j = 1,params%NS

    ! Reset all accumulators:

    p1 = 0. ; p2 = 0. ; p3 = 0. ; p4 = 0.; 

    q1 = 0. ; q2 = 0. ; q3 = 0. ; q4 = 0.;

 

    !$OMP PARALLEL PRIVATE(pcnt1,pcnt2,pcnt3,pcnt4,ecnt1,ecnt2,ecnt3,ecnt4,dresNum,resNum0,resNum1)
@@ -174,7 +180,11 @@ AllTime: do j = 1,params%NS 
    ! Loop over particles:

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

    !$OMP DO SCHEDULE(STATIC)

    AllParticles: DO i = 1,params%NC

    NC_loop1: DO i = 1,params%NC



        ! Initialize plasma: flags and Energy increments

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

        CALL ResetFlags(i,plasma)



        ! Calculate Cyclotron resonance number:

        ! ------------------------------------------------------------------------
@@ -188,11 +198,11 @@ AllTime: do j = 1,params%NS 
        ! ------------------------------------------------------------------------

        IF (plasma%zp(i) .GE. params%zmax) THEN

           plasma%f2(i) = 1

           plasma%dE2(i) = plasma%kep(i)

           plasma%E2(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)

           plasma%E1(i) = plasma%kep(i)

           CALL ReinjectParticles(i,plasma,params,ecnt1,pcnt1)

        END IF
@@ -211,22 +221,44 @@ AllTime: do j = 1,params%NS 
           END IF

        END IF



    END DO AllParticles

    END DO NC_loop1

    !$OMP END DO



    !$OMP CRITICAL

     ecount1(j) = ecount1(j) + ecnt1

     ecount2(j) = ecount2(j) + ecnt2

     ecount3(j) = ecount3(j) + ecnt3

     ecount4(j) = ecount4(j) + ecnt4

     pcount1(j) = pcount1(j) + pcnt1

     pcount2(j) = pcount2(j) + pcnt2

     pcount3(j) = pcount3(j) + pcnt3

     pcount4(j) = pcount4(j) + pcnt4

    !$OMP END CRITICAL

    !!$OMP CRITICAL

     !ecount1(j) = ecount1(j) + ecnt1

     !ecount2(j) = ecount2(j) + ecnt2

     !ecount3(j) = ecount3(j) + ecnt3

     !ecount4(j) = ecount4(j) + ecnt4

     !pcount1(j) = pcount1(j) + pcnt1

     !pcount2(j) = pcount2(j) + pcnt2

     !pcount3(j) = pcount3(j) + pcnt3

     !pcount4(j) = pcount4(j) + pcnt4

    !!$OMP END CRITICAL



    !$OMP DO REDUCTION(+:p1,p2,p3,p4,q1,q2,q3,q4)

    NC_loop2: DO i = 1,params%NC

       p1 = p1 + plasma%f1(i)

       p2 = p2 + plasma%f2(i)

       p3 = p3 + plasma%f3(i)

       p4 = p4 + plasma%f4(i)

       q1 = q1 + plasma%f1(i)*plasma%E1(i)

       q2 = q2 + plasma%f2(i)*plasma%E2(i)

       q3 = q3 + plasma%f3(i)*plasma%E3(i)

       q4 = q4 + plasma%f4(i)*plasma%E4(i)

    END DO NC_loop2

    !$OMP END DO



    !$OMP END PARALLEL



    plasma%Ndot1(j) = p1

    plasma%Ndot2(j) = p2

    plasma%Ndot3(j) = p3

    plasma%Ndot4(j) = p4

    plasma%Edot1(j) = q1

    plasma%Edot2(j) = q2

    plasma%Edot3(j) = q3

    plasma%Edot4(j) = q4

    

    ! Select data to save:

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

    ! Check if data is to be saved
@@ -328,38 +360,38 @@ if (params%iSave) then 
    ! --------------------------------------------------------------------------

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

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

    WRITE(8) pcount1

    WRITE(8) plasma%Ndot1

    CLOSE(unit=8)

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

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

    WRITE(8) pcount2

    WRITE(8) plasma%Ndot2

    CLOSE(unit=8)

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

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

    WRITE(8) pcount3

    WRITE(8) plasma%Ndot3

    CLOSE(unit=8)

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

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

    WRITE(8) pcount4

    WRITE(8) plasma%Ndot4

    CLOSE(unit=8)



    ! Saving ecount to file:

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

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

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

    WRITE(8) ecount1

    WRITE(8) plasma%Edot1

    CLOSE(unit=8)

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

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

    WRITE(8) ecount2

    WRITE(8) plasma%Edot2

    CLOSE(unit=8)

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

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

    WRITE(8) ecount3

    WRITE(8) plasma%Edot3

    CLOSE(unit=8)

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

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

    WRITE(8) ecount4

    WRITE(8) plasma%Edot4

    CLOSE(unit=8)



    ! Write output metadata: