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 / !######################################################################### ¶meters 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/ 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 0Ecut / &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 /