Running your first linear simulation (stellarator)
In this tutorial we set up a linear ion-temperature-gradient (ITG) instability calculation using W7-X stellarator geometry and adiabatic electrons.
Disclaimer: GX is optimized for nonlinear calculations, not necessarily linear ones. Other gyrokinetic codes may be better-suited for some linear studies, which could require sharp velocity-space resolution and/or highly-accurate collision operators.
Contents
Setting up the input file
The input file for this case is included in the GX repository in benchmarks/linear/ITG_w7x/itg_w7x_adiabatic_electrons.in.
All GX input files consist of several groups. For more details about input parameters, see The input file.
Dimensions
The [Dimensions] group controls the number of grid-points/spectral basis functions in each dimension, and the number of evolved kinetic species.
[Dimensions]
ntheta = 256 # number of points along field line (theta)
nperiod = 1 # number of 2pi segments along field line is 2*nperiod-1
nky = 28 # number of (de-aliased) fourier modes in y
nkx = 1 # number of (de-aliased) fourier modes in x
nhermite = 16 # number of hermite moments (v_parallel resolution)
nlaguerre = 8 # number of laguerre moments (mu B resolution)
nspecies = 1 # number of evolved kinetic species (adiabatic electrons don't count towards nspecies)
Here, we have set the number of points along the field line per \(2\pi\) segment to be ntheta = 256. Unlike tokamak linear calculations, for stellarator calculations it is recommended to use nperiod = 1, and control the field line length instead using the npol parameter in [Geometry]] below.
The parameters nky and nkx control the number of (de-aliased) Fourier modes in the perpendicular dimensions, with \(x\) the radial coordinate and \(y\) the binormal coordinate. Specifying nky = 28 means we will evolve 28 \(k_y\) Fourier modes in this calculation. Typically, for linear instability calculations the fastest growing modes will have \(k_x = 0\), so we use nkx = 1.
The parameters nhermite and nlaguerre control the velocity-space resolution for the calculation. Since GX uses a spectral velocity space formulation, the \(v_\parallel\) resolution is controlled by the number of Hermite moments (nhermite), and the \(\mu B\) resolution is controlled by the number of Laguerre moments (nlaguerre).
The parameter nspecies controls the number of evolved kinetic species. Here we are using kinetic ions and adiabatic electrons, so we set nspecies = 1 (adiabatic electrons aren’t a kinetic species and don’t count towards nspecies).
Domain
The [Domain] group controls the size of the domain and the boundary conditions.
[Domain]
y0 = 10.0 # controls box length in y (in units of rho_ref) and minimum ky, so that ky_min*rho_ref = 1/y0
boundary = "linked" # use generalized twist-shift boundary conditions along field line
The parameter y0 is related to the box length in the binormal coordinate \(y\) via \(L_y = 2\pi y_0\), where these lengths are normalized to the reference gyroradius. This effectively sets the minimum \(k_y\) in the box, given by \(k_{y\, \mathrm{min}} \rho_\mathrm{ref} = 1/y_0 = 0.05\). Together with nky = 16, this means we will have a range of Fourier modes with \(k_y = [0, 0.05, 0.1, ..., 0.5, 0.55]\) and \(k_x = 0\).
We use twist-shift boundary conditions along the field line by specifying boundary = "linked". Generalizations by Martin et al for low-global-shear devices are included.
Physics
The [Physics] group controls what physics is included in the simulation.
[Physics]
beta = 0.0 # reference normalized pressure, beta = n_ref T_ref / ( B_ref^2 / (8 pi))
nonlinear_mode = false # this is a linear calculation
Since this is an electrostatic calculation with adiabatic electrons we set the plasma reference beta beta = 0.0, which turns off electromagnetic effects. We also make this a linear calculation by setting nonlinear_mode = false.
Time
The [Time] group controls the timestepping scheme and parameters.
[Time]
t_max = 200.0 # end time (in units of L_ref/vt_ref)
scheme = "rk3" # use RK3 timestepping scheme (with adaptive timestepping)
GX uses explicit timestepping methods, and the timestep size dt will be automatically set based on an estimate of the fastest frequencies in the system. Here we specify an end time with t_max, which means the code will take as many timesteps as it needs to reach this time. We will use the 3rd order Runge-Kutta (rk3) timestepper.
Initialization
The [Initialization] group controls the initial conditions.
[Initialization]
ikpar_init = 0 # parallel wavenumber of initial perturbation
init_field = "density" # initial condition set in density
init_amp = 1.0e-5 # amplitude of initial condition
Here we set up an initial condition that is constant along the field line (ikpar_init = 0) with initial perturbation amplitude init_amp = 1.0e-5 in the density moment (init_field = "density"). By default, a random initial condition satisfying these constraints is used.
Geometry
The [Geometry] group controls the simulation geometry.
[Geometry]
geo_option = "vmec" # use VMEC geometry
vmec_file = "wout_w7x.nc" # name of the vmec file
#Field line label theta - iota*phi. alpha = 0.0 usually corresponds to a
#field line on the outboard side
alpha = 0.0
# Number of poloidal turns.
# The field line goes from (-npol*PI, npol*PI]
npol = 6
# Normalized toroidal flux (or s) is how vmec labels surfaces.
# s goes from [0,1] where 0 is the magnetic axis and 1 is the
# last closed flux surface.
torflux = 0.64
Here we generate the geometry using a VMEC equilibrium file from the W7-X experiment. GX uses a local flux-tube geometry, which means we must pick a particular value of \(\alpha = \theta - \iota \phi\) that labels the field line on which our flux tube is centered. In the non-axisymmetric geometry of a stellarator, the geometry will be different for different values of \(\alpha\). We extend the flux tube for npol = 6 poloidal turns, resulting in a poloidal domain that ranges from \([-6\pi, 6\pi]\). Note, however, that this extended grid is rescaled to an internal poloidal grid \([-\pi, \pi]\); this is relevant for some diagnostic outputs. The flux surface where our flux tube is centered is specified by torflux, which is \(s = \Psi/\Psi_{LCFS}\), the normalized toroidal flux.
Species
The [species] group specifies parameters like charge, mass, and gradients of each species.
# it is okay to have extra species data here; only the first nspecies elements of each item are used
[species]
z = [ 1.0, -1.0 ] # charge (normalized to Z_ref)
mass = [ 1.0, 2.7e-4 ] # mass (normalized to m_ref)
dens = [ 1.0, 1.0 ] # density (normalized to dens_ref)
temp = [ 1.0, 1.0 ] # temperature (normalized to T_ref)
tprim = [ 3.0, 0.0 ] # temperature gradient, L_ref/L_T
fprim = [ 1.0, 0.0 ] # density gradient, L_ref/L_n
vnewk = [ 0.0, 0.0 ] # collision frequency
type = [ "ion", "electron" ] # species type
Note that species parameters z, mass, dens, temp are normalized to the corresponding value of the reference species, which can be chosen arbitrarily. Here we only have a single species, so we simply choose the ions as the reference species, meaning \(z_i = Z_i/Z_\mathrm{ref} = 1.0\) etc. Also, note that the species tables can have extra species data; only the first nspecies elements will be used.
Boltzmann
The [Boltzmann] group sets up a Boltzmann species. Here we use Boltzmann (adiabatic) electrons.
[Boltzmann]
add_Boltzmann_species = true # use a Boltzmann species
Boltzmann_type = "electrons" # the Boltzmann species will be electrons
tau_fac = 1.0 # temperature ratio, T_i/T_e
The parameter tau_fac is defined to be \(T_\mathrm{non-adiabatic}/T_\mathrm{adiabatic}\). Thus when using adiabatic electrons, tau_fac = T_i/T_e.
Dissipation
The [Dissipation] group controls numerical dissipation parameters.
[Dissipation]
closure_model = "none" # no closure assumptions (just truncation)
hypercollisions = true # use hypercollision model
nu_hyper_m = 0.1 # coefficient of hermite hypercollisions
p_hyper_m = 6 # power of hermite hypercollisions
nu_hyper_l = 0.1 # coefficient of laguerre hypercollisions
p_hyper_l = 6 # power of laguerre hypercollisions
We do not make any closure assumptions (closure_model = "none"), instead opting for a simple truncation of the moment series. In lieu of closures, hypercollisions can provide the necessary dissipation at small scales in velocity space (large hermite and laguerre index).
.. Here we use a hypercollision operator of the form
Diagnostics
The [Diagnostics] group controls the diagnostic quantities that are computed and written to the NetCDF output file.
[Diagnostics]
nwrite = 1000 # write diagnostics every 1000 timesteps
omega = true # compute and write growth rates and frequencies
# spectra of W = |G_lm|**2
[Wspectra]
species = false
hermite = false
laguerre = false
hermite_laguerre = true # W(l,m) (summed over kx, ky, z)
kx = false
ky = true # W(ky) (summed over kx, z, l, m)
kxky = false
z = false
# spectra of P = ( 1 - Gamma_0(b_s) ) |Phi|**2
[Pspectra]
species = false
kx = false
ky = true # P(ky) (summed over kx, z)
kxky = false
z = true # P(z) (summed over kx, ky)
For linear calculations, we can request to compute and write growth rates and real frequencies by using omega = true. Additionally, various spectra can be computed as specified by the [Wspectra] and [Pspectra] groups.
Running the simulation
To run the simulation from the benchmarks/linear/ITG_w7x directory, we can use
../../../gx itg_w7x_adiabatic_electrons.in
This will generate an output file in NetCDF format called itg_w7x_adiabatic_electrons.nc (in the same directory).
Plotting the results
Growth Rates & Frequencies
We can plot the growth rates and real frequencies using the following python script:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator
from netCDF4 import Dataset
fig, (ax1, ax2) = plt.subplots(2)
# read gx data
stem = "itg_w7x_adiabatic_electrons"
data = Dataset("%s.nc" % stem, mode='r')
t = data.variables['time'][:]
ky = data.variables['ky'][1:]
omegas = data.groups['Special'].variables['omega_v_time'][:,1:,0,0]
gams = data.groups['Special'].variables['omega_v_time'][:,1:,0,1]
omavg = np.mean(omegas[int(len(t)/2):, :], axis=0)
gamavg = np.mean(gams[int(len(t)/2):, :], axis=0)
# plot growth rates and frequencies
ax1.plot(ky, gamavg, 'o', fillstyle='none')
ax2.plot(ky, omavg, 'o', fillstyle='none', label='GX')
ax1.xaxis.set_minor_locator(AutoMinorLocator())
ax2.xaxis.set_minor_locator(AutoMinorLocator())
ax1.yaxis.set_minor_locator(AutoMinorLocator())
ax2.yaxis.set_minor_locator(AutoMinorLocator())
ax1.set_xlim(left=0)
ax2.set_xlim(left=0)
ax1.set_ylabel(r'$\gamma\ a / v_{ti}$')
ax2.set_ylabel(r'$\omega\ a / v_{ti}$')
ax2.set_xlabel(r'$k_y \rho_i$')
ax2.legend()
plt.tight_layout()
plt.show()