Simulation Class

This class contains data and methods relating to the simulation.

Instance Variables of Simulation Class

class Simulation(config_file)

Attributes

file_ext	(str, ‘.out.nc’) File extension for NetCDF output file.
run_folder	(str, ‘../..’) Path to run folder.
cdf_file	(str, None) Path (relative or absolute) and name of input NetCDF file. If None, the directory is searched for a file ending in ‘.cdf’ and the name is appended to the run_folder path.
g_file	(str, None) Path to the ‘.g’ file. If None, run_folder will be searched and the first returned file will be used.
geometry	(array_like) Array containing entire ‘.g’ file.
input_file	(dict) Dictionary containing all namelist variables from the ‘.inp’ file produced by GS2.
in_field	(str) Name of the field to be read in from NetCDF file.
analysis	(str) Type of analysis to be done. Options are ‘all’, ‘perp’, ‘par’, ‘time’, ‘write_field’, ‘write_field_full’.
out_dir	(str, ‘analysis’) Output directory for analysis.
time_interpolate_bool	(bool, True) Interpolate in time onto a regular grid. Specify as interpolate in configuration file.
time_interp_fac	(int, 1) Sets the time interpolation multiple.
zero_bes_scales_bool	(bool, False) Zero out scales which are larger than the BES. Specify as zero_bes_scales in configuration file.
zero_zf_scales_bool	(bool, False) Zero out the zonal flow (ky = 0) modes. Specify as zero_zf_scales in configuration file.
lab_frame	(bool, False) Transform from rotating to lab frame.
domain	(str, ‘full’) Specifies whether to analyze the full real space domain, or only the middle part of size box_size.
time_slice	(int, 49) Size of time window for averaging
perp_guess_x	(array_like, 0.05) Initial guess for radial correlation length in metres.
perp_guess_y	(array_like, 0.1) Initial guess for poloidal correlation length in metres.
perp_guess_ky	(array_like, 1) Initial guess for poloidal wavenumber in metres^-1.
ky_free	(bool, False) Determines whether ky is free during the poloidal fitting procedure.
time_guess	(array_like, [1e-5,100]) Initial guess for the correlation time and wavenumber in seconds read in from the configuration file.
time_guess_dec	(float) Guess for the time correlation estimated by the decaying exponential.
time_guess_grow	(float) Guess for the time correlation estimated by the growing exponential.
time_guess_osc	(float) Guess for the time correlation and wavenumber estimated by the oscillating exponential.
time_max	(float, 1) Maximum correlation time in seconds. Values above time_max will be excluded.
box_size	(array_like, [0.2,0.2]) When running correlation analysis in the middle of the full GS2 domain, this sets the approximate [radial, poloidal] size of this box in m. This variable is only used when domain = ‘middle’
time_range	(array_like, [0,-1]) Time index range for which analysis is done. Default is entire range. -1 for the final time step is interpreted as up to the final time step, inclusively.
npeaks_fit	(int) Number of peaks to fit when calculating the correlation time.
species_index	(int) Specied index to be read from NetCDF file. GS2 convention is to use 0 for ion and 1 for electron in a two species simulation.
theta_index	(int or None) Parallel index at which to do analysis. If no theta index in array set to None.
amin	(float) Minor radius of device in m.
vth	(float) Thermal velocity of the reference species in m/s
bref	(float) Reference magnetic field at the centre of the LCFS.
rho_ref	(float) Larmor radius of the reference species in m.
rho_star	(float) The expansion parameter defined as rho_ref/amin.
rho_tor	(float) Radial location of the flux tube in terms of the square root of the normalized toroidal magnetic flux.
pitch_angle	(float) Pitch angle of the magnetic field lines in rad.
rmaj	(float, 0) Major radius of the outboard midplane. Used when writing out the field to the NetCDF file. This is not the rmaj value from GS2.
nref	(float, 1) Density of the reference species in m^-3.
tref	(float, 1) Temperature of the reference species in eV.
omega	(float, 0) Angular frequency of the plasma at the radial location of the flux tube.
dpsi_da	(float, 0) Relationship between psi_n and rho_miller (a_n = diameter/diameter LCFS)
drho_dpsi: float, 1	Gradient of flux surface label with respect to psi.
gradpar	(array_like) Value of the parallel gradient as a function of theta.
r_prime	(array_like) Value of Rprime as a function of theta.
seaborn_context	(str) Context for plot output: paper, notebook, talk, poster. See: http://stanford.edu/~mwaskom/software/seaborn/tutorial/aesthetics.html
field	(array_like) Field read in from the NetCDF file. Automatically converted to a complex array.
field_real_space	(array_like) Field in real space coordinates (x,y)
field_real_space_norm_x	(array_like) Field in real space coordinates (x,y) normalized in the x dimension.
field_real_space_norm_y	(array_like) Field in real space coordinates (x,y) normalized in the y dimension.
perp_corr_x	(array_like) Radial correlation function calculated from field_real_space_norm_x.
perp_corr_y	(array_like) Poloidal correlation function calculated from field_real_space_norm_y.
perp_fit_len_x	(array_like) Radial correlation length obtained from perp fitting procedure. Of size (nt_slices).
perp_fit_len_err_x	(array_like) Error in the radial correlation length obtained from perp fitting procedure, calculated from the covariance matrix. Of size (nt_slices).
perp_fit_len_y	(array_like) Poloidal correlation length obtained from perp fitting procedure. Of size (nt_slices).
perp_fit_len_err_y	(array_like) Error in the poloidal correlation length obtained from perp fitting procedure, calculated from the covariance matrix. Of size (nt_slices).
time_corr	(array_like) Correlation function used to calculate the correlation function. It is of size (nt_slices, time_slice, nx, ny), owing to the 2D correlation calculated in the t and y directions.
corr_time	(array_like) Parameters obtained from time fitting procedure. Of size (nt_slices, nx).
kx	(array_like) Values of the kx grid in the following order: 0,...,kx_max,-kx_max,... kx_min.
ky	(array_like) Values of the ky grid.
t	(array_like) Values of the time grid.
dt	(array_like) Values of the time separations. Min and max values will depend on time_slice.
nt_slices	(int) Number of time slices. nt_slices = nt/time_slice
x	(array_like) Values of the real space x (radial) grid.
x_box_size	(float) Real space size of the box in the x-direction.
dx	(array_like) Values of the dx (radial separation) grid.
fit_dx	(array_like) dx values for section that will be fitted
fit_dx_mesh	(array_like) dx grid (in fitting region) as a 2D mesh.
y	(array_like) Values of the real space y (poloidal) grid. This has been transformed from the toroidal plane to the poloidal plane by using the pitch-angle of the magnetic field lines.
dy	(array_like) Values of the dy (poloidal separation) grid.
fit_dy	(array_like) dy values for section that will be fitted
fit_dy_mesh	(array_like) dy grid (in fitting region) as a 2D mesh.
nkx	(int) Number of kx values.
nky	(int) Number of ky values.
nx	(int) Number of real space x points.
ny	(int) Number of real space y points. This is ny = 2*(nky - 1).
nt	(int) Number of time points.
ncfile	(object) SciPy NetCDF object referencing all arrays in the NetCDF file. Arrays are not loaded into memory due to the large amount of memory needed, but simply read from the NetCDF file. Since the field is manipulated, it is copied into a NumPy array, however.
field_max	(float) Maximum value of the real space field.
field_min	(float) Minimum value of the real space field.
write_field_interp_x	(bool, True) Determines whether the radial coordinate is interpolated to match the BES resolution when writing to a NetCDF file.
r	(array_like) Radial coordinate x, centered at the major radius rmaj.
z	(array_like) Poloidal coordinate z centered at 0.
R	(array_like) Radial location of centre of the flux tube as read in from geometry file.
Z	(array_like) Poloidal location of centre of the flux tube as read in from geometry file.
phi_tor	(array_like) Toroidal angle of location of centre of the flux tube as read in from geometry file. It is equal to -alpha4 in geometry file.
dR_drho	(array_like) Derivative of R with respect to the radial coordinate rho.
dZ_drho	(array_like) Derivative of Z with respect to the radial coordinate rho.
dalpha_drho	(array_like) Derivative of alpha with respect to the radial coordinate rho.
bpol	(array_like) Poloidal magnetic field as a function of theta as printed out in the geometry file.
btor	(float) Toroidal magnetic field at the radial location of the flux tube. btor = [r_geo/R(theta=0)]*bref
bmag	(float) Magnitude of the magnetic field at the location of the flux tube.
r_geo	(float) Radial location of the reference magnetic field.
par_corr	(array_like) Parallel correlation function C(t,x,y,dtheta)
l_par	(array_like) Regular real space parallel grid.
dl_par	(array_like) Values of parallel separation in real space.
par_guess	(array_like, [1,0.1]) Initial guess for the parallel fitting in SI units in the form [l_par, k_par].
par_fit_params	(array_like) Array which stores parallel correlation fitting parameters. It is a function of time slice and contains both l and k which define an oscillating Gaussian.
par_fit_params_err	(array_like) Stores errors associated with fitting the parallel correlation function. These are calculated by taking the square root of the diagonal terms in the covariance matrix and give the error in each fitting parameter.

Configuration Variables

read_config()

The following parameters can be set in the configuration file ‘<name>.ini’.

Parameters:

run_folder : str, ‘../..’

Path to run folder.

cdf_file : str, None

Path (relative or absolute) and name of input NetCDF file. If None, the directory is searched for a file ending in ‘.cdf’ and the name is appended to the run_folder path.

g_file : str, None

Path to the ‘.g’ file. If None, run_folder will be searched and the first returned file will be used.

field : str

Name of the field to be read in from NetCDF file.

analysis : str

Type of analysis to be done. Options are ‘all’, ‘perp’, ‘par’, ‘time’, ‘write_field’, ‘write_field_full’.

out_dir : str, ‘analysis’

Output directory for analysis.

time_interpolate : bool, True

Interpolate in time onto a regular grid.

time_interp_fac : int, 1

Sets the time interpolation multiple.

zero_bes_scales : bool, False

Zero out scales which are larger than the BES.

zero_zf_scales : bool, False

Zero out the zonal flow (ky = 0) modes.

lab_frame : bool, False

Transform from rotating to lab frame.

domain : str, ‘full’

Specifies whether to analyze the full real space domain, or only the middle part of size box_size.

time_slice : int, 49

Size of time window for averaging

perp_guess : array_like, [0.02,0.1,1]

Initial guess for the radial and poloidal correlation lengths. Of the form [lx, ly] in metres. The third parameter is optional and only used when the flag ky_free = True.

ky_free : bool, False

Determines whether ky is free during the poloidal fitting procedure.

time_guess : array_like, [1e-5,100]

Initial guess for the correlation time and wavenumber in seconds read in from the configuration file.

time_max : float, 1

Maximum correlation time in seconds. Values above time_max will be excluded.

box_size : array_like, [0.2,0.2]

When running correlation analysis in the middle of the full GS2 domain, this sets the approximate [radial, poloidal] size of this box in m. This variable is only used when domain = ‘middle’

time_range : array_like, [0,-1]

Time index range for which analysis is done. Default is entire range. -1 for the final time step is interpreted as up to the final time step, inclusively.

npeaks_fit : int, 5

Number of peaks to fit when calculating the correlation time.

species_index : int or None

Specied index to be read from NetCDF file. GS2 convention is to use 0 for ion and 1 for electron in a two species simulation. If reading phi or simulation with only one species, set to None.

theta_index : int or None

Parallel index at which to do analysis. If no theta index in array set to None.

geom_file : str

Location of the geometry file. By default searches the run folder for a ‘.g’ file and loads the first one found.

par_guess : array_like, [1,0.1]

Initial guess for the parallel fitting in SI units in the form [l_par, k_par].

amin : float

Minor radius of device in m.

vth : float

Thermal velocity of the reference species in m/s

bref : float

Reference magnetic field at the centre of the LCFS.

rho_ref : float

Larmor radius of the reference species in m.

rho_tor : float

Radial location of the flux tube in terms of the square root of the normalized toroidal magnetic flux.

nref : float, 1

Density of the reference species in m^-3.

tref : float, 1

Temperature of the reference species in eV.

omega : float, 0

Angular frequency of the plasma at the radial location of the flux tube.

dpsi_da : float, 0

Relationship between psi_n and rho_miller (a_n = diameter/diameter LCFS)

seaborn_context : str, ‘talk’

Context for plot output: paper, notebook, talk, poster. See: http://stanford.edu/~mwaskom/software/seaborn/tutorial/aesthetics.html

write_field_interp_x : bool, True

Determines whether the radial coordinate is interpolated to match the BES resolution when writing to a NetCDF file.

Method Documentation

class gs2_correlation.simulation.Simulation(config_file)

Class containing all simulation information.

The class mainly reads from the simulation NetCDF file and operates on the field specified in the configuration file, such as performing correlations, FFTs, plotting, etc.

An exhaustive list of instance variables is provided separately in the documentation since it is very long.

Methods

calculate_l_par
calculate_par_corr
calculate_perp_corr
calculate_time_corr
config_checks
domain_reduce
extract_input_file
field_normalize_perp
field_normalize_time
field_odd_pts
field_to_complex
field_to_real_space
find_file_with_ext
fourier_correction
par_analysis
par_analysis_summary
par_corr_fit
par_plot
perp_analysis
perp_analysis_summary
perp_corr_fit
perp_norm_mask
perp_plots_x
perp_plots_y
read_config
read_geometry_file
read_input_file
read_netcdf
time_analysis
time_analysis_summary
time_corr_fit
time_interpolate
time_norm_mask
time_plot
to_lab_frame
write_field
write_field_full
write_results
zero_bes_scales
zero_zf_scales

__init__(config_file)

Initialized by object using information from configuration file.

Initialization does the following, depending on the parameters in the configuration file:

• Calculates sizes of kx, ky, x, and y.
• Reads data from the NetCDF file.
• Interpolates onto a regular time grid.
• Zeros out BES scales.
• Zeros out ZF scales.
• Transforms to lab frame.
• Calculates the real space field.
• Reduce domain size to according to box_size if domain is specified as ‘middle’, otherwise do nothing.
• Ensures real space field has odd points.

Parameters:

config_file : str

Filename of configuration file and path if not in the same directory.

Notes

The configuration file should contain the following namelists:

• ‘general’: information such as the analysis to be performed and where to find the NetCDF file.
• ‘perp’ : parameters relating to the perpendicular correlation analysis.
• ‘time’ : parameters relating to the time analysis.
• ‘par’ : parameters relating to the parallel analysis.
• ‘normalization’: normalization parameters for the simulation/experiment.
• ‘output’: parameters which relate to the output produced by the code.

The exact parameters read in are documented in the documentation.

calculate_l_par()

Calculates the real space parallel grid.

calculate_par_corr()

Calculate the parallel correlation function and apply normalization mask.

Interpolation onto a regular parallel grid, correlation calculation, and normalization are done in one function to avoid unnecessary looping over x, y, and theta.

calculate_perp_corr()

Calculates the perpendicular correlation function from the real space field.

calculate_time_corr(it)

Calculate the time correlation for a given time window at each x.

Parameters:

it : int

This is the index of the time slice currently being calculated.

config_checks()

This function contains consistency checks of configurations parameters.

domain_reduce()

Initialization consists of:

• Calculating radial and poloidal coordinates r, z.
• Using input parameter box_size to determine the index range to perform the correlation analysis on.
• Reduce extent of real space field using this index.
• Recalculate some real space arrays such as x, y, dx, dy, etc.

extract_input_file()

Extract input file from cdf_file.

field_normalize_perp()

Defines normalized field for the perpandicular correlation by subtracting the mean and dividing by the RMS value.

field_normalize_time()

Defines normalized field for the time correlation by subtracting the mean and dividing by the RMS value.

field_odd_pts()

Ensures real space field has odd number of points in x and y.

This is done so that when sig.fftconvolve is called on the field with the ‘same’ option, the resulting correlation function also has an odd number of points. If this wasn’t done, one axis might have an even number of points meaning the correlation function will not have a convenient middle point to use as the zero separation point, which would then require more general fitting functions. This way we can force dx=0, dy=0 at the self.nx/2, self.ny/2 location and avoid extra parameters in the fitting functions.

field_to_complex()

Converts field to a complex array.

Field is in the following format: field[t, kx, ky, ri] where ri represents a dimension of length 2.

• ri = 0 - Real part of the field.
• ri = 1 - Imaginary part of the field.

field_to_real_space()

Converts field from (kx, ky) to (x, y) and saves as new array attribute.

Notes

• Since python defines x = IFFT[FFT(x)] need to undo the implicit normalization by multiplying by the size of the arrays.
• GS2 fluctuations are O(rho_star) and must be multiplied by rho_star to get their true values.
• In order to avoid memory overloads, the fourier space field is cleared after the real space field is calculated.

find_file_with_ext(ext)

Find a file in the run_folder with the extension ext

Parameters:

ext : str

Extension of the file to be searched for. Of the form ‘.ext’.

fourier_correction()

Correction to GS2s Fourier components to regular Fourier components.

Notes

GS2 fourier components are NOT simple fourier components obtained from regular FFT functions. Due to GS2’s legacy as a linear code, instead of just getting fourier components g_k the output is:

G_k = {g_k for ky = 0, 2g_k for ky > 0}

Therfore converting to regular fourier components simply means dividing all non-zonal components by 2.

par_analysis()

Calculates the parallel correlation function and fits with a Gaussian to find the parallel correlation length.

par_analysis_summary()

Summarize parallel correlation analysis by plotting parallel correlation fitting parameters along with associated errors and writing to results file.

par_corr_fit(it)

Fit the parallel correlation function with an oscillatory Gaussian function for time slice it.

Parameters:

it : int

Time slice to average over and fit.

Before fitting average over time, x, and y.

par_plot(it, corr, corr_std)

Plots and saves the parallel correlation function and its fit for each time window.

Parameters:

it : int

Time window being plotted.

corr : array_like

Parallel correlation averaged in t, x, and y.

corr_std : array_like

Standard deviation in the correlation function as a function of dl_par.

perp_analysis()

Performs a perpendicular correlation analysis on the field.

Notes

• Remove theta dimension before doing analysis to avoid pointless code which accounts for a theta information. Will produce error if more than one element during config_checks.
• Uses sig.fftconvolve to calculate the 2D perp correlation function.
• Splits correlation function into time slices and fits each time slice with a tilted Gaussian using the perp_fit function.
• The fit parameters for the previous time slice is used as the initial guess for the next time slice.
• Also writes information on the mean fluctuation levels

perp_analysis_summary()

Prints out a summary of the perpendicular analysis.

• Plots fitting parameters as a function of time window.
• Averages them in time and calculates a standard deviation.
• Writes summary to a text file.

perp_corr_fit(it)

Fits the appropriate Gaussian to the radial and poloidal correlation functions.

Parameters:

it : int

This is the index of the time slice currently being fitted.

Notes

• The radial correlation function is fitted with a Gaussian.
• The poloidal correlation function is fitted with an oscillating Gaussian.

perp_norm_mask()

Applies the appropriate normalization to the perpendicular correlation function.

Notes

After calling sig.correlate to calculate perp_corr, we are left with an unnormalized correlation function as a function of dx or dy. This function applies a nomalization mask to perp_corr_x and perp_corr_y which is dependent on the number of points that field_real_space_norm has in common with itself for a given dx or dy. field_real_space_norm is already normalized to the standard deviation of the time signal, so the only difference between correlate and np.corrcoef is the number of points in common in the convolution (that aren’t the zero padded values and after averaging over many time steps).

perp_plots_x(it, corr_fn, corr_std, corr_fit)

Plot radial correlation function and fitted Gaussian.

perp_plots_y(it, corr_fn, corr_std, corr_fit)

Plot radial correlation function and fitted Gaussian.

read_config()

Reads analysis and normalization parameters from self.config_file.

The full list of possible configuration parameters is provided separately in the documentation since it is quite long and are technically not parameters of this function.

read_geometry_file()

Read the geometry file for the GS2 run.

read_input_file()

Read the input ‘.inp’ file for the GS2 run.

read_netcdf()

Read array from NetCDF file.

Read array specified in configuration file as ‘in_field’. Function uses information from the configuration object passed to it.

Notes

• An additional axis is added in the case of no theta_idx. This allows usual loops in theta to be kept but have no effect. If only one element in dimension initialization will be performed regardless, however dimension will be removed for perp and time analysis.

time_analysis()

Performs a time correlation analysis on the field.

Notes

• Split into time windows and perform correlation analysis on each window separately.

time_analysis_summary()

Prints out a summary of the time analysis.

• Plots average correlation time as a function of radius.
• Calculates standard deviation.
• Writes summary to a text file.

time_corr_fit(it)

Fit the time correlation function with exponentials to find correlation time.

Parameters:

it : int

This is the index of the time slice currently being fitted.

Notes

The fitting procedure consists of the following steps:

• Loop over radial points
• Identify the time correlation function peaks
• Determine what type of function to fit to the peaks
• Fitting the appropriate function

The peaks can be fitted with a growing/decaying exponential depending on the direction of decrease of the peaks, or an oscillating exponential if the peaks don’t monotically increase or decrease.

time_interpolate()

Interpolates in time onto a regular grid

Depending on whether the user specified time_interpolate, the time grid is interpolated into a regular grid. This is required in order to do FFTs in time. Interpolation is done by default if not specified. time_interp_fac sets the multiple of interpolation.

time_norm_mask(it)

Applies the appropriate normalization to the time correlation function.

Parameters:

it : int

This is the index of the time slice currently being calculated.

Notes

After calling sig.fftconvolve to calculate time_corr, we are left with an unnormalized correlation function as a function of dt and dy. This function applies a 2D nomalization mask to time_corr which is dependent on the number of points that field_real_space_norm has in common with itself for a given dt, dy, and time window. field_real_space_norm is already normalized to the standard deviation of the time signal, so the only difference between fftconvolve and np.corrcoef is the number of points in common in the convolution (that aren’t the zero padded values and after averaging over many time steps).

time_plot(it, ix, max_index, peaks, plot_type, **kwargs)

Plots the time correlation peaks as well as the apprpriate fitting function.

Parameters:

it : int

Time slice currently being fitted.

ix : int

Radial index currently being fitted.

max_index : array_like

Array containing radial indices of the peaks. Size: (nx, npeaks_fit)

peaks : array_like

Array of peak values at the max_index values. Size: (nx, npeaks_fit)

plot_type : str

Type of fitting function to plot. One of: ‘growing’/’decaying’/ ‘oscillating’

to_lab_frame()

Transforms from the rotating frame to the lab frame.

Notes

The important thing here is that ky is NOT the toroidal wavenumber n0. It is related to n0 by [R1]:

ky_gs2 = n0*(rho_ref/a_min)*dpsi_da

[R1]	C. M. Roach, “Equilibrium flow shear implementation in GS2”, http://gyrokinetics.sourceforge.net/wiki/index.php/Documents, http://svn.code.sf.net/p/gyrokinetics/code/wikifiles/CMR/ExB_GS2.pdf

write_field()

Outputs the field to NetCDF in real space.

Notes

• The radial and poloidal coordinates are centered at 0.
• The radial coordinate is interpolated if neccessary to ensure a 0.5cm resolution consistent with the BES.

write_field_full()

Outputs the full 3D field to NetCDF in real space.

Notes

• The radial, poloidal, and parallel coordinates are centered at 0.
• The radial coordinate is interpolated if neccessary to ensure a 0.5cm resolution consistent with the BES.

write_results(analysis, result_dict)

Write results to the main results file.

Parameters:

result_dict : dict

Dictionary containing results from a given analysis. Will overwrite any existing results for that analysis.

zero_bes_scales()

Sets modes larger than the BES to zero.

The BES is approximately 160x80mm(rad x pol), so we would set kx < 0.25 and ky < 0.5 to zero, since k = 2 pi / L.

zero_zf_scales()

Sets zonal flow (ky = 0) modes to zero.