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

Reorganized the main for loop by having a single OMP parallel construct for each time step

parent ddc7b424

src/Main.f90

+169 −190
Original line number Diff line number Diff line
@@ -36,7 +36,7 @@ INTEGER(i4), DIMENSION(:), ALLOCATABLE :: jrng 
INTEGER(i4) :: id

! CPU time at start and end of computation:

REAL(r8) :: tstart, tend

! openMP computational time:

! OPENMP computational time:

DOUBLE PRECISION :: ostart, oend, oend_estimate

! Cyclotron resonance number change:

REAL(r8) :: df
@@ -56,9 +56,9 @@ REAL(r8), DIMENSION(:) , ALLOCATABLE :: pcount1, pcount2, pcount3, pcount4 
REAL(r8), DIMENSION(:)  , ALLOCATABLE :: ecount1, ecount2, ecount3, ecount4

! Cyclotron resonance number:

REAL(r8), DIMENSION(:)  , ALLOCATABLE :: fcurr, fnew

! Local simulation diagnostics:

REAL(r8) :: ecnt, ecnt1, ecnt2

REAL(r8) :: pcnt, pcnt1, pcnt2

! Local simulation diagnostic counters:

REAL(r8) :: ecnt1, ecnt2, ecnt3, ecnt4

REAL(r8) :: pcnt1, pcnt2, pcnt3, pcnt4

! To store system commands and fileNames:

CHARACTER*300 :: command, mpwd

INTEGER(i4) :: n_mpwd, STATUS
@@ -67,7 +67,6 @@ CHARACTER*300 :: inputFileDir, inputFile, fileName, xpSelector, rootDir, dir0, d 
! Create input namelist from the user-defined structures:

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

namelist/in_nml/in

!namelist/xp_nml/xpSelector



! Get root directory:

! ==============================================================================
@@ -80,9 +79,9 @@ CALL GET_ENVIRONMENT_VARIABLE('INPUT_FILE_DIR',inputFileDir) 


! Read input data into "in" structure:

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

open(unit=4,file=inputFileDir,status='old',form='formatted')

OPEN(unit=4,file=inputFileDir,status='old',form='formatted')

read(4,in_nml)

close(unit=4)

CLOSE(unit=4)



! Populate the in%rootDir field:

! =============================================================================
@@ -100,36 +99,36 @@ end if 
! ==============================================================================

WRITE(*,*) ''

WRITE(*,*) '*********************************************************************'

print *, 'Input file:         ', TRIM(inputFile)

print *, 'fileDescriptor:     ', TRIM(in%fileDescriptor)

print *, 'Number of particles:', in%Nparts

print *, 'Number of steps:    ', in%Nsteps

print *, 'Particle BC:        ', in%particleBC

print *, 'dt [ns]:            ', in%dt*1E+9

print *, 'iPush:              ', in%iPush

print *, 'iDrag:              ', in%iDrag

print *, 'iColl:              ', in%iColl

print *, 'iHeat:              ', in%iHeat

print *, 'iSave:              ', in%iSave

print *, 'elevel:             ', in%elevel

print *, 'zTarget [m]:        ', in%zmax

print *, 'zDump [m]:          ', in%zmin

print *, 'zp_init [m]:        ', in%zp_init

print *, 'B field file:       ', TRIM(in%BFieldFile)

print *, 'Ew:                 ', in%Ew

print *, 'Te0:                ', in%Te0

print *, 'ne0:                ', in%ne0

if (in%CollOperType .EQ. 1) print *, 'Boozer-Only collision operator'

if (in%CollOperType .EQ. 2) print *, 'Boozer-Kim collision operator'

WRITE(*,*) 'Input file:         ', TRIM(inputFile)

WRITE(*,*) 'fileDescriptor:     ', TRIM(in%fileDescriptor)

WRITE(*,*) 'Number of particles:', in%Nparts

WRITE(*,*) 'Number of steps:    ', in%Nsteps

WRITE(*,*) 'Particle BC:        ', in%particleBC

WRITE(*,*) 'dt [ns]:            ', in%dt*1E+9

WRITE(*,*) 'iPush:              ', in%iPush

WRITE(*,*) 'iDrag:              ', in%iDrag

WRITE(*,*) 'iColl:              ', in%iColl

WRITE(*,*) 'iHeat:              ', in%iHeat

WRITE(*,*) 'iSave:              ', in%iSave

WRITE(*,*) 'elevel:             ', in%elevel

WRITE(*,*) 'zTarget [m]:        ', in%zmax

WRITE(*,*) 'zDump [m]:          ', in%zmin

WRITE(*,*) 'zp_init [m]:        ', in%zp_init

WRITE(*,*) 'B field file:       ', TRIM(in%BFieldFile)

WRITE(*,*) 'Ew:                 ', in%Ew

WRITE(*,*) 'Te0:                ', in%Te0

WRITE(*,*) 'ne0:                ', in%ne0

if (in%CollOperType .EQ. 1) WRITE(*,*) 'Boozer-Only collision operator'

if (in%CollOperType .EQ. 2) WRITE(*,*) 'Boozer-Kim collision operator'

WRITE(*,*) '*********************************************************************'

WRITE(*,*) ''



! Allocate memory to splines:

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

call InitSpline(spline_Bz  ,in%nz,0._8,0._8,1,0._8)

call InitSpline(spline_ddBz,in%nz,0._8,0._8,1,0._8)

call InitSpline(spline_Phi ,in%nz,0._8,0._8,1,0._8)

call InitSplineTest(spline_Test,in%nz)

CALL InitSpline(spline_Bz  ,in%nz,0._8,0._8,1,0._8)

CALL InitSpline(spline_ddBz,in%nz,0._8,0._8,1,0._8)

CALL InitSpline(spline_Phi ,in%nz,0._8,0._8,1,0._8)

CALL InitSplineTest(spline_Test,in%nz)



! Allocate memory to main simulation variables:

! ==============================================================================
@@ -168,11 +167,12 @@ CALL ComputeSpline(spline_ddBz) 
spline_Phi%x = spline_Bz%x

spline_Phi%y = 0

if (in%iPotential) then

  call PotentialProfile(spline_Phi,in)

  CALL PotentialProfile(spline_Phi,in)

end if

CALL ComputeSpline(spline_Phi)



! Test the splines:

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

if (.true.) then

    do i=1,in%nz

        !spline_Test%x(i)  = in%zmin + i*0.03
@@ -185,12 +185,12 @@ if (.true.) then 
        spline_Test%y6(i) = BESSEL_JN(1,spline_Test%x(i))

    end do

    fileName = "B_spline.dat"

    open(unit=8,file=fileName,form="formatted",status="unknown")

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

    do j = 1,in%nz

        write(8,*) spline_Test%x(j), spline_Test%y1(j), spline_Test%y2(j),&

        WRITE(8,*) spline_Test%x(j), spline_Test%y1(j), spline_Test%y2(j),&

         spline_Test%y3(j), spline_Test%y4(j), spline_Test%y5(j), spline_Test%y6(j)

    end do

    close(unit=8)

    CLOSE(unit=8)

end if



! Record start time:
@@ -199,24 +199,25 @@ ostart = OMP_GET_WTIME() 


! Initialize pseudo random number generator:

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

call random_seed(size=seed_size)

CALL random_seed(size=seed_size)

ALLOCATE(seed(seed_size))

call random_seed(get=seed)

CALL random_seed(get=seed)

seed = 314159565

call random_seed(put=seed)

CALL random_seed(put=seed)



! Inititalize zp, kep, xip

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

kep = 0.; xip = 0.; zp = 0.;

call loadParticles(zp,kep,xip,in)

CALL loadParticles(zp,kep,xip,in)



! Test initial distribution:

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

fileName = "LoadParticles.dat"

open(unit=8,file=fileName,form="formatted",status="unknown")

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

do i = 1,in%Nparts

    write(8,*) zp(i), kep(i), xip(i)

    WRITE(8,*) zp(i), kep(i), xip(i)

end do

close(unit=8)

CLOSE(unit=8)



! Initialize simulation diagnostics:

! ==============================================================================
@@ -230,115 +231,93 @@ ecount1 = 0; ecount2 = 0; ecount3 = 0; ecount4 = 0 
! OMP setup:

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

! Set the number of threads:

call OMP_SET_NUM_THREADS(in%threads_request)

!$OMP PARALLEL PRIVATE(id)

id = OMP_GET_THREAD_NUM()

CALL OMP_SET_NUM_THREADS(in%threads_request)

in%threads_given = OMP_GET_NUM_THREADS()

if (id .EQ. 0) then

WRITE(*,*) ''

WRITE(*,*) '*********************************************************************'

	write(*,*) "number of threads given: ", in%threads_given

WRITE(*,*) "number of threads given: ", in%threads_given

WRITE(*,*) '*********************************************************************'

WRITE(*,*) ''

end if



!$OMP PARALLEL PRIVATE(id)

    !id = OMP_GET_THREAD_NUM()

    !if (id .EQ. 0) then

    !	WRITE(*,*) ''

    !	WRITE(*,*) '*********************************************************************'

    !	WRITE(*,*) "number of threads given: ", in%threads_given

    !	WRITE(*,*) '*********************************************************************'

    !	WRITE(*,*) ''

  !  end if

!$OMP END PARALLEL



! Start of simulation:

! Loop over time:

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

! Begin time stepping:

TimeStepping: do j = 1,in%Nsteps



    !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(pcnt1,pcnt2,pcnt3,pcnt4,ecnt1,ecnt2,ecnt3,ecnt4,i,df)



          ! Initialize private particle counters:

          pcnt1 = 0; pcnt2 = 0; pcnt3 = 0; pcnt4 = 0

          ! Initialize private energy counters:

          ecnt1 = 0; ecnt2 = 0; ecnt3 = 0; ecnt4 = 0

          ! Initialize resonance number difference:

          df = 0



          ! Loop over particles:

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

          !$OMP DO

              AllParticles: do i = 1,in%Nparts



                ! Calculate Cyclotron resonance number:

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

    if (in%iHeat) then

        do i = 1,in%Nparts

            call CyclotronResonanceNumber(zp(i),kep(i),xip(i),fcurr(i),in,spline_Bz)

        end do

    end if

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

                if (in%iHeat) CALL CyclotronResonanceNumber(zp(i),kep(i),xip(i),fcurr(i),in,spline_Bz)

                ! fcurr and fnew could be declared private



                ! Push particles adiabatically:

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

    if (in%iPush) then

      !$OMP PARALLEL DO PRIVATE(i)

        do i = 1,in%Nparts

            call MoveParticle(zp(i),kep(i),xip(i),in,spline_Bz,spline_Phi)

        end do

      !$OMP END PARALLEL DO

    end if

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

                if (in%iPush) CALL MoveParticle(zp(i),kep(i),xip(i),in,spline_Bz,spline_Phi)



                ! Re-inject particles:

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

    if (.true.) then

      !$OMP PARALLEL PRIVATE(ecnt1, ecnt2, pcnt1, pcnt2)

          ecnt1 = 0; ecnt2 = 0; pcnt1 = 0; pcnt2 = 0

          !$OMP DO

              do i = 1,in%Nparts

                  if (zp(i) .GE. in%zmax) then

                    call ReinjectParticles(zp(i),kep(i),xip(i),in,ecnt2,pcnt2)

                  else if (zp(i) .LE. in%zmin) then

                    call ReinjectParticles(zp(i),kep(i),xip(i),in,ecnt1,pcnt1)

                  end if

              end do

          !$OMP END DO

          !$OMP CRITICAL

              ecount1(j) = ecount1(j) + ecnt1

              pcount1(j) = pcount1(j) + pcnt1

              ecount2(j) = ecount2(j) + ecnt2

              pcount2(j) = pcount2(j) + pcnt2

          !$OMP END CRITICAL

      !$OMP END PARALLEL

    end if

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

                if (zp(i) .GE. in%zmax) CALL ReinjectParticles(zp(i),kep(i),xip(i),in,ecnt2,pcnt2)

                if (zp(i) .LE. in%zmin) CALL ReinjectParticles(zp(i),kep(i),xip(i),in,ecnt1,pcnt1)



                ! Apply Coulomb collision operator:

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

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

                if (in%iColl) then

    		!$OMP PARALLEL PRIVATE(i, id, ecnt, pcnt)

              ecnt = 0; pcnt = 0

              id = OMP_GET_THREAD_NUM()



                    ! "in" needs to be private to avoid race condition. This can be

                    ! fixed by looping over species inside the subroutine "collisionOperator"

                    in%species_b = 1

          		!$OMP DO

              		do i = 1,in%Nparts

                    !if (id .EQ. 0) write(*,*) "Thread", id, " at i = ", i

              			call collisionOperator(zp(i),kep(i),xip(i),ecnt,pcnt,in)

              		end do

          		!$OMP END DO



                    CALL collisionOperator(zp(i),kep(i),xip(i),ecnt4,pcnt4,in)

                    in%species_b = 2

          		!$OMP DO

              		do i = 1,in%Nparts

              		    call collisionOperator(zp(i),kep(i),xip(i),ecnt,pcnt,in)

              		end do

          		!$OMP END DO



          		!$OMP CRITICAL

              		ecount4(j) = ecount4(j) + ecnt

              		pcount4(j) = pcount4(j) + pcnt

          		!$OMP END CRITICAL



    		!$OMP END PARALLEL

                    CALL collisionOperator(zp(i),kep(i),xip(i),ecnt4,pcnt4,in)

                end if



                ! Apply RF heating operator:

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

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

                if (in%iHeat) then

      !$OMP PARALLEL PRIVATE(i, ecnt, pcnt, df, fnew)

          ecnt = 0; pcnt = 0; df = 0;

          !$OMP DO

              do i = 1,in%Nparts

                      call CyclotronResonanceNumber(zp(i),kep(i),xip(i),fnew(i),in,spline_Bz)

                  CALL CyclotronResonanceNumber(zp(i),kep(i),xip(i),fnew(i),in,spline_Bz)

                  df = dsign(1.d0,fcurr(i)*fnew(i))

                  if (df .LT. 0 .AND. zp(i) .GT. in%zRes1 .AND. zp(i) .LT. in%zRes2)  then

                        call RFHeatingOperator(zp(i),kep(i),xip(i),ecnt,pcnt,in,spline_Bz,spline_ddBz,spline_Phi)

                    CALL RFHeatingOperator(zp(i),kep(i),xip(i),ecnt3,pcnt3,in,spline_Bz,spline_ddBz,spline_Phi)

                  end if

              end do

                end if



              end do AllParticles

          !$OMP END DO

          !$OMP CRITICAL

              ecount3(j) = ecount3(j) + ecnt

              pcount3(j) = pcount3(j) + pcnt

          !$OMP END CRITICAL



          !$OMP CRITICAL AccumulateCounters

          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 AccumulateCounters



    !$OMP END PARALLEL

    end if



    ! Update time array:

    ! =========================================================================
@@ -385,7 +364,7 @@ oend = OMP_GET_WTIME() 
in%tComputeTime = oend-ostart

WRITE(*,*) ''

WRITE(*,*) '*********************************************************************'

print *, 'Reached End of Program, Computational time [s] = ', in%tComputeTime

WRITE(*,*) 'Reached End of Program, Computational time [s] = ', in%tComputeTime

WRITE(*,*) '*********************************************************************'

WRITE(*,*) ''
@@ -396,7 +375,7 @@ if (in%iSave) then 
    ! --------------------------------------------------------------------------

    dir1 = trim(in%rootDir)//'/OutputFiles'

    command = 'mkdir '//dir1

    call system(command,STATUS)

    CALL system(command,STATUS)

    WRITE(*,*) 'Status: ', STATUS



    ! Create subdirectory with input file name:
@@ -406,99 +385,99 @@ if (in%iSave) then 
    dir0 = dir0(1:n_mpwd)

    dir1 = trim(in%rootDir)//'/OutputFiles/'//trim(dir0)

    command = 'mkdir '//dir1

    call system(command,STATUS)

    CALL system(command,STATUS)



    ! Create subdirectory based on "fileDescriptor" inside the input file:

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

    dir1 = trim(dir1)//'/'//trim(in%fileDescriptor)

    command = 'mkdir '// dir1

    call system(command,STATUS)

    call getcwd(mpwd)

    CALL system(command,STATUS)

    CALL getcwd(mpwd)



    ! Saving zp_hist to file:

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

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

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

    write(8) zp_hist

    close(unit=8)

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

    WRITE(8) zp_hist

    CLOSE(unit=8)



    ! Saving kep_hist to file:

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

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

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

    write(8) kep_hist

    close(unit=8)

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

    WRITE(8) kep_hist

    CLOSE(unit=8)



    ! Saving xip_hist to file:

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

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

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

    write(8) xip_hist

    close(unit=8)

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

    WRITE(8) xip_hist

    CLOSE(unit=8)



    ! Saving t_hist to file:

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

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

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

    write(8) t_hist

    close(unit=8)

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

    WRITE(8) t_hist

    CLOSE(unit=8)



    ! Saving pcount to file:

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

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

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

    write(8) pcount1

    close(unit=8)

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

    WRITE(8) pcount1

    CLOSE(unit=8)

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

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

    write(8) pcount2

    close(unit=8)

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

    WRITE(8) pcount2

    CLOSE(unit=8)

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

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

    write(8) pcount3

    close(unit=8)

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

    WRITE(8) pcount3

    CLOSE(unit=8)

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

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

    write(8) pcount4

    close(unit=8)

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

    WRITE(8) pcount4

    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

    close(unit=8)

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

    WRITE(8) ecount1

    CLOSE(unit=8)

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

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

    write(8) ecount2

    close(unit=8)

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

    WRITE(8) ecount2

    CLOSE(unit=8)

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

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

    write(8) ecount3

    close(unit=8)

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

    WRITE(8) ecount3

    CLOSE(unit=8)

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

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

    write(8) ecount4

    close(unit=8)

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

    WRITE(8) ecount4

    CLOSE(unit=8)



    ! Write output metadata:

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

    fileName = trim(trim(dir1)//'/'//'metadata.txt')

    open(unit=8,file=fileName,form="formatted",status="unknown")

    write(8,NML = in_nml)

    close(unit=8)

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

    WRITE(8,NML = in_nml)

    CLOSE(unit=8)



    ! Copy input file to output directory:

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

    dir0 = trim(in%rootDir)//'/InputFiles/'//trim(inputFile)

    command = 'cp '//trim(trim(dir0)//' '//trim(dir1))

    call system(command)

    CALL system(command)



    ! Copy magnetic field data:

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

    dir0 = trim(in%rootDir)//'/BfieldData'//trim(in%BFieldFile)

    command = 'cp '//trim(trim(dir0)//' '//trim(dir1))//'/Bfield.txt'

    call system(command)

    CALL system(command)



    ! Create text file with commit Hash:

    ! --------------------------------------------------------------------------
@@ -508,7 +487,7 @@ if (in%iSave) then 
    dir0 = trim(in%rootDir)//'/OutputFiles/'//trim(dir0)//'/'//trim(in%fileDescriptor)

    fileName = trim(dir0)//'/commitHash.txt'

    command = 'git log --oneline -1 > '//trim(fileName)

    call system(command)

    CALL system(command)



end if

src/MoveParticlePack.f90

+1 −13
Original line number Diff line number Diff line
@@ -191,18 +191,6 @@ else if (in0%particleBC .EQ. 3) then 
  end if

end if



! Particle velocity standard deviation:

!sigma_u0 = sqrt(e_c*T0/m_t)



! Re-inject particle at source with new zp, kep, xip

!call random_number(Rm6)

!zp0 = in0%zp_init_std*sqrt(-2.*log(Rm6(1)))*cos(2.*pi*Rm6(2)) + in0%zp_init

!uper = sigma_u0*sqrt(-2.*log(Rm6(3)))

!upar = sigma_u0*sqrt(-2.*log(Rm6(4)))*cos(2.*pi*Rm6(5))

!u    = sqrt( uper**2 + upar**2 )

!kep0 = (m_t*u**2.)/(2.*e_c)

!xip0 = upar/u



return

END SUBROUTINE ReinjectParticles
@@ -231,7 +219,7 @@ q_t = in0%q_t 
upar = sqrt(2.*e_c*kep0/m_t)*xip0

! Magnetic field at location zp0 of test particle:

Bf = Interp1(zp0,spline0)

! Cyclotron frequenc of test particle:

! Cyclotron frequency of test particle:

Omega = abs(q_t)*Bf/m_t

! RF frequency in rad/s:

Omega_RF = 2*pi*in0%f_RF