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

succesfully implemented a working power contrained RF operator. Results show... 

succesfully implemented a working power contrained RF operator. Results show that time average follows desired RF power but its spread is quite large
parent 4cfb8815

InputFiles/xp_TestDrag.in

+3 −3
Original line number Diff line number Diff line 
&params_nml

! Simulation name:

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

params%fileDescriptor = '2021_02_19d',

params%fileDescriptor = '2021_02_22d',



! Magnetic field input data:

! =========================
@@ -75,8 +75,8 @@ params%zRes1 = 0. 
params%zRes2      = +1.5                                                        ! Defines interval where RF is present

params%kper       = 125.                                                        ! Perpendicular wave number of EBW

params%kpar       = 30.                                                         ! Parallel wave number of EBW

params%Ew         = 15000                                                       ! Wave electric field magnitude [V/m]

params%Prf        = 300000                                                      ! Mean RF power absorbed via cyclotron interaction

params%Ew         = 10000                                                       ! Wave electric field magnitude [V/m]

params%Prf        = 2000000                                                      ! Mean RF power absorbed via cyclotron interaction

params%n_harmonic = 2                                                           ! Cyclotron harmonic



! Electric potential conditions:

src/MoveParticlePack.f90

+6 −17
Original line number Diff line number Diff line
@@ -158,10 +158,6 @@ zp0 = plasma%zp(i) 
kep0 = plasma%kep(i)

xip0 = plasma%xip(i)



! Record event:

!ecnt = ecnt + kep0

!pcnt = pcnt + 1



! Choose particle BC:

! 1: Isotropic plasma source

! 2: NBI
@@ -287,9 +283,6 @@ REAL(r8) :: Bf, Omega, dOmega, ddOmega 
REAL(r8) :: Omega_dot, Omega_ddot, tau_rf

REAL(r8) :: rl, flr, besselterm

REAL(r8) :: mean_dkep_per

!REAL(r8) :: dkep_par, dkep, kep1

!REAL(r8) :: kep_per1, kep_par1

!REAL(r8) :: upar1, u1

REAL(r8) :: Ma, qa



! Input variables:
@@ -338,7 +331,7 @@ rl = uper0/(abs(qa)*Bf/Ma) 
flr      = params%kper*rl

besselterm = BESSEL_JN(params%n_harmonic-1,flr)



!Calculate the mean RF kick per unit electric field:

!Calculate the mean RF kick per unit electric field squared:

mean_dkep_per = 0.5*(e_c/Ma)*(besselterm*tau_rf)**2.



! Populate plasma structure:
@@ -349,7 +342,7 @@ RETURN 
END SUBROUTINE RFoperatorTerms



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

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

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

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

USE local

USE spline_fits
@@ -358,7 +351,7 @@ USE dataTYP 


IMPLICIT NONE

! Define interface variables:

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

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

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

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

TYPE(paramsTYP)        , INTENT(IN)    :: params
@@ -366,12 +359,7 @@ TYPE(fieldSplineTYP) , INTENT(IN) :: fieldspline 


! Define local variables:

REAL(r8) :: kep0, xip0, Ew2

!REAL(r8) :: u0, upar0, uper0

REAL(r8) :: kep_par0, kep_per0

!REAL(r8) :: dB, ddB, dV

!REAL(r8) :: Bf, Omega, dOmega, ddOmega

!REAL(r8) :: Omega_dot, Omega_ddot, tau_rf

!REAL(r8) :: rl, flr, besselterm

REAL(r8) :: mean_dkep_per, dkep_per, Rm1

REAL(r8) :: dkep_par, dkep, kep1

REAL(r8) :: kep_per1, kep_par1
@@ -390,7 +378,7 @@ kep_par0 = kep0*xip0**2. 
kep_per0 = kep0*(1. - xip0**2.)



! Calculate the mean RF kick:

Ew2 = params%Ew**2.

Ew2 = (params%Prf/Prf_0)

mean_dkep_per = plasma%Erf_hat(i)*Ew2



! Calculate the change in perp, parallel and total energy:
@@ -415,6 +403,7 @@ kep_par1 = kep_par0 + dkep_par 
kep1     = kep_per1 + kep_par1



if (kep1 .LT. 0) then

  WRITE(*,*) 'final kep < 0'

  print *, 'kep0', kep0

  print *, 'dkep', dkep

  print *, 'kep1', kep1

src/linFP.f90

+26 −11
Original line number Diff line number Diff line
@@ -16,10 +16,10 @@ TYPE(plasmaTYP) :: plasma 
TYPE(fieldSplineTYP) :: fieldspline

TYPE(outputTYP)      :: output

! Accumulators:

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

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

REAL(r8) :: p1,p2,p3,p4,SP2RP

REAL(r8) :: q1,q2,q3,q4,Prf_0

! DO loop indices:

INTEGER(i4) :: i,j,k

INTEGER(i4) :: i,j,k,N3

! Thread ID:

INTEGER(i4) :: id

! OPENMP computational time:
@@ -68,6 +68,7 @@ WRITE(*,*) 'zDump [m]: ', params%zmin 
WRITE(*,*) 'BC_zp_mean [m]:     ', params%BC_zp_mean

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

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

WRITE(*,*) 'Prf:                ', params%Prf

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

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

IF (params%CollOperType .EQ. 1) WRITE(*,*) 'Boozer-Only collision operator'
@@ -109,10 +110,6 @@ fieldspline%U%y = 0. 
! Complete setting up the spline data :

CALL ComputeFieldSpline(fieldspline)



! Record start time:

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

ostart = OMP_GET_WTIME()



! Initialize simulation variables:

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

CALL InitializePlasma(plasma,params)
@@ -126,6 +123,10 @@ WRITE(*,*) "Number of threads given: ", params%threads_given 
WRITE(*,*) '*********************************************************************'

WRITE(*,*) ''



! Record start time:

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

ostart = OMP_GET_WTIME()



! Initialize simulation diagnostics:

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

! Record time index:
@@ -143,6 +144,9 @@ NS_loop: DO j = 1,params%NS 
    ! Reset resonance number:

    dresNum = 0.; resNum0 = 0.; resNum1 = 0.



    ! Reset reference RF power:

    Prf_0 = 0.



    ! Calculate scale factor: super-particle to real-particle

    SP2RP = plasma%NR(j)/plasma%NSP(j)
@@ -181,7 +185,7 @@ NS_loop: DO j = 1,params%NS 
           dresNum = dsign(1.d0,resNum0*resNum1)

           IF (dresNum .LT. 0 .AND. plasma%zp(i) .GT. params%zRes1 .AND. plasma%zp(i) .LT. params%zRes2)  THEN

              plasma%f3(i) = 1

              CALL RFoperatorTerms(i,plasma,params,fieldspline)

              CALL RFoperatorTerms(i,plasma,fieldspline,params)

           END IF

        END IF
@@ -189,12 +193,21 @@ NS_loop: DO j = 1,params%NS 
    !$OMP END DO



    ! Calculate RF power per unit electric field:

    !$OMP DO REDUCTION(+:Edot3_hat)

    !$OMP DO REDUCTION(+:Prf_0)

    NC_loop2: DO i = 1,params%NC

       Edot3_hat = Edot3_hat + (SP2RP/params%dt)*plasma%a(i)*plasma%f3(i)*plasma%Erf_hat(i)*(1 + plasma%doppler(i))

       Prf_0 = Prf_0 + e_c*(SP2RP/params%dt)*plasma%a(i)*plasma%f3(i)*plasma%Erf_hat(i)*(1 + plasma%doppler(i))

    END DO NC_loop2

    !$OMP END DO

  

    !IF (OMP_GET_THREAD_NUM() .EQ. 0) THEN

       !IF (ISNAN(Prf_0)) THEN

        !   WRITE(*,*) 'N3:', N3

      ! ELSE 

        !   WRITE(*,*) 'Prf_0:', Prf_0

       !    WRITE(*,*) 'N3:', N3          

      ! END IF

   ! END IF



    ! Apply RF operator and particle re-injection:

    !$OMP DO

    NC_loop3: DO i = 1,params%NC
@@ -202,7 +215,7 @@ NS_loop: DO j = 1,params%NS 
           CALL ReinjectParticles(i,plasma,params)

       END IF

       IF (plasma%f3(i) .EQ. 1) THEN

           CALL RFOperator(i,Edot3_hat,plasma,fieldspline,params)

           CALL RFOperator(i,Prf_0,plasma,fieldspline,params)

       END IF

    END DO NC_loop3

    !$OMP END DO
@@ -222,6 +235,8 @@ NS_loop: DO j = 1,params%NS 


    !$OMP END PARALLEL



!    WRITE(*,*) 'Ew: ', sqrt(params%Prf/(e_c*Edot3_hat))



    ! Calculate particle and energy rates in physical units [P/s] and [J/s]

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

    plasma%Ndot1(j) = SP2RP*p1/params%dt