gvec.scripts.quasr#

The pyGVEC script for converting a QUASR configuration to a G-Frame for use with GVEC.

QUASR is the QUAsi-symmetric Stellarator Repository: https://quasr.flatironinstitute.org/. The algorithm used to construct a G-Frame is described in [1] and is as follows:

Evaluate surface in cartesian space

First use the simsopt interface to load a json file from QUASR, which can also be retrieved automatically with a given ID. Then evaluate the surface cartesian position \((x,y,z)(\vartheta_i,\zeta_j)\) at a meshgrid on the full torus:

\(\vartheta_i=2\pi \frac{i}{n_t},i=0\dots,n_t-1,\quad \zeta_j=2\pi\frac{j}{n_z},j=0,\dots,n_z-1\)

where the angles \(\vartheta,\zeta\) are just the parametrization of the given surface. (For the quasr surfaces, its a boozer angle parameterization!). The number \(n_z\) is chosen as a multiple of the number of field periods \(n_{FP}\), to be able to reduce the discrete dataset exactly to one field period.
Project to surface with elliptical cross-sections

Project \({\bm x}_m(\zeta)=\frac{1}{2\pi}\int_{\vartheta=0}^{2\pi}{\bm x}(\vartheta,\zeta)\sigma_m(\vartheta)d\vartheta\) with \(\sigma_0(\vartheta)=1,\sigma_{s}(\vartheta)=2\sin(\vartheta),\sigma_{c}(\vartheta)=2\cos(\vartheta)\), leading to a surface

\({\bm x}_m(\vartheta,\zeta)={\bm x}_0(\zeta) + {\bm x}_s(\zeta)\sin(\vartheta)+ {\bm x}_c(\zeta)\cos(\vartheta)\)

where cross-sections of \(\zeta=\text{const.}\) are planar elliptical curves.
Compute the plane of the ellipse cross-sections

First choice is to set the first basis unit vector \(N\) from the center point of the ellipse to a point on the boundary at \(\vartheta=0\) position. Then use a first guess for the second basis unit vector \(B\) from the center point to \(\theta=\frac{\pi}{2}\) position to span the unit normal of the plane \(K=(N \times B)\), and then set the second unit vector \(B=K\times N\), such that \(N\) and \(B\) are orthonormal and describe the plane of the ellipse.
compute fourier coefficients of the ellipse

The ellipse in a single \(N,B\) plane is defined as \(X^k(\vartheta)= x^k_c\cos(\vartheta)+x^k_s\sin(\vartheta),k=1,2\). We can deduce from the four coefficients the shift \(\vartheta_0\) and the rotation angle \(\Gamma\).
Final frame

The final frame is obtained by rotating the \(N\) and \(B\) vectors by \(-(\Gamma-\vartheta_0)\), which yields constant rotation speed along \(\zeta\). Thus, in the final frame, the rotating ellipse is represented by a single poloidal and a single toroidal Fourier mode.
cut the original surface with the planes of the frame

For each discrete \(N,B\) plane, we compute the intersection of the all curves \(\bm x(\vartheta_i,\zeta)\) and compute its position \(X^1,X^2\) in the \(N,B\) plane. This gives the final surface.

See also

gvec.gframe: module for constructing and manipulating G-Frames
QUASR and G-Frame: high-level documentation for this script

References

gvec.scripts.quasr.check_args(parser, args)#

gvec.scripts.quasr.generate_quasr_case(ID: int, save_path: str | Path, totalIter: int = 100000, tol: float = 1e-08, time: float = 1200.0, **kwargs)#

Generate a gvec quasr equilibrium and compare it to fieldline tracing.

Parameters:

ID (int) – Quasr ID
save_path (str | Path) – Directory where the case is saved
totalIter (int, optional) – GVEC totalIter, by default 100000
tol (float, optional) – Tolerance for determining minimal necessary (M, N) for the output Fourier modes of X1,X2. default is 1e-8
time (float, optional) – Time to trace fieldlines, by default 1200.0

gvec.scripts.quasr.get_coils_from_json_file(filename: Path | str)#

Get coils from a quasr json file and translate them into a CoilSet

Parameters:: filename (Path | str) – json filename/path
Returns:: Coils as GVEC CoilSet
Return type:: CoilSet

gvec.scripts.quasr.get_json_from_quasr(configuration: int, filename: str | Path | None = None)#: Retrieve a simsopt-compatible JSON for a given QUASR configuration.

gvec.scripts.quasr.get_surface_from_json_file(filename: Path | str)#: Get the boundary surface as a SIMSOPT Surface object from a QUASR JSON file.

gvec.scripts.quasr.get_xyz_from_surface(nt: int, nz: int, surface)#

Sample a SIMSOPT Surface object in cartesian coordinates.

Sample surface at nt,nz*nfp point positions on the full torus. Gives cartesian positions xyz[0:nz*nfp, 0:nt, 0:2].

gvec.scripts.quasr.load_quasr(configuration: int | str | Path, filetype: Literal['netcdf', 'json'] | None = None, tol: float = 1e-08, nt: int = 81, nz: int = 81, quiet: bool = True, verbose: int = 0, name: str | None = None, save_boundary: bool = True, cutoff: int = -1, param_type: Literal['toml', 'yaml'] = 'toml', clean: int = 0, stellsym: bool = False, boundary_coefficients: bool = False, logger: Logger | None = None)#

Load a QUASR configuration and convert it to a G-Frame and boundary for use with GVEC.

Parameters:

configuration (int | str | Path) – QUASR configuration ID or file.
filetype (Literal["netcd","json"] | None, optional) – If configuration is a string or Path, specifies the boundary file format. Either a netCDF file containing boundary data or a SIMSOPT JSON file of the boundary (e.g. QUASR configuration).
tol (float, optional) – _description_, by default 1e-8
nt (int, optional) – _description_, by default 81
nz (int, optional) – _description_, by default 81
quiet (bool, optional) – _description_, by default True
verbose (int, optional) – verbosity level: 1 for info, 2 for debug, by default 0
name (str, optional) – Name for outputfiles. If not set, the name of the input boundary is used. By default None
save_boundary (bool, optional) – Save the boundary points to a netCDF file, by default True
cutoff (int, optional) – Cutoff toroidal mode number only for G-Frame construction, reduces the number of Fourier modes in the G-Frame. Default is -1, which means no cutoff.
param_type (Literal["toml";,"yaml"], optional) – Parameterfile format, by default “toml”
clean (int, optional) – Tolerance to reduce necessary Fourier modes (M, N) for the input surface. Default is 0., which means no cleaning.
stellsym (bool, optional) – If set, imposes stellarator symmetry for the input surface. Use this with great care! By default False
boundary_coefficients (bool, optional) – Write the boundary data as fourier coefficients in the parameter file instead of points in the netCDF file.
logger (Logger | None, optional) – If None a new logger will be created.

Returns:

parameters (dict) – Dictionary containing the GVEC parameters.
dict_out (dict) – Dictionary containing the G-frame data.

gvec.scripts.quasr.load_xyz(filename: Path | str)#

gvec.scripts.quasr.main(args: Sequence[str] | Namespace | None = None)#

Perform fieldline tracing with initial conditions on equilibrium flux surfaces

Parameters:

state (State) – GVEC state used for initialization
coil_set (CoilSet) – Coils for evaluating the magnetic field
n_fieldlines (int, optional) – Number of initial conditions, by default 20
max_step (float, optional) – Maximum stepsize as used by solve_ivp, by default 0.1
t (int, optional) – Time to trace fieldlines, by default 1200
zetas (ArrayLike | None | int, optional) – Poloidal cut positions, by default None
axs (pyplot axes, optional) – Axes to draw the flux-surface contours / Poincaré sections into, default None

Returns:

Returns either the fieldlines and newly generated pyplot axes or only the fieldlines.

Return type:

tuple(dt, axs) or dt

gvec.scripts.quasr.save_xyz(xyz: ndarray, nfp: int, filename: Path | str, attrs: dict = {})#

Save cartesian surface points xyz[0:nz*nfp,0:nt,0:2] to a netcdf file.

The surface cartesian position \((x,y,z)(\vartheta_i,\zeta_j)\) must be evaluated at a meshgrid on the full torus, excluding the periodic endpoint: \(\vartheta_i=2\pi \frac{i}{n_t},i=0\dots,n_t-1,\quad \zeta_j=2\pi\frac{j}{n_z},j=0,\dots,n_z-1\)

Parameters:

xyz (3D np.ndarray) – cartesian surface points of shape (nz*nfp, nt, 3)
nfp (int) – number of field periods
filename (Path | str) – path to the output netcdf file (with .nc extension)
attrs (dict, optional) – additional attributes to add to the netcdf file

Examples

Given a function eval_surface that evaluates the surface cartesian position \(\vec{x}(\vartheta,\zeta)\):

>>> theta = np.linspace(0, 2 * np.pi, ntheta, endpoint=False)
>>> zeta = np.linspace(0, 2 * np.pi, nzeta * nfp, endpoint=False)
>>> xyz = np.zeros((nzeta * nfp, ntheta, 3))
>>> for j in range(nzeta * nfp):
...     for i in range(ntheta):
...         xyz[j,i,:] = eval_surface(theta[i], zeta[j])
>>> save_xyz(xyz, nfp, "example-boundary.nc")