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

succesfully added fueled constrained re-injection scheme

parent c38651bb

InputFiles/xp_TestDrag.in

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

! Simulation name:

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

params%fileDescriptor = '2021_02_22e',

params%fileDescriptor = '2021_02_25a',



! Magnetic field input data:

! =========================
@@ -15,6 +15,7 @@ params%nz = 501, 
params%NC = 70000,                                                              ! Total number of particles

params%NS = 10000,                                                              ! Total number of time steps

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

params%G  = 15.0E+21,                                                           ! Real particle injection rate 



! Time steps to record:

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

src/CoulombCollisions.f90

+28 −28
Original line number Diff line number Diff line
@@ -41,29 +41,29 @@ xip0 = plasma%xip(i) 
Ma = params%Ma

qa = params%qa

species_a = params%species_a

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

IF (species_a .EQ. 1) THEN      ! Test electrons

   Za = -1.

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

ELSE IF (species_a .EQ. 1) THEN ! Test ions

   Za = params%Zion

end if

END IF



species_b_loop: do ss = 1,2

species_b_loop: DO ss = 1,2

! Background species properties:

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

if (ss .EQ. 1) species_b = 1

if (ss .EQ. 2) species_b = 2

IF (ss .EQ. 1) species_b = 1

IF (ss .EQ. 2) species_b = 2



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

IF (species_b .EQ. 1) THEN       ! Background electron

  Mb = m_e

  Tb = params%Te0

  Zb = -1.

  nb = params%ne0

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

ELSE IF (species_b .EQ. 2) THEN  ! Background ion

  Mb = params%Aion*m_p

  Tb = params%Ti0

  Zb = params%Zion

  nb = params%ne0

end if

END IF



! Coordinate systems:

! -----------------------------------------------------------------------------------------------------------------
@@ -80,15 +80,15 @@ vpar = xip0*v 
vper = sqrt(1 - (xip0**2.) )*v



! Plasma drift in the lab frame:

if (params%iDrag) then

   if (zp0 .GE. params%BC_zp_mean) then

IF (params%iDrag) THEN

   IF (zp0 .GE. params%BC_zp_mean) THEN

      u = +1.*sqrt( e_c*(params%Te0 + params%Ti0)/(params%Aion*m_p) )

   else

   ELSE

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

   end if

else

   END IF

ELSE

   u = 0.

end if

END IF



! Plasma frame:

! ----------------------------------------------------------------------------------------------------------------
@@ -102,9 +102,9 @@ xip_pf_0 = wpar/w 
kep_pf_0 = (0.5*Ma/e_c)*w**2.



! Diagnostic:

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

IF (xip_pf_0**2. .GT. 1.) THEN

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

end if

END IF



! Test particle to background thermal velocity ratio:

wTb = sqrt(2.*e_c*Tb/Mb)
@@ -127,8 +127,8 @@ xip_pf_1 = S1 + S3 
! Apply energy scattering in plasma frame:

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

! Select mass term

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

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

IF (params%CollOperType .EQ. 1) mass_term = 1.                 ! Boozer-Only term

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



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

E0 = (1.5 + E_nuE_d_nu_E_dE(xab))*Tb
@@ -138,24 +138,24 @@ 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

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

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. params%elevel*Tb) then

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

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

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

	if(species_a .EQ. species_b) then

	IF (species_a .EQ. species_b) THEN

		!pcnt = pcnt + 1

                plasma%f4(i) = 1

	end if

end if

	END IF

END IF



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

w    = sqrt(2.*e_c*kep_pf_1/Ma)
@@ -173,7 +173,7 @@ v = sqrt( (vpar**2.) + (vper**2.) ) 
xip0 = vpar/v

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



end do species_b_loop

END DO species_b_loop



! Output data from calculation:

plasma%xip(i) = vpar/v

src/Modules.f90

+74 −69
Original line number Diff line number Diff line
@@ -227,6 +227,7 @@ TYPE paramsTYP 
  ! Simulation conditions:

  INTEGER(i4) :: NC, NS

  REAL(r8)    :: dt

  REAL(r8)    :: G

  ! Time steps to record:

  INTEGER(i4) :: jstart, jend, jincr

  ! Domain geometry:
@@ -258,12 +259,12 @@ 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 :: 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

 REAL(r8)   , DIMENSION(:), ALLOCATABLE :: Edot1, Edot2, Edot3, Edot4, Edot3_hat

 REAL(r8)   , DIMENSION(:), ALLOCATABLE :: dE1, dE2, dE3, dE4, dE5

 REAL(r8)   , DIMENSION(:), ALLOCATABLE :: udE3, udErf, doppler

 REAL(r8) :: NR , NSP, alpha

 REAL(r8) :: Eplus, Eminus

 REAL(r8) :: Ndot1, Ndot2, Ndot3, Ndot4

 REAL(r8) :: Edot1, Edot2, Edot3, Edot4, uEdot3

END TYPE plasmaTYP



! -----------------------------------------------------------------------------
@@ -275,7 +276,10 @@ END TYPE fieldSplineTYP 
TYPE outputTYP

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

 REAL(r8) , DIMENSION(:)  , ALLOCATABLE :: tp, jrng

 INTEGER(i4) :: jsize

 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

END TYPE outputTYP



CONTAINS
@@ -287,27 +291,18 @@ SUBROUTINE AllocatePlasma(plasma,params) 
   TYPE(paramsTYP), INTENT(IN)    :: params



   ! Declare local variables:

   INTEGER(i4) :: NC, NS

   INTEGER(i4) :: NC



   ! NC: Number of computational particles

   NC = params%NC

   ! NS: Number of time steps

   NS = params%NS

 

   ! 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%E1(NC) ,plasma%E2(NC)  ,plasma%E3(NC)  ,plasma%E4(NC))

   ALLOCATE(plasma%Erf_hat(NC))

   ALLOCATE(plasma%dE1(NC) ,plasma%dE2(NC) ,plasma%dE3(NC) ,plasma%dE4(NC), plasma%dE5(NC))

   ALLOCATE(plasma%udErf(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(plasma%Edot3_hat(NS))

   ALLOCATE(plasma%udE3(NC))



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

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

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

   INTEGER(i4) :: i, j

   

   ! Declare local variables:

   INTEGER(i4) :: i



   ! Derived parameters: Reference cross sectional area

   params%Area0 = 0.5*params%dtheta*( params%r2**2. - params%r1**2.)



   ! Derived parameters: Select test species

   if (params%species_a .EQ. 1) then

   IF (params%species_a .EQ. 1) THEN

       params%qa = -e_c

       params%Ma = m_e

   else

   ELSE

       params%qa = +params%Zion*e_c

       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(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 DO

   ! Initialize plasma: For all computational particles

   END IF

   

   ! Initialize plasma: scalar quantities:

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

   plasma%NSP     = params%NC

   plasma%alpha   = plasma%NR/plasma%NSP

   plasma%Eplus   = 0.

   plasma%Eminus  = 0.

   plasma%Ndot1   = 0.

   plasma%Ndot2   = 0.

   plasma%Ndot3   = 0.

   plasma%Ndot4   = 0.

   plasma%Edot1   = 0.

   plasma%Edot2   = 0.

   plasma%Edot3   = 0.

   plasma%Edot4   = 0.

   plasma%uEdot3  = 0.



   !$OMP PARALLEL DO

   ! Initialize plasma: zp, kep, xip and a:

   DO i = 1,params%NC

     plasma%a(i)  = 1.

     CALL loadParticles(i,plasma,params)

   END DO

   !$OMP END DO

   !$OMP END PARALLEL

   !$OMP END PARALLEL DO



END SUBROUTINE InitializePlasma
@@ -375,9 +364,17 @@ SUBROUTINE ResetFlags(i,plasma) 
   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%f1(i)      = 0. 

   plasma%f2(i)      = 0. 

   plasma%f3(i)      = 0. 

   plasma%f4(i)      = 0.

   plasma%dE1(i)     = 0. 

   plasma%dE2(i)     = 0. 

   plasma%dE3(i)     = 0. 

   plasma%dE4(i)     = 0.

   plasma%dE5(i)     = 0.

   plasma%udE3(i)    = 0.

   plasma%udErf(i)   = 0.

   plasma%doppler(i) = 0.;



END SUBROUTINE ResetFlags
@@ -430,22 +427,30 @@ SUBROUTINE AllocateOutput(output,params) 
   ! Declare interface variables:

   TYPE(outputTYP), INTENT(INOUT) :: output

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

   INTEGER(i4) :: j

   

   ! Declare local variables:

   INTEGER(i4) :: j, jsize, NS



   ! Determine size of temporal snapshots to record:

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

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

   

   ! Allocate memory:

   ALLOCATE(output%jrng(output%jsize))

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

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

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

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

   ALLOCATE(output%tp(output%jsize))

   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%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

   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))

  

END SUBROUTINE AllocateOutput



END MODULE dataTYP

src/MoveParticlePack.f90

+30 −24
Original line number Diff line number Diff line
@@ -162,15 +162,15 @@ xip0 = plasma%xip(i) 
! 1: Isotropic plasma source

! 2: NBI

! 3: Periodic boundary

if (params%BC_Type .EQ. 1 .OR. params%BC_Type .EQ. 2) then

IF (params%BC_Type .EQ. 1 .OR. params%BC_Type .EQ. 2) THEN



  if (params%BC_Type .EQ. 1) then

  IF (params%BC_Type .EQ. 1) THEN

    T = params%Te0

    E = 0

  else if (params%BC_Type .EQ. 2) then

  ELSE IF (params%BC_Type .EQ. 2) THEN

    T = params%BC_Tp

    E = params%BC_Ep

  end if

  END IF



  ! Velocity distribution:

  Ma = params%Ma
@@ -204,19 +204,19 @@ if (params%BC_Type .EQ. 1 .OR. params%BC_Type .EQ. 2) then 
  ! Position distribution:

  plasma%zp(i) = params%BC_zp_std*sqrt(-2.*log(Rm6(5)))*cos(2.*pi*Rm6(6)) + params%BC_zp_mean



else if (params%BC_Type .EQ. 3) then

  if (zp0 .GE. params%zmax) then

ELSE if (params%BC_Type .EQ. 3) THEN

  IF (zp0 .GE. params%zmax) THEN

    plasma%zp(i) = params%zmin

  else

  ELSE

    plasma%zp(i) = params%zmax

  end if

end if

  END IF

END IF



RETURN

END SUBROUTINE ReinjectParticles



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

SUBROUTINE CyclotronResonanceNumber(i,plasma,resNum0,fieldspline,params)

SUBROUTINE CyclotronResonanceNumber(i,plasma,fieldspline,params,resNum0)

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



USE local
@@ -296,7 +296,7 @@ Ma = params%Ma 
! Test particle charge:

qa = params%qa



! Calculate derived quantities

! Calculate derived quantities:

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

upar0    = u0*xip0

uper0    = u0*(1. - xip0**2)**0.5
@@ -319,12 +319,12 @@ Omega_dot = upar0*dOmega 
Omega_ddot = (upar0**2.)*ddOmega  - (uper0**2.)*dOmega*dOmega/(2.*Omega) - qa*dV*dOmega/Ma



! Calculate the interaction time (tau_RF):

if ( (Omega_ddot**2.) .GT. 4.8175*ABS(Omega_dot**3.) )  then

IF ( (Omega_ddot**2.) .GT. 4.8175*ABS(Omega_dot**3.) )  then

        ! Approximate Ai(x) ~ 0.3833

        tau_rf = (2.*pi)*(ABS(2./Omega_ddot)**(1/3.))*0.3833

else

ELSE

        tau_rf = sqrt(2.*pi/ABS(Omega_dot))

end if

END IF



! Calculate Bessel term:

rl       = uper0/(abs(qa)*Bf/Ma)
@@ -335,24 +335,25 @@ besselterm = BESSEL_JN(params%n_harmonic-1,flr) 
mean_dkep_per = 0.5*(e_c/Ma)*(besselterm*tau_rf)**2.



! Populate plasma structure:

plasma%Erf_hat(i) = mean_dkep_per

plasma%doppler(i) = params%kpar*(upar0)/Omega

plasma%E3_hat(i)  = plasma%Erf_hat(i)*(1. + plasma%doppler(i))

plasma%udErf(i)   = mean_dkep_per

plasma%doppler(i) = params%kpar*upar0/Omega

plasma%udE3(i)    = plasma%udErf(i)*(1. + plasma%doppler(i))



RETURN

END SUBROUTINE RFoperatorTerms



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

SUBROUTINE RFOperator(i,plasma,fieldspline,params,q3_hat)

SUBROUTINE RFOperator(i,plasma,fieldspline,params,uE3)

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

USE local

USE spline_fits

USE PhysicalConstants

USE dataTYP

USE OMP_LIB



IMPLICIT NONE

! Define interface variables:

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

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

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

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

TYPE(paramsTYP)        , INTENT(IN)    :: params
@@ -379,8 +380,13 @@ kep_par0 = kep0*xip0**2. 
kep_per0 = kep0*(1. - xip0**2.)



! Calculate the mean RF kick:

Ew2 = (params%Prf/q3_hat)

mean_dkep_per = plasma%Erf_hat(i)*Ew2

Ew2 = params%Prf/uE3

mean_dkep_per = plasma%udErf(i)*Ew2



IF (OMP_GET_THREAD_NUM() .EQ. 0) THEN

!   WRITE(*,*) 'Ew', sqrt(Ew2)

END IF





! Calculate the change in perp, parallel and total energy:

CALL RANDOM_NUMBER(Rm1)
@@ -421,7 +427,7 @@ plasma%xip(i) = upar1/u1 
! Record resonance event:

plasma%f3(i) = 1

! Record energy kick:

plasma%E3(i) = dkep

plasma%dE3(i) = dkep



RETURN

END SUBROUTINE RFOperator
@@ -453,8 +459,8 @@ SUBROUTINE loadParticles(i,plasma,params) 
  xip0 = plasma%xip(i)



  ! Particle position:

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

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

  zmin = params%zmin 

  zmax = params%zmax 

  if (params%IC_Type .EQ. 1) then

      ! Uniform load

      call random_number(X1)

src/linFP.f90

+125 −50

File changed.

Preview size limit exceeded, changes collapsed.