Growth rate calculation functionality added to gslayer.f

docs/examples/DIIID_ideal_example/README

@@ -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

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
+/

input/slayer.in

@@ -1,44 +1,62 @@
 &SLAYER_INPUT
-
-    !!! FOR READING IN ASCII TABLE, SET INPUT_FLAG = TRUE
     input_flag=f         ! reads profile quantities (n_e, t_e, etc.) from an ascii table
-    infile=''  ! Path to ascii table of profile quantities read when using input_flag
-
-    !!! FOR MANUALLY ENTERING KINETIC QUANTITIES @ RATIONAL SURFACE, SET PARAMS_FLAG = TRUE
+    infile=''            ! Path to the kinetic profile file (read_eq=t), or to an ascii table of surface quantities (input_flag=t)
+    ncfile=''            ! 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=7.77E+17         ! electron density [m^-3]
-    t_e=25.8             ! electron temperature [eV]
-    t_i=25.8             ! ion temperature [eV]
-    sval=2.0             ! magnetic shear at the layer
-    bt=1.0               ! toroidal field [T]
-    rs=0.17              ! minor radius of resonant surface [m]
-    R0=2.0               ! major radius of magnetic axis [m]
-    omega=5.45E+04       ! ExB frequency [1/s]
+    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 ??
-    zeff=2.0             ! plasma Z_effective
-    inpr=10.0            ! Prantdl number
-    inpr_prof=10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0 ! Prandtl number profile
-    inpe=0.0             ! Electron viscosity (if > 0.0, Pe_flag=t will be overridden and this value used instead)
-
-    !!! IF INPUT_FLAG & PARAMS_FLAG ARE FALSE, MANUALLY ENTER DIMENSIONLESS QUANTITIES @ RATIONAL SURFACE HERE
-    inQ=23.0             ! normalized ExB frequency
-    inQ_e=2.0            ! normalized electron diamagnetic frequency
-    inQ_i=-2.0           ! normalized ion diamagnetic frequency
-    inc_beta=0.7         ! dimensionless measure of the plasma pressure
-    inds=6.0             ! normalized ion sound radius
-    intau=1.0            ! ion temperature divided by electron temperature
-    Q0=4.0               ! unnecessary placeholder variable for inQ?
-    delta_n_p=(0,0)      ! complex delta offset used for jxb calculations. Default is (1e-2,1e-2).
+    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 alont with the other.
+    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=2          ! 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=f            ! 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.
@@ -49,18 +67,27 @@
     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            ! Not operational yet. If true, it will calculate classical inpe value from inpr and include in SLAYER run
+    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=t         ! writes results to ascii files
+    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=t     ! 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
-    params_check=t       ! prints the primary parameters to the terminal when they are calculated
+    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
 /

install/DEPENDENCIES.inc

@@ -183,10 +183,13 @@
 
 # SLAYER
 ../slayer/delta.o: ../slayer/sglobal.o
-../slayer/gslayer.o: ../slayer/delta.o ../slayer/sglobal.o
+../slayer/growthrates.o: ../slayer/delta.o ../slayer/layerinputs.o ../slayer/params.o ../slayer/sglobal.o ../slayer/slayer_netcdf.o
+../slayer/gslayer.o: ../slayer/delta.o ../slayer/layerinputs.o ../slayer/params.o ../slayer/sglobal.o ../slayer/slayer_netcdf.o
+../slayer/layerinputs.o: ../equil/bicube.o ../equil/equil.o ../equil/spline.o ../pentrc/inputs.o ../slayer/params.o ../slayer/sglobal.o ../slayer/slayer_netcdf.o
 ../slayer/params.o: ../slayer/sglobal.o
 ../slayer/sglobal.o: ../equil/local.o
-../slayer/slayer.o: ../slayer/delta.o ../slayer/params.o ../slayer/sglobal.o
+../slayer/slayer.o: ../slayer/delta.o ../slayer/growthrates.o ../slayer/gslayer.o ../slayer/layerinputs.o ../slayer/sglobal.o
+../slayer/slayer_netcdf.o: ../slayer/sglobal.o ../slayer/version.inc
 
 # COIL
 ../coil/coil.o: ../coil/version.inc ../equil/bicube.o ../equil/fspline.o

install/RULES.inc

@@ -125,7 +125,7 @@ FCMULTI = $(FC) $(FFLAGS) -I../equil
 	$(FCMULTI) $(MODFLAG) ../multi -c $< -o $@
 
 # SLAYER
-FCSLAYER = $(FC) $(FFLAGS) $(OMPFLAG) -I$(MATHINC) -I$(NETCDFINC) -I../equil
+FCSLAYER = $(FC) $(FFLAGS) $(OMPFLAG) -I$(MATHINC) -I$(NETCDFINC) -I../equil -I../slayer -I../pentrc
 ../slayer/%.o: ../slayer/%.f
 	$(FCSLAYER) $(MODFLAG) ../slayer -c $< -o $@
 ../slayer/%.o: ../slayer/%.F

install/TARGETS.inc

@@ -298,7 +298,18 @@ rmatch: v neededdeps $(RMATCH_OBJS) equil lsode
 	$(call print_success,rmatch -> $(BINDIR)/rmatch)
 
 # SLAYER
+../slayer/version.inc:
+	@ver="`git describe --tags --match 'v[0-9]*'`"; \
+	if [ ! -f $@ ] || ! grep -F "'$$ver'" $@ >/dev/null 2>&1; then \
+		echo ">>> Updating version file"; \
+		echo "      CHARACTER(len=*), PARAMETER :: version ='$$ver'" | tee $@; \
+	else \
+		echo ">>> Git version unchanged!"; \
+	fi
+	$(call print_success,version file updated: $@)
+
 SLAYER_LIBS = \
+	-lpentrc \
 	-llsode \
 	-lequil \
 	$(MATHLIBS) \
@@ -308,7 +319,7 @@ SLAYER_LIBDIRS = \
 	-L$(LIBDIR) \
 	-L$(MATHDIR) \
 	-L$(NETCDFDIR)
-slayer: v neededdeps $(SLAYER_OBJS) lsode equil
+slayer: v neededdeps $(SLAYER_OBJS) lsode equil pentrc
 	mkdir -p $(BINDIR)
 	rm -f ../slayer/slayer
 	$(FCSLAYER) -o ../slayer/slayer $(SLAYER_OBJS) $(SLAYER_LIBDIRS) $(SLAYER_LIBS) $(LDFLAGS)

install/gen_deps.py

@@ -39,7 +39,7 @@
     "rmatch": ["equil"],
     "sum":    ["equil"],
     "multi":  ["equil"],
-    "slayer": ["equil"],
+    "slayer": ["equil", "pentrc"],
     "coil":   ["equil"],
     "gpec":   ["equil", "vacuum", "pentrc", "coil", "slayer"],
 }