Input file for Kotschenreuther-Dorland-Liu Gyrokinetic Code GS2

GWH 99/08/10: With comments to document input variables.
(You may insert text between namelists.)

miller12 based on parameters in Miller's paper, and in Table I of Waltz
and Miller, PoP (1999).

More documentation is in the comments near beginning of the geometry.f90
routine. 

This case uses a Miller local equilibrium.
Should be able to recover s-alpha (at r/R=0) by making rhoc very small, 
and turning off elongation, etc.

Also note, that for electromagnetic runs it is important that beta
beta_prime, and R/LT etc. be consistent.

!#########################################################################

Macroscale lengths are normalized to some macroscopic length "L".  The way
I usually run the code (see the discussion in docs/units.txt), k_theta is
normalized to rho_i, and frequencies are normalized to v_t/L, where v_t is
a reference ion thermal speed sqrt(T_ref/m_ref) at a reference temperature
and a reference mass, and rho_i=v_t/omega_c_ref.

The macroscopic scale length "L" does not need to be the plasma edge minor
radius a_minor. You can play games by setting Rmaj=R_geo=1 (this
corresponds to setting the macroscopic scale length "L" to be the major
radius) to get frequencies normalized to v_t/R instead of v_t/a_minor (but
then rhoc must be set to r_mid/R, not r_mid/a_minor, tprim becomes R/LT not
a_minor/LT, etc).  

But you must be careful: for irho=2, rho is defined by rho=r_midplane/L, so
quantities involving d/drho will change if the choice of L changes.

There is some freedom in the choice of the reference ion parameters.  If
one uses hydrogen as the reference ion species, then in the
species_parameters namelists below, one would set mass=2 for deuterium,
mass=12 for Carbon, etc., and the reported frequencies and growth rates
will be in units of v_t/a for hydrogen, However, one can use deuterium as
the reference ion species.  Then in the species_parameters namelists below,
one would set mass=1 for deuterium, mass=6 for Carbon (Carbon's Z is still
12, because the reference Z is 1), etc, and the reported frequencies will
be in units of v_t/a for deuterium.

!###########################################################################

&theta_grid_knobs
!equilibrium_option='grid.out' ! Read equilibrium from gridgen file
!equilibrium_option='s-alpha'
 equilibrium_option='eik' ! local or numerical equil. (depending on iflux)
/

&theta_grid_parameters

!###########################################################################
! If equilibrium_option='s-alpha', then use:
!
! eps =0.1666667   ! r/R, used in B = B_0 / (1 + eps * cos(theta))
!                  ! eps=r/R=0 allowed, to turn off trapped particles.
! epsl = 2.0       ! 2.0 * a/R
! pk = 1.333333    ! 2.0 * a/(q*R)
! shat = 1.0       ! r/q * (dq/dr)
! shift = 0.0      ! =alpha in the s-alpha model (comment out Miller shift)
         ! (s-alpha shift has opposite sign from Miller shift defined below) 
!
 ! Enters into the curvature, grad-B and grad_perp**2 operators as follows:
 ! grad-B and curvature: 
 !
 ! epsl*(cos(theta) + (shat*theta-shift*sin(theta))*sin(theta))

!###########################################################################
! If equilibrium_option='eik', then set most of the Miller/local shape 
! parameters here: 
! (also use s_hat_input and beta_prime_input given later).

 rhoc = 0.849   ! r_mid/a, where r_mid is the half-width of the flux surface
 Rmaj = 2.69    ! R at the center of this flux surface (in units of a)
 R_geo= 2.69    ! R  where B_T is specified in (R_geo also in units of a)
                ! I(psi) = R*B_T is constant on a flux surface.

 qinp = 3.03     ! safety factor q
 shift = -0.354  ! Shafranov shift gradient dRmaj/drho (Rmaj normed to a)
 akappa = 1.66   ! Elongation kappa
 akappri = 1.37  ! dkappa/drho  (rho should be half diameter)
 tri = 0.429     ! triangularity tri = asin[(R(max(Z))-R_major)/r_mid]
                 ! Dorland's tri = asin(delta_Miller)
 tripri = 1.61   ! d(tri)/drho

 ! For local/Miller equilibria, the flux surface shape (in physical units)
 ! is given by: 
 ! R(theta) = [Rmaj + rhoc*cos(theta+tri*sin(theta))] *a
 ! Z(theta) = [ akappa*rhoc*sin(theta)              ] *a

!###########################################################################
! Used for both s-alpha and Miller/local equilibrium:

 ntheta= 32   ! # of theta grid points per 2*pi (needs to be even)
 nperiod= 2   ! # of 2*pi cells = 2*nperiod-1

! For numerically generated equilibria, ntheta is set elewhere (in
! gridgen?).  Actually ntheta controls # of theta grid points for calculating
! omega_d(theta), etc.)

! Initialization CPU time scales as N_theta_tot**3=(ntheta*(2*nperiod-1))**3
! After initialization, CPU time / time step scales as 
!       C1*N_theta_tot**2 + C2*N_theta_tot*negrid*ngauss 

/

!#########################################################################
&parameters
 beta = 0.0142   ! Not the total beta, this beta is defined as
                 !    403. * n_e_19 * T_ref(kev) / 1.e5 / B_T0**2
                 ! (this beta parameter needed for electromagnetic runs)
 ! where n_e_19 is the electron density in units of 1.e19/m^3 and
 ! B_T0 is the toroidal field at R_geo in Tesla, such that the local 
 ! toroidal field B_T(theta) on the flux surface is:
 !
 !    B_T(theta) = B_T0*R_geo/R(theta) = I(Psi)/R
 !
 ! since R*B_T is constant on a flux surface.  B_T is often reported in
 ! experiments at the R at the center of the last-closed flux surface
 ! (LCFS), so some times it is convenient to allow R_geo (where B_T0 is
 ! defined) and Rmaj (of this flux surface) to differ.  Dia/paramagnetic
 ! effects can make R_geo*B_T0 on the LCFS differ a bit from the value
 ! on this flux-surface.

 zeff =   1.00000    ! only used for electron collisionality
 ! TiTe = 1.0   ! Ti/Te, used only if electrons are assumed to be adiabatic
 ! adiabatic electrons assumed if there is no species of type='electron'
 ! default value is TiTe=1.0
/

&collisions_knobs
 conserve_momentum=.true.   ! conserves momentum of each ion species
                            ! (electron momentum not conserved).
 ! See collision_notes.ps/pdf for more info.

 collision_model='default'  ! Lorentz collision operator
 !conserve_number=.true.    ! applicable only to Krook operator
 ! use_shmem = .true.  ! is this ever needed ??
/

!#########################################################################
&theta_grid_eik_knobs     ! Only read if equilibrium_option='eik'
 itor = 1   ! =1 only working option of relevance

!iflux = 1  ! =1 to use numerical equilbrium
 iflux = 0  ! =0 to use Miller/local equilibrium.

!irho = 1   ! rho = sqrt(Phi_tor/Phi_tor_edge)
!irho = 3   ! rho = Phi_pol/Phi_pol_edge
 irho = 2   ! =2 for rho to be defined as r_mid/a, where r_mid is the
            ! midplane minor radius, i.e., the half-width of the
            ! flux-surface measured at the height of the magnetic axis. 
            ! irho=2 recommended, required for Miller equilibrium (iflux=0)

 ! Options for various types of equilibrium
 ppl_eq = .false.   ! PPPL equilibrium
 gen_eq =  .false.  ! generic
 vmom_eq = .false.  ! VMOMS
 efit_eq = .false.  ! EFIT
 local_eq = .true.  ! .true. required when iflux=0 for Miller equilibrium
                    ! local_eq and iflux=0 a little redundant

 eqfile = 'g081499.04000.beer'  ! EFIT results file name
                 ! JSOLVER can also make a file in the EQDSK (EFIT) format

 equal_arc = .true.  ! equal poloidal arc lengths, same as NCLASS choice
                     ! =.true. usually best choice

 ! See geometry.f90 for defns of bishop parameter.
 bishop = 0 ! do not use Bishop relations.
 bishop = 1 ! use Bishop relations with actual equilibrium s_hat, p'
 bishop = 4 ! use Bishop relations with s_hat_input and beta_prime_input:

 ! used only for bishop>1, or for local 
 s_hat_input = 2.85   ! magnetic shear s_hat = rho/q*dq/drho
 beta_prime_input = -0.113  ! total dbeta/drho used in equilibrium.  

 ! in principle, beta_prime_input can be determined from the beta parameter
 ! and all of the species parameters entered below.  I.e.
 !
 ! beta_prime_input = beta * sum_s dens * temp * (fprim + tprim))
 !
 ! where sum_s is a sum over all species.  The beta_prime_input entered
 ! separately here is used only in determining the plasma equilibrium shape
 ! and the resulting k_r(theta) etc. in the ballooning coordinates.  The
 ! parameter beta (and the species gradients a/LT and a/Ln) separately
 ! effects the strength of the magnetic fluctuations in the gyrokinetic Eq.

 ! alpha_input = 0  ! not used for bishop=4
 ! invLp_input = 2.7  ! not used for bishop=4

 delrho = 1.e-3  ! step size for radial derivatives. 1.e-5 good except
                 ! (1) very close to edge, or (2) in a local equilibrium
                 ! with rhoc very small to turn off trapped particles.
                 ! need delrho < rhoc and delrho < 1-rhoc.
 isym = 1  ! 0 allows up-down asymetry, 1 to symmetrize up-down equilibrium
 writelots = .false.
/ 
!#########################################################################

&fields_knobs
 field_option='implicit'
/
&gs2_diagnostics_knobs
 print_line=.true.  ! to terminal
 write_line=.true.  ! to file
 write_phi=.true.                                            
 write_apar=.false.                                            
 write_aperp=.false.                                            
 write_omega=.true.  ! complex omega
 write_omavg=.true.  ! time-averaged omega
 write_qheat=.false. ! quasilinear heat flux
 write_pflux=.false. ! quasilinear particle flux
 write_vflux=.false. ! v_par flux ??
 write_qmheat=.false. ! Magnetic flutter part of heat flux
 write_pmflux=.false. ! Magnetic flutter part of particle flux
 write_vmflux=.false. ! Magnetic flutter part of v_par flux
 write_dmix=.false.     ! calculate D_mixing = gamma/<k_perp^2>
 write_kperpnorm=.false.                                            
 write_phitot=.false.                                            
 write_eigenfunc=.true.  ! write normalized eigenfunctions phi(theta)
 write_final_fields=.true.
 write_fieldline_avg_phi=.false.
 write_neoclassical_flux=.false.
 dump_neoclassical_flux=.false.
 dump_check1=.false.   ! .true. to write some zonal flow diagnostics
 write_fieldcheck=.false.
 write_fcheck=.false.

!print_old_units=.t.  ! to get units used in KRT paper (with sqrt(2) factors)?
 print_old_units=.f.  ! to get units in v_t/a (without sqrt(2), etc)

 nwrite=           50  ! # of time steps between writing
 navg=          20     ! time averaging for omega
 omegatol=  1.0e-4     ! run till omega(t) doesn't vary much
 omegatinst = 500.0    ! give up if omega > omegatinst (numerical instability?)
/

!#########################################################################
&le_grids_knobs
 ngauss=10  ! # of passing pitch angles. ngauss=1-6, 8, 10, or 12
            ! # of trapped particle pitch angles = ntheta/2
 negrid=8   ! # of energy grid points.  negrid = 3-12, 14, or 18.

 Ecut=   2.500000    ! negrid-2 energy grid points cover 0<E<Ecut
                     ! 2 more energy grid points exist for E>Ecut

/
&dist_fn_knobs
 gridfac=   1.0    ! no reason to change this.
                   ! fields at theta boundary reduced by 1/gridfac hack
 apfac=   1.000000     ! related to delta B_par ?
 driftknob=   1.0000000    
 !adiabatic_option = "iphi00=2"  ! uncomment for zonal flow test
 !boundary_option ="default"     ! uncomment for zonal flow test
/
&init_g_knobs   ! "g" refers to non-adiabatic part of distribution function 
 !ginit_option= "zero"  ! for zonal flow test
 !ginit_option= "default"
 ginit_option= "noise"

 phiinit=   1.e-5     ! initial amplitude

 ! If ginit_option="default", then the next 3 parameters are used:
 width0=   3.00000    ! initial width in theta of g
 chop_side = .false.  ! initialize even parity only
 chop_side = .true.   ! initialize even and odd parities

 !next 3 are values for zonal flow test:
 !width0=5.e4
 !phiinit=1.0
 !chop_side=.false.

/

!#########################################################################
!
! k_theta (and theta_0) grid selection/layout
!
! aky is normalized to the reference rho_i = v_t_ref/omega_ref
! where v_t_ref=sqrt(T_ref/m_ref)
!
!#########################################################################

&kt_grids_knobs
!grid_option='specified'
!grid_option='single'    ! for zonal flows
 grid_option='range'
!norm_option='default'   ! aky defined as k_theta*rho_i*sqrt(2)  (KRT paper)
 norm_option='t_over_m'  ! aky defined as k_theta*rho_i
/

! Waltz's k_theta and rho_s are different than Dorland's k_theta and rho_i,
! But Waltz's k_theta*rho_s is the same as Dorland's ky=k_theta*rho_i (at
! Ti=Te).

&kt_grids_range_parameters
 naky=1       ! # of k_thetas, spanning aky_min to aky_max:
 ntheta0=1    ! # of theta_0's, spanning theta0_min to theta0_max:
 aky_min=0.3
 aky_max=0.3
 theta0_min=0.
 theta0_max=0.
/
&kt_grids_single_parameters
 theta0=0.
 aky=0.4
 !aky=0.      ! zonal flow test
 !akx=0.013   ! zonal flow test (specify akx since theta0=infinity at ky=0)
/ 
&kt_grids_specified_parameters
 naky=2
 ntheta0=1
/ 
&kt_grids_specified_element_1
 aky=0.2828
 theta0=0.
/
&kt_grids_specified_element_2
 aky=0.5
 theta0=0.
/

&knobs
 fphi=   1.000000     ! multiplier on electrostatic potential phi.
 fapar=   1.000000    ! 1 to keep A_par, 0 for electrostatic runs
 faperp=  1.000000    ! 1 to keep A_perp (delta B_par), 0 for electrostatic

 delt=0.05  ! time_step in code is delt*a/v_t/sqrt(2) in NL norm

!delt=1.06066017177982E-02 ! time_step = delt*a/v_t/aky*sqrt(2) in KRT norm,
      ! using aky in the KRT norm.  In KRT norm, time_step*omega_*a = delt,
      ! where omega_*a=(cT/eB)*k_theta/a is omega_* evaluated at at L_n=a.

     ! Implicit time-centered solutions of a simple equation usually require
     !       dt*Re(omega) < 0.5 and dt*Im(omega) < 0.5 
     ! to determine growth rates to 5% accuracy.

 nstep=        1000  ! # of time steps to take
!wstar_units = .true.   ! normalize omega to omega_star_e
 wstar_units = .false.  ! frequency normalized to v_t/a = sqrt(T_ref/m_ref)/a
/

!#########################################################################
!
! Multiple species definitions
!
!#########################################################################

&species_knobs
 nspec=          2   ! number of species
/
&species_parameters_1                                                 
 z=      1.000000    ! Z_s charge of species s
 mass=   1.000000    ! mass_s relative to reference mass (i.e. H or D)
 dens=   1.00        ! n_s/n_e
 temp=   1.000000    ! T_s/T_ref
 tprim=  3.0         ! a/LT_s
 fprim=  1.0         ! a/L_n_s
 ! More accurately for any choice of irho, tprim = -dlog(T_s)/drho
 !                                         fprim = -dlog(n_s)/drho
 uprim=  0.0000000E+00  ! u_parallel gradient
 vnewk=  0.000          ! collisionality for this species, see below
 type='ion'
/
&dist_fn_species_knobs_1
 fexpr= 0.5000000      ! 0 = fully implicit, 0.5 = standard centered implicit
 fexpi= 0.0000000E+00  ! Imaginary part of implicit parameter, always =0
 bakdif=  0.05         ! upwind parameter (=0 for 2cd order method in KRT)
                       ! sometimes bakdif=0.05 helps smooth the eigenfunction
/

! Note that to preserve quasineutrality radially, should satisfy the
! following sum over species (where z, dens, fprim are functions of the
! species index s, and z=-1 for electrons):
! 
!          sum_s z dens fprim = 0

&species_parameters_2                                              
 z=   -1.000000    
 mass=   2.72e-4
 dens=   1.000000
 temp=   1.0
 tprim=  3.0
 fprim=  1.0
 uprim=   0.00000E+00
 vnewk=   0.0    
 type='electron'
/

 vnewk_e= 0.00279 *coulog *ne_19 /Te(keV)**1.5 *a *(mass_ref/ T_ref(kev))**0.5
 coulog = 15.94 - 0.5 * log(ne_19 / Te(kev)**2)
 see below for more details

&dist_fn_species_knobs_2
 fexpr= 0.5000000
 fexpi= 0.0000000E+00
 bakdif=  0.05
/
&species_parameters_3                                              
 z=      6.000000    
 mass=   6.000000    
 dens=   0.0
 temp=   1.0
 tprim=  2.482
 fprim=  0.795
 uprim=  0.0000000E+00
 vnewk=  0.000
 type='ion'
/
&dist_fn_species_knobs_3
 fexpr= 0.5000000
 fexpi= 0.0000000E+00
 bakdif=  0.05
/
&species_parameters_4                                              
 z=      1.000000    
 mass=   1.000000    
 dens=   0.0
 temp=   32.6  ! If type='beam', then temp=E_inj/T_ref
 tprim=  0.5   ! a/L_"Tbeam"  (? Ecrit or Einj varied radially)
 fprim=  1.872
 uprim=  0.0000000E+00
 vnewk=  0.000
 type='ion'
!type='beam'  ! ='beam' for a rough model of a slowing down distribution
    ! GWH: I recommend treating beams as a hot Maxwellian with type='ion',
    ! see docs/beams.txt.
/
&dist_fn_species_knobs_4
 fexpr= 0.5000000
 fexpi= 0.0000000E+00
 bakdif=  0.05
/

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Collisionality parameters

vnewk for the electrons is

 vnewk_e= 0.00279 *coulog *ne_19 /Te(keV)**1.5 *a *(mass_ref/ T_ref(kev))**0.5

 coulog = 15.94 - 0.5 * log(ne_19 / Te(kev)**2)

where ne_19 is the electron density in 1.e19/m^3, the reference macroscale
a is in meters, the electron temperature Te(keV) and reference temperature
T_ref(keV} are in keV, and mass_ref is the reference mass in AMU (i.e.,
relative to the proton mass).

vnewk for ion species i, approximately accounting for collisions with all
other species b (including b=i), is:

vnewk_i = vnewk_e *sqrt(m_e/m_i) *(T_e/T_i)**1.5 *Z_i**2 *
          sum_b (n_b/n_e) *Z_b**2 *2 /(1+(mass_i/mass_b*T_b/T_i)**0.5)

For more info, see docs/collisions_notes.ps (or *.pdf).

!#########################################################################

&theta_grid_file_knobs
 gridout_file='grid.out'  ! name of the rungridgen output file to read
                          ! to get geometry info from
/
&theta_grid_gridgen_knobs
 npadd = 0
 alknob = 0.0
 epsknob = 1.e-5
 extrknob = 0.0
 tension = 1.0
 thetamax = 0.0
 deltaw = 0.0
 widthw = 1.0
/
&source_knobs
 ! These next 3 options are for the Rosenbluth-Hinton zonal flow tests.
 !
 !source_option="phiext_full"	! turn on external phi
 !phi_ext=-1.0			! magnitude of external phi
 !t0=-10.0                      ! sudden turn on if t0<0.

 ! rarely used:
 !source_option="test2_full"
 !thetas=1000.0
 !omega0=0.0
 !gamma0=0.0
 !gamma_damp=0.0
 !gamma_damp_e=0.0
/
&nonlinear_terms_knobs
/
&additional_linear_terms_knobs
/
&reinit_knobs
/
&theta_grid_salpha_knobs  ! empty namelist needed for s-alpha to work
 !model_option='s-alpha'  ! no longer used?
 !alpmhdfac=   0.000000   ! no longer used?
 !alpha1=  0.0000000E+00  ! no longer used?
/