Procedures

Add the adiabatic eletron contribution for globally radial simulations. This actually entails solving for the whole ky = 0 slice of phi at once (not really adding!)

add vM . grad y d/dy or vM . grad x d/dx (or equivalents with g) or omega_* * d/dy term to RHS of GK equation

adds the contributions to the GKE RHS that comes from switching from gbar^{n+1} = g^{n+1} + (Ze/T)(vpa/c)*F0 to g^{n+1} = in the time derivative;

Add radial variation of the Jacobian and gyroaveraing in the velocity integration of , needed for radially global simulations

Layouts Parameters Grids Calculations

switch the vpa integration weights to ensure correct integration by parts

compute phase factor needed when running with equilibrium flow shear compute ikyg FFT to get dg/dy in (y,x) space compute ikx zero out the zonal contribution to d/dx if requested if running with equilibrium flow shear, make adjustment to the term multiplying dg/dy FFT to get d/dx in (y,x) space multiply by the geometric factor appearing in the Poisson bracket; i.e., (dx/dpsidy/dalpha)0.5 compute the contribution to the Poisson bracket from dg/dy*d/dx

advance_explicit takes as input the guiding centre distribution function in k-space and updates it to account for all of the terms in the GKE that are advanced explicitly in time

advance_explicit_euler uses forward Euler to advance one time step

advance_expliciit_rk2 uses strong stability-preserving RK2 to advance one time step

strong stability-preserving RK3

standard RK4

This calls the appropriate routines needed to all fields in the main code. This routine calls the appropriate update depending on the effects included in the simulation (e.g. Electrostatic, Full Flux surface effects or Radiatl Variation effects).

This calls the appropriate routines needed to advance phi in the main code when using fluxtube stella, depending on the distribution (i.e. if the information is parallelised over (kx,ky,z) or (vpa,mu) ). Note that Apar and Bpar are only advanced when using EM so these are in fields_electromagnetic.fpp

avoid spatially dependent kperp add in hyper-dissipation of form dg/dt = -D(k/kmax)^4g add in hyper-dissipation of form dg/dt = -D(k/kmax)^4g

computes the fourth derivative of g in vpa and returns this in dgdvpa multiplied by the vpa diffusion coefficient

computes the fourth derivative of g in z and returns this in dgdz multiplied by the z hyper diffusion coefficient

dist_choice indicates whether the non-Boltzmann part of the pdf (h) is evolved in parallel streaming or if the guiding centre distribution (g = ) is evolved

advance_mirror_explicit calculates the contribution to the RHS of the gyrokinetic equation due to the mirror force term; it treats all terms explicitly in time

check estimated cfl_dt to see if the time step size needs to be changed

if flux tube simulation parallel streaming stays in ky,kx,z space with ky,kx,z local if full flux surface (flux annulus), will need to calculate in y space start the timer for the parallel streaming part of the time advance

unless running in multibox mode, no need to worry about mb_communicate calls as the subroutine is immediately exited if not in multibox mode. save value of phi & apar for use in diagnostics (to obtain frequency) reverse the order of operations every time step as part of alternating direction operator splitting this is needed to ensure 2nd order accuracy in time

advance_wdriftx_explicit subroutine calculates and adds the x-component of the magnetic drift term to the RHS of the GK equation

advance_wdrifty_explicit subroutine calculates and adds the y-component of the magnetic drift term to the RHS of the GK equation

start timing the time advance due to the driving gradients

Allocate arrays needed for solving fields for all versions of stella

Allocate arrays needed for solving electromagnetic fields This includes Apar and Bpar

set up offset_apar and offset_bpar consistently so that the array slices below are consistent with the size of the response matrix

Calculate the total particle, momentum and heat fluxes (pflux_vs_s, vflux_vs_s, qflux_vs_s) and the contributions from a given (kx,ky,z) location (pflux_kxkyz, vflux_kxkyz, qflux_kxkyz) inputs are the particle density (dens), parallel flow (upar) and pressure (pres)

For momentum flux species-dependent factor by which velocity moments must be multiplied to get density, pressure, etc. set species-dependent factors needed for density, parallel flow and pressure Already allocated arrays (allocated in time_advance.f90). Set to zero just in case

Use R(ialpha,izeta) and Z(ialpha,izeta), to compute X = R * cos(zeta) Y = R * sin(zeta)

center_zed_segment_complex takes as arguments the vpa index (iv) the z-depenendent conplex function f, and the starting iz index for the array f (llim), and overwrites f with the cell-centered version;

center_zed_segment_real takes as arguments the vpa index (iv) the z-depenendent real function f, and the starting iz index for the array f (llim), and overwrites f with the cell-centered version

check_transforms checks the various physics flag choices to determine if FFTs are needed for the simulation

subroutine takes a set of Fourier coefficients (ft) and returns the minimum number of coeffients that must be retained (idx) to ensure that the relative error in the total spectral energy is below a specified tolerance (tol_floor)

Parameters Arrays TODO-GA: move apar stuff to EM fields Routines for deallocating arrays fields depending on the physics being simulated TODO-GA: REMOVE TODO-GA: move the above deallocations into 'finish_fields_electromagnetic' when EM is decoupled

TODO-GA:

arrays only allocated/used if simulating a full flux surface

Finish a simulation, call the finialisation routines of all modules

the Fourier components of the guiding centre distribution function normalized by the equilibrium Maxwellian is passed in as g, along with the Fourier components of the electrostatic potential, phi. g_to_f calculates the Maxwellian-normalized distribution function f, which is related to g via f = g + (Ze/T)*(_R - phi)

compute phi * maxwellian in real space and transform back to k-space phiy = maxwellian * phi Put this back into k-space for future calculations Note that adjust here is (_R - phi) * maxwellian calculate the normalized f, given (_R - phi) f = g + Z/T * ( - phi) * maxwellian f = g + J0 * phi * maxwellian

adjust bpar part of Zs / Ts

adjust bpar part of Zs / Ts

adjust apar part of Zs / Ts

adjust apar part of Zs / Ts

adjust apar part of

Get_apar solves pre-factor * Apar = beta_ref * sum_s Z_s n_s vth_s int d3v vpa * J0 * pdf for apar, with pdf being either g or gbar (specified by dist input). the input apar is the RHS of the above equation and is overwritten by the true apar the pre-factor depends on whether g or gbar is used (kperp2 in former case, with additional term appearing in latter case)

get_contributions_from_apar takes as input the appropriately averaged parallel component of the vector potential, apar, and returns in rhs the sum of the source terms involving apar that appear on the RHS of the GK equation when g is the pdf

get_contributions_from_bpar takes as input the appropriately averaged electrostatic potential bpar and returns in rhs the sum of the source terms involving bpar that appear on the RHS of the GK equation when g is the pdf

get_contributions_from_fields takes as input the appropriately averaged electrostatic potential phi and magnetic vector potential components apar and returns in rhs the sum of the source terms involving phi and apar that appear on the RHS of the GK equation when g is the pdf

get_contributions_from_pdf takes as an argument the evolved pdf (either guiding centre distribution g= or maxwellian-normlized, non-Boltzmann distribution h/F0=f/F0+(Ze*phi/T)) and the scratch array rhs, and returns the source terms that depend on the pdf in rhs

get_contributions_from_phi takes as input the appropriately averaged electrostatic potential phi and returns in rhs the sum of the source terms involving phi that appear on the RHS of the GK equation when g is the pdf

Compute d/dx in (ky,kx) space

TODO-GA: MOVE

Compute d/dy in (ky,kx) space

Compute d/dy and d/dx in (ky,kx) space where <.> is a gyroaverage d/dy = i * ky * J0 * chi d/dx = i * kx * J0 * chi chi = phi - Z/T * vpa * apar There are different routines depending on the size of the input array TODO-GA: maybe separate for EM and electrostatic Compute d/dy in (ky,kx,z,tube) space

compute dg/dx in k-space accepts g(ky,kx)

compute dg/dx in k-space accepts g(ky,kx,z,tube)

compute dg/dx in k-space accepts g(ky,kx,z,tube,(vpa,mu,spec))

compute dg/dy in k-space accepts g(ky,kx)

compute dg/dy in k-space accepts g(ky,kx,z,tube)

compute dg/dy in k-space accepts g(ky,kx,z,tube,(vpa,mu,spec))

TOGO-GA: check division rather than multiplication -- kept division for now to be consistent with parallel_streaming phase shift

Note that these advance fields routines are only needed when advancing the collision operators This is used in coll_dougherty.f90

This is used in coll_fokkerplanck.f90 Note that is looks identical to the routine above - we don't know why they are separated

Layouts Arrays Parameters Grids Calculations

If we are parallelising over (vpa,mu) then this subroutine is called This is the more common version used compared with parallelising over (kx,ky,z) and is the default for stella. This advances the fields when Electromagnetic effects are included, so we advance , , and .

get_fields_ffs accepts as input the guiding centre distribution function g and calculates/returns the electronstatic potential phi for full_flux_surface simulations

If we are parallelising over (kx,ky,z) then this subroutine is called This is the less common version used compared with parallelising over (vpa, mu). This is NOT the default routine for stella, and the flag must be set to use this routine. Here we calculate: sum_s int dv (J0 * g) and then call get_phi which divides this by the appropriate gamtot factor TODO-GA: remove apar from this and make it only needed for EM stella

If we are parallelising over (vpa,mu) then this subroutine is called This is the more common version used compared with parallelising over (kx,ky,z) and is the default for stella. Here we calculate: sum_s int dv (J0 * g) and then call get_phi which divides this by the appropriate gamtot factor TODO-GA: remove apar from this and make it only needed for EM stella

Layouts Parameters Arrays Grids Calculations

Layouts Parameters Grids Calculations assume there is only a single flux surface being simulated integrate over velocity space and sum over species within each processor as v-space and species possibly spread over processors, wlil need to gather sums from each proceessor and sum them all together below gather sub-sums from each processor and add them together store result in phi, which will be further modified below to account for polarization term no longer need , so deallocate

Calculate = dkx/dky * jtwist Minus its sign gives the direction of the shift in kx for the twist-and-shift BC Note that the default BC are the unconnected BC in the module

b_dot_grad_z is the alpha-dependent b . grad z, and gradpar is the constant-in-alpha part of it. for a z-pinch, b_dot_grad_z is independent of alpha.

get_gke_rhs calculates the RHS of the GK equation. as the response matrix approach requires separate solution of the 'inhomogeneous' GKE, the homogeneous GKE (to obtain the response matrix itself), and the full GKE, which RHS is obtained depends on the input values for 'pdf', 'phi', 'apar', 'aparnew' and 'aparold'

is_end is either nspec for fluxes, or ntubes when calcualting |phi|^2

Get the index of the time dimension in the netCDF file that corresponds to a time no larger than tstart

Multiply input moment (e.g. density, momentum, or energy) but Jacobian (called flxfac) Need to do this because Jacobian has y-dependence then Fourier transform back to (ky,kx)-space pflx_vs_kxkyz is the particle flux before summing over (kx,ky) and integrating over z calculate the volume average of the particle flux note that the factor of 1/B that appears in the Jacobian has already been taken into account in the numerator of the flux surface average

The 'phi' variable passed in is: sum_s int dv (J0 * g) This routine divides by the appropriate gamtot factor depending on if we have kinetic or adiabatic electrons, and also on whether we are using 'g' or 'h' as our distribution function that we are evolving. Note that this routine is only called in the Electrostatic, Fluxtube case.

Arrays Parameters Grids

change from rhs defined on grid with ky >=0 and kx from 0,...,kxmax,-kxmax,...,-dkx to rhs_swap defined on grid with ky = -kymax,...,kymax and kx >= 0 solve sum_s Z_s int d^3v = gam0*phi where sum_s Z_s int d^3v is initially passed in as rhs_swap and then rhs_swap is over-written with the solution to the linear system

Non-perturbative approach to solving quasineutrality for radially global simulations

the following routine gets the correction in phi both from gyroaveraging and quasineutrality

calculate_vmec_geometry returns psitor/psitor_lcfs as rhoc stella uses rhoc = rho = sqrt(psitor/psitor_lcfs) = rhotor

if maxwwellian_normalization = .true., the pdf is normalized by F0 (which is not the case otherwise) unless reading in g from a restart file, normalise g by F0 for a full flux surface simulation

Initialise the distribution function with random noise. This is the default

if simulating a full flux surface, the alpha dependence present in kperp makes gyro-averaging non-local in k-space if simulating a flux tube, a gyro-average is local in k-space

if simulating a full flux surface, the alpha dependence present in kperp makes gyro-averaging non-local in k-space if simulating a flux tube, a gyro-average is local in k-space

if simulating a full flux surface, the alpha dependence present in kperp makes gyro-averaging non-local in k-space if simulating a flux tube, a gyro-average is local in k-space

if simulating a full flux surface, the alpha dependence present in kperp makes gyro-averaging non-local in k-space if simulating a flux tube, a gyro-average is local in k-space

gyro_average_local takes a field at a given ky, kx, z and (vpa, mu, s) value and returns the gyro-average of that field;

solves Delta * phi_hom = -delta_{ky,0} * ne/Te for phi_hom this is the vector describing the response of phi_hom to a unit impulse in phi_fsa it is the sum over ky and integral over kx of this that is needed, and this is stored in adiabatic_response_factor

aj0_alpha will contain J_0 as a function of k_alpha and alpha j0_B will contain J_0Bexp(-v^2) as a function of k_alpha and alpha for each value of alpha, take kperp^2 calculated on domain kx = [-kx_max, kx_max] and ky = [0, ky_max] and use symmetry to obtain kperp^2 on domain kx = [0, kx_max] and ky = [-ky_max, ky_max] this makes later convolutions involving sums over all ky more straightforward calculate the argument of the Bessel function, which depends on both alpha and k_alpha

TODO:GA- add correct CFL condition

allocate and initialise kperp2 and dkperp2dr allocate and initialise vperp2 init_bessel sets up arrays needed for gyro-averaging; for a flux tube simulation, this is j0 and j1; for a flux annulus simulation, gyro-averaging is non-local in ky and so more effort is required

init_dkperp2dr allocates and initialises the dkperp2dr array, needed for radial variation

phase shift due to the twist-and-shift boundary condition Usually set to zero for standard local simulation, but can have an effect for global simulations and simulations with low magnetic shear that use periodic boundary conditions everywhere

Parameters Routined needed to initialise the different field arrays depending on the physics being simulated Allocate arrays such as phi that are needed throughout the simulation

Fill arrays needed for the electromagnetic calculations

Grids

Initialise arrays that are needed during the main time advance loop in the field solve for the flux tube simulations only These are initialised once and then used throughout the rest of the simulation gamtot = 1 - gamma0 = Z^2*n/T * int e^(-v^2) * (1 - J0^2) dv If using adiabatic electrons then this factor is modified and we use gamtot3 which includes the Boltmann response

calculate the Fourier coefficients in y of the Jacobian this is needed in the computation of the flux surface average of phi

calculate and LU factorise the matrix multiplying the electrostatic potential in quasineutrality this involves the factor 1-Gamma_0(kperp(alpha))

determine if iky corresponds to zonal mode

avoid spatially dependent kperp (through the geometric coefficients) still allowed to vary along zed with global magnetic shear useful for full_flux_surface and radial_variation runs

mirror_int_fac = exp(vpa^2 * (mudB/dz)/(mudB/dz + Zedpihnc/dz)) is the integrating factor needed to turn the dg/dvpa part of the GKE advance into an advection equation a, b and c contain the sub-, main- and super-diagonal terms, respectively if running in full-flux-surface mode, solve mirror advance in y-space rather than ky-space due to alpha-dependence of coefficients corresponds to sign of mirror term positive on RHS of equation must treat boundary carefully assumes fully upwinded at outgoing boundary corresponds to sign of mirror term negative on RHS of equation must treat boundary carefully assumes fully upwinded at outgoing boundary time_upwind = 0.0 corresponds to centered in time time_upwind = 1.0 corresponds to fully implicit (upwinded) account for fact that we have expanded d(gnorm)/dvpa, where gnorm = g/exp(-v^s); this gives rise to d(gnormexp(-vpa^2))/dvpa + 2vpagnormexp(-vpa^2) term we solve for gnorm*exp(-vpa^2) and later multiply by exp(vpa^2) to get gnorm multiply by mirror coefficient

init_kperp2 allocates and initialises the kperp2 array

mirror has sign consistent with being on RHS of GKE; it is the factor multiplying dg/dvpa in the mirror term mirror_sign set to +/- 1 depending on the sign of the mirror term. NB: mirror_sign = -1 corresponds to positive advection velocity set up the tridiagonal matrix that must be inverted for the implicit treatment of the mirror operator

allocate arrays and initialize to zero dvpe * vpe = d(2muB0) * B/2B0 first get equally spaced grid in mu with max value mu_max = vperp_max2/(2max(bmag)) want first grid point at dmu/2 to avoid mu=0 special point dmu/2 + (nmu-1)dmu = mu_max so dmu = mu_max/(nmu-1/2) do simplest thing to start leave dmu(nmu) uninitialized. should never be used, so want valgrind or similar to return error if it is

call all the multibox communication subroutines to make sure all the jobs have the appropriate information

stream_sign set to +/- 1 depending on the sign of the parallel streaming term. NB: stream_sign = -1 corresponds to positive advection velocity only need to consider ia=1, iz=0 and is=1 because alpha, z and species dependences do not lead to change in sign of the streaming pre-factor get gradpar centred in zed for negative vpa (affects upwinding) get gradpar centred in zed for positive vpa (affects upwinding)

Initialise stella

read time_advance_knobs namelist from the input file; sets the explicit time advance option, as well as allows for scaling of the x and y components of the magnetic drifts and of the drive term allocate distribution function sized arrays needed, e.g., for Runge-Kutta time advance set up neoclassical corrections to the equilibrium Maxwellian; only calculated/needed when simulating higher order terms in rhostar for intrinsic rotation calculate the term multiplying dg/dvpa in the mirror term and set up either the semi-Lagrange machinery or the tridiagonal matrix to be inverted if solving implicitly calculate the term multiplying dg/dz in the parallel streaming term and set up the tridiagonal matrix to be inverted if solving implicitly allocate and calculate the factors multiplying dg/dx, dg/dy, dphi/dx and dphi/dy in the magnetic drift terms allocate and calculate the factor multiplying dphi/dy in the gradient drive term

vpa is the parallel velocity at grid points wgts_vpa are the integration weights assigned to the parallel velocity grid points this is the Maxwellian in vpa parallel velocity grid goes from -vpa_max to vpa_max, with no point at vpa = 0; the lack of a point at vpa=0 avoids treating the vpa=z=0 phase space location, which is isolated from all other phase space points in the absence of collisions equal grid spacing in vpa

set up the vpa grid points and integration weights set up the mu grid points and integration weights

allocate wdriftx_phi, the factor multiplying dphi/dx in the magnetic drift term allocate wdrifty_phi, the factor multiplying dphi/dy in the magnetic drift term allocate wdriftx_bpar, the factor multiplying dbpar/dx in the magnetic drift term allocate wdrifty_bpar, the factor multiplying dbpar/dy in the magnetic drift term allocate wdriftx_g, the factor multiplying dg/dx in the magnetic drift term allocate wdrifty_g, the factor multiplying dg/dy in the magnetic drift term this is the curvature drift piece of wdrifty with missing factor of vpa vpa factor is missing to avoid singularity when including non-Maxwellian corrections to equilibrium this is the grad-B drift piece of wdrifty if including neoclassical correction to equilibrium Maxwellian, then add in v_E^{nc} . grad y dg/dy coefficient here if maxwwellian_normalization = .true., evolved distribution function is normalised by a Maxwellian otherwise, it is not; a Maxwellian weighting factor must thus be included assign wdrifty_bpar, neoclassical terms not supported if including neoclassical corrections to equilibrium, add in -(Ze/m) * v_curv/vpa . grad y d/dy * dF^{nc}/dvpa term and v_E . grad z dF^{nc}/dz (here get the dphi/dy part of v_E) NB: the below neoclassical correction needs to be divided by an equilibrium Maxwellian if maxwellian_normalization = .true. this is the curvature drift piece of wdriftx with missing factor of vpa vpa factor is missing to avoid singularity when including non-Maxwellian corrections to equilibrium this is the grad-B drift piece of wdriftx if including neoclassical correction to equilibrium Maxwellian, then add in v_E^{nc} . grad x dg/dx coefficient here if maxwellian_normalizatiion = .true., evolved distribution function is normalised by a Maxwellian otherwise, it is not; a Maxwellian weighting factor must thus be included assign wdriftx_bpar, neoclassical terms not supported if including neoclassical corrections to equilibrium, add in (Ze/m) * v_curv/vpa . grad x d/dx * dF^{nc}/dvpa term and v_E . grad z dF^{nc}/dz (here get the dphi/dx part of v_E) and v_E . grad alpha dF^{nc}/dalpha (dphi/dx part of v_E) NB: the below neoclassical correction needs to be divided by an equilibrium Maxwellian if running with maxwellian_normalzation = .true.

NB: for FFS, assume that there is only one flux annulus the inclusion of the Maxwellian term below is due to the fact that g/F is evolved for FFS

use the LU-decomposed response matrix and the contributions from the 'inhomogeneous' fields (phi, apar) to solve for (phi^{n+1}, apar^{n+1})

gvmu is initialised with a Maxwellian weighting for flux tube simulations, with the Maxwellian evaluated at ia = 1 we are undoing that weighting here, so also need to use ia = 1

Parse some basic command line arguments. Currently just 'version' and 'help'.

If not specified in the input file these are the default options that will be set for all parameters under the namelist &debug_flags'. The default here is that all debug flags are set to .false.

Write the time traces by default, while we will not write any higher dimensional data by default to have a small netCDF file.

Read which option to select for the kxky grid layouts

If not specified in the input file these are the default options that will be set for all parameters under the namelist &numerical'.

If not specified in the input file these are the default options that will be set for all parameters under the namelist &numerical'.

If not specified in the input file these are the default options that will be set for all parameters under the namelist &parameters_physics'.

Rescale fields, including the distribution function

Save the input file in the NetCDF file

solve_gke accepts as argument pdf, the guiding centre distribution function in k-space, and returns rhs_ky, the right-hand side of the gyrokinetic equation in k-space; i.e., if dg/dt = r, then rhs_ky = rdt; note that if include_apar = T, then the input pdf is actually gbar = g + (Ze/T)(vpa/c)F0

split n tasks over current communicator. Returns the low and high indices for a given processor. Assumes indices start at 1

Flush netCDF file to disk

input galph array is real and contains values on the padded alpha grid gkalph is output array; it contains the Fourier coefficients of galph for positive ky values only (reality can be used to obtain the negative ky coefs) the highest 1/3 of the ky modes from the FFT have been discarded to avoid de-aliasing

transform routines start here

Write time trace of electromagnetic field A|| to netCDF

Write time trace of electromagnetic field B|| to netCDF

TODO-GA: use g to f routine rather than g to h -- but check first TODO-GA: use g to f routine rather than g to h -- but check first

Note that to obtain , , to write the fluxes to the txt files, we already calculate fluxes(kx,ky,z,s)