Skip to content

Replace existing .geqdsk and .gpeckf files with TokaMaker DIII-D-like synthetic equilibrium #245

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
d-burg wants to merge 16 commits into
base: develop
Choose a base branch
from
Open

Replace existing .geqdsk and .gpeckf files with TokaMaker DIII-D-like synthetic equilibrium #245

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
253a6e5
WIP: DCON DIII-D-like example
d-burg Feb 21, 2026
92e82be
WIP: DCON DIII-D-like example - add output figures
d-burg Feb 21, 2026
fe963a1
replace all geqdsk and gpeckf files, update equil and pentrc namelists
d-burg Feb 27, 2026
b499c02
DOCS - MINOR - Makes the DIIID resistive example have it's own self c…
logan-nc Feb 27, 2026
934b97b
DOCS - UPDATE - Adds readme files for all DIIID-like examples
logan-nc Feb 27, 2026
048222b
DOCS - MINOR - Updates etol for the new DIIID-like case in a run_* ex…
logan-nc Feb 27, 2026
ec2a172
rename DIIID-like equilibrium
d-burg Feb 28, 2026
7e0f0e0
Merge remote-tracking branch 'origin/develop' into D3D-like_examples
d-burg Mar 31, 2026
370cae9
DOCS - UPDATE - add updated geqdsks and fix stale READMEs
d-burg Mar 31, 2026
078d239
DOCS - MINOR - Improves DIIID ideal example README
logan-nc Apr 1, 2026
af21aaa
fix example geqdsk header spacing
d-burg Apr 1, 2026
687ffd8
Merge branch 'D3D-like_examples' of github.com:PrincetonUniversity/GP…
d-burg Apr 1, 2026
294d8ad
Update D3D-like equilibrium to Ip=1.15 MA + rotation-corrected kineti…
d-burg Jun 22, 2026
869ec87
Update D3D-like example equilibrium to edge-current-matched variant
d-burg Jun 23, 2026
678d977
Update D3D-like kinetic profiles: analytic ExB rotation
d-burg Jun 24, 2026
53d0488
DOCS - MINOR - Fix plotting typo in DIIID ideal example README
d-burg Jun 25, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
241 changes: 28 additions & 213 deletions docs/examples/DIIID_ideal_example/README
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,227 +1,42 @@
GPEC DIII-D Ideal Example
==========================

This example uses a baseline H-mode Error Field Correction (EFC) plasma from a n=1 EFC study done in 2013. The equilibrium used is a kinetic EFIT from DIII-D shot 147131 at 2300ms, and the profiles are self consistent with the equilibrium. Pre-GPEC results using a combination of GPEC and PENT for this calculation are published in `[C. Paz-Soldan, M. J. Lanctot, N. C. Logan, et al., Phys. Plasmas 21, 072503 (2014).] <http://scitation.aip.org/content/aip/journal/pop/21/7/10.1063/1.4886795>`_.
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 and pentrc input namelists use relative paths to the package's data files. If you would like to run the example in a different location, you must change the input appropriately.
Note that the coil and pentrc input namelists use relative paths to the package's data files. If you would like to run the example in a different location, you must update those paths appropriately.

Following a typical workflow we first cd to the example's directory, then run DCON::
Following a typical workflow we first ``cd`` to the example's directory, then run DCON::

dcon

DCON START => GPEC version 1.0.0
__________________________________________
psihigh = 9.930E-01
Equilibrium: g147131.02300_DIIID_KEFIT, Type: efit
Jac_type = hamada, power_bp = 0, power_b = 0, power_r = 0
Diagnosing Grad-Shafranov solution
Evaluating Mercier criterion
Evaluating ballooning criterion
Fourier analysis of metric tensor components
q0 = 1.062E+00, qmin = 1.062E+00, qmax = 5.321E+00, q95 = 4.228E+00
sas_flag = T, dmlim = 2.000E-01, qlim = 5.200E+00, psilim = 9.915E-01
betat = 2.184E-02, betan = 1.905E+00, betap1 = 7.977E-01
nn = 1, mlow = -12, mhigh = 21, mpert = 34, mband = 33
Computing F, G, and K Matrices
Starting integration of ODE's
psi = 1.000E-04, q = 1.062
psi = 5.936E-01, q = 2.000
psi = 8.186E-01, q = 3.000
psi = 9.282E-01, q = 4.000
psi = 9.884E-01, q = 5.000
psi = 9.915E-01, q = 5.200
Computing free boundary energies
Energies: plasma = -1.754E+00, vacuum = 2.485E+00, total = 7.309E-01
All free-boundary modes stable for nn = 1.
Total cpu time = 8.328E+00 seconds
DCON STOP => Normal termination.


From the positive total energy, we see that the plasma is reasonably stable for n=1 perturbations. Totals are often of this order 1e0-1e-2, while numbers of order 1e-3 or smaller might be considered approching instability and require some investigation. It is
always a good idea to check any solution, using xdraw or just skimming dcon.out.

Satisfied with the DCON result, we run GPEC::

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.

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 version 1.0.0
__________________________________________
Reading dcon eigenfuctions
Counting and reading dcon solutions
mlow = -12 , mhigh = 21 , mpert = 34
mstep = 1574 , mfix = 17 , msing = 4
Recontructing flux functions and metric tensors
Reading vacuum energy matrices
Calculating field on the boundary from coils
Building free boundary solutions
eigenmode = 1, dw = -2.364E-03, error = 3.608E-09
eigenmode = 2, dw = 3.918E-01, error = 1.151E-12
eigenmode = 3, dw = 1.413E+00, error = 1.776E-11
eigenmode = 4, dw = 2.932E+00, error = 4.966E-12
eigenmode = 5, dw = 3.617E+00, error = 6.209E-12
eigenmode = 6, dw = 4.939E+00, error = 5.212E-12
eigenmode = 7, dw = 7.459E+00, error = 4.384E-12
eigenmode = 8, dw = 1.033E+01, error = 3.098E-12
eigenmode = 9, dw = 1.050E+01, error = 7.427E-12
eigenmode = 10, dw = 1.301E+01, error = 4.758E-12
eigenmode = 11, dw = 1.602E+01, error = 4.668E-12
eigenmode = 12, dw = 1.760E+01, error = 2.422E-11
eigenmode = 13, dw = 1.913E+01, error = 5.817E-12
eigenmode = 14, dw = 2.227E+01, error = 6.625E-12
eigenmode = 15, dw = 2.553E+01, error = 6.178E-12
eigenmode = 16, dw = 2.650E+01, error = 4.967E-11
eigenmode = 17, dw = 2.876E+01, error = 6.067E-12
eigenmode = 18, dw = 3.129E+01, error = 3.658E-11
eigenmode = 19, dw = 3.229E+01, error = 5.655E-13
eigenmode = 20, dw = 3.677E+01, error = 1.296E-11
eigenmode = 21, dw = 3.718E+01, error = 3.762E-11
eigenmode = 22, dw = 4.168E+01, error = 1.153E-11
eigenmode = 23, dw = 4.429E+01, error = 2.995E-11
eigenmode = 24, dw = 5.306E+01, error = 2.081E-11
eigenmode = 25, dw = 6.506E+01, error = 1.280E-11
eigenmode = 26, dw = 8.191E+01, error = 8.032E-12
eigenmode = 27, dw = 9.998E+01, error = 7.700E-12
eigenmode = 28, dw = 1.069E+02, error = 5.822E-12
eigenmode = 29, dw = 1.389E+02, error = 1.135E-12
eigenmode = 30, dw = 1.347E+02, error = 1.387E-11
eigenmode = 31, dw = 1.374E+02, error = 1.473E-11
eigenmode = 32, dw = 2.171E+02, error = 6.004E-12
eigenmode = 33, dw = 2.329E+02, error = 9.530E-12
eigenmode = 34, dw = 6.934E+02, error = 2.111E-12
Calculating inductrances and permeability
Single mode permeability = 7.923E+00
Required energy to perturb vacuum = 9.309E+00
Required energy to perturb plasma = 3.862E+00
Amplification factor = 2.411E+00
Computing total resonant fields
psi = 5.936E-01, q = 2.000, total resonant field = 3.705E-04
psi = 8.186E-01, q = 3.000, total resonant field = 2.805E-04
psi = 9.282E-01, q = 4.000, total resonant field = 5.450E-05
psi = 9.884E-01, q = 5.000, total resonant field = 4.211E-04
Computing Clebsch displacements
volume = 10% Clebsch decomposition
volume = 20% Clebsch decomposition
volume = 30% Clebsch decomposition
volume = 40% Clebsch decomposition
volume = 50% Clebsch decomposition
volume = 60% Clebsch decomposition
volume = 70% Clebsch decomposition
volume = 80% Clebsch decomposition
volume = 90% Clebsch decomposition
volume = 100% Clebsch decomposition
Total cpu time for GPEC = 27 seconds
GPEC STOP => Normal termination.


Here we see the plasma has a single mode permeability ~8, meaning that the application of the first permeability eigenmode on the plasma control surface would result in a total flux (external + plasma response) 8 times the applied mode. The full permeability eigenmode and eigevector information can be found in gpec_response_n1.out. We also see a preview of the resonant field information, a more detailed collection of which is in gpec_singfld_n1.out. Finally, we have told GPEC to calculate and output the clebsch coordinate displacments in gpec_xclebsch_n1.out. This becomes an input for PENT.

Finally, we run PENT::

pentrc

PENTRC START => GPEC version 1.0.0
__________________________________________
clearing working directory
rm: cannot remove `pentrc_*.out': No such file or directory
Set idconfile:
euler.bin
Reading dcon eigenfuctions
Counting and reading dcon solutions
mlow = -12 , mhigh = 21 , mpert = 34
mstep = 1574 , mfix = 17 , msing = 4
Recontructing flux functions and metric tensors
Computing perturbed b field for gpec
Reading table from file:
g147131.02300_DIIID_KEFIT.kin
Reading table from file:
gpec_xclebsch_n1.out
-> WARNING: Assuming DCON hamada coordinates
-> calculating deltaB/B, divxi_prp
|----------| 0% iterations complete
|oo--------| 20% iterations complete
|oooo------| 40% iterations complete
|oooooo----| 60% iterations complete
|oooooooo--| 80% iterations complete
|oooooooooo| 100% iterations complete
PENTRC - full general-aspect-ratio calculation
psi = 0.0 -> T_phi = 7.63E-05 1.36E-04j
psi = 0.1 -> T_phi = 1.42E-02 5.57E-03j
psi = 0.2 -> T_phi = 2.26E-02 5.39E-03j
psi = 0.3 -> T_phi = 2.50E-02 9.84E-03j
psi = 0.4 -> T_phi = 2.80E-02 1.49E-02j
WARNING: vpar zero crossing internal to magnetic well at psi 4.783E-01
5.000E-01 <= 5.006E-01 <= 1.500E+00
-> Lambda, t/p boundry = 0.7737391481753563 0.7737464409937628
-> consider changing mtheta in equil.in
psi = 0.5 -> T_phi = 3.56E-02 2.03E-02j
psi = 0.6 -> T_phi = 4.20E-02 2.21E-02j
WARNING: vpar zero crossing internal to magnetic well at psi 6.333E-01
5.000E-01 <= 1.497E+00 <= 1.500E+00
-> Lambda, t/p boundry = 0.7299058478670797 0.7299127275381094
-> consider changing mtheta in equil.in
psi = 0.7 -> T_phi = 4.32E-02 2.33E-02j
psi = 0.8 -> T_phi = 5.03E-02 2.12E-02j
WARNING: vpar zero crossing internal to magnetic well at psi 8.265E-01
4.844E-01 <= 1.481E+00 <= 1.484E+00
-> Lambda, t/p boundry = 0.6744367173720330 0.6744430742230700
-> consider changing mtheta in equil.in
psi = 0.9 -> T_phi = 5.55E-02 3.58E-02j
WARNING: vpar zero crossing internal to magnetic well at psi 9.326E-01
4.492E-01 <= 4.495E-01 <= 1.449E+00
-> Lambda, t/p boundry = 0.6424701444514599 0.6424762000040400
-> consider changing mtheta in equil.in
---------------------------------------------
Total torque = 4.263E-001
Total Kinetic Energy = 2.657E-002
alpha/s = -8.022E+000
---------------------------------------------
WARNING: vpar zero crossing internal to magnetic well at psi 8.000E-01
4.844E-01 <= 4.847E-01 <= 1.484E+00
-> Lambda, t/p boundry = 0.6821990656221174 0.6822054956365666
-> consider changing mtheta in equil.in
WARNING: vpar zero crossing internal to magnetic well at psi 4.575E-01
5.039E-01 <= 1.502E+00 <= 1.504E+00
-> Lambda, t/p boundry = 0.7797137521775525 0.7797211013091262
-> consider changing mtheta in equil.in
PENTRC - trapped particle general-aspect-ratio calculation
psi = 0.0 -> T_phi = 9.05E-05 1.38E-04j
psi = 0.1 -> T_phi = 1.41E-02 3.46E-03j
psi = 0.2 -> T_phi = 2.25E-02 -9.46E-05j
psi = 0.3 -> T_phi = 2.49E-02 1.45E-03j
psi = 0.4 -> T_phi = 2.74E-02 4.33E-03j
psi = 0.5 -> T_phi = 3.56E-02 7.62E-03j
psi = 0.6 -> T_phi = 4.18E-02 7.99E-03j
psi = 0.7 -> T_phi = 4.31E-02 8.85E-03j
psi = 0.8 -> T_phi = 5.07E-02 5.19E-03j
psi = 0.9 -> T_phi = 5.42E-02 1.97E-02j
---------------------------------------------
Total torque = 4.027E-001
Total Kinetic Energy = 1.002E-002
alpha/s = -2.010E+001
---------------------------------------------
PENTRC - reduced large-aspect-ratio calculation
Reading F^-1/2_mnl from file:
../../../pentrc/fkmnl.dat
psi = 0.0 -> T_phi = 7.47E-06 1.01E-04j
psi = 0.1 -> T_phi = 9.18E-01 2.02E-01j
psi = 0.2 -> T_phi = 1.40E+00 1.76E+00j
psi = 0.3 -> T_phi = 1.66E+00 5.13E+00j
psi = 0.4 -> T_phi = 1.94E+00 1.07E+01j
psi = 0.5 -> T_phi = 2.33E+00 1.74E+01j
psi = 0.6 -> T_phi = 2.98E+00 2.32E+01j
psi = 0.7 -> T_phi = 3.62E+00 3.68E+01j
psi = 0.8 -> T_phi = 6.27E+00 5.75E+01j
psi = 0.9 -> T_phi = 1.47E+01 3.04E+02j
---------------------------------------------
Total torque = 2.039E+003
Total Kinetic Energy = -2.748E+002
alpha/s = 3.710E+000
---------------------------------------------
total cpu time = 8 minutes, 0 seconds
PENTRC STOP => Normal termination.
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``.

.. note:: The warnings about magnetic wells are ok in this case.
Finally, we run PENTRC to compute the neoclassical toroidal viscosity (NTV) torque from the perturbed field::

We see that the full and trapped-only calculations of the torque agree reasonably well, which is consistent with the understanding that particles on banana orbits dominate the nonambipolar transport. The RLAR result in this case is orders of magnitude off, but does have generaly the correct qualitative behavior.
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 typically differs significantly from the general-aspect-ratio calculations for shaped, high-beta H-mode plasmas.
Loading