GPEC - UPDATE - Synchronize the production fork with upstream develop 54991242 (1.6.0)

docs/conf.py

@@ -61,9 +61,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = '1.5'
+version = '1.6'
 # The full version, including alpha/beta/rc tags.
-release = '1.5.7'
+release = '1.6.0'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

docs/examples.rst

@@ -3,3 +3,11 @@
 .. _DIIID_ideal_example:
 
 .. include:: examples/DIIID_ideal_example/README
+
+.. _DIIID_kinetic_example:
+
+.. include:: examples/DIIID_kinetic_example/README
+
+.. _DIIID_resistive_example:
+
+.. include:: examples/DIIID_resistive_example/README

docs/examples/DIIID_ideal_example/README

@@ -3,7 +3,7 @@ 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 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.
+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::
 
@@ -23,6 +23,20 @@ A ``stride.in`` is also included, allowing STRIDE to be run as an alternative st
 
 .. 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
@@ -39,4 +53,4 @@ Finally, we run PENTRC to compute the neoclassical toroidal viscosity (NTV) torq
   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.
+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.

docs/examples/DIIID_ideal_example/slayer.in

@@ -0,0 +1,93 @@
+&SLAYER_INPUT
+    input_flag=f         ! reads profile quantities (n_e, t_e, etc.) from an ascii table
+    infile='TkMkr_D3Dlike_Hmode.gpeckf' ! Path to the kinetic profile file (read_eq=t), or to an ascii table of surface quantities (input_flag=t)
+    ncfile='stride_output_n1.nc' ! Path to netCDF file of STRIDE outputs
+    params_flag=f        ! calculates normalized parameters from profile quantities
+    ! >>>
+    ! >>> Input rational surface quantities
+    ! >>>
+    mm=2                 ! poloidal mode number
+    nn=1                 ! toroidal mode number
+    n_e=5.50E+19         ! electron density [m^-3]
+    t_e=1000.0           ! electron temperature [eV]
+    t_i=1000.0           ! ion temperature [eV]
+    sval=1.5             ! magnetic shear at the layer
+    bt=2.0               ! toroidal field [T]
+    rs=0.5               ! minor radius of resonant surface [m]
+    R0=1.67              ! major radius of magnetic axis [m]
+    omega=5.00E+03       ! ExB rotation frequency [rad/s]
+    l_t=0.11             ! temperature gradient scale length
+    l_n=0.25             ! density gradient scale length
+    qval=2.0             ! q of resonant surface
+    mu_i=2.0             ! ion mass ratio to proton mass, 2 for deuterium plasma. Also used when read_eq=t
+    zeff=2.0             ! plasma Z_effective. Also used when read_eq=t (the kinetic file's Z_eff is not used)
+    dr_val=-0.1          ! GGJ resistive interchange index D_R (taken from STRIDE when read_eq=t)
+    dgeo_val=0.0         ! unused (the toroidal critical Delta was removed); accepted so older inputs still read
+    chi_p_prof=0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2 ! Rational surface anomalous perpendicular energy diffusivity profile
+    chi_t_prof=0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2 ! Rational surface anomalous perpendicular ion momentum diffusivity profile
+    kappa_prof=0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 ! Rational surface thermal conductivity profile (NOTE: set > 0 to override chi_p_prof)
+    ! >>>
+    ! >>> Input normalized layer parameters. If read_eq=f, each will individually override internal calcuation using above inputs
+    ! >>>
+    inpr=0.0             ! Prantdl number. Set /= 0.0 to override internal calculation
+    inpe=0.0             ! electron viscosity
+    inQ=1.0              ! normalized ExB frequency. TO DO: update for mode rotation frequency calculation! 
+    inQ_e=0.0            ! normalized electron diamagnetic frequency. Set /= 0.0 to override internal calculation
+    inQ_i=0.0            ! normalized ion diamagnetic frequency. Set /= 0.0 to override internal calculation
+    inc_beta=0.0         ! dimensionless measure of the plasma pressure. Set /= 0.0 to override internal calculation
+    inds=0.0             ! normalized ion sound radius, or normalized ion skin depth. Set /= 0.0 to override internal calculation
+    intau=1.0            ! ratio of ion to electron temperature
+    Q0=1.0               ! unnecessary placeholder variable for inQ?
+    delta_prime=(0.0,0.0)  ! input Delta' for the growth rate calculation when read_eq=f, referenced to the minor radius (r_s*Delta')
+    delta_n_p=(0.01,0.01)  ! delta offset used for jxb calculations. Default for latter is (1e-2,1e-2).
+    ingamma=(0.0,0.0)     ! input (Re(Q),Im(Q)) to calculate inner layer Delta. Set /= 0.0 to override internal root finding
+/
+
+&SLAYER_CONTROL
+    inum=400             ! resolution to find error field thresholds.
+    jnum=100             ! resolution for 2d scan along with Q,omega.
+    knum=100             ! resolution for 2d scan along with the other.
+    Q_num=100            ! number of points per axis of the initial Re(Q),Im(Q) scan grid
+    scan_width=1.5       ! half-width of the Re(Q),Im(Q) scan window, in normalized Q
+    AMR_passes=4         ! number of adaptive refinement passes of the scan grid
+    !scan_center=0.0     ! Re(Q) center of the scan window. Unset: 0, or the middle of the Doppler-shifted branches of a coupled scan
+    msing_max=3          ! number of rational surfaces (counted from the core) to scan, and to couple when coupling_flag=t
+    dc_type="rfitzp"     ! critical Delta model: "rfitzp" (default) or "none"
+    read_eq=t            ! read the equilibrium and Delta' from STRIDE (ncfile) and profiles from the kinetic file (infile)
+    Pperp_Ptor_flag=t    ! Use R. Fitzpatrick formalism (inc. anomalous diffusion) for growth rate calculation
+    coupling_flag=f      ! solve the coupled multi-surface determinant instead of each surface separately. Requires read_eq=t
+    rotation_flag=t      ! Doppler-shift each surface by its ExB rotation in the coupled determinant (root reported in the lab frame)
+    QPscan_flag=f        ! scan (Q,P) space for delta and torque.
+    Qscan_flag=f         ! scan Q space
+    QPescan_flag=f       ! scan (Q,Pe) space for delta and torque.
+    Qbscan_flag=f        ! scan (Q,beta) space for delta and torque.
+    onscan_flag=f        ! scan (omega,n) space for error fields.
+    otscan_flag=f        ! scan (omega,t) space for error fields.
+    ntscan_flag=f        ! scan (n,te) space for error fields.
+    nbtscan_flag=f       ! scan (n,bt) space for error fields.
+    parflow_flag=f       ! set parallel flow on
+    PeOhmOnly_flag=t     ! only include Pe from Ohm's law
+    Pe_flag=f            ! if true, calculate classical inpe value from inpr and include in SLAYER run
+    layfac=0.02          ! layfac*EXP(ifac*ATAN2(AIMAG(Q-Q_e),REAL(Q-Q_e)) is added to Q_e if ABS(Q-Q_e)<layfac). Default is 0.02.
+/
+
+&SLAYER_OUTPUT
+    verbose=t            ! prints progress to terminal
+    ascii_flag=f         ! writes results to ascii files
+    bin_flag=f           ! writes results to binary files
+    netcdf_flag=f        ! writes results to netcdf files
+    est_gamma_flag=f     ! Calculate and output growth rate estimate
+    match_gamma_flag=t   ! match the layer Delta to Delta' and solve for the growth rate
+    stability_flag=f     ! calculate delta dependence on complex Q
+    stabscan_flag=f      ! Generate Re(Q),Im(Q) scan and output as ASCII
+    coupled_stabscan_flag=f      ! uniform-grid scan of the coupled determinant, output as ASCII
+    amr_flag=t           ! adaptive Re(Q),Im(Q) scan and root extraction for match_gamma_flag (results in the netcdf file)
+    br_th_flag=f         ! Calculate br_th
+/
+
+&SLAYER_DIAGNOSE
+    riccati_out=f        ! writes LSDOE Riccati integration to an ascii file
+    riccatiscan_flag=f
+    params_check=f       ! prints the primary parameters to the terminal when they are calculated
+    bal_flag=f           ! calculate the resonant field penetration threshold from torque balance
+/

docs/examples/DIIID_kinetic_example/README

@@ -30,4 +30,4 @@ Finally, we run PENTRC for the full neoclassical torque calculation::
   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 serves as a cross-check.
+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.

docs/index.rst

@@ -67,3 +67,5 @@ The GPEC repository contains a number of examples intended to provide both consi
 Example console commands and output are given below.
 
 - :ref:`DIII-D Ideal Example <DIIID_ideal_example>`
+- :ref:`DIII-D Kinetic Example <DIIID_kinetic_example>`
+- :ref:`DIII-D Resistive Example <DIIID_resistive_example>`

docs/installation.rst

@@ -50,7 +50,7 @@ Important make targets
 ----------------------
 
 * all (default): build all main executables and libraries
-* allf: same as all but Fortran-only modules
+* allf: currently builds the same modules as all
 * clean: remove built objects and module files
 * clear: remove example outputs in docs/examples
 * realclean: deep clean including dependencies built into deps
@@ -82,7 +82,7 @@ netCDF libraries:
 * NETCDF_FORTRAN_HOME: preferred when C and Fortran installs are split
 * NETCDFHOME or NETCDF_DIR: netCDF installation prefix
 * NETCDFINC: netCDF include path
-* NETCDF_C_HOME or NETCDF_C_DIR: netCDF-C location if separate
+* NETCDF_C_HOME: netCDF-C location if separate
 
 Dependency install location:
 
@@ -180,8 +180,8 @@ Then build with::
 X11 (xdraw)
 ^^^^^^^^^^^
 
-If xdraw fails to link on macOS, ensure XQuartz is installed and add include
-and library paths in xdraw/makefile, for example::
+If xdraw fails to link on macOS, ensure XQuartz is installed and set
+X11_HOME (read by install/DEFAULTS.inc), for example::
 
     brew install libx11
     export X11_HOME=/opt/homebrew

docs/namelists.rst

@@ -1,5 +1,5 @@
 
-Fully annotated namelist input files can be found in the $GPECHOME/input directory and any of the examples in $GPECHOME/docs/examples directories of any `public installation <installation>`.
+Fully annotated namelist input files can be found in the $GPECHOME/input directory and any of the examples in $GPECHOME/docs/examples directories of any :ref:`public installation <installation>`.
 
 The major input files requiring user interaction are as follows.
 
@@ -40,7 +40,7 @@ This namelist controls the MATCH code, including which DCON eigenfunction (sol_n
 The PENTRC Namelist
 ====================
 
-This namelist sets the inputs, controls the behavior, and determines the outputs of the GPEC code.
+This namelist sets the inputs, controls the behavior, and determines the outputs of the PENTRC code.
 
 .. include:: ../input/pentrc.in
    :literal: