GPEC DIII-D Ideal Example¶
This example uses a baseline DIII-D-like H-mode. The equilibrium is an artificial one made with TokaMaker to have 0D parameters matching typical DIII-D H-mode values, and the kinetic profiles are self-consistent with the equilibrium.
Note that the coil input namelist uses a relative path to the package’s data files. If you would like to run the example in a different location, you must update that path appropriately.
Following a typical workflow we first cd to the example’s directory, then run DCON:
dcon
DCON START => ...
DCON STOP => Normal termination.
From the positive total free-boundary energy, we see that the plasma is stable for n=1 perturbations. Total energies of order 1e0–1e-2 indicate clear stability, while values of order 1e-3 or smaller may warrant further investigation. It is always a good idea to verify the solution using xdraw, OMFIT, or your own plotting routines. One can also skim the human readable dcon.out ascii output.
A stride.in is also included, allowing STRIDE to be run as an alternative stability solver. STRIDE parallelizes the ODE integration across radial intervals and should give consistent n=1 stability results:
stride
STRIDE START => ...
STRIDE STOP => Normal termination.
Note
STRIDE does not write euler.bin (its bin_euler flag is a dummy variable), so DCON must be run when GPEC will be executed in subsequent steps.
STRIDE also writes stride_output_n1.nc, which holds the outer-region tearing stability matrix Delta’ and the rational-surface data that SLAYER needs to compute resistive-layer tearing growth rates. With read_eq=t, the included slayer.in reads that file together with the kinetic profiles in TkMkr_D3Dlike_Hmode.gpeckf. It then matches each layer response to Delta’ on the first three rational surfaces:
slayer
Calculating growth rate on m = 2 rational surface:
Delta_prime (r_s ref.): 5.94 ...
Delta_crit: 33.3
Dispersion roots found: 1, poles: 1
root 1: Q= -1.162293E+00 -4.163214E-02 omega[rad/s]= -1.059284E+04 gamma[1/s]= -3.794248E+02 polished=T
Delta’ is converted from STRIDE’s normalized-flux reference length to the minor-radius reference that the layer quantities use before matching. With the default dc_type="rfitzp" critical Delta, the 2/1, 3/1 and 4/1 tearing modes are all damped here, because Delta_crit exceeds Delta’ on each surface. Setting dc_type="none" drops this stabilization and leaves the 2/1 weakly unstable (gamma of about 80 /s). Each root is reported in its surface’s ExB frame. The results, including every root and pole found in the scan, are written to slayer_output_n1.nc.
Setting coupling_flag=t instead solves the coupled multi-surface dispersion relation for the first msing_max surfaces. Each surface is then Doppler-shifted by its ExB rotation from the kinetic file, so the root is reported in the lab frame. The Doppler shifts spread the branches apart in Re(Q): with the shipped msing_max=3 and scan_width=1.5 the 2/1 branch falls outside the scan window (SLAYER prints a warning), so for the coupled run set msing_max=2 or widen scan_width.
Satisfied with the DCON result, we run GPEC, which computes the ideal plasma response to the external 3D field from the coils defined in coil.in:
gpec
GPEC START => ...
GPEC STOP => Normal termination.
The single mode permeability quantifies how much the plasma amplifies the applied coil field. A value of ~8 means the total plasma-boundary response is 8 times the applied perturbation for the dominant permeability eigenmode. Detailed permeability eigenvalues and eigenvectors are in gpec_response_n1.out. Resonant field amplitudes at each rational surface are summarized on screen and stored in gpec_singfld_n1.out. The Clebsch coordinate displacements needed by PENTRC are written to gpec_xclebsch_n1.out.
Finally, we run PENTRC to compute the neoclassical toroidal viscosity (NTV) torque from the perturbed field:
pentrc
PENTRC START => ...
PENTRC STOP => Normal termination.
The full (fgar) and trapped-only (tgar) torque calculations agree reasonably well, consistent with the expectation that banana-orbit particles dominate the nonambipolar transport. The reduced large-aspect-ratio (rlar) result (enabled with rlar_flag=t; off in this example) typically differs significantly from the general-aspect-ratio calculations for shaped, high-beta H-mode plasmas.
GPEC DIII-D Kinetic Example¶
This example demonstrates the kinetic extension to the ideal MHD workflow using the same artificial DIII-D-like equilibrium as the ideal example. Here DCON is run with the kinetic Euler-Lagrange equation enabled (kin_flag=t), which adds finite-Larmor-radius and wave-particle resonance contributions to the MHD potential energy. The resulting kinetically-modified eigenfunctions are then used by GPEC to compute self-consistent energy and torque profiles, and by PENTRC for the full neoclassical torque calculation.
Note that the pentrc input namelist uses a relative path to the package’s pre-formed fkmnl.dat data file. If you would like to run the example in a different location, you must update the data_dir path in pentrc.in appropriately.
Following a typical workflow we first cd to the example’s directory, then run DCON with the kinetic flag enabled:
dcon
DCON START => ...
DCON STOP => Normal termination.
With kin_flag=t and con_flag=t in dcon.in, the Euler-Lagrange equation incorporates kinetic corrections from the profiles in TkMkr_D3Dlike_Hmode.gpeckf. The con_flag=t setting allows the integration to continue through singular (rational) surfaces, which is required for the kinetic calculation. The resulting euler.bin contains kinetically-modified eigenfunctions. Checking dcon.out or the output NetCDF file is recommended to confirm the solution looks reasonable.
Satisfied with the DCON result, we run GPEC. With dw_flag=t in gpec.in, GPEC outputs self-consistent kinetic energy and torque profiles in addition to the standard perturbed field quantities:
gpec
GPEC START => ...
GPEC STOP => Normal termination.
The kinetic energy and torque profiles are written to gpec_dw_profile_n1.out. The Clebsch coordinate displacements in gpec_xclebsch_n1.out carry the kinetically-modified field structure needed by PENTRC. Resonant surface field amplitudes are in gpec_singfld_n1.out.
Finally, we run PENTRC for the full neoclassical torque calculation:
pentrc
PENTRC START => ...
PENTRC STOP => Normal termination.
The full general-aspect-ratio (fgar) and trapped-only (tgar) torque results can be compared to assess the trapped-particle dominance. The data_dir in pentrc.in points to the precomputed special function matrix used by the reduced large-aspect-ratio (rlar) method, which can be enabled with rlar_flag=t as a cross-check.
GPEC DIII-D Resistive Example¶
This example demonstrates the resistive MHD workflow using the same baseline H-mode DIII-D-like equilibrium as the ideal example. The resistive workflow extends the ideal analysis by computing tearing-stable (resistive) plasma responses near rational surfaces. RDCON generates both the ideal outer-regime solutions and the resistive inner-layer Galerkin solutions; RMATCH then matches them across each singular surface to produce a complete set of resistive eigenfunctions; and finally GPEC computes the resistive plasma response to the applied coil field.
Following a typical workflow we first cd to the example’s directory.
An optional first step is to run DCON (or STRIDE, for which a stride.in is also provided) to quickly assess the ideal MHD stability before proceeding with the full resistive workflow:
dcon
DCON START => ...
DCON STOP => Normal termination.
Note
When running the full resistive workflow, RDCON will overwrite euler.bin with its own ideal outer solutions, so the DCON step is optional and serves primarily as a quick stability check.
Next, run RDCON. With rpec_flag=t and gal_flag=t set in rdcon.in, RDCON performs both the ideal outer-regime integration and the resonant Galerkin inner-layer computation at each rational surface. The outputs gal_solution.dat and delta_gw.dat contain the Galerkin solutions needed by RMATCH, and euler.bin is updated with the ideal outer solutions for GPEC:
rdcon
RDCON START => ...
RDCON STOP => Normal termination.
Then run RMATCH to match the resistive inner-layer solutions to the ideal outer solutions. This produces globalsol.bin, which contains the complete set of resistive plasma eigenfunctions:
rmatch
RMATCH START => ...
RMATCH STOP => Normal termination.
The resistive matching combines the mpert individual resistive solutions into a single binary file. Key parameters in rmatch.in include the resistivity (eta) and mass density (massden) at each rational surface.
Finally, run GPEC with gal_flag=t to use the Galerkin resistive eigenfunctions for the plasma response calculation. The resp_index=1 setting in gpec.in selects the inductance matrix formulation appropriate for resistive solutions:
gpec
GPEC START => ...
GPEC STOP => Normal termination.
The resistive plasma response accounts for field penetration through the tearing layers at each rational surface. Resonant field amplitudes are in gpec_singfld_n1.out and vacuum resonant fields in gpec_vsingfld_n1.out. Together these allow comparison of screened (plasma) and unscreened (vacuum) resonant fields.