Commit 7b55e9e8 authored by Juan Caneses Marin (nfc)'s avatar Juan Caneses Marin (nfc)
Browse files

Succesfully optmized openMP up to x20 with 48 threads with collisions and... 

Succesfully optmized openMP up to x20 with 48 threads with collisions and verified against analytic theory
parent e4ef9d6c

InputFiles/xp_IonSlowingDown.in

+1 −2
Original line number Diff line number Diff line
@@ -12,7 +12,7 @@ in%nz = 501, ! N 


! Simulation conditions:

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

in%Nparts = 5000,                                                              ! Total number of particles

in%Nparts = 50000,                                                              ! Total number of particles

in%Nsteps = 5000,                                                               ! Total number of time steps

in%dt     = 0.5E-7,                                                             ! Time step in [s]

in%zmax = +5.0,
@@ -37,7 +37,6 @@ in%Te0 = 10., ! 
in%Ti0 = 10.,                                                                  ! eV

in%ne0 = 1.0E+19,	                                                             ! m**-3

in%Aion = 1.,                                                                  ! Ion mass

in%Zeff = 1.,                                                                  ! Zeff for ion/impurity mixture

in%Zion = 1.,                                                                  ! Charge number for main ion species

in%species_a = 2,                                                              ! =1 for electrons, =2 for ions

in%elevel = 10,                                                                ! Energy at which slowing down power is calculated E = elevel*Te0

InputFiles/xp_Test.in

+0 −1
Original line number Diff line number Diff line
@@ -37,7 +37,6 @@ in%Te0 = 100., 
in%Ti0 = 100.,                                                                  ! eV

in%ne0 = 5.0E+19,	                                                             ! m**-3

in%Aion = 1.,                                                                  ! Ion mass

in%Zeff = 1.,                                                                  ! Zeff for ion/impurity mixture

in%Zion = 1.,                                                                  ! Charge number for main ion species

in%species_a = 2,                                                              ! =1 for electrons, =2 for ions

in%elevel = 10,                                                                ! Energy at which slowing down power is calculated E = elevel*Te0

InputFiles/xp_TestCollisions.in

0 → 100644
+77 −0
Original line number Diff line number Diff line 
&in_nml

! Simulation name:

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

in%fileDescriptor = '2021_02_09b',



! Magnetic field input data:

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

in%repoDir ="/home/nfc/myRepos/LinearFokkerPlanck_Axisymmetric",

in%BFieldFileDir = "/BfieldData",

in%BFieldFile = "/Bfield_b.txt",

in%nz = 501,			                                                ! Number of points in BFieldFile



! Simulation conditions:

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

in%Nparts = 100000,                                                               ! Total number of particles

in%Nsteps = 5000,                                                               ! Total number of time steps

in%dt     = 0.5E-7,                                                             ! Time step in [s]

in%zmax   = +5.0,

in%zmin   = -5.0,



in%iSave      = .true.,

in%iPush      = .false.,                                                       ! Enable RK4 particle pusher

in%iColl      = .true.,                                                        ! Collisions or no collisions

in%iHeat      = .false.,                                                       ! Enable RF heating operator

in%iPotential = .false.,                                                       ! Enable electric potential

in%iDrag      = .false.,                                                       ! Enable plasma drag



! Time steps to record:

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

in%jstart = 1,                                                                 ! start frame

in%jend   = 5000,                                                              ! end frame

in%jincr  = 50,                                                                ! increment



! Collision operator conditions:

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

in%Te0 = 10.,	                                                               ! eV

in%Ti0 = 10.,                                                                  ! eV

in%ne0 = 1.0E+19,	                                                             ! m**-3

in%Aion = 1.,                                                                  ! Ion mass

in%Zion = 1.,                                                                  ! Charge number for main ion species

in%species_a = 2,                                                              ! =1 for electrons, =2 for ions

in%elevel = 10,                                                                ! Energy at which slowing down power is calculated E = elevel*Te0

in%CollOperType = 2,                                                           ! 1 for Boozer-Only, 2 for Boozer-Kim



! Particle boundary conditions:

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

in%particleBC = 1,                                                             ! 1: Isotropic plasma source, 2: NBI, 3: Periodic



! Initial conditions:

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

in%zp_InitType = 2,                                                            ! 1: uniform load, 2: gaussian load

in%zp_init = 0.0,                                                              ! Mean particle injection

in%zp_init_std = 0.3,                                                          ! STD of particle injection

in%Ep_init = 1000.,                                                            ! Drift kinetic energy

in%Tp_init = 0.1,                                                              ! Thermal kinetic energy

in%xip_init = 0.9240,                                                          ! Mean pitch angle where xip = cos(theta) = vpar/v



! RF heating operator conditions:

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

in%f_RF = 28.0E+9                                                              ! RF frequency

in%zRes1 = 2.4                                                                 ! Defines interval where RF is present

in%zRes2 = 3.0                                                                 ! Defines interval where RF is present

in%kper = 29920.                                                               ! Perpendicular wave number of EBW

in%kpar = 5864.                                                                ! Parallel wave number of EBW

in%Ew   = 10000                                                                ! EBW E-minus wave electric field magnitude

in%n_harmonic = 4                                                              ! Cyclotron harmonic



! Electric potential conditions:

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

in%s1 = 0,

in%s2 = 1.3,

in%s3 = 4.0,

in%phi1 = 0.,

in%phi2 = 0,

in%phi3 = 0.,



/

src/CoulombCollisions.f90

+250 −141
Original line number Diff line number Diff line
@@ -7,137 +7,103 @@ USE PhysicalConstants 
USE dataTYP



IMPLICIT NONE

TYPE(inTYP)  :: in0

REAL(r8) :: xip0, kep0, zp0, ecnt, pcnt

REAL(r8) :: xip1, kep1

REAL(r8) :: u, xb

REAL(r8) :: ln_lambda, Tb

REAL(r8) :: phi, phip, phi2p, Gb

REAL(r8) :: nuab0, nu_s, nu_D, nu_E, nu_prll

REAL(r8) :: E_nuE_d_nu_E_dE

REAL(r8) :: xsi1, xsi2, ran_num

REAL(r8) :: d1, d2, S1, S2, E0, E1, E2, S3, mass_term



REAL(r8) :: upar, uper

REAL(r8) :: vpar, vper, v, Cs, vthb, vthb3

! Define interface variables:

TYPE(inTYP), INTENT(IN) :: in0

REAL(r8), INTENT(IN) :: zp0

REAL(r8), INTENT(INOUT) :: xip0, kep0, ecnt, pcnt



! Define local variables:

INTEGER(i4) :: species_a, species_b, ss

REAL(r8) :: Za, Ma, qa

REAL(r8) :: Zb, Mb, Tb, nb

REAL(r8) :: u, w, wpar, wper, v, vpar, vper

REAL(r8) :: xip_pf_0, kep_pf_0

REAL(r8) :: xip_pf_1, kep_pf_1

REAL(r8) :: xab, wTb

REAL(r8) :: nu_ab0, nu_s, nu_D, nu_E, nu_prll

REAL(r8) :: xsi1, xsi2, ran_num

REAL(r8) :: d1, d2, S1, S2, E0, E1, E2, S3, mass_term

REAL(r8) :: dE_pf, dE_lf

REAL(r8) :: m_t, q_t

INTEGER(i4) :: species_a, species_b

REAL(r8) :: Za, Zb, Ma, Mb



! Test particle mass:

m_t = in0%m_t



! Test particle charge:

q_t = in0%q_t

! Define functions:

REAL(r8) :: lnA, phi, phip, phi2p, Gb, E_nuE_d_nu_E_dE



! Species to use:

! Test particle properties:

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

Ma = in0%Ma

qa = in0%qa

species_a = in0%species_a

species_b = in0%species_b

if (species_a .EQ. 1) then      ! Test electrons

   Za = -1.

else if (species_a .EQ. 1) then ! Test ions

   Za = in0%Zion

end if



species_b_loop: do ss = 1,2

! Background species properties:

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

if (ss .EQ. 1) species_b = 1

if (ss .EQ. 2) species_b = 2



! In the following, Za needs to be specified in "data.in"

if(species_b .EQ. 1) then       ! Field electron

if(species_b .EQ. 1) then       ! Background electron

  Mb = m_e

  Tb = in0%Te0

	Za = 1.; Zb = 1.

else if(species_b .EQ. 2) then  ! Field ion

  Zb = -1.

  nb = in0%ne0

else if(species_b .EQ. 2) then  ! Background ion

  Mb = in0%Aion*m_p

  Tb = in0%Ti0

	Za = in0%Zion; Zb = in0%Zion

  Zb = in0%Zion

  nb = in0%ne0

end if



! Particle velocities in lab frame:

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

! Input: kep(i), xip(i)

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

upar = xip0*u

uper = sqrt(1 - (xip0**2.) )*u

! Coordinate systems:

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

! We make use of two frames: (1) laboratory frame and (2) plasma frame.

! In the non-relativistic limit, the particle velocity in the lab frame is the following:

! v_vec = u_vec + w_vec

! Where u_vec is the mean drift velocit vector and w_vec is the thermal motion relative

! to the mean drift



! Lab frame: 

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

v = sqrt(2.*e_c*kep0/Ma)

vpar = xip0*v

vper = sqrt(1 - (xip0**2.) )*v



! Plasma drift in the lab frame:

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

if (in0%iDrag) then

   if (zp0 .GE. in0%zp_init) then

		Cs = +1.*sqrt(e_c*(in0%Te0 + in0%Ti0)/(in0%Aion*m_p))

      u = +1.*sqrt(e_c*(in0%Te0 + in0%Ti0)/Mb)

   else

		Cs = -1.*sqrt(e_c*(in0%Te0 + in0%Ti0)/(in0%Aion*m_p))

      u = -1.*sqrt(e_c*(in0%Te0 + in0%Ti0)/Mb)

   end if

else

	Cs = 0.

   u = 0.

end if



! Plasma frame:

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

! Convert quantities to plasma frame:

vper = uper

vpar = upar - Cs

v = sqrt( (vpar**2.) + (vper**2.) )

xip_pf_0 = vpar/v

kep_pf_0 = (0.5*m_t/e_c)*v**2.



! Define initial pitch and energy in plasma frame:

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

!xip0 = xip_pf

!kep0 = kep_pf



! Diagnostic

if (xip_pf_0**2 .GT. 1) then

	!print *,'Start Sub: xip > 1', xip(i)

! Thermal velocities in the plasma frame:

wper = vper

wpar = vpar - u

w = sqrt( (wpar**2.) + (wper**2.) )



! Initial energy and pitch angle in the plasma frame:

xip_pf_0 = wpar/w

kep_pf_0 = (0.5*Ma/e_c)*w**2.



! Diagnostic:

if (xip_pf_0**2. .GT. 1.) then

	print *,'xip_pf_0^2 > 1', xip_pf_0

end if



! Test particle to background thermal velocity ratio:

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

vthb = sqrt(2.*e_c*Tb/Mb)

xb = v/vthb



! Define Chandrasekhar functions:

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

phi = erf(xb)

phip = (2./sqrt(pi))*exp(-xb*xb)

phi2p = -(4.*xb/sqrt(pi))*exp(-xb*xb)

Gb = (phi - xb*phip)/(2.*xb*xb)



! Fundamental collision rate:

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

vthb3 = vthb**3.

ln_lambda = 30. - log(sqrt(in0%ne0)*in0%Te0**(-3/2))

nuab0 = in0%ne0*(e_c**4.)*((Za*Zb)**2.)*ln_lambda/(2.*pi*m_t*m_t*e_0*e_0*vthb3)



! Perpendicular deflection rate:

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

nu_D = nuab0*(phi - Gb)/(xb**3)



! Slowing down rate:

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

nu_s = nuab0*(1. + (m_t/Mb))*Gb/xb



! Parallel diffusion rate:

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

nu_prll = nuab0*Gb/(xb**3)



! Energy loss rate:

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

if (.false.) then

	! From Hinton 1983 EQ 92 and L. Chen 1988 EQ 50

	nu_E = nuab0*( (2*(m_t/Mb)*Gb/xb) - (phip/(xb**2)) )

else

	! From L. Chen 1983 EQ 57 commonly used for NBI

	nu_E = nuab0*(2.*(m_t/Mb))*(Gb/xb)

end if



! Energy loss derivative:

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

E_nuE_d_nu_E_dE = 0.5*(  ( 3*xb*phip - 3*phi - (xb**2)*phi2p )/( phi - (xb*phip) )  )



! Special case for electron-ion-impurity:need only pitch angle scatt.:

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

if(species_a .eq. 1 .and. species_b .eq. 2) then

	nu_D = nuab0*in0%Zeff*(phi - Gb)/(xb**3)

end if

wTb = sqrt(2.*e_c*Tb/Mb)

xab = w/wTb



! Generate random numbers:

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

call random_number(ran_num)

xsi1 = sign(1._r8, ran_num - 0.5_r8)

call random_number(ran_num)
@@ -145,41 +111,35 @@ xsi2 = sign(1._r8, ran_num - 0.5_r8) 


! Apply pitch angle scattering in plasma frame:

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

d1 = nu_D*in0%dt

d1 = nu_D(xab,nb,Tb,Mb,Zb,Za,Ma)*in0%dt

S1 = xip_pf_0*(1. - d1)

S2 = ( 1. - xip_pf_0*xip_pf_0 )*d1

S3 = xsi1*sqrt(S2)



xip_pf_1 = S1 + S3



! Diagnotics

if (isnan(xip_pf_1)) xip_pf_1 = xip_pf_0

if  (xip_pf_1**2 .GT. 1) then

	call random_number(ran_num)

	xip_pf_1 = 2.*ran_num - 1.

end if



! Apply energy scattering in plasma frame:

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

! Select mass term

if (in0%CollOperType .EQ. 1) mass_term = 1.                 ! Boozer-Only term

if (in0%CollOperType .EQ. 2) mass_term = 1. + (Mb/m_t)   ! Boozer-Kim term

if (in0%CollOperType .EQ. 2) mass_term = 1. + (Mb/Ma)   ! Boozer-Kim term



d2 = nu_E*in0%dt/mass_term

E0 = (1.5 + E_nuE_d_nu_E_dE)*Tb

d2 = nu_E(xab,nb,Tb,Mb,Zb,Za,Ma)*in0%dt/mass_term

E0 = (1.5 + E_nuE_d_nu_E_dE(xab))*Tb

E1 = d2*(kep_pf_0 - E0)

E2 = sqrt(Tb*kep_pf_0*d2)



dE_pf = 2.*E1



kep_pf_1 = kep_pf_0 - dE_pf + (2.*xsi2*E2)



! Diagnotics

if (isnan(xip_pf_1)) xip_pf_1 = xip_pf_0

if  (xip_pf_1**2. .GT. 1.) then

	call random_number(ran_num)

	xip_pf_1 = 2.*ran_num - 1.

end if

if (kep_pf_1 .le. 0.) kep_pf_1 = kep_pf_0



! Record energy loss due to collisions in the plasma frame:

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

if (kep_pf_1 .GT. in0%elevel*in0%Te0 ) then

if (kep_pf_1 .GT. in0%elevel*Tb) then

	! Record slowing down energy during time step dt

	ecnt = ecnt + dE_pf

	! Count how many particles are involved in the slowing down power calculation
@@ -188,23 +148,172 @@ if (kep_pf_1 .GT. in0%elevel*in0%Te0 ) then 
	end if

end if



! End of plasma frame calculation

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

! Based on the modifed xip_pf and kep_pf, we get:

w    = sqrt(2.*e_c*kep_pf_1/Ma)

wpar = xip_pf_1*w

wper = sqrt(1. - (xip_pf_1**2.) )*w



! Convert back to lab frame:

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

! Lab frame:

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

dE_lf = dE_pf

v = sqrt(2.*e_c*kep_pf_1/m_t)

vpar = xip_pf_1*v

vper = sqrt( 1 - (xip_pf_1**2.) )*v

upar = vpar + Cs

uper = vper

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

vpar  = u + wpar

vper  = wper

v     = sqrt( (vpar**2.) + (vper**2.) )



! New pitch and energy in the lab frame:

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

xip0 = upar/u

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

xip0 = vpar/v

kep0 = 0.5*(Ma/e_c)*v**2.



return

end do species_b_loop 



RETURN

END SUBROUTINE collisionOperator





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

REAL(r8) FUNCTION phi(x)

USE local



IMPLICIT NONE

REAL(r8), INTENT(IN) :: x



phi = erf(x)



END FUNCTION phi



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

REAL(r8) FUNCTION phip(x)

USE local

USE PhysicalConstants



IMPLICIT NONE

REAL(r8), INTENT(IN) :: x



phip = (2./sqrt(pi))*exp(-x*x)



END FUNCTION phip



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

REAL(r8) FUNCTION phi2p(x)

USE local

USE PhysicalConstants



IMPLICIT NONE

REAL(r8), INTENT(IN) :: x



phi2p = -(4.*x/sqrt(pi))*exp(-x*x)



END FUNCTION phi2p



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

REAL(r8) FUNCTION Gb(x)

USE local



IMPLICIT NONE

REAL(r8), INTENT(IN) :: x

REAL(r8) :: phi, phip



Gb = (phi(x) - x*phip(x))/(2.*x*x)



END FUNCTION Gb



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

REAL(r8) FUNCTION lnA(nb,Tb)

USE local



IMPLICIT NONE

REAL(r8), INTENT(IN) :: nb,Tb



lnA = 30. - log( sqrt( nb*Tb**(-3./2.) ) )



END FUNCTION lnA



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

REAL(r8) FUNCTION nu_ab0(nb,Tb,Mb,Zb,Za,Ma)

! Fundamental collision rate:

USE local

USE PhysicalConstants



IMPLICIT NONE

REAL(r8), INTENT(IN) :: nb,Tb,Mb,Zb,Za,Ma

REAL(r8) :: lnA

REAL(r8) :: wTb



wTb = sqrt(2.*e_c*Tb/Mb)



nu_ab0 = nb*(e_c**4.)*((Za*Zb)**2.)*lnA(nb,Tb)/(2.*pi*Ma*Ma*e_0*e_0*wTb**3.)



END FUNCTION nu_ab0



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

REAL(r8) FUNCTION E_nuE_d_nu_E_dE(x)

USE local



IMPLICIT NONE

REAL(r8), INTENT(IN) :: x

REAL(r8) :: phi, phip, phi2p



E_nuE_d_nu_E_dE = 0.5*(  ( 3.*x*phip(x) - 3.*phi(x) - (x**2.)*phi2p(x) )/( phi(x) - (x*phip(x)) )  )



END FUNCTION E_nuE_d_nu_E_dE



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

REAL(r8) FUNCTION nu_D(x,nb,Tb,Mb,Zb,Za,Ma)

USE local



IMPLICIT NONE

REAL(r8), INTENT(IN) :: x

REAL(r8), INTENT(IN) :: nb,Tb,Mb,Zb,Za,Ma

REAL(r8) :: nu_ab0, phi, Gb



nu_D = nu_ab0(nb,Tb,Mb,Zb,Za,Ma)*(phi(x) - Gb(x))/(x**3.)



END FUNCTION nu_D



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

REAL(r8) FUNCTION nu_s(x,nb,Tb,Mb,Zb,Za,Ma)

USE local



IMPLICIT NONE

REAL(r8), INTENT(IN) :: x

REAL(r8), INTENT(IN) :: nb,Tb,Mb,Zb,Za,Ma

REAL(r8) :: nu_ab0, Gb



nu_s = nu_ab0(nb,Tb,Mb,Zb,Za,Ma)*(1. + (Ma/Mb))*Gb(x)/x



END FUNCTION nu_s



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

REAL(r8) FUNCTION nu_prll(x,nb,Tb,Mb,Zb,Za,Ma)

USE local



IMPLICIT NONE

REAL(r8), INTENT(IN) :: x

REAL(r8), INTENT(IN) :: nb,Tb,Mb,Zb,Za,Ma

REAL(r8) :: nu_ab0, Gb



nu_prll = nu_ab0(nb,Tb,Mb,Zb,Za,Ma)*Gb(x)/(x**3.)



END FUNCTION nu_prll



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

REAL(r8) FUNCTION nu_E(x,nb,Tb,Mb,Zb,Za,Ma)

USE local



IMPLICIT NONE

REAL(r8), INTENT(IN) :: x

REAL(r8), INTENT(IN) :: nb,Tb,Mb,Zb,Za,Ma

REAL(r8) :: nu_ab0, Gb, phip



! Energy loss rate:

if (.false.) then

	! From Hinton 1983 EQ 92 and L. Chen 1988 EQ 50

	nu_E = nu_ab0(nb,Tb,Mb,Zb,Za,Ma)*( (2.*(Ma/Mb)*Gb(x)/x) - (phip(x)/(x**2.)) )

else

	! From L. Chen 1983 EQ 57 commonly used for NBI

	nu_E = nu_ab0(nb,Tb,Mb,Zb,Za,Ma)*(2.*(Ma/Mb))*(Gb(x)/x)

end if



END FUNCTION nu_E

src/Modules.f90

+2 −2
Original line number Diff line number Diff line
@@ -37,7 +37,7 @@ TYPE inTYP 
  CHARACTER*150 :: fileDescriptor, repoDir

  CHARACTER*150 :: BFieldFile, BFieldFileDir

  REAL(r8) :: Ti0, Te0, ne0, dt

  REAL(r8) :: Aion, Zeff, Zion

  REAL(r8) :: Aion, Zion

  INTEGER(i4) :: Nparts, Nsteps, nz, species_a, species_b

  INTEGER(i4) :: jstart, jend, jincr

  INTEGER(i4) :: threads_given, particleBC
@@ -51,7 +51,7 @@ TYPE inTYP 
  REAL(r8) :: f_RF, kpar, kper, Ew, zRes1, zRes2

  INTEGER(i4) :: n_harmonic

  REAL(r8) :: tComputeTime, tSimTime                              ! Variables to hold cpu time at start and end of computation

  REAL(r8) :: m_t, q_t

  REAL(r8) :: Ma, qa

END TYPE inTYP



END MODULE dataTYP