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

Modified how particles are loaded, removed kep_initType, xip_initType. Loading... 

Modified how particles are loaded, removed kep_initType, xip_initType. Loading function is more general
parent 6f540b1c

InputFiles/xp_IonSlowingDown.in

+4 −8
Original line number Diff line number Diff line 
&in_nml

in%fileDescriptor = '2021_01_27k',

in%fileDescriptor = '2021_01_29a',



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

in%BFieldFileDir = "/BfieldData",                                              ! Location of file
@@ -39,13 +39,9 @@ in%CollOperType = 2, ! 
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%kep_InitType = 2,                                                           ! 1: Maxwellian, 2: Beam at E = kep_init

in%kep_init = 1000.,                                                           ! Initial energy of EEDF



in%xip_InitType = 2,                                                           ! 1: Isotropic, 2: Beam at xip = xip_init

in%xip_init = 0.707,                                                           ! Mean value of xip = upar/u

in%xip_init_std = 0.98,                                                        ! STD in xip

in%Ep_init = 1000.,                                                            ! Drift kinetic energy

in%Tp_init_std = 10.;                                                          ! Thermal kinetic energy

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



in%f_RF = 28.0E+9                                                              ! RF frequency

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

src/Modules.f90

+2 −2
Original line number Diff line number Diff line
@@ -28,8 +28,8 @@ TYPE inTYP 
  INTEGER(i4) :: jstart, jend, jincr

  INTEGER(i4) :: threads_request, threads_given, particleBC

  LOGICAL:: iDrag, iPotential, iSave, iPush, iHeat, iColl

  INTEGER(i4) :: zp_InitType,kep_InitType, xip_InitType

  REAL(r8) :: zp_init, kep_init, xip_init, zp_init_std, xip_init_std

  INTEGER(i4) :: zp_InitType

  REAL(r8) :: zp_init, Ep_init, xip_init, zp_init_std, Tp_init

  REAL(r8) :: zmin, zmax

  INTEGER(i4) :: CollOperType

  REAL(r8) :: elevel

src/MoveParticlePack.f90

+46 −37
Original line number Diff line number Diff line
@@ -332,53 +332,62 @@ SUBROUTINE loadParticles(zp,kep,xip,in0) 
  ! Declare internal variables:

  TYPE(inTYP)  :: in0

  REAL(r8), DIMENSION(in0%Nparts) :: zp, kep, xip

  REAL(r8), DIMENSION(in0%Nparts) :: RmArray1, RmArray2, RmArray3

  REAL(r8), DIMENSION(in0%Nparts) :: uperArray, uparArray, uArray

  REAL(r8) :: zmin, zmax, sigma_u_init, m_t

  REAL(r8), DIMENSION(in0%Nparts) :: X1, X2, X3, X4

  REAL(r8), DIMENSION(in0%Nparts) :: R1, R2, R3, R4, t2, t4

  REAL(r8), DIMENSION(in0%Nparts) :: wx, wy, wz, vx, vy, vz, v

  REAL(r8) :: zmin, zmax, sigma_v, m_t

  REAL(r8) :: Ux, Uy, Uz, vT, U, T, E



  ! Particle position:

  zmin = in0%zmin + .01*(in0%zmax-in0%zmin)

  zmax = in0%zmax - .01*(in0%zmax-in0%zmin)

  if (in0%zp_InitType .EQ. 1) then

      ! Uniform load

      call random_number(RmArray1)

      zp = zmin + (zmax - zmin)*RmArray1

      call random_number(X1)

      zp = zmin + (zmax - zmin)*X1

  elseif (in0%zp_InitType .EQ. 2) then

      ! Gaussian load

      call random_number(RmArray1)

      call random_number(RmArray2)

      zp = in0%zp_init_std*sqrt(-2.*log(RmArray1))*cos(2.*pi*RmArray2)  +  in0%zp_init

      call random_number(X1)

      call random_number(X2)

      zp = in0%zp_init_std*sqrt(-2.*log(X1))*cos(2.*pi*X2)  +  in0%zp_init

  end if



  ! Particle kinetic energy:

  call random_number(RmArray1)

  call random_number(RmArray2)

  call random_number(RmArray3)

  ! Particle kinetic energy and pitch angle:

  ! Generate random numbers:

  call random_number(X1)

  call random_number(X2)

  call random_number(X3)

  call random_number(X4)



  ! Derived quantities:

  m_t = in0%m_t

  sigma_u_init = sqrt(e_c*in0%kep_init/m_t)

  uperArray = sigma_u_init*sqrt(-2.*log(RmArray1))

  uparArray = sigma_u_init*sqrt(-2.*log(RmArray2))*cos(2.*pi*RmArray3)

  uArray    = sqrt( uperArray**2 + uparArray**2 )



  if (in0%kep_InitType .EQ. 1) then

      ! Maxwellian EEDF:

      kep = (m_t*uArray**2.)/(2.*e_c)

  elseif (in0%kep_InitType .EQ. 2) then

      ! Beam EEDF:

      kep = in0%kep_init

  end if

  T = in0%Tp_init

  E = in0%Ep_init

  vT = sqrt(2.*e_c*T/m_t)

  U  = sqrt(2.*e_c*E/m_t)

  Ux = 0

  Uy = U*sqrt(1. - in0%xip_init**2.)

  Uz = U*in0%xip_init

  sigma_v = vT/sqrt(2.)



  ! Box-muller:

  R1 = sigma_v*sqrt(-2.*log(X1))

  t2 = 2.*pi*X2

  R3 = sigma_v*sqrt(-2.*log(X3))

  t4 = 2.*pi*X4



  wx = R1*cos(t2)

  wy = R1*cos(t2)

  wz = R3*cos(t4)



  vx = Ux + wx

  vy = Uy + wy

  vz = Uz + wz

  v = sqrt( vx**2. + vy**2. + vz**2. )



  ! Populate output variables:

  kep = 0.5*(m_t/e_c)*v**2

  xip = vz/v



  ! Particle pitch angle:

  if (in0%xip_InitType .EQ. 1) then

      ! Isotropic pitch angle

      xip = uparArray/uArray

  elseif (in0%xip_InitType .EQ. 2) then

      ! Beam pitch angle

      !xip = in0%xip_init

      ! Beam pitch angle with a Gaussian distribution in pitch angle space:

      call random_number(RmArray1)

      call random_number(RmArray2)

      xip = in0%xip_init_std*sqrt(-2.*log(RmArray1))*cos(2.*pi*RmArray2)  +  in0%xip_init

  end if

return

END SUBROUTINE loadParticles