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

**********************************************************************
test for general geometry from efit, 2/27/01
**********************************************************************

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 a_minor. You can play games by setting Rmaj=R_geo=1 (this
corresponds to setting the macroscopic scale length "a" 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).

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='file' ! Read equilibrium from gridgen file
 !equilibrium_option='s-alpha'
 !equilibrium_option='eik' ! local (miller) or numerical equil. (depending on iflux)
/

&theta_grid_parameters
 !nothing needed because reading from a file
/

!#########################################################################
&parameters
 beta = 0.001072696! 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.62286612 ! 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'
 ! nothing because not needed
/ 
!#########################################################################

&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=.true.                                            
 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), a might be R, 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
 !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.282842712 
 aky_max=0.282842712
 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.00000000  ! T_s/T_ref
 tprim=  8.20185527  ! a/LT_s
 fprim=  1.479384438
 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.0000E+00   ! upwind parameter (=0 for 2cd order method in KRT)
/

! 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.7196653e-4 ! electron mass / D mass (= el mass / ref mass)
 dens=   1.000000
 temp=   1.076150849
 tprim=  6.927957324
 fprim=  1.479384438
 uprim=   0.00000E+00
 vnewk=   0.2324629    
 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.0000E+00
/
&species_parameters_3                                              
 z=      6.000000    
 mass=   6.000000    
 dens=   0.020762204
 temp=   1.0
 tprim=  8.20185527
 fprim=  1.479384438
 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.07617951
 temp=   24.411  ! If type='beam', then temp=E_inj/T_ref
 tprim=  2.206830862 ! a/L_"Tbeam"  (? Ecrit or Einj varied radially)
 fprim=  9.242109312
 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.0000E+00
/

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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?
/