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 "a".  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/a, 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 "a" does not need to be the plasma edge minor
radius. You can play games by setting Rmaj=R_geo=1 to get frequencies
normalized to v_t/R instead of v_t/a (but then rhoc must be set to r_mid/R,
not r_mid/a).

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.  In the collisionless limit, either of
these choices are okay.  But at present the collision operator has some
normalizations that actually assume that the reference mass is 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))
! 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 of the 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(Z=0))/r]
                 ! Dorland's tri = asin(delta_Miller)
 tripri = 1.61   ! d(tri)/drho

!###########################################################################
! 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, so 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 collisionality
 ! TiTe = 1.0   ! Ti/Te, used only if there is no species of type='electron'
                ! default value is TiTe=1.0
/

&collisions_knobs
 conserve_momentum=.true.   ! Only for e-e and not e-i presumably?
 collision_model='default'  ! Lorentz collision 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 = sqrt(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'  !not used
 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 assymetry, 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.
 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    
/
&init_g_knobs   ! "g" refers to non-adiabatic part of distribution function 
 ginit_option= "zero"
 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
/

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

&kt_grids_knobs
 grid_option='specified'
 grid_option='single'
 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
/ 
&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 in KRT norm

     ! 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
 uprim=  0.0000000E+00  ! u_parallel gradient
 vnewk=  0.000          ! collisionality for this species
 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.0000E+00   ! upwind parameter (=0 for 2cd order method in KRT)
/
&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    ! 0.00394 * coulog * ne_19 / Te(keV)**1.5 / Ti(kev)**0.5 * a
 ! for light impurities, vnewk = vnewke * 0.291 * (Zeff - 1) * (Te/Ti)**1.5
 ! ne_19 is the electron density in 1.e19/m^3, a is in meters.
 type='electron'
/

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

&dist_fn_species_knobs_2
 fexpr= 0.5000000
 fexpi= 0.0000000E+00
 bakdif=  0.0000E+00
/
&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.0000E+00
/
&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 slowing down distribution
/
&dist_fn_species_knobs_4
 fexpr= 0.5000000
 fexpi= 0.0000000E+00
 bakdif=  0.0000E+00
/

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

&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
/
&nonlinear_terms_knobs
/
&additional_linear_terms_knobs
/
&reinit_knobs
/
&theta_grid_salpha_knobs  ! empty namelist needed for s-alpha to work
/