diff --git a/docs/conf.py b/docs/conf.py index 046fe5f3f..12c00f0a2 100644 --- a/docs/conf.py +++ b/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. diff --git a/docs/examples.rst b/docs/examples.rst index e9f79fa29..f83b9830e 100644 --- a/docs/examples.rst +++ b/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 diff --git a/docs/examples/DIIID_ideal_example/README b/docs/examples/DIIID_ideal_example/README index bc4e65dcb..d27014464 100644 --- a/docs/examples/DIIID_ideal_example/README +++ b/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. diff --git a/docs/examples/DIIID_ideal_example/slayer.in b/docs/examples/DIIID_ideal_example/slayer.in new file mode 100644 index 000000000..418990132 --- /dev/null +++ b/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) ... 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. diff --git a/docs/index.rst b/docs/index.rst index aa1df6a0c..09be708fc 100644 --- a/docs/index.rst +++ b/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 ` +- :ref:`DIII-D Kinetic Example ` +- :ref:`DIII-D Resistive Example ` diff --git a/docs/installation.rst b/docs/installation.rst index f7b488392..43813b04e 100644 --- a/docs/installation.rst +++ b/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 diff --git a/docs/namelists.rst b/docs/namelists.rst index 4629dade7..491cd44d2 100644 --- a/docs/namelists.rst +++ b/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 `. +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 `. 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: diff --git a/docs/outputs.rst b/docs/outputs.rst index a0d80b2ad..22d35140b 100644 --- a/docs/outputs.rst +++ b/docs/outputs.rst @@ -125,46 +125,46 @@ The binary GPEC outputs can be viewed using the commandxdraw filenamewhere filen Output Files ============ -The GPEC outputs are entirely specified by flags (bool types t or f) set in the GPEC_output section of GPEC.in. All outputs are ASCII files, and can be grouped into two major categories. +The GPEC outputs are entirely specified by flags (bool types t or f) set in the gpec_output section of gpec.in. Outputs are written as ASCII tables (ascii_flag), xdraw binary files (bin_flag, bin_2d_flag) and netcdf files (netcdf_flag), and can be grouped into the categories below. ASCII File Outputs ------------------ -A number of the flag options in GPEC.in instruct GPEC to output ASCII file data. Some of these outputs are always available. Some, however, require a input error field instead of a hard coded harmonic_flag call. Both groups are listed in detail here. +A number of the flag options in gpec.in instruct GPEC to output ASCII file data. Some of these outputs are always available. Some, however, require an external 3D field input (data_flag, coil_flag, harmonic_flag) or a single eigenmode (mode_flag). Both groups are listed in detail here. Outputs Always Available -These GPEC outputs can always be obtained from a equil.bin file output from DCON. +These GPEC outputs can always be obtained from a euler.bin file output from DCON. -GPEC_response_n#.out +gpec_response_n#.out **Flag** resp_flag **Info** Energy for DCON eigenmodes and stability indices. Eigenvalues and eigenvectors for vacuum and plasma inductance (virtual casing currents to fields), plasma permeability (external fields to total fields), and plasma reluctance. -GPEC_singcoup_matrix_n#.out +gpec_singcoup_matrix_n#.out **Flag** singcoup_flag **Info** The coupling matrix to resonant fields, coupling matrix to singular currents, and to island half-widths is given for each rational surface within the plasma (q=2, 3, etc) for each surface the real and imaginary coupling constants are given for each poloidal mode number on the control surface. -GPEC_singcoup_svd_n#.out +gpec_singcoup_svd_n#.out **Flag** singcoup_flag - **Info** The SVD singular values (s) and eigen vectors for each coupling matrix in GPEC_singcoup_matrix_n#.out. Large s corresponds to large amplification, with the largest (most important mode) listed at the top. The results should be dotted with the unweighted normal field spectrum to give physical meaning. + **Info** The SVD singular values (s) and eigen vectors for each coupling matrix in gpec_singcoup_matrix_n#.out. Large s corresponds to large amplification, with the largest (most important mode) listed at the top. The results should be dotted with the unweighted normal field spectrum to give physical meaning. Outputs Available When Error Field is Provided -These outputs are only available when an external error field file is provided as an input to GPEC. This means GPEC.in must have the data_flag turned on and a infile specified. +These outputs are only available when an external 3D field is provided as an input to GPEC. This means gpec.in must have data_flag (with an infile specified), coil_flag, harmonic_flag, or mode_flag turned on. -GPEC_control_n#.out +gpec_control_n#.out - **Flag** + **Flag** none (always written with ascii_flag) **Info** The Plasma response for an external perturbation on the control surface. This includes the vacuum energy, surface energy, plasma energy, real and imaginary vacuum input \mathbf{B}_{in} and total field on plasma boundary\mathbf{B}_{out}as a function of poloidal mode number. -GPEC_singfld_n#.out +gpec_singfld_n#.out **Flag** singfld_flag @@ -174,21 +174,21 @@ GPEC_singfld_n#.out **Info** Additional section showing the overlap field (Tesla) and overlap percentage for each eigenmode in the singcoup_svd output. See :ref:`coupling_overlap` for definitions. -GPEC_pmod_n#.out +gpec_pmodb_n#.out **Flag** pmodb_flag - **Info** Eulerian and Lagrangian \left|\mathbf{B}\right|(real and imaginary) for each poloidal mode number at each value of \Psi_{N} output. This output is necessary for NTV post processing. + **Info** Perturbed field strength (real and imaginary, Tesla) for each poloidal mode number at each value of \Psi_{N} output: the Eulerian \delta B (eulb), the Lagrangian \delta B of the displacement with zero \xi^\zeta in the working (jac_type) coordinates (lagb), B\nabla\cdot\xi_\perp (Bdivxprp) and the force-balance form of B\kappa\cdot\xi_\perp (Bkxprp). lagb is not the \xi_\perp-only Lagrangian \delta B; that is -(Bdivxprp+Bkxprp) up to the equilibrium force-balance residual. See :doc:`pmodb_gauge` before using these columns for NTV post processing. -GPEC_xbnormal_n#.out +gpec_xbnormal_n#.out **Flag** xbnormal_flag - **Info** The normal components of the displacement, magnetic field without the plasma response, and magnetic field with the plasma response included for each poloidal mode number at each value of \Psi_{N} output.?? + **Info** The normal displacement (xno), the normal perturbed magnetic field (bno), and the Jacobian-weighted contravariant psi component of the perturbed field J\mathbf{B}\cdot\nabla\psi (bwp) for each poloidal mode number at each value of \Psi_{N} output. -GPEC_*rzphi_n#.out +gpec_*rzphi_n#.out -A number of output files have a similar structure. Here the * in the file name is replaced by the appropriate leading letters of the corresponding flag. For example the xrzphi_flag for n=1 creates a GPEC_xrzphi_n1.out file. some common properties of these files are: +A number of output files have a similar structure. Here the * in the file name is replaced by the appropriate leading letters of the corresponding flag. For example the xrzphi_flag for n=1 creates a gpec_xrzphi_n1.out file. some common properties of these files are: • real and imaginary components: Output files contain two dimensional data on an \left(r,z\right) grid for a single toroidal harmonic. To translate into three dimensions, perform the transformationB\left(r,z,\phi\right)=B_{real}\left(r,z\right)\cos\left(n\phi\right)+B_{imag}\left(r,z\right)\sin\left(n\phi\right) @@ -211,7 +211,7 @@ A number of output files have a similar structure. Here the * in the file name i **Info** The false perturbed magnetic field in the vacuum region on the \left(r,z,\phi\right) grid calculated using the GPEC boundary surface current composed of both the vacuum component and the plasma response. - **Flag** vpbrzphi_flag + **Flag** pbrzphi_flag (writes gpec_vpbrzphi_n#.out) **Info** The true perturbed magnetic field in the vacuum region on the \left(r,z,\phi\right) grid due to the plasma response alone calculated from the plasma response surface condition. @@ -219,9 +219,9 @@ A number of output files have a similar structure. Here the * in the file name i **Info** The false perturbed magnetic field in the vacuum region on the \left(r,z,\phi\right) grid calculated using the GPEC boundary surface current from the external fields alone. - **Flag** ssbrzphi_flag + **Flag** vsbrzphi_flag (with ss_flag(i) selecting surfaces; requires singfld_flag) - **Info** The false perturbed magnetic field in the vacuum region on the \left(r,z,\phi\right) grid calculated using the GPEC boundary surface current from the external fields alone. + **Info** The vacuum field on the \left(r,z,\phi\right) grid due to the currents on the selected resonant surfaces. Energy Decomposition Outputs ---------------------------- @@ -281,7 +281,7 @@ pflux_re(im)_2d.bin bnormal_spectrum.bin - **Flag** bin_flag + **Flag** bin_2d_flag **Info** Surfmn type contours of the normal perturbed magnetic fields as a function of poloidal harmonic number and \Psi_{N}. diff --git a/docs/pmodb_gauge.rst b/docs/pmodb_gauge.rst new file mode 100644 index 000000000..5da0559f9 --- /dev/null +++ b/docs/pmodb_gauge.rst @@ -0,0 +1,247 @@ +.. _pmodb_gauge: + +****************************************** +Perturbed field strength (pmodb) and gauge +****************************************** + +This page states what the columns of ``gpec_pmodb_n#.out`` are, how they are +related, and which of them a field-strength consumer needs. + +.. contents:: On this page + :local: + :depth: 1 + + +Coordinates and displacement +============================ + +Assume smooth nested surfaces with B>0 and nonzero J, chi1 and q. +Fourier inversions below additionally exclude nonzero resonant sources. + +GPEC works in the DCON chart :math:`(\psi,\theta,\zeta)` with angles in turns, +any straight-field-line Jacobian :math:`\mathcal{J}(\psi,\theta)` selected by +``jac_type``, :math:`\mathcal{J}B^\theta=\chi_1`, +:math:`\mathcal{J}B^\zeta=\chi_1 q` and Fourier factor +:math:`\exp[2\pi i(m\theta-n\zeta)]`. The Euler-Lagrange solution supplies +:math:`\xi^\psi`, its radial derivative and the surface displacement +(Glasser 2016, Eqs. 9 and 18) + +.. math:: + + \xi_s = (\boldsymbol\xi\times\mathbf B)_\psi = \chi_1(q\xi^\theta-\xi^\zeta), + +the in-surface displacement across :math:`\mathbf B` (``xms``; +``gpec_xclebsch`` writes :math:`\xi^\alpha=\xi_s/\chi_1`). A field-aligned +displacement changes neither :math:`\xi^\psi` nor :math:`\xi_s`, so +:math:`(\xi^\psi,\xi_s)` determine :math:`\boldsymbol\xi_\perp` completely. +The parallel part used elsewhere in GPEC is fixed in ``gpeq_contra`` by +:math:`\nabla\cdot\boldsymbol\xi=0`, which gives ``xmt`` +:math:`=\mathcal J\xi^\theta` and ``xmz`` :math:`=\mathcal J\xi^\zeta`. +The pmodb columns do not use ``xmt`` or ``xmz``. + +For every ``jac_type`` the Jacobian-weighted perturbed field of +``gpeq_sol`` is :math:`\mathcal J\,\nabla\times(\boldsymbol\xi\times\mathbf B)`. + + +The pmodb columns +================= + +``eulb`` + :math:`\delta B_E=\hat{\mathbf b}\cdot\delta\mathbf B`, the Eulerian + field strength. + +``lagb`` + :math:`\delta B_E+\xi^\psi\partial_\psi B+\xi_s/(\chi_1 q)\,\partial_\theta B`. + This is the Lagrangian field strength of the displacement + :math:`\boldsymbol\xi_g=\boldsymbol\xi_\perp+\lambda\mathbf B` with + :math:`\xi_g^\zeta=0` in the *working* chart (``jac_type``). + It is neither the :math:`\boldsymbol\xi_\perp`-only value nor Park et + al. (2009) Eq. 7 with the divergence-free :math:`\boldsymbol\xi` that + GPEC constructs (that value is ``lagb`` + :math:`+\,\xi^\zeta/q\,\partial_\theta B`). ``jac_out`` only + relabels the harmonics; it does not change this gauge. + +``Bdivxprp`` + :math:`B\,\nabla\cdot\boldsymbol\xi_\perp`. ``gpout_pmodb`` forms + :math:`\mathcal J\boldsymbol\xi_\perp` directly from :math:`\xi^\psi` and + :math:`\xi_s`, + + .. math:: + + \mathcal J\xi_\perp^\theta=\chi_1^2\left(-f_2\,\xi^\psi+f_3\,\xi_s/\chi_1\right), + \qquad + \mathcal J\xi_\perp^\zeta=q\,\mathcal J\xi_\perp^\theta-\mathcal J\xi_s/\chi_1, + + with the equilibrium coefficients + :math:`f_2=(g_{\psi\theta}+qg_{\psi\zeta})/(\mathcal JB^2)` and + :math:`f_3=(g_{\theta\zeta}+qg_{\zeta\zeta})/(\mathcal JB^2)` (``eqfun`` + columns 2 and 3, covariant metric). The same projection gives the + ``xmat``-``zmat`` matrices of ``idcon_build``. + +``Bkxprp`` + :math:`\boldsymbol\xi_\perp\cdot\nabla B+\mu_0p'\xi^\psi/B`. This is + :math:`B\boldsymbol\kappa\cdot\boldsymbol\xi_\perp` only through + equilibrium force balance: exactly, + :math:`B\boldsymbol\kappa\cdot\boldsymbol\xi_\perp=\boldsymbol\xi_\perp\cdot\nabla B + +\boldsymbol\xi_\perp\cdot[(\nabla\times\mathbf B)\times\mathbf B]/B`. + The PENTRC matrices ``smat``, ``tmat`` (``idcon_build``) form the + same quantity times :math:`\mathcal J B^2` from :math:`\xi^\psi` and + :math:`\xi^\alpha`. + +``llagb`` (``gpec_diagnostics_pmodb``, ``pmodbmn_flag``) + :math:`\delta B_L(\boldsymbol\xi_\perp)=\delta B_E+\boldsymbol\xi_\perp\cdot\nabla B` + ``lllagb`` is :math:`-(`\ ``Bdivxprp``\ :math:`+`\ ``Bkxprp``\ :math:`)`, + and ``Bkxprp`` :math:`=` ``llagb`` :math:`-` ``eulb`` :math:`+\mu_0p'\xi^\psi/B`. + + +Relations between the columns +============================= + +Gauge of ``lagb``: + +.. math:: + + \mathrm{lagb}=\delta B_L(\boldsymbol\xi_\perp)-\frac{\xi_\perp^\zeta}{q}\partial_\theta B, + \qquad + \xi_\perp^\zeta=-\frac{B_\theta\xi_s/\chi_1+qB_\psi\xi^\psi}{B_\theta+qB_\zeta}. + +Kinematic identity (any :math:`\boldsymbol\xi_\perp`, no force balance) +and its pmodb form: + +.. math:: + + \delta B_L(\boldsymbol\xi_\perp)=-B(\nabla\cdot\boldsymbol\xi_\perp + +\boldsymbol\kappa\cdot\boldsymbol\xi_\perp) + =-(\mathrm{Bdivxprp}+\mathrm{Bkxprp})-\boldsymbol\xi_\perp\cdot\mathbf F/B, + \qquad \mathbf F=(\nabla\times\mathbf B)\times\mathbf B-\mu_0\nabla p. + +Hence + +.. math:: + + \mathrm{lagb}+\mathrm{Bdivxprp}+\mathrm{Bkxprp} + =-\frac{\xi_\perp^\zeta}{q}\partial_\theta B-\frac{\boldsymbol\xi_\perp\cdot\mathbf F}{B}. + +``lagb`` :math:`=-(`\ ``Bdivxprp``\ :math:`+`\ ``Bkxprp``\ :math:`)` is +therefore not an identity. The first term is a gauge term of first order in +the perturbation; its size depends on :math:`B_\theta/(qB_\zeta)`, +:math:`B_\psi`, :math:`\partial_\theta B` and the displacement spectrum. +The second term vanishes for an exact equilibrium. In the output, finite +Fourier truncation and spline error add to both. + + +Regularisation +============== + +With ``reg_flag`` (Park et al. 2009 Eq. 10), ``gpeq_sol`` replaces the +radial derivative by :math:`R\,\partial_\psi\xi^\psi`, +:math:`R=s^2/(s^2+\sigma^2)`, :math:`s=m-nq`, in the field and in ``xms``. +The regularised field is then not divergence free: + +.. math:: + + \mathcal J\nabla\cdot\delta\mathbf B=2\pi i\chi_1 s(1-R)\,\partial_\psi\xi^\psi. + +``Bdivxprp`` uses the same derivative ``xmp1`` as the field, and +:math:`\boldsymbol\xi_\perp` is built from the same ``xms``, so the +regularisation cancels between the columns and the relation above holds +unchanged with the regularised quantities. +Near rational surfaces the regularised field and displacement do not form +an exact ideal-induction pair. This is a limitation of the smoothing +policy, not evidence of resistive physics. With ``reg_flag=f`` the +columns are unregularised and singular at rational surfaces. + + +What consumers need +=================== + +PENTRC + The drive of Logan et al. (2013) Eqs. 4-5 needs :math:`\delta B` and + :math:`\nabla\cdot\boldsymbol\xi` of the *same* displacement; the + combination is then invariant under + :math:`\boldsymbol\xi\to\boldsymbol\xi+\lambda\mathbf B` up to a total + derivative along the field line. Both PENTRC input paths use the + :math:`\boldsymbol\xi_\perp` pair + :math:`\delta B/B=-(\nabla\cdot\boldsymbol\xi_\perp+\boldsymbol\kappa\cdot\boldsymbol\xi_\perp)` + and :math:`\nabla\cdot\boldsymbol\xi_\perp`: the default path + (``peq_file`` = ``gpec_xclebsch``, ``set_peq``) from the ``smat``-``zmat`` + matrices, and the ``pmodb_file`` path (``read_pmodb``) from + :math:`-(\mathrm{Bdivxprp}+\mathrm{Bkxprp})` and ``Bdivxprp``. The + ``lagb`` column is not read. Regenerate ``pmodb_file`` inputs with GPEC + 1.6 or later. + +Bounce-averaged codes in perturbed Boozer or Hamada coordinates (NEO-RT, NEO-2, Shaing-type) + need :math:`\delta B` in the perturbed straight-field-line coordinates. + That follows from :math:`\delta B_L(\boldsymbol\xi_\perp)` (``llagb``, or + :math:`-(\mathrm{Bdivxprp}+\mathrm{Bkxprp})` away from rational surfaces) + and :math:`\nabla\cdot\boldsymbol\xi_\perp` by a magnetic differential + equation for the parallel gauge (next section). ``lagb`` must not be + used as :math:`\delta B_L(\boldsymbol\xi_\perp)`. + + +Target gauges in GPEC units +=========================== + +Let :math:`D=\partial_\theta+q\partial_\zeta`, so +:math:`\mathbf B\cdot\nabla=(\chi_1/\mathcal J)D`. +For a matching equilibrium chart :math:`\mathcal J B^a=H(\psi)`, +:math:`a=2` is Boozer and :math:`a=0` is Hamada. A displacement with +:math:`b_X=\delta B_L(\boldsymbol\xi_X)` and +:math:`d_X=\nabla\cdot\boldsymbol\xi_X` has target-gauge source + +.. math:: + + T_a=d_X+a b_X/B,\qquad + \boldsymbol\xi_{target}=\boldsymbol\xi_X+\lambda\mathbf B,\qquad + T_a(\boldsymbol\xi_{target})=T_a(\boldsymbol\xi_X) + +\frac{\mathbf B\cdot\nabla(\lambda B^a)}{B^a}. + +For nonaxisymmetric modes the target flux-function perturbation is zero. +With :math:`w=\chi_1\lambda/\mathcal J`, this gives + +.. math:: + + Dw=-T_a,\qquad w_m=\frac{i T_{a,m}}{2\pi(m-nq)},\qquad + b_{target}=b_X+w\partial_\theta B. + +The source is formed in real space before Fourier transformation. A +nonzero source at :math:`m=nq` has no periodic solution. A resonant zero +source leaves an arbitrary homogeneous shift, set to zero here. +The factor :math:`2\pi` is from angles in turns and cancels the same +factor in :math:`\partial_\theta B`. Relabeling equilibrium harmonics +with ``jac_out`` alone does not solve this equation. + +For a PEST chart :math:`\mathcal J/R^2` is a flux function; its +linear perturbation is :math:`\nabla\cdot\boldsymbol\xi-2\xi_R/R`. +PEST output is not automatically either target gauge. Transform the +equilibrium chart first, then solve the target equation. + + +History +======= + +Before GPEC 1.6.0, ``Bdivxprp`` and ``Bkxprp`` subtracted a +Fourier-truncated parallel displacement, ``gpout_pmodb`` regularised the +divergence even with ``reg_flag=f``, ``gpdiag_pmodb`` used a stale Jacobian, +PENTRC's ``read_pmodb`` paired ``lagb`` with an inconsistent divergence, and +both PENTRC CGL paths halved the curvature term. These were fixed in +1.6.0 (#340-#343). + + +Source map and literature +========================= + +Code: ``gpeq_sol`` and ``gpeq_contra`` (``gpec/gpeq.f``, regularisation and +the divergence-free closure), ``gpout_pmodb`` (``gpec/gpout.f``, the pmodb +columns), ``gpdiag_pmodb`` (``gpec/gpdiag.f``, ``llagb`` and ``lllagb``), +``idcon_build`` (``gpec/idcon.f``, ``smat``-``zmat``), ``read_pmodb`` and +``set_peq`` (``pentrc/inputs.f90``), and ``direct_run`` (``equil/direct.f``, +``eqfun``). + +The physical definitions were checked against Park 2009, PoP 16 056115, +Eqs. 7 and 10; Park 2011, PoP 18 110702, Eqs. 4, 10 and 11; +Park and Logan 2017, PoP 24 032505, Eq. 13 and p. 6; Logan and Park 2013, +PoP 20 122507, Eqs. 4 and 5; and Boozer and Nuehrenberg 2006, +PoP 13 102501, Eqs. 42 and 43. The latter's Delta includes twice the +curvature; its Eq. 43 is Eulerian and is not the perpendicular Lagrangian +scalar. diff --git a/docs/releases.rst b/docs/releases.rst index f24f36b92..17d70c578 100644 --- a/docs/releases.rst +++ b/docs/releases.rst @@ -10,6 +10,129 @@ Note, backwards compatibility is defined as the ability to return to a previous The `github release notes `_ are reproduced below. +GPEC v1.6.0 +=========== + +This is planned to be the last minor release of the Fortran GPEC 1.x series. +New feature development continues in the Julia GPEC 2.x code; the 1.6.x line will receive bug-fix patch releases only. +GPEC is now distributed under the MIT license. + +Read the Upgrade notes first: the build system, several binary interface files, and a few namelist and output names changed, and many numerical results change because long-standing bugs were fixed (see Fixes). + +Upgrade notes +------------- +- INSTALL - The build system was rewritten. Build only from ``install/`` with ``make`` (``make -j N`` is supported); the per-directory makefiles are gone. ``FC`` and ``CC`` must be set, git is required at build time, and missing netCDF/HDF5/OpenBLAS are built from bundled git submodules (``git submodule update --init``). ``DEPSINSTALLDIR`` installs them to a shared location; ``make realclean`` is required before switching compilers. +- ALL - Rerun the full chain after upgrading. DCON ``euler.bin`` has a new record (``wv_farwall_flag``), and the resistive interface files ``delta.bin``, ``gal_solution.bin``, ``globalsol.bin`` and ``vmat.bin`` changed format, so RDCON, RMATCH and GPEC must be rerun together. +- DCON, RDCON - ``vacuum.bin`` is no longer written by default (set ``bin_vac=t``). GPEC always recomputes the vacuum matrices and ignores ``ivacuumfile``. ``ahg2msc.out`` is no longer written; vacuum data are passed in memory (``out_ahg2msc`` in ``equil.in`` restores the files). +- RDCON - ``coil%rpec_flag`` moved from ``&GAL_OUTPUT`` to ``&RDCON_CONTROL`` in ``rdcon.in``. Setting it in ``&GAL_OUTPUT`` still works but prints a deprecation warning. +- RMATCH - ``ntor`` is no longer an input (n is read from ``delta.bin``); remove it from any ``rmatch.in`` that sets it. The default inner-layer solver is now ``inps_type="inps"`` with ``nx=512``, and the eigenvalue search uses the ``2*pi*n*Omega`` rotation convention of RPEC. +- DCON, PENTRC, COIL, STRIDE - Thread counts default to ``OMP_NUM_THREADS`` (inputs default to 0). ``parallel_threads`` is renamed ``dcon_kin_threads`` and ``openmp_threads`` is renamed ``pentrc_threads``; the old names still work with a deprecation warning. +- EQUIL - ``newq0 < 0`` (including the undocumented ``newq0=-1`` sign flip) now stops with an error. ``eq_type`` ``marklin_direct``, ``marklin_inverse`` and ``hansen_inverse`` were removed; ``chease2`` is renamed ``chease_ascii`` (the old name still works with a warning). +- GPEC - Renamed outputs: netCDF ``w_isl_crit`` is now ``w_isl_v_crit``; cylindrical ``vcbr``, ``vcbz``, ``vcbp`` are now ``b_r_vacuum``, ``b_z_vacuum``, ``b_p_vacuum``. The ``gpec_singfld_n#.out`` ASCII columns changed (new threshold, island-width and pressure columns). ``I_overlap_norm`` and the DCON/RDCON ``W_v_eigenvalue`` are now 1-D real arrays. +- RDCON - The netCDF output is ``rdcon_output_n#.nc`` (it no longer overwrites DCON's ``dcon_output_n#.nc``). +- PENTRC - With ``pmodb_file``, PENTRC now reads the supplied ``Bdivxi_perp`` column. Regenerate ``gpec_pmodb_n#.out`` with GPEC 1.6 before using it. +- DOCS - All DIII-D examples now use a public TokaMaker DIII-D-like H-mode equilibrium and kinetic profiles in place of the KEFIT 147131 files. + +Adds +------ +- SLAYER - New slab-layer linear drift-MHD module (Park, PoP 2022) with a standalone ``slayer`` executable, ``input/slayer.in`` and an a10 example. +- SLAYER - Tearing-mode growth rates at the STRIDE rational surfaces by matching the inner-layer and outer Delta' (optional critical Delta, ``dc_type``), per surface or coupled, with the E x B Doppler shift (``slayer_output_n#.nc``; see the DIII-D ideal example). +- GPEC - SLAYER resonant-field penetration threshold at each rational surface (``singthresh_slayer_flag``; ``singthresh_callen_flag`` selects the Callen threshold and ``singthresh_flag`` enables both). New netCDF ``Phi_res_crit``. +- GPEC - Prandtl number for the SLAYER threshold set by ``singthresh_slayer_inpr`` (default 5.0) or per rational surface by ``singthresh_slayer_inpr_prof``. +- GPEC - Callen saturated and minimum island widths (``w_isl_sat``, ``w_isl_min``) and the Callen critical resonant field ``Phi_res_crit_callen``. +- GPEC - Estimate of the pressure lost to Callen-saturated islands (``dP_res``, ``P_res`` and the modified pressure profile ``p_mod``). +- GPEC - Bernstein/Chance decomposition of the eigenmode potential energy (``effective_field_flag``, ``standard_form_flag``, ``effective_field_terms_flag``, ``curl_check_flag``, ``recon_out``, ``recon_int``) with ``gpec_recon*`` outputs and a derivation note in ``docs/tex/recon``. +- GPEC - Mutual inductance between the plasma surface and a conformal wall and its inverse (``M_w``, ``M_w_inv``); ``wv_farwall_flag`` in ``dcon.in`` keeps no-wall eigenmodes for GPEC while the wall inductance is computed; ``mutual_test_flag`` diagnostics. +- GPEC - Resistive perturbed equilibria (RPEC) carry the analytic Galerkin derivative d(xi)/dpsi and the RMATCH field through to GPEC. Per-surface resistive sweet spots from the Lundquist number (``use_res_spot``, ``sing_resspot``, ``sing_interpspot``) and new outputs ``sing_spots``, ``Sfac``, ``tauA``, ``tauR``, ``Jbgradpsi_rmatch`` and ``b_pen``. +- GPEC - Dominant-mode overlap outputs ``dominant_mode_overlap`` and ``dominant_mode_local_overlap``. +- GPEC - Full set of DCON plasma attributes and local SLAYER inputs in the profile and control netCDF outputs. +- COIL - Shift and tilt coils from ``coil.in`` (``coil_shift[xyz]``, ``coil_tilt[xyz]``, ``coil_[xyz]nom``, ``cnpert``), with optional netCDF output of the modified geometry (``coil_out``). +- COIL - OpenMP-parallel Biot-Savart fields with a new ``coil_threads`` input; the maximum number of coils is raised from 50 to 500; MAST-U, JET ex-vessel and extra ITER coil sets. +- RDCON - netCDF output ``rdcon_output_n#.nc`` (``netcdf_out``, default true) with the ``dcon.out`` content, the Delta, Delta', Delta_gw and Delta_coil matrices and the Galerkin cell layout. +- RDCON - Optional modified Rutherford equation terms (``MRE_flag``, ``geom_flag``) with a radial Zeff profile (``Zeff``, ``psi_N_Zeff``). +- RDCON - STRIDE-style truncation controls ``qlow``, ``qhigh`` and ``reform_eq_with_psilim``; ``gnstep`` sets the maximum LSODE steps for the Galerkin basis; ``interp_np_res`` refines output in resonant cells. +- MATCH - ``sol_num`` selects which DCON eigenfunction MATCH reconstructs. +- PENTRC, DCON - ``kinetic_format="rates-v1"`` reads a 9-column kinetic file that supplies logLambda, nu_i and nu_e directly (default ``"legacy"`` unchanged). +- PENTRC - Separate energy-integration tolerances ``atol_x``, ``rtol_x`` (tighter by default, so runs take longer); Zeff and ``omega_p`` in the netCDF profile output. +- EQUIL - ``profile_source="integrate"`` builds F and p from the g-file's FF' and p' (recommended only for single-precision g-files such as TokaMaker's); default ``"values"`` is unchanged. +- EQUIL - ``ldp_i`` (TokaMaker i-file) reader accepts F < 0 and an optional ``[psi_axis, psi_lcfs]`` record. +- EQUIL - ``enstep`` sets the maximum field-line integration steps in ``direct_run`` (default 16384, previously fixed at 2048). +- EQUIL - LAR ``sigma_type='wesson'`` current profile; ``chease_ascii`` eq_type. +- VACUUM - New Green's function evaluation (Chance 2007), more accurate at high n; ``use_legacy_greens_function`` restores the old one. +- INSTALL - gfortran (with OpenMP), ifx and ARM macOS support; ``OPENBLASHOME``; ``make deps`` / ``make checkdeps``; build instructions for NERSC, SNU 3DST, PLARE and WSL. +- CI - GitHub Actions builds with gfortran and ifort and runs Solov'ev ideal, kinetic, resistive and energy-decomposition examples on every pull request. +- DOCS - Installation page, sign-conventions reference, resonant coupling, energy-decomposition and pmodb gauge output docs, and the ``docs/tex`` notes (moved from ``tex/``). + +Changes +-------- +- ALL - GPEC is released under the MIT license (``install/LICENSE``), replacing the previous user agreement. +- ALL - Floating-point kinds are standardized on ``REAL(r8)``; LSODE updated to Fortran 2018; duplicate utility, grid, spline and jacobi code consolidated; obsolete files removed. +- ALL - Consistent ``!! WARNING`` terminal messages; GPEC reports CPU and wall time separately for OpenMP runs. +- DCON - Stops with a clear error when ``nn = 0``. +- GPEC - Singular-coupling overlap outputs have correct units (Tesla, or % for ``*_overlap_norm``) and consistent descriptions. +- GPEC - ``singthresh_callen_flag`` forces ``vsingfld_flag`` when coils are present; ``gal_flag`` runs always use the resistive sweet spot. +- RDCON - Free-boundary energies use the DCON vacuum-call pattern and a complex eigen-solve; ``euler.bin`` stores the complex eigenvalues. +- RDCON - ``solver="cholesky"`` switches to LU with a warning when ``rpec_flag=t``. +- RDCON, STRIDE - Modules renamed with ``rdcon_``/``stride_`` prefixes; STRIDE builds without OpenMP and uses ``OMP_SET_MAX_ACTIVE_LEVELS``. +- STRIDE - Banded solve of the Delta' boundary-value problem (faster at high n). +- RDCON - Lower peak memory in the LU Galerkin solve. +- PENTRC - Exact-zero ``omega_E`` knots are replaced by 1e-9 rad/s with a warning. +- PENTRC, STRIDE - LSODE failures in the PENTRC energy and radial integrations and ZVODE failures in STRIDE stop with a nonzero exit status. +- EQUIL - Tabulated Zeff and kinetic profiles use monotone (pchip) fits. +- DOCS - Template and example decks document the new inputs (including ``enstep`` in ``equil.in``), ``rmatch.in`` is fully commented, and the a10 example writes the PENTRC ASCII tables. + +Fixes +------ +- PENTRC - ``fcgl`` recovers all of kappa.xi_perp from dB/B and div(xi_perp), not half. +- PENTRC - Removed a spurious factor of the major radius from the general-aspect-ratio ``omega_D`` (``fgar``, ``tgar``, ``pgar``, ``*mm``); recompute results from these methods. +- PENTRC - Removed spurious factors of ``sqrt(B0)`` from the ``clar`` ``omega_b`` and ``B0`` from its ``omega_D``. +- PENTRC - Fixed OpenMP data races (energy contour flag, ``tpsi``, shared splines), so results no longer depend on the thread count. +- PENTRC - Fixed an out-of-bounds ``bpts`` array and a kinetic-spline length mismatch in ``tpsi``. +- PENTRC - ``rlar`` uses the evaluated ``fnml_f`` in the pitch integrand, so it no longer returns zero torque. +- PENTRC, DCON - Large automatic arrays in the pitch and energy integrations are heap allocated, so OpenMP kinetic runs no longer overflow the worker-thread stacks. +- PENTRC, DCON - The kinetic Euler-Lagrange matrices (``kin_flag``) conjugate the bounce-averaged W vectors, so they act on the actual perturbation instead of its up-down mirror image (up-down asymmetric equilibria only). +- PENTRC - Collision frequencies use a monotone (pchip) fit with temperature and logLambda floors, so they stay positive near a cold edge; non-positive densities and negative temperatures stop with an error. +- PENTRC - Parallel-velocity bounce-root splines use periodic end conditions. +- PENTRC - The ``omega_D`` integrand is no longer padded with ``omega_b`` values near the trapped/passing boundary. +- PENTRC - Fixed anti-damping with ``nutype="small"`` for negative ``omega_E``. +- PENTRC - Radial torque-density export writes the last surface's own value. +- DCON - Fixed the kinetic matrix ``kinfac`` for ``method > 0`` and race conditions in the kinetic matrix loop. +- DCON - Plasmas with q0 >= qa are no longer truncated in the deep core. +- DCON - Fixed out-of-bounds singular-surface access when the q range has no rational surface. +- DCON, RDCON, STRIDE - With ``sas_flag``, ``dmlim`` no longer carries ``qlim`` past ``qhigh``. +- DCON, RDCON, STRIDE - The free-boundary vacuum matrix and eigenvector norm are put in the plasma's Fourier frame, so delta-W no longer depends on the theta=0 origin. +- DCON, RDCON - The singular reduction at a rational surface is kept when a norm reset falls just before it. +- DCON - ``W_v_eigenvalue`` holds the correct real eigenvalues. +- DCON, RDCON, STRIDE - Frobenius exponents at rational surfaces include the trace of the resonant block (no spurious imaginary parts on the Delta' diagonal); Mercier alpha is computed in double precision. +- DCON, RDCON - Fixed Taylor-series slot typos in ``sing_mmat`` (F 6th order; DCON G 3rd order for ``sing_order >= 3``). +- DCON, RDCON, STRIDE, GPEC - Fixed LAPACK and ZVODE argument types, sizes and work arrays (``W_xe_*`` used a stale matrix; ``s_P`` precision; uninitialized ``out_evals`` columns). +- RDCON - Added missing zero initializations, safe deallocations and the correct ``restore_us`` condition in the Galerkin solver; ``vmat.bin`` is written correctly. +- RDCON - Galerkin F is built from the same factored spline as the asymptotics. +- RDCON - Galerkin vacuum matrix keeps its imaginary part; the LU solve is diagonally scaled; a failed Cholesky falls back to LU. +- RDCON - Galerkin elements below a singular surface are sized with ``|n q1|`` instead of ``|q1|``. +- RDCON - Fixed the ``delta_gw.out`` crash with 25 or more singular surfaces and the netCDF failure with ``rpec_flag=f``; ODE integration checks that it reaches psilim. +- RDCON - The netCDF Galerkin cell grid is initialized when the Galerkin output flags are off. +- RDCON, GPEC - Removed spurious shielding at the RPEC control surface; fixed the ``globalsol.bin`` point-count mismatch that blocked field penetration at rational surfaces. +- RMATCH - Removed spurious "imap error" messages with the default ``inps_type="inps"``. +- RMATCH - Fixed an out-of-bounds read in ``match_solution`` when ``msing`` is less than RDCON's surface count, an uninitialized ``delta`` and an undefined-behavior logical in ``inps_ua``. +- GPEC - Control-surface normals ``R_n`` and ``z_n`` use the current geometry derivatives. +- GPEC - ``xbnormal`` normalizes each surface with its own major radius. +- GPEC - pmodb netCDF functions (``b_eul_fun``, ``b_lag_fun``, ``Bdivxi_perp_fun``, ``Bkappaxi_perp_fun``) have correct imaginary parts and include the theta=1 row. +- GPEC - pmodb divergence output honours ``reg_flag=f``; the pmodb diagnostic projection uses the current Jacobian. +- GPEC, PENTRC - pmodb ``Bdivxi_perp`` and ``Bkappaxi_perp`` are formed directly from xi_perp instead of subtracting a Fourier-truncated parallel displacement; PENTRC's ``pmodb_file`` reader uses them instead of ``b_lag``. +- GPEC - Fixed a use-after-free of the local coupling vectors; exported singular vectors use the adjoint area weighting; ``C_l_x_out`` holds the local coupling matrix. +- GPEC - ``q1`` at resonant surfaces is computed instead of set to 0; fixed out-of-bounds access with a single rational surface, a netCDF buffer error when ``coil_flag`` is set without coils, the ``gpout_xbnormal`` array length and long_name attributes written to the wrong files. +- STRIDE - Fixed an OpenMP data race in ZLANGE, stack overflow at large ``mpert``, an out-of-bounds write and a double allocation. +- STRIDE - The interval from the axis to the first rational surface is always split into at least 8 pieces. +- VACUUM - Fixed an off-by-one between ``mthsurf`` and ``mthvac``, a divide by zero in ``defglo`` and out-of-bounds writes; observer-grid coordinates are double precision. +- COIL - Fixed undefined spline memory and an out-of-range write in ``field.F``; zero-amplitude shifts no longer distort the coils. +- EQUIL - Fixed spline "extrap" end derivatives on clustered grids and real-valued classic complex spline coefficients; ``spline_roots`` finds roots near cell endpoints. +- EQUIL - Fixed uninitialized ``ro``/``zo`` for EFIT input, ``eq_type="dump"`` for direct equilibria, an out-of-bounds write for g-files without a q column, a crash re-splining inverse equilibria, and the CHEASE ASCII reader; the field-line integration failure message is printed in full. +- EQUIL - Direct equilibria are now sampled exactly at the common straight-fieldline angles instead of interpolated from each surface's own field-line trace. +- EQUIL - ``grid_type="original"`` honours ``psilow`` and ``psihigh``. +- INSTALL - The embedded version string uses release tags only; bundled dependencies honour ``FC``/``CC`` and embed an rpath; xdraw compiles with strict C99/C11 compilers; added missing build dependencies. + + GPEC v1.5.7 =========== diff --git a/docs/sign_conventions.rst b/docs/sign_conventions.rst index 40dfccb45..25290474a 100644 --- a/docs/sign_conventions.rst +++ b/docs/sign_conventions.rst @@ -57,8 +57,8 @@ Toroidal Coordinate :math:`\zeta` and :math:`\phi` Working Coordinate Options -------------------------- -Controlled by ``jac_type`` in ``equil.in``. The Jacobian is -:math:`J \propto B_p^{p_{bp}} \, B^{p_b} \, R^{-p_r}`: +Controlled by ``jac_type`` in ``equil.in``. The Jacobian is set by +:math:`\mathbf{B}\cdot\nabla\theta \propto 1/J \propto B_p^{p_{bp}} \, B^{p_b} \, R^{-p_r}`: .. list-table:: :header-rows: 1 @@ -84,6 +84,14 @@ Controlled by ``jac_type`` in ``equil.in``. The Jacobian is - 1 - 0 - 0 + * - Park + - 0 + - 1 + - 0 + * - Other (powers set explicitly in ``equil.in``) + - + - + - See also the coordinate discussion in :doc:`outputs`. diff --git a/docs/site.rst b/docs/site.rst index 59e5c59c0..7265bca24 100644 --- a/docs/site.rst +++ b/docs/site.rst @@ -11,6 +11,7 @@ GPEC documentation website site map contact.rst outputs.rst sign_conventions.rst + pmodb_gauge.rst dcon_readme.rst examples.rst releases.rst diff --git a/docs/source_documentation.rst b/docs/source_documentation.rst index 038d69188..5c3767513 100644 --- a/docs/source_documentation.rst +++ b/docs/source_documentation.rst @@ -1,6 +1,6 @@ .. _source_documentation: -The recommended post processing for GPEC is now to use the devoted `OMFIT module `. +The recommended post processing for GPEC is now to use the devoted `OMFIT module `_. The deprecated code below has been ported there and improved, but is included here for reference. diff --git a/gpec/gpdiag.f b/gpec/gpdiag.f index 9f1c7e3c9..57b503048 100644 --- a/gpec/gpdiag.f +++ b/gpec/gpdiag.f @@ -1090,8 +1090,8 @@ SUBROUTINE gpdiag_pmodb(egnum,xspmn) $ mthsurf,curv2_mn) eulbparmns(istep,:)=eulbpar_mn ! p'xsp_mn+bxcdeltamn - lagbparmns(istep,:)=lagbpar_mn ! with parallel components - llagbparmns(istep,:)=llagbpar_mn ! without parallel components + lagbparmns(istep,:)=lagbpar_mn ! xi^zeta=0 gauge (pmodb lagb) + llagbparmns(istep,:)=llagbpar_mn ! xi_perp only divxmns(istep,:)=-divx_mn ! Bdivxprp curvmns(istep,:)=-curv_mn ! Bkxprp divx2mns(istep,:)=divx2_mn ! jb-weighting diff --git a/gpec/gpout.f b/gpec/gpout.f index 9f4d5cbb1..4250fa70d 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -3053,6 +3053,9 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) eulbparfun(istep,itheta)= $ chi1*(bvt_fun(itheta)+sq%f(4)*bvz_fun(itheta)) $ /(rzphi%f(4)*eqfun%f(1)) +c lagb: Lagrangian dB of xi_perp + lambda B with xi^zeta=0 +c in the working coordinates, not the xi_perp-only value; see +c docs/pmodb_gauge.rst. lagbparfun(istep,itheta)= $ eulbparfun(istep,itheta)+ $ xsp_fun(itheta)*eqfun%fx(1)+ @@ -3491,6 +3494,7 @@ SUBROUTINE gpout_xbnormal(egnum,xspmn,spots,interpspot,nspot) eta=twopi*(theta(itheta)+rzphi%f(2)) rs(istep,itheta)=ro+rfac*COS(eta) zs(istep,itheta)=zo+rfac*SIN(eta) + r(itheta)=rs(istep,itheta) jac=rzphi%f(4) jacs(itheta)=jac w(1,1)=(1+rzphi%fy(2))*twopi**2*rfac*r(itheta)/jac diff --git a/input/dcon.in b/input/dcon.in index 8db588137..57ec14a12 100644 --- a/input/dcon.in +++ b/input/dcon.in @@ -52,7 +52,7 @@ ahb_flag=f ! Output normal magnetic field eigenvalues and eigenfunctons at plasma-vacuum interface (must be false for GPEC) msol_ahb=1 ! Number of eigenfunctions output by ahb_flag=t ? - mthsurf0=1 ! Linear multiplier on number of boundary points used to display surface eigenfunctions for ahgb_flag=t + mthsurf0=1 ! Linear multiplier on number of boundary points used to display surface eigenfunctions for ahb_flag=t bin_euler=t ! Output M psi-by-M euler-lagrange solutions to binary file euler.bin bin_vac=f ! Output deprecated binary vacuum.bin containing vacuum energy matrices diff --git a/input/equil.in b/input/equil.in index 86f3b7f3c..6189d78c0 100644 --- a/input/equil.in +++ b/input/equil.in @@ -2,7 +2,7 @@ eq_type="efit" ! Type of the input 2D equilibrium file. Accepts efit, chease, fluxgrid, transp, jsolver, lar, sol, etc. eq_filename="mypath" ! Path to input file (geqdsk, etc.) Must specified (no default). - jac_type="hamada" ! Working coordinate system for all DCON and GPEC calculations. Overrides individual powers. Accepts hamada, pest, boozer, equal_arc + jac_type="hamada" ! Working coordinate system for all DCON and GPEC calculations. Overrides individual powers. Accepts hamada, pest, boozer, equal_arc, park, other power_bp=0 ! del.B ~ B_p**power_bp * B**power_b / R**power_r power_b=0 ! del.B ~ B_p**power_bp * B**power_b / R**power_r power_r=0 ! del.B ~ B_p**power_bp * B**power_b / R**power_r diff --git a/input/slayer.in b/input/slayer.in index 06e3c06e1..6e03c57a9 100644 --- a/input/slayer.in +++ b/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)&1 | tee make.log + make 2>&1 | tee make.log where the logging part is merely a useful suggestion. Appropriate compiler details for the major fusion clusters on which public installations are maintained is recorded in the $GPECHOME/install/DEFAULTS.inc file annotations. \ No newline at end of file diff --git a/install/RULES.inc b/install/RULES.inc index 5cf187176..392ed4818 100644 --- a/install/RULES.inc +++ b/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 diff --git a/install/TARGETS.inc b/install/TARGETS.inc index 6046fb414..8c9a5f1d0 100644 --- a/install/TARGETS.inc +++ b/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) diff --git a/install/gen_deps.py b/install/gen_deps.py index 44700f1b6..b60cbfb78 100644 --- a/install/gen_deps.py +++ b/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"], } diff --git a/pentrc/torque.F90 b/pentrc/torque.F90 index 7e594bd71..e9354632c 100644 --- a/pentrc/torque.F90 +++ b/pentrc/torque.F90 @@ -786,8 +786,8 @@ subroutine tpsi(tpsi_var,psi,n,l,zi,mi,wdfac,divxfac,electron, & enddo ! build complete action mxm matrices W_X, W_Y, W_Z - wmmt(1,:) = bwspl(1)%fsi(bwspl(1)%mx,:) - wemt(1,:) = bwspl(2)%fsi(bwspl(2)%mx,:) + wmmt(1,:) = conjg(bwspl(1)%fsi(bwspl(1)%mx,:)) + wemt(1,:) = conjg(bwspl(2)%fsi(bwspl(2)%mx,:)) wxmt = matmul(wmmt,xmat) wymt = matmul(wmmt,3*smat+ymat)-2.0*matmul(wemt,smat) wzmt = matmul(wmmt,3*tmat+zmat)-2.0*matmul(wemt,tmat) diff --git a/slayer/delta.f b/slayer/delta.f index 8b5c085f2..7eb4500a3 100644 --- a/slayer/delta.f +++ b/slayer/delta.f @@ -1,32 +1,81 @@ MODULE delta_mod +c----------------------------------------------------------------------- +c delta_mod: Riccati-based tearing-mode layer Delta solvers. +c +c Contains three Riccati formulations: +c riccati - standard formulation (non-stiff, lsode mf=10) +c riccati_del_s - del_s formulation (stiff, lsode mf=21) +c riccati_f - Fitzpatrick P_perp/P_tor formulation (mf=21) +c +c Each solver integrates a Riccati ODE for W(x) from a large-|x| +c asymptotic boundary condition inward to x ~ 0, then extracts +c Delta = pi / W'(0). +c +c Associated ODE subroutines (w_der, w_der_del_s, w_der_f) and +c Jacobian subroutines (jac_del_s, jac_f) follow below. +c +c Module-level flags: +c riccati_out - write W(x) profile to binary + text file +c parflow_flag - include parallel electron flow terms in w_der +c PeOhmOnly_flag - retain only Pe-Ohm coupling in w_der +c----------------------------------------------------------------------- USE sglobal_mod + USE ieee_arithmetic, ONLY: ieee_value, ieee_quiet_nan IMPLICIT NONE - LOGICAL :: riccati_out,parflow_flag,PeOhmOnly_flag +c --- module-level control flags + LOGICAL :: riccati_out ! write W(x) profile to binary + text + LOGICAL :: parflow_flag ! enable parallel-flow terms in w_der + LOGICAL :: PeOhmOnly_flag ! Pe-Ohm-only coupling in w_der CONTAINS c----------------------------------------------------------------------- -c calculate delta based on riccati w_der formulation. +c riccati: compute Delta via Riccati integration of w_der. +c +c Integrates W(x) from x_start inward to x_min using lsode +c (non-stiff, mf=10). Optionally applies the layfac singularity +c guard when Q is near Q_e. Returns Delta = pi / W'(x_min). +c +c If riccati_out = .TRUE., writes the W(x) profile to +c slayer_riccati_profile_n.{bin,out}. +c c----------------------------------------------------------------------- FUNCTION riccati(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe, $ iinQ,inx,iny) - REAL(r8),INTENT(IN) :: inQ,inQ_e,inQ_i,inpr,inpe,inc_beta,inds - REAL(r8),INTENT(IN) :: intau - REAL(r8),INTENT(IN),OPTIONAL :: iinQ,inx - COMPLEX(r8), INTENT(IN), OPTIONAL :: iny +c --- input arguments: physical parameters for this surface + REAL(r8),INTENT(IN) :: inQ ! real part of Q + REAL(r8),INTENT(IN) :: inQ_e ! electron diamagnetic freq + REAL(r8),INTENT(IN) :: inQ_i ! ion diamagnetic freq + REAL(r8),INTENT(IN) :: inpr ! Prandtl number + REAL(r8),INTENT(IN) :: inpe ! electron Prandtl number + REAL(r8),INTENT(IN) :: inc_beta ! c_beta parameter + REAL(r8),INTENT(IN) :: inds ! ds, normalised rho_s + REAL(r8),INTENT(IN) :: intau ! ion-to-electron temp ratio +c --- optional arguments + REAL(r8),INTENT(IN),OPTIONAL :: iinQ ! imaginary part of Q + REAL(r8),INTENT(IN),OPTIONAL :: inx ! override starting x + COMPLEX(r8),INTENT(IN),OPTIONAL :: iny ! override starting W +c --- function result COMPLEX(r8) :: riccati - INTEGER :: istep,neq,itol,itask,istate,liw,lrw,iopt,mf - - REAL(r8) :: xintv,x,xout,rtol,jac,xmin - COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: y,dy - - INTEGER, DIMENSION(:), ALLOCATABLE :: iwork - REAL(r8), DIMENSION(:), ALLOCATABLE :: xfac,atol,rwork +c --- lsode solver control + INTEGER :: istep ! integration step counter + INTEGER :: neq ! number of equations (=2) + INTEGER :: itol,itask ! lsode tolerance/task flags + INTEGER :: istate,iopt,mf ! lsode state/option/method flags + INTEGER :: liw,lrw ! lsode work array sizes + REAL(r8) :: x,xout ! current and target x + REAL(r8) :: xmin ! inner integration bound + REAL(r8) :: rtol ! relative tolerance +c --- work arrays + COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: y,dy ! W and dW + INTEGER, DIMENSION(:), ALLOCATABLE :: iwork ! lsode int work + REAL(r8), DIMENSION(:), ALLOCATABLE :: atol,rwork ! lsode real work +c --- copy input arguments to module-level globals for w_der Q=inQ IF(present(iinQ)) Q=inQ+ifac*iinQ Q_e=inQ_e @@ -37,75 +86,430 @@ FUNCTION riccati(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe, ds=inds tau=intau +c --- singularity guard: displace Q away from Q_e when too close IF ((layfac>0).AND.(ABS(Q-Q_e)xout) istep=istep+1 CALL lsode(w_der,neq,y,x,xout,itol,rtol,atol, - $ itask,istate,iopt,rwork,lrw,iwork,liw,jac,mf) + $ itask,istate,iopt,rwork,lrw,iwork,liw,jac_dummy,mf) WRITE(bin_unit)REAL(x,4),REAL(REAL(y),4),REAL(AIMAG(y),4) WRITE(out2_unit,'(1x,3(es17.8e3))') x,REAL(y),AIMAG(y) ENDDO CLOSE(bin_unit) CLOSE(out2_unit) ELSE +c single-shot integration to x_min istep = 1 itask = 1 CALL lsode(w_der,neq,y,x,xout,itol,rtol,atol, - $ itask,istate,iopt,rwork,lrw,iwork,liw,jac,mf) - + $ itask,istate,iopt,rwork,lrw,iwork,liw,jac_dummy,mf) ENDIF - ! w=0 when Q=Q_e. Why? - +c --- extract Delta from final W derivative at x_min +c NOTE: W -> 0 when Q -> Q_e (see layfac guard above). CALL w_der(neq,x,y,dy) riccati=pi/dy(1) DEALLOCATE(atol,y,dy,iwork,rwork) END FUNCTION riccati c----------------------------------------------------------------------- -c riccati integration. +c riccati_del_s: compute Delta via the del_s Riccati formulation. +c +c Uses a stiff solver (lsode mf=21) with user-supplied Jacobian +c (jac_del_s). Integrates W(q) from large q inward to q_min. +c Returns Delta = -(pi / sqrt(1+1/tau)) * W'(q_min). +c +c Module-level variables tau, D_norm, and P_tor must be valid +c before entry; callers are responsible for setting them. +c Q_e, Q_i, and P_perp are set from arguments inside this +c function. Q, c_beta, and d_beta are not referenced here. +c Note: slayer.f passes Q_e_arr (not Q_arr) as the first +c argument -- this is intentional. +c----------------------------------------------------------------------- + FUNCTION riccati_del_s(inQ_e,inQ_i,inpr,inx,iny) + +c --- input arguments + REAL(r8),INTENT(IN) :: inQ_e ! electron diamagnetic freq + REAL(r8),INTENT(IN) :: inQ_i ! ion diamagnetic freq + REAL(r8),INTENT(IN) :: inpr ! mapped to P_perp (see below) +c --- optional arguments + REAL(r8),INTENT(IN) :: inx ! starting q for integration + COMPLEX(r8),INTENT(IN),OPTIONAL :: iny ! override starting W +c --- function result + COMPLEX(r8) :: riccati_del_s + +c --- lsode solver control + INTEGER :: istep ! integration step counter + INTEGER :: neq ! number of equations (=2) + INTEGER :: itol,itask ! lsode tolerance/task flags + INTEGER :: istate,iopt,mf ! lsode state/option/method flags + INTEGER :: liw,lrw ! lsode work array sizes + REAL(r8) :: x ! secondary x variable (set from inx) + REAL(r8) :: xout ! target integration endpoint + REAL(r8) :: xmin ! inner integration bound + REAL(r8) :: rtol ! relative tolerance + REAL(r8) :: my_q ! integration variable (large -> small) + REAL(r8) :: P_hat ! normalized P_perp + REAL(r8) :: alpha ! boundary condition coefficient +c --- work arrays + COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: W,dW_dq ! W and dW/dq + INTEGER, DIMENSION(:), ALLOCATABLE :: iwork ! lsode int work + REAL(r8), DIMENSION(:), ALLOCATABLE :: atol,rwork ! lsode real work + +c --- configure lsode: stiff BDF method with user Jacobian (mf=21) + neq = 2 + itol = 2 + rtol = 1e-10 + ALLOCATE(atol(neq),W(1),dW_dq(1)) + atol(:) = 1e-10 + itask = 2 + istate = 1 + iopt = 1 ! enable optional inputs (iwork(6)) + mf = 21 ! stiff, user-supplied Jacobian (jac_del_s) + liw = 20*2 + lrw = 22+9*neq+neq**2 ! stiff work array size + ALLOCATE(iwork(liw+neq),rwork(lrw)) + +c --- set maximum internal steps + iwork=0 + iwork(6)=50000 ! MXSTEP: max internal steps per call + rwork=0 + +c --- set starting integration point + my_q=inx ! start backwards integration at large q + xmin=1e-5 + x=inx + xout=xmin + +c --- copy input arguments to module-level globals for w_der_del_s + Q_e = inQ_e + Q_i = inQ_i + P_perp = inpr + + P_hat = P_perp / D_norm**6.0 + +c --- asymptotic boundary condition at large q + alpha = (P_hat/(1+1/tau))**0.5 + W(1) = -alpha*my_q**2 - 0.5 + IF(present(iny)) W(1)=iny + +c --- integrate W(q) from q_start inward to q_min via lsode + IF (riccati_out) THEN +c profile output: step-by-step integration with file writes + istep = 1 + itask = 2 + OPEN(UNIT=bin_unit,FILE='slayer_riccati_profile_n'// + $ TRIM(sn_str)//'.bin',STATUS='UNKNOWN', + $ POSITION='REWIND',FORM='UNFORMATTED') + + OPEN(UNIT=out2_unit,FILE='slayer_riccati_profile_n'// + $ TRIM(sn_str)//'.out',STATUS='UNKNOWN') + WRITE(out2_unit,'(1x,3(a17))'),"x","RE(y)","IM(y)" + DO WHILE (my_q>xout) + istep=istep+1 + CALL lsode(w_der_del_s,neq,W,my_q,xout,itol,rtol,atol, + $ itask,istate,iopt,rwork,lrw,iwork,liw,jac_del_s,mf) + WRITE(bin_unit)REAL(my_q,4),REAL(REAL(W),4),REAL(AIMAG(W),4) + WRITE(out2_unit,'(1x,3(es17.8e3))')my_q,REAL(W),AIMAG(W) + ENDDO + CLOSE(bin_unit) + CLOSE(out2_unit) + ELSE +c single-shot integration to q_min + istep = 1 + itask = 1 + CALL lsode(w_der_del_s,neq,W,my_q,xout,itol,rtol,atol, + $ itask,istate,iopt,rwork,lrw,iwork,liw,jac_del_s,mf) + ENDIF + +c --- extract Delta from final W derivative at q_min + CALL w_der_del_s(neq,my_q,W,dW_dq) + riccati_del_s=-( pi/((1+1/tau)**0.5) )*dW_dq(1) + DEALLOCATE(atol,W,dW_dq,iwork,rwork) + + END FUNCTION riccati_del_s +c----------------------------------------------------------------------- +c jac_dummy: placeholder Jacobian for lsode calls whose method flag +c does not use a user-supplied Jacobian (mf=10, mf=22). lsode never +c calls it; it exists so that every lsode call in this module passes +c a subroutine in the JAC position. +c----------------------------------------------------------------------- + SUBROUTINE jac_dummy(neq, x, y, ml, mu, pd, nrpd) + INTEGER, INTENT(IN) :: neq, ml, mu, nrpd + REAL(r8), INTENT(IN) :: x + COMPLEX(r8), DIMENSION(neq), INTENT(IN) :: y + COMPLEX(r8), DIMENSION(nrpd,neq), INTENT(INOUT) :: pd + END SUBROUTINE jac_dummy +c----------------------------------------------------------------------- +c jacobian for riccati_del_s(): pd = dF/dW for stiff lsode. +c----------------------------------------------------------------------- + SUBROUTINE jac_del_s(neq, my_q, W, ml, mu, pd, nrpd) + INTEGER, INTENT(IN) :: neq, ml, mu, nrpd + REAL(r8), INTENT(IN) :: my_q + COMPLEX(r8), DIMENSION(neq), INTENT(IN) :: W + COMPLEX(r8), DIMENSION(nrpd,neq), INTENT(INOUT) :: pd + pd(1,1) = 1.0/my_q - 2.0d0*W(1)/my_q + END SUBROUTINE jac_del_s +c----------------------------------------------------------------------- +c w_der_del_s: ODE right-hand side dW/dq for riccati_del_s. +c Implements the del_s dispersion relation using normalised +c quantities Q_hat, P_perp_hat, P_tor_hat. +c----------------------------------------------------------------------- + SUBROUTINE w_der_del_s(neq,my_q,W,dW_dq) + + INTEGER, INTENT(IN) :: neq + REAL(r8), INTENT(IN) :: my_q + COMPLEX(r8), DIMENSION(neq), INTENT(IN) :: W + COMPLEX(r8), DIMENSION(neq), INTENT(OUT) :: dW_dq + REAL(r8) :: Q_hat, P_tor_hat, P_perp_hat + COMPLEX(r8) :: E,F + +c --- normalise physical quantities + Q_hat = (Q_e*(1+tau)/tau) / D_norm**4.0 + P_perp_hat = P_perp / D_norm**6.0 + P_tor_hat = P_tor / D_norm**6.0 +c --- build the E and F dispersion coefficients + E = (-(Q_hat**2)/(1+1/tau)) - ifac*Q_hat*(P_perp_hat+ + $ P_tor_hat)*(my_q**2) + P_perp_hat*P_tor_hat*(my_q**4) + F = P_perp_hat - ifac*Q_hat + (1+1/tau)*P_tor_hat*my_q**2 + +c --- Riccati ODE: dW/dq = W/q - W^2/q + q*E/F + dW_dq(1)=W(1)/my_q - (W(1)**2)/my_q + (my_q*E)/F + RETURN + END SUBROUTINE w_der_del_s +c----------------------------------------------------------------------- +c riccati_f: compute Delta via Fitzpatrick P_perp / P_tor +c Riccati formulation. +c +c Uses a stiff solver (lsode mf=21) with user-supplied Jacobian +c (jac_f). Boundary conditions are set analytically from the +c large-p asymptotic behaviour; the branch depends on whether +c D_norm^2 exceeds iota_e * P_perp / P_tor^(2/3). +c +c----------------------------------------------------------------------- + FUNCTION riccati_f() + +c --- function result + COMPLEX(r8) :: riccati_f + +c --- lsode solver control + INTEGER :: istep ! integration step counter + INTEGER, PARAMETER :: neq = 2 ! number of equations + INTEGER :: itol,itask ! lsode tolerance/task flags + INTEGER :: istate,iopt,mf ! lsode state/option/method flags + INTEGER, PARAMETER :: liw = 42 ! 20*2 + neq + INTEGER, PARAMETER :: lrw = 44 ! 22 + 9*neq + neq**2 + REAL(r8) :: xout ! target integration endpoint + REAL(r8) :: xmin ! inner integration bound + REAL(r8) :: rtol ! relative tolerance + REAL(r8) :: my_p ! integration variable p (large -> small) + REAL(r8) :: bk ! asymptotic coefficient b_k +c --- boundary-condition intermediates + COMPLEX(r8) :: ak ! asymptotic coefficient a_k + COMPLEX(r8) :: ck ! asymptotic coefficient c_k + COMPLEX(r8) :: xk ! asymptotic coefficient x_k + COMPLEX(r8) :: W_bound ! boundary value for W(p_start) +c --- work arrays (fixed-size, stack-allocated to avoid heap overhead) + COMPLEX(r8) :: W(1), dWdp(1) ! W and dW/dp + REAL(r8) :: atol(neq) ! absolute tolerance + INTEGER :: iwork(liw) ! lsode int work + REAL(r8) :: rwork(lrw) ! lsode real work + +c --- configure lsode: stiff BDF method with user Jacobian (mf=21) + itol = 2 + rtol = 1.0d-10 + atol(:) = 1.0d-10 + itask = 2 + istate = 1 + iopt = 1 ! enable optional inputs (iwork(6)) + mf = 21 ! stiff, user-supplied Jacobian (jac_f) + +c --- set maximum internal steps (only iwork(5:6) matter for lsode) + iwork(:) = 0 + iwork(6) = 50000 ! MXSTEP: max internal steps per call + rwork(:) = 0.0d0 + + xmin=1e-6 + xout=xmin + +c --- compute starting p and W boundary condition +c Branch on asymptotic regime: D_norm^2 vs iota_e*P_perp/P_tor^(2/3) +c --- branch 1: D_norm^2 > iota_e * P_perp / P_tor^(2/3) +c large-D_norm regime: p scales with (P_tor*D_norm^2/(iota_e*...)) + IF (D_norm**2 > iota_e*P_perp / P_tor**(2.0d0/3.0d0)) THEN + my_p = ( (P_tor*D_norm**2)/(iota_e*P_tor*P_perp) )**0.25d0 + my_p = MAX(my_p, 6.0d0) + + ak = -(g_tmp + ifac*Q_e) + bk = (iota_e*P_perp*P_tor)/(P_tor*D_norm**2) + + ck = bk*(1 + (g_tmp+ifac*Q_i) + $ *((P_tor+P_perp)/(P_tor*P_perp)) + $ - (P_perp + (g_tmp+ifac*Q_i)*D_norm**2) + $ *(iota_e/(P_tor*D_norm**2))) + + xk = (ck - SQRT(bk)*(1 - SQRT(bk)*ak)) / (2.0d0*SQRT(bk)) + + W_bound = xk - SQRT(bk)*my_p + ELSE +c --- branch 2: D_norm^2 <= iota_e * P_perp / P_tor^(2/3) +c small-D_norm regime: p scales with 1/P_tor^(1/6) + my_p = 1.0d0 / P_tor**(1.0d0/6.0d0) + my_p = MAX(my_p, 6.0d0) + + ak = -(g_tmp + ifac*Q_e) + bk = P_tor + ck = -ifac*(Q_e - Q_i)*(P_tor/P_perp) + (g_tmp + ifac*Q_i) + xk = (ak*bk - ck)/(2.0d0*SQRT(bk)) + + W_bound = -1.0d0 + xk*my_p - SQRT(bk)*my_p**3 + END IF + + W(1) = W_bound + +c --- integrate W(p) from p_start inward to p_min via lsode + IF (riccati_out) THEN +c profile output: step-by-step integration with file writes + istep = 1 + itask = 2 + OPEN(UNIT=bin_unit,FILE='slayer_riccati_profile_n'// + $ TRIM(sn_str)//'.bin',STATUS='UNKNOWN', + $ POSITION='REWIND',FORM='UNFORMATTED') + + OPEN(UNIT=out2_unit,FILE='slayer_riccati_profile_n'// + $ TRIM(sn_str)//'.out',STATUS='UNKNOWN') + WRITE(out2_unit,'(1x,3(a17))'),"x","RE(y)","IM(y)" + DO WHILE (my_p>xout) + istep=istep+1 + CALL lsode(w_der_f,neq,W,my_p,xout,itol,rtol,atol, + $ itask,istate,iopt,rwork,lrw,iwork,liw,jac_f,mf) + WRITE(bin_unit)REAL(my_p,4),REAL(REAL(W),4), + $ REAL(AIMAG(W),4) + WRITE(out2_unit,'(1x,3(es17.8e3))')my_p,REAL(W),AIMAG(W) + ENDDO + CLOSE(bin_unit) + CLOSE(out2_unit) + ELSE +c single-shot integration to p_min + istep = 1 + itask = 1 + CALL lsode(w_der_f,neq,W,my_p,xout,itol,rtol,atol, + $ itask,istate,iopt,rwork,lrw,iwork,liw,jac_f,mf) + ENDIF + +c --- extract Delta from final W derivative at p_min + CALL w_der_f(neq,my_p,W,dWdp) + riccati_f = pi / dWdp(1) +c A failed integration (istate < 0) leaves W short of p_min. Return +c NaN so that scans and the root search skip the point instead of +c using a Delta built from an unconverged W. + IF (istate < 0) riccati_f = CMPLX( + $ ieee_value(1.0_r8, ieee_quiet_nan), + $ ieee_value(1.0_r8, ieee_quiet_nan), KIND=r8) + + END FUNCTION riccati_f +c----------------------------------------------------------------------- +c jacobian for riccati_f(): pd = dF/dW for stiff lsode. +c----------------------------------------------------------------------- + SUBROUTINE jac_f(neq, my_p, W, ml, mu, pd, nrpd) + INTEGER, INTENT(IN) :: neq, ml, mu, nrpd + REAL(r8), INTENT(IN) :: my_p + COMPLEX(r8) :: fA_p, denom + REAL(r8) :: p2 + COMPLEX(r8), DIMENSION(neq), INTENT(IN) :: W + COMPLEX(r8), DIMENSION(nrpd,neq), INTENT(INOUT) :: pd + + p2 = my_p * my_p + denom = g_tmp + ifac*Q_e + p2 + fA_p = (denom - 2.0d0*p2) / denom + + pd(1,1) = (-fA_p/my_p) - (2.0d0*W(1))/my_p + END SUBROUTINE jac_f +c----------------------------------------------------------------------- +c w_der_f: ODE right-hand side dW/dp for riccati_f. +c Implements the Fitzpatrick P_perp / P_tor dispersion relation. +c Coefficients fA, fB, fC are evaluated at the current p. +c----------------------------------------------------------------------- + SUBROUTINE w_der_f(neq,my_p,W,dWdp) + + INTEGER, INTENT(IN) :: neq + REAL(r8), INTENT(IN) :: my_p + COMPLEX(r8), DIMENSION(neq), INTENT(IN) :: W + COMPLEX(r8), DIMENSION(neq), INTENT(OUT) :: dWdp + COMPLEX(r8) :: fA, fA_prime, fB, fC + COMPLEX(r8) :: denom ! cached g_tmp + i*Q_e + p^2 + REAL(r8) :: p2, p4 ! cached p^2, p^4 + REAL(r8) :: D2 ! cached D_norm^2 + +c cache powers and shared denominator + p2 = my_p * my_p + p4 = p2 * p2 + D2 = D_norm * D_norm + denom = g_tmp + ifac*Q_e + p2 + +c evaluate coefficients at the current p + fA = p2 / denom + fA_prime = (denom - 2.0d0*p2) / denom ! (g+iQe - p^2)/(g+iQe + p^2) + + fB = g_tmp*(g_tmp + ifac*Q_i) + $ + (g_tmp + ifac*Q_i)*(P_perp + P_tor)*p2 + $ + (P_perp*P_tor)*p4 + + fC = g_tmp + ifac*Q_e + $ + (P_perp + (g_tmp + ifac*Q_i)*D2)*p2 + $ + (1.0d0/iota_e)*P_tor*D2*p4 + + dWdp(1) = -(fA_prime/my_p)*W(1) - W(1)*W(1)/my_p + $ + (fB/(fA*fC))*(p2*my_p) + + RETURN + END SUBROUTINE w_der_f +c----------------------------------------------------------------------- +c W derivative for riccati() c----------------------------------------------------------------------- SUBROUTINE w_der(neq,x,y,dy) @@ -122,7 +526,7 @@ SUBROUTINE w_der(neq,x,y,dy) COMPLEX(r8) :: C2p COMPLEX(r8) :: A1 COMPLEX(r8) :: A2 - COMPLEX(r8), PARAMETER :: ifac=(0,1) + !COMPLEX(r8), PARAMETER :: ifac=(0,1) IF (parflow_flag) THEN C1=((1 + tau)*x**2*pe* @@ -586,7 +990,7 @@ SUBROUTINE w_der_temp(neq,x,y,dy) REAL(r8), INTENT(IN) :: x COMPLEX(r8), DIMENSION(neq), INTENT(IN) :: y COMPLEX(r8), DIMENSION(neq), INTENT(OUT) :: dy - COMPLEX(r8), PARAMETER :: ifac=(0,1) + !COMPLEX(r8), PARAMETER :: ifac=(0,1) dy(1)=(2.0*x/(ifac*(Q-Q_e)+x**2.0)-1.0/x)*y(1)-y(1)*y(1)/x $ +x*(ifac*(Q-Q_e)+x**2.0) @@ -614,7 +1018,7 @@ FUNCTION directred(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe) INTEGER :: istep,neq,itol,itask,istate,liw,lrw,iopt,mf - REAL(r8) :: xintv,x,xout,rtol,jac,xmax + REAL(r8) :: xintv,x,xout,rtol,xmax COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: y,dy INTEGER, DIMENSION(:), ALLOCATABLE :: iwork @@ -654,7 +1058,7 @@ FUNCTION directred(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe) DO istep=istep+1 CALL lsode(phi_der,neq,y,x,xout,itol,rtol,atol, - $ itask,istate,iopt,rwork,lrw,iwork,liw,jac,mf) + $ itask,istate,iopt,rwork,lrw,iwork,liw,jac_dummy,mf) xout=xout+xintv IF (xmax-xout 0 flags an exactly singular matrix, whose zero + ! pivot already makes the product zero. + lu_mat = matk + CALL zgetrf(nk, nk, lu_mat, nk, ipiv, info) + IF (info < 0) THEN + detk = CMPLX(0.0_r8, 0.0_r8, KIND=r8) + status = -1 + ELSE + detk = CMPLX(1.0_r8, 0.0_r8, KIND=r8) + DO i = 1, nk + detk = detk * lu_mat(i,i) + IF (ipiv(i) /= i) detk = -detk + END DO + END IF + END IF + + END SELECT + RETURN + END SUBROUTINE calc_determinant +c----------------------------------------------------------------------- +c subprogram 7. dispersion_det. +c Coupled dispersion determinant for the first msing_max rational +c surfaces: +c +c det( Dp - diag(Delta_crit_k) - diag(S_k^(1/3)*Delta_k(g_k)) ) +c +c Dp is the r_s-referenced Deltaprime matrix (complex), of which +c the leading msing_max x msing_max block is used, so a subset of +c the surfaces STRIDE found can be coupled. Delta_crit is a local +c layer offset and enters the diagonal only. Each layer response +c is evaluated at that surface's own normalised eigenvalue g_k. +c For msing_max = 1 this is the single-surface matching condition. +c----------------------------------------------------------------------- + FUNCTION dispersion_det(g_in,n_k,sl_in,msing_max) + +c --- arguments + COMPLEX(r8), INTENT(IN) :: g_in ! complex growth rate + INTEGER, INTENT(IN) :: n_k ! number of surfaces + INTEGER, INTENT(IN) :: msing_max ! surfaces to include + TYPE(slayer_inputs_type), INTENT(IN) :: sl_in +c --- function result and locals + COMPLEX(r8) :: dispersion_det ! returned determinant + COMPLEX(r8) :: det_val ! intermediate determinant + COMPLEX(r8), ALLOCATABLE :: result_matrix(:,:) ! Dp - layer terms + INTEGER :: k ! surface loop index + INTEGER :: det_status ! calc_determinant error status + + IF (msing_max < 1 .OR. msing_max > n_k) THEN + WRITE(*,*) 'ERROR: dispersion_det needs 1 <= msing_max <= ', + $ 'number of surfaces; got msing_max=',msing_max, + $ ' with',n_k,' surfaces' + STOP 'dispersion_det: msing_max out of range' + END IF + IF (.NOT. ALLOCATED(sl_in%dp_matrix)) THEN + WRITE(*,*) 'ERROR: the coupled determinant needs the STRIDE' + WRITE(*,*) ' Deltaprime matrix; run with read_eq=t.' + STOP 'dispersion_det: no Deltaprime matrix' + END IF + + ALLOCATE(result_matrix(msing_max,msing_max)) + result_matrix = sl_in%dp_matrix(1:msing_max,1:msing_max) + + DO k=1,msing_max +c set module-level variables for this surface + Q_e = sl_in%Q_e_arr(k) + Q_i = sl_in%Q_i_arr(k) + P_perp = sl_in%P_perp_arr(k) + P_tor = sl_in%P_tor_arr(k) + tau = sl_in%tau_arr(k) + D_norm = sl_in%D_norm_arr(k) + c_beta = sl_in%c_beta_arr(k) + tauk = sl_in%Qconv_arr(k) + iota_e = Q_e / (Q_e - Q_i) + +c g_in is normalised by the first surface's tau_k, i.e. +c g_in = tau_1*(gamma + i*omega). The same physical +c eigenvalue reaches surface k as tau_k*(gamma + i*omega) +c = g_in*tau_k/tau_1, the normalisation its Q_e and Q_i use. + g_tmp = g_in*(tauk/sl_in%Qconv_arr(1)) +c The layer's Q is the mode frequency in the local ExB +c frame, while the coupled scan runs in the lab frame: +c Q_layer = tau_k*omega_lab - tau_k*n*Omega_E. q_shift_arr +c holds that real offset (zero unless the driver sets it), +c and g = i*Q. + g_tmp = g_tmp + ifac*sl_in%q_shift_arr(k) + +c layer response, de-normalised by S^(1/3), and Delta_crit + result_matrix(k,k) = result_matrix(k,k) + $ - riccati_f()*sl_in%lu_arr(k)**(1.0_r8/3.0_r8) + $ - CMPLX(sl_in%d_crit_arr(k), 0.0_r8, KIND=r8) + END DO + + CALL calc_determinant(result_matrix, msing_max, det_val, + $ det_status) + IF (det_status /= 0) THEN + WRITE(*,*) 'ERROR: calc_determinant failed for rank=', + $ msing_max + STOP 'dispersion_det: calc_determinant failed' + END IF + DEALLOCATE(result_matrix) + dispersion_det = det_val + END FUNCTION dispersion_det + +c----------------------------------------------------------------------- +c get_or_compute_v2: hash-cached dispersion evaluation for AMR v2. +c Applies ifac Wick rotation: g_tmp = q_in * ifac. +c----------------------------------------------------------------------- + SUBROUTINE get_or_compute_v2(q_in, idx_out, n_k, + $ sl_in, msing_max, + $ coupling_flag, full) + + IMPLICIT NONE + +c --- arguments + COMPLEX(r8), INTENT(IN) :: q_in + INTEGER, INTENT(OUT) :: idx_out + INTEGER, INTENT(IN) :: n_k, msing_max + TYPE(slayer_inputs_type), INTENT(IN) :: sl_in + LOGICAL, INTENT(IN) :: coupling_flag + LOGICAL, INTENT(OUT) :: full + +c --- locals + INTEGER :: h, curr + COMPLEX(r8) :: delta_val + INTEGER(8) :: ix8, iy8, h8 + + full = .FALSE. +c --- 1. compute hash bucket + ix8 = NINT(REAL(q_in) * HASH_SCALE, KIND=8) + iy8 = NINT(AIMAG(q_in) * HASH_SCALE, KIND=8) + h8 = MOD(ABS(ix8 * 73856093_8 + iy8 * 19349663_8), + $ INT(HASH_SZ, 8)) + 1_8 + h = INT(h8) + +c --- 2. search hash chain for existing point + curr = hash_head(h) + DO WHILE (curr /= 0) + IF (ABS(Q_store(curr) - q_in) < 1.0d-8) THEN + idx_out = curr + RETURN + END IF + curr = hash_next(curr) + END DO + +c --- 3. not found: evaluate with ifac rotation and store + n_pts = n_pts + 1 + IF (n_pts > MAX_PTS) THEN + n_pts = n_pts - 1 + full = .TRUE. + idx_out = -1 + RETURN + END IF + + idx_out = n_pts + Q_store(idx_out) = q_in + + delta_val = eval_dispersion(q_in, n_k, sl_in, msing_max, + $ coupling_flag) + D_store(idx_out) = delta_val + +c --- 4. insert into hash chain (prepend) + hash_next(idx_out) = hash_head(h) + hash_head(h) = idx_out + + END SUBROUTINE get_or_compute_v2 + +c----------------------------------------------------------------------- +c dispersion_AMR_v2: cell-based adaptive mesh refinement scanner. +c Unlike v1 (hash-based point deduplication), v2 stores complete +c cells (TYPE amr_cell_type), each carrying 4 corner Q- and D- +c values. Refinement subdivides cells that contain a zero in +c Re(D) or Im(D) and re-evaluates the dispersion relation at the +c 5 new midpoints. +c +c All dispersion evaluations go through get_or_compute_v2, which +c caches results in Q_store / D_store via a hash table. This +c eliminates redundant evaluations for shared corners (initial grid) +c and shared edge-midpoints (refinement). At completion, Q_store +c and D_store are trimmed to n_pts unique output points. +c +c----------------------------------------------------------------------- + SUBROUTINE dispersion_AMR_v2(n_k, sl_in, msing_max, + $ scan_width, Q_num, AMR_passes, + $ coupling_flag, scan_center) + + IMPLICIT NONE + +c --- arguments + INTEGER, INTENT(IN) :: n_k ! number of rational surfaces + INTEGER, INTENT(IN) :: msing_max ! max surfaces for coupling + INTEGER, INTENT(IN) :: Q_num ! grid points per axis + INTEGER, INTENT(IN) :: AMR_passes ! refinement passes + REAL(r8), INTENT(IN) :: scan_width ! half-width of Re/Im scan window + REAL(r8), INTENT(IN) :: scan_center ! Re(Q) centre of the window + TYPE(slayer_inputs_type), INTENT(IN) :: sl_in + LOGICAL, INTENT(IN) :: coupling_flag ! coupled dispersion_det? +c --- locals + TYPE(amr_cell_type), ALLOCATABLE :: new_cells(:) + TYPE(amr_cell_type), ALLOCATABLE :: swap_tmp(:) ! for pointer swap + INTEGER :: i, j, c, corner, pass ! loop counters + REAL(r8) :: step ! grid spacing + REAL(r8) :: x, y ! real / imag grid coords + LOGICAL :: cross_real, cross_imag ! zero-crossing flags + LOGICAL :: pts_full ! MAX_PTS reached flag + INTEGER :: n_new_cells ! count during refinement + INTEGER :: cells_to_refine ! cells flagged per pass + INTEGER :: cells_kept ! cells kept per pass + INTEGER :: idx_tmp ! hash-cache index + COMPLEX(r8), ALLOCATABLE :: temp_Q(:) ! for trimming output + COMPLEX(r8), ALLOCATABLE :: temp_D(:) ! for trimming output + +c --- 1. initialise cell storage and hash cache + + IF (ALLOCATED(amr_cells)) DEALLOCATE(amr_cells) + IF (ALLOCATED(Q_store)) DEALLOCATE(Q_store) + IF (ALLOCATED(D_store)) DEALLOCATE(D_store) + IF (ALLOCATED(hash_head)) DEALLOCATE(hash_head) + IF (ALLOCATED(hash_next)) DEALLOCATE(hash_next) + + ALLOCATE(amr_cells(MAX_CELLS)) + ALLOCATE(new_cells(MAX_CELLS)) + ALLOCATE(Q_store(MAX_PTS)) + ALLOCATE(D_store(MAX_PTS)) + ALLOCATE(hash_head(HASH_SZ)) + ALLOCATE(hash_next(MAX_PTS)) + + hash_head = 0 + hash_next = 0 + n_pts = 0 + n_amr_cells = 0 + step = (2.0d0 * scan_width) / DBLE(Q_num - 1) + +c --- 2. build initial coarse grid of (Q_num-1)^2 cells + + DO i = 1, Q_num - 1 + DO j = 1, Q_num - 1 + x = scan_center - scan_width + DBLE(i-1) * step + y = -scan_width + DBLE(j-1) * step + + n_amr_cells = n_amr_cells + 1 + + IF (n_amr_cells > MAX_CELLS) THEN + WRITE(*,*) 'ERROR: Exceeded MAX_CELLS in init' + STOP 'dispersion_AMR_v2: MAX_CELLS in init' + END IF + +c corner order: BL=1, BR=2, TL=3, TR=4 + amr_cells(n_amr_cells)%Q(1) = CMPLX(x, y, KIND=r8) + amr_cells(n_amr_cells)%Q(2) = CMPLX(x+step, y, KIND=r8) + amr_cells(n_amr_cells)%Q(3) = CMPLX(x, y+step, KIND=r8) + amr_cells(n_amr_cells)%Q(4) = CMPLX(x+step, y+step, + $ KIND=r8) + +c evaluate dispersion at each corner (hash-cached) + DO corner = 1, 4 + CALL get_or_compute_v2( + $ amr_cells(n_amr_cells)%Q(corner), + $ idx_tmp, n_k, sl_in, msing_max, + $ coupling_flag, pts_full) + IF (pts_full) GOTO 800 + amr_cells(n_amr_cells)%D(corner) = + $ D_store(idx_tmp) + END DO + + amr_cells(n_amr_cells)%needs_refine = .FALSE. + END DO + END DO + +c --- 3. refinement passes: subdivide cells with zero crossings + DO pass = 1, AMR_passes + WRITE(*,'(A,I2,A,I7,A)') ' Pass ', pass, + $ ': Processing ', n_amr_cells, ' cells' + +c flag cells that span a zero in Re(D) or Im(D) + cells_to_refine = 0 + DO c = 1, n_amr_cells + CALL check_cell_crossing_sub(amr_cells(c), + $ cross_real, cross_imag) + amr_cells(c)%needs_refine = (cross_real .OR. cross_imag) + IF (amr_cells(c)%needs_refine) THEN + cells_to_refine = cells_to_refine + 1 + END IF + END DO + +c build new cell list: subdivide flagged, keep the rest + n_new_cells = 0 + cells_kept = 0 + + DO c = 1, n_amr_cells + IF (amr_cells(c)%needs_refine) THEN + CALL subdivide_cell_sub( + $ amr_cells(c), + $ new_cells, n_new_cells, + $ MAX_CELLS, n_k, sl_in, + $ msing_max, coupling_flag, + $ pts_full) + IF (pts_full) GOTO 800 + ELSE + n_new_cells = n_new_cells + 1 + IF (n_new_cells > MAX_CELLS) THEN + WRITE(*,*) 'ERROR: Exceeded MAX_CELLS in refine' + STOP 'dispersion_AMR_v2: MAX_CELLS in refine' + END IF + new_cells(n_new_cells) = amr_cells(c) + cells_kept = cells_kept + 1 + END IF + END DO + +c swap arrays for next pass (pointer swap, no element copy) + CALL MOVE_ALLOC(new_cells, swap_tmp) + CALL MOVE_ALLOC(amr_cells, new_cells) ! old amr_cells becomes new_cells + CALL MOVE_ALLOC(swap_tmp, amr_cells) ! filled array becomes amr_cells + n_amr_cells = n_new_cells + + END DO + GOTO 810 + + 800 CONTINUE + WRITE(*,'(A,I8,A)') + $ ' WARNING: MAX_PTS (', MAX_PTS, + $ ') reached during AMR v2.' + WRITE(*,'(A)') + $ ' Saving existing results.' + + 810 CONTINUE +c --- 4. output: Q_store/D_store already populated by hash cache. +c Trim to exact size n_pts and deallocate hash infrastructure. + + ALLOCATE(temp_Q(n_pts)) + ALLOCATE(temp_D(n_pts)) + temp_Q(1:n_pts) = Q_store(1:n_pts) + temp_D(1:n_pts) = D_store(1:n_pts) + CALL MOVE_ALLOC(temp_Q, Q_store) + CALL MOVE_ALLOC(temp_D, D_store) + + IF (ALLOCATED(hash_head)) DEALLOCATE(hash_head) + IF (ALLOCATED(hash_next)) DEALLOCATE(hash_next) + IF (ALLOCATED(new_cells)) DEALLOCATE(new_cells) +c keep amr_cells allocated for potential post-run inspection + + WRITE(*,*) 'AMR v2 Complete. Unique output points:', n_pts + WRITE(*,'(A,2ES14.6)') ' D_store checksum (Re,Im):', + $ SUM(REAL(D_store(1:n_pts))), + $ SUM(AIMAG(D_store(1:n_pts))) + WRITE(*,'(A,2ES14.6)') ' D_store(1) sample:', + $ REAL(D_store(1)), AIMAG(D_store(1)) + + RETURN + END SUBROUTINE dispersion_AMR_v2 + +c----------------------------------------------------------------------- +c check_cell_crossing_sub: test whether a cell’s 4 corner D-values +c span a zero crossing in Re(D) and/or Im(D). Used by +c dispersion_AMR_v2 to decide which cells to refine. +c----------------------------------------------------------------------- + SUBROUTINE check_cell_crossing_sub(cell, cross_real, cross_imag) + + IMPLICIT NONE + + TYPE(amr_cell_type), INTENT(IN) :: cell + LOGICAL, INTENT(OUT) :: cross_real, cross_imag + + REAL(r8) :: r_vals(4), i_vals(4) ! corner Re/Im values + REAL(r8) :: r_min, r_max, i_min, i_max + INTEGER :: k + + DO k = 1, 4 + r_vals(k) = REAL(cell%D(k), KIND=r8) + i_vals(k) = AIMAG(cell%D(k)) + END DO + + r_min = MINVAL(r_vals) + r_max = MAXVAL(r_vals) + cross_real = (r_min * r_max <= 0.0d0) + + i_min = MINVAL(i_vals) + i_max = MAXVAL(i_vals) + cross_imag = (i_min * i_max <= 0.0d0) + + RETURN + END SUBROUTINE check_cell_crossing_sub + + +c----------------------------------------------------------------------- +c subdivide_cell_sub: split a parent cell into 4 child cells by +c computing 5 midpoints (bottom-mid, top-mid, left-mid, right-mid, +c centre) and evaluating the dispersion relation at each. The 4 +c resulting child cells are appended to new_cells. +c----------------------------------------------------------------------- + SUBROUTINE subdivide_cell_sub(parent, + $ new_cells, n_new, max_cells, n_k, + $ sl_in, msing_max, coupling_flag, + $ pts_full) + + IMPLICIT NONE + +c --- arguments + TYPE(amr_cell_type), INTENT(IN) :: parent + TYPE(amr_cell_type), INTENT(INOUT) :: new_cells(*) + INTEGER, INTENT(INOUT) :: n_new ! running count of new cells + INTEGER, INTENT(IN) :: max_cells ! capacity of new_cells + INTEGER, INTENT(IN) :: n_k + INTEGER, INTENT(IN) :: msing_max + TYPE(slayer_inputs_type), INTENT(IN) :: sl_in + LOGICAL, INTENT(IN) :: coupling_flag + LOGICAL, INTENT(OUT) :: pts_full ! MAX_PTS flag +c --- corner coordinates and D-values from parent + COMPLEX(r8) :: q_bl, q_br, q_tl, q_tr + COMPLEX(r8) :: d_bl, d_br, d_tl, d_tr +c --- midpoint coordinates and D-values (cached via hash) + COMPLEX(r8) :: q_bm, q_tm, q_lm, q_rm, q_mm + COMPLEX(r8) :: d_bm, d_tm, d_lm, d_rm, d_mm + INTEGER :: idx_tmp ! hash-cache index + +c --- extract parent corners (BL=1, BR=2, TL=3, TR=4) + q_bl = parent%Q(1) + q_br = parent%Q(2) + q_tl = parent%Q(3) + q_tr = parent%Q(4) + + d_bl = parent%D(1) + d_br = parent%D(2) + d_tl = parent%D(3) + d_tr = parent%D(4) + +c --- compute 5 midpoint coordinates + q_bm = 0.5d0 * (q_bl + q_br) + q_tm = 0.5d0 * (q_tl + q_tr) + q_lm = 0.5d0 * (q_bl + q_tl) + q_rm = 0.5d0 * (q_br + q_tr) + q_mm = 0.25d0 * (q_bl + q_br + q_tl + q_tr) + +c --- evaluate dispersion at new midpoints (hash-cached) + pts_full = .FALSE. + CALL get_or_compute_v2(q_bm, idx_tmp, + $ n_k, sl_in, msing_max, + $ coupling_flag, pts_full) + IF (pts_full) RETURN + d_bm = D_store(idx_tmp) + CALL get_or_compute_v2(q_tm, idx_tmp, + $ n_k, sl_in, msing_max, + $ coupling_flag, pts_full) + IF (pts_full) RETURN + d_tm = D_store(idx_tmp) + CALL get_or_compute_v2(q_lm, idx_tmp, + $ n_k, sl_in, msing_max, + $ coupling_flag, pts_full) + IF (pts_full) RETURN + d_lm = D_store(idx_tmp) + CALL get_or_compute_v2(q_rm, idx_tmp, + $ n_k, sl_in, msing_max, + $ coupling_flag, pts_full) + IF (pts_full) RETURN + d_rm = D_store(idx_tmp) + CALL get_or_compute_v2(q_mm, idx_tmp, + $ n_k, sl_in, msing_max, + $ coupling_flag, pts_full) + IF (pts_full) RETURN + d_mm = D_store(idx_tmp) + +c --- check space for 4 new cells + IF (n_new + 4 > max_cells) THEN + WRITE(*,*) 'ERROR: Would exceed MAX_CELLS in subdivide' + STOP 'subdivide_cell_sub: MAX_CELLS exceeded' + END IF + +c --- child 1: bottom-left quadrant (BL, BM, LM, MM) + n_new = n_new + 1 + new_cells(n_new)%Q(1) = q_bl + new_cells(n_new)%Q(2) = q_bm + new_cells(n_new)%Q(3) = q_lm + new_cells(n_new)%Q(4) = q_mm + new_cells(n_new)%D(1) = d_bl + new_cells(n_new)%D(2) = d_bm + new_cells(n_new)%D(3) = d_lm + new_cells(n_new)%D(4) = d_mm + new_cells(n_new)%needs_refine = .FALSE. + +c --- child 2: bottom-right quadrant (BM, BR, MM, RM) + n_new = n_new + 1 + new_cells(n_new)%Q(1) = q_bm + new_cells(n_new)%Q(2) = q_br + new_cells(n_new)%Q(3) = q_mm + new_cells(n_new)%Q(4) = q_rm + new_cells(n_new)%D(1) = d_bm + new_cells(n_new)%D(2) = d_br + new_cells(n_new)%D(3) = d_mm + new_cells(n_new)%D(4) = d_rm + new_cells(n_new)%needs_refine = .FALSE. + +c --- child 3: top-left quadrant (LM, MM, TL, TM) + n_new = n_new + 1 + new_cells(n_new)%Q(1) = q_lm + new_cells(n_new)%Q(2) = q_mm + new_cells(n_new)%Q(3) = q_tl + new_cells(n_new)%Q(4) = q_tm + new_cells(n_new)%D(1) = d_lm + new_cells(n_new)%D(2) = d_mm + new_cells(n_new)%D(3) = d_tl + new_cells(n_new)%D(4) = d_tm + new_cells(n_new)%needs_refine = .FALSE. + +c --- child 4: top-right quadrant (MM, RM, TM, TR) + n_new = n_new + 1 + new_cells(n_new)%Q(1) = q_mm + new_cells(n_new)%Q(2) = q_rm + new_cells(n_new)%Q(3) = q_tm + new_cells(n_new)%Q(4) = q_tr + new_cells(n_new)%D(1) = d_mm + new_cells(n_new)%D(2) = d_rm + new_cells(n_new)%D(3) = d_tm + new_cells(n_new)%D(4) = d_tr + new_cells(n_new)%needs_refine = .FALSE. + + RETURN + END SUBROUTINE subdivide_cell_sub + +c----------------------------------------------------------------------- +c subprogram 12. eval_dispersion. +c Dispersion function at one scan-plane point q_in, with the ifac +c rotation g = i*q_in applied. Uncoupled: riccati_f() - delta_eff +c using the module-level layer parameters of the current surface. +c Coupled: the determinant from dispersion_det. +c----------------------------------------------------------------------- + FUNCTION eval_dispersion(q_in, n_k, sl_in, msing_max, + $ coupling_flag) + + IMPLICIT NONE + + COMPLEX(r8), INTENT(IN) :: q_in + INTEGER, INTENT(IN) :: n_k, msing_max + TYPE(slayer_inputs_type), INTENT(IN) :: sl_in + LOGICAL, INTENT(IN) :: coupling_flag + COMPLEX(r8) :: eval_dispersion + + IF (coupling_flag) THEN +c dispersion_det sets g_tmp per-surface internally; +c pass q_in*ifac directly as g_in argument. + eval_dispersion = dispersion_det(q_in * ifac, n_k, sl_in, + $ msing_max) + ELSE + g_tmp = q_in * ifac + eval_dispersion = riccati_f() - delta_eff + END IF + + END FUNCTION eval_dispersion + +c----------------------------------------------------------------------- +c subprogram 13. extract_roots. +c Locate the zeros and poles of the dispersion function from the +c leaf cells left in amr_cells by dispersion_AMR_v2. +c +c The dispersion function is analytic in the scan variable, so the +c phase of D winds by +2*pi counter-clockwise around a cell that +c holds a simple zero and by -2*pi around one that holds a simple +c pole. The winding number is taken from the four corner values. +c Each zero is then refined by a secant iteration confined to the +c neighbourhood of its cell; if that leaves the neighbourhood or +c does not converge, the best point evaluated is kept and the root +c is marked unpolished. +c +c Roots are returned in the scan plane (q = omega_hat - i*gamma_hat) +c sorted by descending growth rate. No physical filtering of the +c roots is attempted here. +c----------------------------------------------------------------------- + SUBROUTINE extract_roots(n_k, sl_in, msing_max, coupling_flag, + $ tauk_ref, roots) + + IMPLICIT NONE + +c --- arguments + INTEGER, INTENT(IN) :: n_k, msing_max + TYPE(slayer_inputs_type), INTENT(IN) :: sl_in + LOGICAL, INTENT(IN) :: coupling_flag + REAL(r8), INTENT(IN) :: tauk_ref ! tau_k of the scanned Q + TYPE(roots_outputs_type), INTENT(OUT) :: roots +c --- parameters + INTEGER, PARAMETER :: ccw(5) = (/1,2,4,3,1/) ! BL,BR,TR,TL,BL + INTEGER, PARAMETER :: MAXIT_POLISH = 30 + REAL(r8), PARAMETER :: POLISH_TOL = 1.0e-9_r8 +c --- locals + INTEGER :: c, j, it, nw, i, ins + REAL(r8) :: wsum, h, nan + COMPLEX(r8) :: d(4), z, q_c, d_c, slope + COMPLEX(r8) :: qa, qb, qn, da, db, dn, q_best, d_best + LOGICAL :: conv, dup + + nan = ieee_value(1.0_r8, ieee_quiet_nan) + roots%n_root = 0 + roots%n_pole = 0 + roots%q_root = CMPLX(nan, nan, KIND=r8) + roots%q_pole = CMPLX(nan, nan, KIND=r8) + roots%polished = .FALSE. + roots%tauk = tauk_ref + + IF (.NOT. ALLOCATED(amr_cells)) RETURN + + DO c = 1, n_amr_cells + d = amr_cells(c)%D +c skip cells with a failed (non-finite) or exactly zero corner + IF (ANY(d /= d) .OR. ANY(ABS(d) > HUGE(1.0_r8)) + $ .OR. ANY(ABS(d) == 0.0_r8)) CYCLE + +c winding number of D around the cell, counter-clockwise + wsum = 0.0_r8 + DO j = 1, 4 + z = (d(ccw(j+1))/ABS(d(ccw(j+1)))) + $ * CONJG(d(ccw(j))/ABS(d(ccw(j)))) + wsum = wsum + ATAN2(AIMAG(z), REAL(z)) + END DO + nw = NINT(wsum / (2.0_r8*pi)) + IF (nw == 0) CYCLE + + q_c = 0.25_r8 * SUM(amr_cells(c)%Q) + h = REAL(amr_cells(c)%Q(2) - amr_cells(c)%Q(1)) + +c pole: record the cell centre + IF (nw < 0) THEN + IF (roots%n_pole < MAX_ROOTS) THEN + roots%n_pole = roots%n_pole + 1 + roots%q_pole(roots%n_pole) = q_c + END IF + CYCLE + END IF + +c zero: linear estimate from the corner values, then secant + d_c = 0.25_r8 * SUM(d) + slope = ((d(2) - d(1)) + (d(4) - d(3))) / (2.0_r8*h) + qa = q_c + da = eval_dispersion(qa, n_k, sl_in, msing_max, + $ coupling_flag) + qb = q_c + 0.25_r8*h + IF (ABS(slope) > 0.0_r8) THEN + qn = q_c - d_c/slope + IF (ABS(REAL(qn - q_c)) <= 0.5_r8*h .AND. + $ ABS(AIMAG(qn - q_c)) <= 0.5_r8*h .AND. + $ ABS(qn - q_c) > 0.0_r8) qb = qn + END IF + db = eval_dispersion(qb, n_k, sl_in, msing_max, + $ coupling_flag) + IF (ABS(db) < ABS(da) .OR. da /= da) THEN + q_best = qb + d_best = db + ELSE + q_best = qa + d_best = da + END IF + conv = .FALSE. + DO it = 1, MAXIT_POLISH + IF (da /= da .OR. db /= db) EXIT + IF (ABS(db - da) == 0.0_r8) EXIT + qn = qb - db*(qb - qa)/(db - da) + IF (qn /= qn) EXIT +c stay within 1.5 cell widths of the cell centre + IF (ABS(REAL(qn - q_c)) > 1.5_r8*h .OR. + $ ABS(AIMAG(qn - q_c)) > 1.5_r8*h) EXIT + dn = eval_dispersion(qn, n_k, sl_in, msing_max, + $ coupling_flag) + IF (dn /= dn) EXIT + qa = qb + da = db + qb = qn + db = dn + IF (ABS(db) < ABS(d_best)) THEN + q_best = qb + d_best = db + END IF + IF (ABS(qb - qa) <= + $ POLISH_TOL*MAX(1.0_r8, ABS(qb))) THEN + conv = .TRUE. + q_best = qb + d_best = db + EXIT + END IF + END DO + +c drop a root already found from a neighbouring cell + dup = .FALSE. + DO i = 1, roots%n_root + IF (ABS(roots%q_root(i) - q_best) < 0.5_r8*ABS(h)) + $ dup = .TRUE. + END DO + IF (dup .OR. roots%n_root >= MAX_ROOTS) CYCLE + +c insert, keeping descending growth rate (ascending AIMAG(q)) + ins = roots%n_root + 1 + DO i = 1, roots%n_root + IF (AIMAG(q_best) < AIMAG(roots%q_root(i))) THEN + ins = i + EXIT + END IF + END DO + DO i = roots%n_root, ins, -1 + roots%q_root(i+1) = roots%q_root(i) + roots%polished(i+1) = roots%polished(i) + END DO + roots%q_root(ins) = q_best + roots%polished(ins) = conv + roots%n_root = roots%n_root + 1 + END DO + +c --- report + WRITE(*,'(A,I0,A,I0)') ' Dispersion roots found: ', + $ roots%n_root, ', poles: ', roots%n_pole + DO i = 1, roots%n_root + WRITE(*,'(A,I2,A,2ES14.6,A,ES14.6,A,ES14.6,A,L1)') + $ ' root', i, ': Q=', REAL(roots%q_root(i)), + $ -AIMAG(roots%q_root(i)), + $ ' omega[rad/s]=', REAL(roots%q_root(i))/tauk_ref, + $ ' gamma[1/s]=', -AIMAG(roots%q_root(i))/tauk_ref, + $ ' polished=', roots%polished(i) + END DO + IF (roots%n_root == 0) WRITE(*,'(A)') + $ ' WARNING: no dispersion root in the scan window.' + + RETURN + END SUBROUTINE extract_roots + + END MODULE growthrates_mod diff --git a/slayer/gslayer.f b/slayer/gslayer.f index 1e51f569c..76653d162 100644 --- a/slayer/gslayer.f +++ b/slayer/gslayer.f @@ -1,110 +1,200 @@ MODULE gslayer_mod +c----------------------------------------------------------------------- +c gslayer_mod: Single-surface SLAYER driver called by GPEC. +c +c This module contains the gpec_slayer subroutine, which computes +c Delta, torque, and critical field threshold for one rational +c surface. It is the primary interface used by gpec/gpout.f. +c +c Growth-rate scanning, AMR dispersion solvers, and I/O helpers +c have been split out into growthrates_mod (growthrates.f). +c +c Subprograms: +c 1. gpec_slayer - single-surface delta/torque driver +c----------------------------------------------------------------------- + + USE omp_lib - USE sglobal_mod, ONLY: out_unit, r8, mu0, m_p, chag, lnLamb, + USE sglobal_mod, ONLY: out_unit,r8,mu0,m_p,chag,lnLamb, $ Q_e,Q_i,pr,pe,c_beta,ds,tau, $ eta,visc,rho_s,lu,omega_e,omega_i, - $ delta_n, - $ Q - USE delta_mod, ONLY: riccati,riccati_out, + $ delta_n,Q, + $ ifac,g_tmp,pi, ! used by AMR + $ tauk,iota_e,D_norm,P_perp,P_tor,delta_eff, ! used by det + $ amr_cell_type,amr_cells,n_amr_cells, ! v2 types + $ Q_store,D_store,n_pts, ! output arrays + $ MAX_PTS,HASH_SZ,HASH_SCALE, ! v1 constants + $ hash_head,hash_next, ! v1 hash + $ MAX_CELLS, ! v2 limit + $ slayer_inputs_type,slayer_outputs_type, + $ deltas_outputs_type, + $ tau_r,dc_tmp,dc_type, + $ sn_str,sm_str + USE delta_mod, ONLY: riccati,riccati_f,riccati_out, $ parflow_flag,PeOhmOnly_flag + USE params_mod + USE layerinputs_mod + USE slayer_netcdf_mod IMPLICIT NONE +c --- reconnection regulariser used in psi0 / JxB expressions; +c expose via namelist to make user-configurable. + REAL(r8), PARAMETER :: DELTA_N_PERT = 1.0e-2_r8 + CONTAINS c----------------------------------------------------------------------- c subprogram 1. gpec_slayer. -c run slayer to provide b_crit(ising). -c----------------------------------------------------------------------- -c----------------------------------------------------------------------- -c declarations. +c Single-surface SLAYER driver: compute Delta, torque, and b_crit +c for one rational surface characterised by its plasma profiles. +c +c Steps: +c 1. Derive normalised SLAYER parameters (lu, ds, Q, Q_e, etc.) +c from dimensional inputs. +c 2. Compute baseline Delta via riccati(). +c 3. Scan over a range of Q (rotation) to build a torque +c balance curve, then identify the critical threshold br_th. +c +c TODO: `zeff` and `qval` are INTENT(IN) +c arguments accepted in the interface but never referenced in +c the subroutine body. Remove from signature (breaking API +c change) or add a comment documenting their reserved intent. c----------------------------------------------------------------------- SUBROUTINE gpec_slayer(n_e,t_e,n_i,t_i,zeff,omega,omega_e, $ omega_i,qval,sval,bt,rs,R0,mu_i,inpr,mms,nns,ascii_flag, $ delta,psi0,jxb,omega_sol,br_th) - REAL(r8),INTENT(IN) :: n_e,t_e,n_i,t_i,omega,omega_e,omega_i, - $ qval,sval,bt,rs,R0,zeff,inpr - INTEGER, INTENT(IN) :: mms,nns,mu_i - LOGICAL, INTENT(IN) :: ascii_flag - COMPLEX(r8),INTENT(OUT) :: delta,psi0 - REAL(r8),INTENT(OUT) :: jxb,omega_sol,br_th - - INTEGER :: i,inum - INTEGER, DIMENSION(1) :: index - - REAL(r8) :: inQ,inQ_e,inQ_i,inpe,inc_beta,inds,intau,inlu - REAL(r8) :: mrs,nrs,rho,b_l,v_a,Qconv,Q0,delta_n_p, - $ lbeta,tau_i,tau_h,tau_r,tau_v - REAL(r8) :: inQ_min,inQ_max,Q_sol,maxbal - INTEGER :: ipass - INTEGER, PARAMETER :: nref=4 - REAL(r8) :: xpk,rlo,rhi,rdq,rdqc,xq,bloc,jloc - COMPLEX(r8) :: dloc - - REAL(r8), DIMENSION(:), ALLOCATABLE :: inQs,iinQs,jxbl,bal - COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: deltal - CHARACTER(3) :: sn,sm +c --- input arguments: dimensional plasma profiles for this surface + REAL(r8),INTENT(IN) :: n_e ! electron density [m^-3] + REAL(r8),INTENT(IN) :: t_e ! electron temperature [eV] + REAL(r8),INTENT(IN) :: n_i ! ion density [m^-3] + REAL(r8),INTENT(IN) :: t_i ! ion temperature [eV] + REAL(r8),INTENT(IN) :: zeff ! (TODO: unused; see header) + REAL(r8),INTENT(IN) :: omega ! plasma rotation frequency + REAL(r8),INTENT(IN) :: omega_e ! electron diamagnetic freq + REAL(r8),INTENT(IN) :: omega_i ! ion diamagnetic freq + REAL(r8),INTENT(IN) :: qval ! (TODO: unused; see header) + REAL(r8),INTENT(IN) :: sval ! magnetic shear + REAL(r8),INTENT(IN) :: bt ! toroidal field [T] + REAL(r8),INTENT(IN) :: rs ! minor radius of surface [m] + REAL(r8),INTENT(IN) :: R0 ! major radius [m] + REAL(r8),INTENT(IN) :: inpr ! Prandtl number + INTEGER, INTENT(IN) :: mms ! poloidal mode number + INTEGER, INTENT(IN) :: nns ! toroidal mode number + INTEGER, INTENT(IN) :: mu_i ! ion mass number (AMU) + LOGICAL, INTENT(IN) :: ascii_flag ! write ASCII torque-balance file +c --- output arguments + COMPLEX(r8),INTENT(OUT) :: delta ! complex Delta + COMPLEX(r8),INTENT(OUT) :: psi0 ! reconnected flux (a.u.) + REAL(r8),INTENT(OUT) :: jxb ! JxB torque (a.u.) + REAL(r8),INTENT(OUT) :: omega_sol ! rotation at torque-balance + REAL(r8),INTENT(OUT) :: br_th ! critical radial field threshold + +c --- loop / scan control + INTEGER :: i ! loop index + INTEGER :: inum ! number of scan points + INTEGER, DIMENSION(1) :: max_idx ! index of max torque balance +c --- local copies of normalised parameters for riccati() + REAL(r8) :: inQ,inQ_e,inQ_i ! normalised frequencies + REAL(r8) :: inpe ! electron Prandtl (set to 0) + REAL(r8) :: inc_beta,inds ! c_beta, ds copies + REAL(r8) :: intau ! tau copy +c --- derived dimensional quantities + REAL(r8) :: mrs,nrs ! real-valued mode numbers + REAL(r8) :: rho ! mass density [kg/m^3] + REAL(r8) :: b_l ! characteristic magnetic field + REAL(r8) :: Qconv ! Q normalisation factor + REAL(r8) :: Q0 ! initial Q (before scan) + REAL(r8) :: delta_n_p ! Delta_n perturbation + REAL(r8) :: lbeta ! local beta + REAL(r8) :: tau_i ! ion collision time + REAL(r8) :: tau_h ! Alfven time across surface + REAL(r8) :: tau_v ! viscous diffusion time +c --- scan workspace + REAL(r8) :: inQ_min,inQ_max ! scan bounds in Q + REAL(r8) :: Q_sol ! Q at torque-balance threshold + REAL(r8) :: maxbal ! refined maximum of the balance + INTEGER :: ipass ! nose-refinement pass index + INTEGER, PARAMETER :: nref=4 ! number of nose-refinement passes + REAL(r8) :: xpk,rlo,rhi ! running peak and its bracket in Q + REAL(r8) :: rdq,rdqc,xq ! refined/coarse spacing, scan point + REAL(r8) :: bloc,jloc ! local balance and JxB values + COMPLEX(r8) :: dloc ! local delta value + REAL(r8), DIMENSION(:), ALLOCATABLE :: inQs ! Q scan array + REAL(r8), DIMENSION(:), ALLOCATABLE :: jxbl ! JxB scan array + REAL(r8), DIMENSION(:), ALLOCATABLE :: bal ! balance scan array + COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: deltal ! delta scan array + CHARACTER(3) :: l_sn,l_sm ! local mode number strings +c --- set delta_mod control flags for this single-surface call parflow_flag=.FALSE. PeOhmOnly_flag=.TRUE. riccati_out=.FALSE. mrs = real(mms,4) nrs = real(nns,4) - - ! String representations of the m and n mode numbers +c----------------------------------------------------------------------- +c build string representations of m and n for file names. +c----------------------------------------------------------------------- IF (nns<10) THEN - WRITE(UNIT=sn,FMT='(I1)') nns - sn=ADJUSTL(sn) + WRITE(UNIT=l_sn,FMT='(I1)') nns + l_sn=ADJUSTL(l_sn) ELSE - WRITE(UNIT=sn,FMT='(I2)') nns + WRITE(UNIT=l_sn,FMT='(I2)') nns ENDIF IF (mms<10) THEN - WRITE(UNIT=sm,FMT='(I1)') mms - sm=ADJUSTL(sm) + WRITE(UNIT=l_sm,FMT='(I1)') mms + l_sm=ADJUSTL(l_sm) ELSEIF (mms<100) THEN - WRITE(UNIT=sm,FMT='(I2)') mms - sm=ADJUSTL(sm) + WRITE(UNIT=l_sm,FMT='(I2)') mms + l_sm=ADJUSTL(l_sm) ELSE - WRITE(UNIT=sm,FMT='(I3)') mms + WRITE(UNIT=l_sm,FMT='(I3)') mms ENDIF - inpe=0.0 ! Waybright added this + inpe=0.0 ! electron Prandtl not used here - tau= t_i/t_e ! ratio of ion to electron temperature - tau_i = 6.6e17*mu_i**0.5*(t_i/1e3)**1.5/(n_e*lnLamb) ! ion colls. - eta= 1.65e-9*lnLamb/(t_e/1e3)**1.5 ! spitzer resistivity (wesson) +c----------------------------------------------------------------------- +c derive normalised SLAYER parameters from dimensional inputs. +c tau, eta, rho, b_l, v_a are intermediate dimensional quantities; +c lu (Lundquist), ds, Qconv, Q, Q_e, Q_i, c_beta are the +c normalised parameters used by riccati(). +c----------------------------------------------------------------------- + tau= t_i/t_e ! ion-to-electron temp ratio + tau_i = 6.6e17*mu_i**0.5*(t_i/1e3)**1.5/(n_e*lnLamb) ! ion coll. time + eta= 1.65e-9*lnLamb/(t_e/1e3)**1.5 ! Spitzer resistivity (Wesson) rho=(mu_i*m_p)*n_e ! mass density - b_l=(nrs/mrs)*rs*sval*bt/R0 ! characteristic magnetic field - v_a=b_l/(mu0*rho)**0.5 ! alfven velocity - rho_s=1.02e-4*(mu_i*t_e)**0.5/bt ! ion Lamour by elec. Temp. - - tau_h=R0*(mu0*rho)**0.5/(nns*sval*bt) ! alfven time across surface - tau_r=mu0*rs**2.0/eta ! resistive time scale - tau_v=tau_r/inpr ! rho*rs**2.0/visc ! viscous time scale + b_l=(nrs/mrs)*nrs*sval*bt/R0 ! characteristic magnetic field + rho_s=1.02e-4*(mu_i*t_e)**0.5/bt ! ion Larmor radius at T_e - ! this one must be anomalous. calculated back from pr. - visc= rho*rs**2.0/tau_v + tau_h=R0*(mu0*rho)**0.5/(nns*ABS(sval)*bt) ! Alfven time, |s| + tau_r=mu0*rs**2.0/eta ! resistive diffusion time + tau_v=tau_r/inpr ! viscous diffusion time + visc= rho*rs**2.0/tau_v ! back-calculated viscosity - lu=tau_r/tau_h ! Lundquist number + lu=tau_r/tau_h ! Lundquist number + Qconv=lu**(1.0/3.0)*tau_h ! Q normalisation factor - Qconv=lu**(1.0/3.0)*tau_h ! conversion to Qs based on Cole +c --- normalised frequencies. GPEC passes the bare ExB and +c diamagnetic frequencies with no toroidal mode number attached +c (gpout.f). The layer sees the mode at the Doppler frequency +c n*omega_E, and omega_* = k_theta*v_* = n*(dp/dpsi)/(e*n_e), so +c GPEC's n (nns) multiplies all three here, as in build_inputs +c and the coupled determinant's Doppler shift (slayer.f). + Q=Qconv*nns*omega + Q_e=-Qconv*nns*omega_e + Q_i=-Qconv*nns*omega_i - ! note Q depends on Qconv even if omega is fixed. - Q=Qconv*omega - Q_e=-Qconv*omega_e - Q_i=-Qconv*omega_i - - ! This is the most critical parameter - ds=lu**(1.0/3.0)*rho_s/rs ! conversion based on Cole. + ds=lu**(1.0/3.0)*rho_s/rs ! normalised ion sound radius lbeta=(5.0/3.0)*mu0*n_e*chag*(t_e+t_i)/bt**2.0 - c_beta=(lbeta/(1.0+lbeta))**0.5 + c_beta=(lbeta/(1.0+lbeta))**0.5 ! compressibility parameter - delta_n=lu**(1.0/3.0)/rs ! norm factor for delta primes + delta_n=lu**(1.0/3.0)/rs ! Delta normalisation factor +c --- copy normalised values into local variables for riccati() call inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -113,14 +203,24 @@ SUBROUTINE gpec_slayer(n_e,t_e,n_i,t_i,zeff,omega,omega_e, intau=tau Q0=Q c----------------------------------------------------------------------- -c calculate basic delta, torque, balance, error fields. +c compute baseline Delta, reconnected flux, and JxB torque. +c delta_n_p uses the module-level DELTA_N_PERT constant; promote +c that constant to a namelist input to make it user-configurable. c----------------------------------------------------------------------- - delta_n_p=1e-2 + delta_n_p = DELTA_N_PERT delta=riccati(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe) - psi0=1.0/ABS(delta+delta_n_p) ! a.u. - jxb=-AIMAG(1.0/(delta+delta_n_p)) ! a.u. + psi0=1.0/ABS(delta+delta_n_p) ! reconnected flux (a.u.) + jxb=-AIMAG(1.0/(delta+delta_n_p)) ! JxB torque (a.u.) c----------------------------------------------------------------------- -c find solutions based on simple torque balance. +c torque-balance scan: sweep Q over [inQ_min, inQ_max] and +c compute delta(Q), JxB(Q), and the balance parameter. +c The threshold br_th is sqrt(max(bal)/lu * s^2/2). +c +c TODO: physics bounds from Q0/Q_e (computed +c in the IF/ELSE above) are immediately overridden by the two +c fixed assignments below. To use physics bounds, remove the +c inQ_max=10 / inQ_min=-10 lines. To keep fixed bounds, remove +c the dead IF/ELSE block above. c----------------------------------------------------------------------- IF (Q0>inQ_e) THEN inQ_max=2.0*Q0 @@ -134,8 +234,7 @@ SUBROUTINE gpec_slayer(n_e,t_e,n_i,t_i,zeff,omega,omega_e, ENDIF ENDIF - ! Scan of rotation - inQ_max=10.0 + inQ_max=10.0 ! TODO: fixed bound; see header inQ_min=-10.0 inum=200 ALLOCATE(inQs(0:inum),deltal(0:inum),jxbl(0:inum),bal(0:inum)) @@ -147,10 +246,10 @@ SUBROUTINE gpec_slayer(n_e,t_e,n_i,t_i,zeff,omega,omega_e, bal(i)=2.0*inpr*(Q0-inQs(i))/jxbl(i) ENDDO - ! Write torque balance curves to file for diagnostic purposes +c --- optionally write torque balance curve to ASCII file IF(ascii_flag)THEN OPEN(UNIT=out_unit,FILE="gpec_slayer_torque_balance_m"// - $ TRIM(sm)//"_n"//TRIM(sn)//".out", + $ TRIM(l_sm)//"_n"//TRIM(l_sn)//".out", $ STATUS="UNKNOWN") WRITE(out_unit,'(1x,5(a17))') "inQ","RE(delta)", $ "IM(delta)","jxb","bal" @@ -169,9 +268,9 @@ SUBROUTINE gpec_slayer(n_e,t_e,n_i,t_i,zeff,omega,omega_e, ! maximum and iteratively refine the scan around the running peak ! so the nose is resolved regardless of how narrow it is. Exclude ! non-finite entries (NaN from 0/0, Inf from a near-zero jxbl). - index=MAXLOC(bal,MASK=(bal==bal .AND. ABS(bal) r_s) + REAL(r8) :: k_ref ! K = r_s*dpsi_N/dr + REAL(r8) :: alpha_m ! Mercier exponent sqrt(-D_I) + INTEGER :: jsing ! second surface index +c----------------------------------------------------------------------- +c read STRIDE NetCDF: Deltaprime matrix, geometry, equilibrium scalars. +c----------------------------------------------------------------------- + CALL read_stride_netcdf_diagonal(ncfile,msing,dp_mat, + $ Re_dp_diagonal,Im_dp_diagonal,q_rational, + $ psi_n_rational,shear,r_o,my_bt0,my_psio,dr_vals, + $ di_vals,mpsi_arr,nn_arr,resm) + +c The toroidal mode number of the STRIDE run defines the layer +c problem: it sets the module-level nn used by params() (Alfven +c time) and by the output file names, whatever slayer.in holds. + IF (nn /= nn_arr(1) .AND. nn /= 0) THEN + WRITE(*,'(1x,a,i0,a,i0)') 'WARNING: slayer.in nn = ',nn, + $ ' replaced by the STRIDE toroidal mode number ',nn_arr(1) + END IF + nn = nn_arr(1) + IF (nn < 10) THEN + WRITE(UNIT=sn_str,FMT='(I1)') nn + sn_str = ADJUSTL(sn_str) + ELSE + WRITE(UNIT=sn_str,FMT='(I2)') nn + ENDIF + + mpsi = INT(mpsi_arr(1)) + mthsurf = 512 ! poloidal segments for surface integral + +c --- allocate slayer_inputs_type arrays (one entry per surface) + ALLOCATE(sl_in%qval_arr(msing), sl_in%omegas_arr(msing), + $ sl_in%omegas_e_arr(msing), sl_in%dp_matrix(msing,msing), + $ sl_in%omegas_i_arr(msing), + $ sl_in%Q_e_arr(msing), sl_in%Q_i_arr(msing), + $ sl_in%psi_n_arr(msing), + $ sl_in%Re_dp_arr(msing), sl_in%Im_dp_arr(msing), + $ sl_in%d_crit_arr(msing), sl_in%P_tor_arr(msing), + $ sl_in%P_perp_arr(msing), sl_in%tau_arr(msing), + $ sl_in%D_norm_arr(msing), sl_in%D_fac_arr(msing), + $ sl_in%d_beta_arr(msing), sl_in%gammafac_arr(msing), + $ sl_in%c_beta_arr(msing), sl_in%lu_arr(msing), + $ sl_in%Qconv_arr(msing), + $ sl_in%Re_dp_psi_arr(msing), sl_in%Im_dp_psi_arr(msing), + $ sl_in%k_ref_arr(msing), sl_in%alpha_arr(msing), + $ sl_in%q_shift_arr(msing), sl_in%q_arr(msing), + $ sl_in%Qexb_arr(msing), sl_in%shear_arr(msing), + $ sl_in%ds_arr(msing)) + sl_in%q_shift_arr = 0.0_r8 + +c --- allocate local kinetic arrays (for future NetCDF diagnostic output) + ALLOCATE(ne_arr(msing), te_arr(msing), ni_arr(msing), + $ ti_arr(msing), zeff_arr(msing), bt_arr(msing), + $ rs_arr(msing), R0_arr(msing), mu_i_arr(msing), + $ nns_arr(msing), dr_arr(msing), + $ omegas_e_arr(msing), omegas_i_arr(msing)) + + ALLOCATE(jacs(0:mthsurf), delpsi(0:mthsurf), + $ rsurf(0:mthsurf), asurf(0:mthsurf)) +c----------------------------------------------------------------------- +c set up kinetic profiles and equilibrium. +c [read_kin]: external, reads kinetic input file into kin spline. +c [equil_read]: external, reads equilibrium into EQUIL module. +c----------------------------------------------------------------------- + zi = 1 ! main-ion charge + zimp = 6 ! impurity charge (carbon) + mi = 2 ! main-ion mass (deuterium) + mimp = 12 ! impurity mass + nfac = 1.0 + tfac = 1.0 + wefac = 1.0 + wpfac = 1.0 + e = chag ! elementary charge [C] + chi1 = twopi * my_psio(1) ! total poloidal flux (module-level) + + CALL read_kin(infile,zi,zimp,mi,mimp,nfac, + $ tfac,wefac,wpfac,.false.) + + CALL equil_read(out_unit) + +c store full complex Deltaprime matrix in sl_in. It is converted +c to the r_s reference length after the surface loop below. + sl_in%dp_matrix(:,:) = CMPLX(dp_mat(:,:,1), dp_mat(:,:,2), + $ KIND=r8) + + dpsi_h = 0.002_r8 ! finite-difference step for Jacobian + +c----------------------------------------------------------------------- +c loop over singular surfaces: evaluate kinetic/equilibrium +c quantities via spline interpolation, compute derived layer +c parameters via params(), and populate sl_in arrays. +c----------------------------------------------------------------------- + DO ising = 1, msing + + respsi = psi_n_rational(ising) ! normalised psi + firstsurf = .TRUE. + unitfun = 1 + +c compute flux-surface-averaged minor radius + a_surf = issurfint(unitfun,mthsurf,respsi,3,1, + $ fsave,psave,jacs,delpsi,rsurf,asurf,firstsurf) + +c compute Jacobian da_surf/dpsiN by central difference + a_surf_p = issurfint(unitfun,mthsurf, + $ MIN(respsi+dpsi_h, 1.0_r8),3,1, + $ fsave,psave,jacs,delpsi,rsurf,asurf,firstsurf) + a_surf_m = issurfint(unitfun,mthsurf, + $ MAX(respsi-dpsi_h, 0.001_r8),3,1, + $ fsave,psave,jacs,delpsi,rsurf,asurf,firstsurf) + da_dpsiN = (a_surf_p - a_surf_m) + $ / (MIN(respsi+dpsi_h, 1.0_r8) + $ - MAX(respsi-dpsi_h, 0.001_r8)) + +c----------------------------------------------------------------------- +c evaluate kinetic splines at this surface. +c [spline_eval]: external, evaluates kin spline at respsi. +c kin%f(1..5) = n_i, n_e, t_i, t_e, omega (SI units) +c kin%f1(1..5) = d/d(psi_n) of the above +c----------------------------------------------------------------------- + CALL spline_eval(kin,respsi,1) + +c diamagnetic frequencies (rad/s) from GPEC kinetic splines. +c omega_* = k_theta*v_* = (m/r)*(dp/dr)/(e*n*B). In flux +c coordinates dpsi/dr = r*B_phi/q absorbs the q of m = n*q and +c leaves omega_* = n*(dp/dpsi)/(e*n_e): the toroidal mode +c number multiplies the n = 1 expression. + omega_e = twopi*kin%f(4)*kin%f1(2)/(e*chi1*kin%f(2)) + $ +twopi*kin%f1(4)/(e*chi1) + omega_i = -twopi*kin%f(3)*kin%f1(1)/(e*zi*chi1*kin%f(1)) + $ -twopi*kin%f1(3)/(e*zi*chi1) + omega_e = REAL(nn,r8)*omega_e + omega_i = REAL(nn,r8)*omega_i + sl_in%omegas_e_arr(ising) = omega_e + sl_in%omegas_i_arr(ising) = omega_i + +c extract local plasma quantities from spline + n_e = kin%f(2) + t_e = kin%f(4) / e ! convert J -> eV + n_i = kin%f(1) + t_i = kin%f(3) / e + +c Z_eff and the ion mass ratio come from slayer.in (zeff, mu_i), +c not from the kinetic file: its n_i/n_e ratio is often not a +c usable Z_eff. Unset (non-positive) values fall back to 2, a +c deuterium plasma with carbon impurities. + IF (zeff_in > 0.0_r8) THEN + zeff = zeff_in + ELSE + zeff = 2.0_r8 + END IF + + omega = kin%f(5) + my_qval = q_rational(ising) + my_sval = shear(ising) + dgeo_val = 0.0_r8 + my_bt = my_bt0(1) + my_rs = a_surf + R_0 = r_o(1) + IF (mu_i_in > 0.0_r8) THEN + mu_i = mu_i_in + ELSE + mu_i = 2.0_r8 ! deuterium + END IF + dr_val = dr_vals(ising) + +c convert STRIDE shear (psiN-based) to Fitzpatrick shear (r-based). +c STRIDE's shear(ising) = psi_N * d(ln q)/d(psi_N) (psi_N-based). +c Fitzpatrick's rfitzp formula needs s = (r/q) * dq/dr (r-based). +c Chain-rule: s_Fitz = s_psiN * r_s / (psi_N * da/dpsi_N). +c Typical magnitudes: s_psiN ~ 0.3-0.5 mid-plasma, while s_Fitz +c ~ 0.7-1.5 -- r-based shear is usually LARGER than psi_N-based +c because (r/psi_N) > (dr/dpsi_N) inside the mid-radius region. + s_fitz = my_sval * my_rs / (respsi * da_dpsiN) + WRITE(*,'(A,I2,A,F7.4,A,F6.3,A,F8.4,A,F8.4,A,F7.3,A,F7.3)') + $ ' [layerinputs] m=', resm(ising), + $ ' psiN=', respsi, + $ ' r_s=', my_rs, + $ ' da/dpsiN=', da_dpsiN, + $ ' s_psiN=', my_sval, + $ ' -> s_Fitz=', s_fitz + +c transport coefficients from caller-provided arrays. +c guard: arrays may be smaller than msing (e.g. from +c fixed-size namelist); reuse last element if exceeded. + i = MIN(ising, SIZE(sl_in%chi_p_arr)) + chi_s(1) = sl_in%chi_p_arr(i) + i = MIN(ising, SIZE(sl_in%chi_t_arr)) + chi_s(2) = sl_in%chi_t_arr(i) + i = MIN(ising, SIZE(sl_in%kappa_arr)) + chi_s(3) = sl_in%kappa_arr(i) + +c store local kinetic arrays (for future NetCDF diagnostic output) + ne_arr(ising) = n_e + te_arr(ising) = t_e + ni_arr(ising) = n_i + ti_arr(ising) = t_i + zeff_arr(ising) = zeff + bt_arr(ising) = my_bt + rs_arr(ising) = my_rs + R0_arr(ising) = R_0 + mu_i_arr(ising) = mu_i + + mms = resm(ising) + nns = nn + nns_arr(ising) = nn + nr = REAL(nn,r8) ! module-level mode numbers (real) + mr = REAL(mms,r8) + + l_n = 0.0 + l_t = 0.0 + +c----------------------------------------------------------------------- +c compute derived layer parameters using Fitzpatrick (r-based) +c shear. params() sees s_fitz as `sval`, so lu, tauk, D_norm, +c dc_tmp, etc. are all Fitzpatrick-consistent. Gradient lengths +c are zero here; params() skips its own omega_e/omega_i and we +c set Q_e/Q_i below from the spline-derived frequencies. +c----------------------------------------------------------------------- + CALL params(n_e,t_e,t_i,omega,chi_s,dr_val,dgeo_val, + $ l_n,l_t,my_qval,s_fitz,my_bt,my_rs,R_0,mu_i, + $ zeff,.false.) + +c----------------------------------------------------------------------- +c reference-length conversion of the outer Deltaprime. +c STRIDE's Deltaprime is referenced to unit normalised poloidal +c flux, while the layer response S^(1/3)*Delta(Q), S and the +c critical Delta are all referenced to x = (r - r_s)/r_s. +c Near the surface the outer solution scales as +c |x|^(1/2 -+ alpha) with the Mercier exponent +c alpha = sqrt(-D_I) (Glasser, Greene & Johnson 1975), so +c rescaling x_psi = K*x with K = r_s*dpsi_N/dr gives +c Dp_rs(i,j) = K_i^(1/2+alpha_i) * Dp(i,j) * K_j^(alpha_j-1/2) +c whose diagonal factor is K^(2*alpha), reducing to the slab +c r_s*Deltaprime at D_I = -1/4. Mercier-unstable surfaces +c (D_I >= 0) are left unconverted through alpha = 0, and K = 1 +c is used if da/dpsi_N is not a usable positive number. +c----------------------------------------------------------------------- + IF (da_dpsiN > 0.0_r8 .AND. da_dpsiN < HUGE(da_dpsiN)) THEN + k_ref = my_rs / da_dpsiN + ELSE + k_ref = 1.0_r8 + WRITE(*,*) 'WARNING: da/dpsi_N =',da_dpsiN,' at psi_N =', + $ respsi,' is not usable;' + WRITE(*,*) ' Delta_prime is left unconverted there.' + END IF + IF (di_vals(ising) < 0.0_r8) THEN + alpha_m = SQRT(-di_vals(ising)) + ELSE + alpha_m = 0.0_r8 + END IF + sl_in%k_ref_arr(ising) = k_ref + sl_in%alpha_arr(ising) = alpha_m + sl_in%Re_dp_psi_arr(ising) = Re_dp_diagonal(ising) + sl_in%Im_dp_psi_arr(ising) = Im_dp_diagonal(ising) + Re_dp_diagonal(ising) = Re_dp_diagonal(ising) + $ * k_ref**(2.0_r8*alpha_m) + Im_dp_diagonal(ising) = Im_dp_diagonal(ising) + $ * k_ref**(2.0_r8*alpha_m) + +c growth-rate conversion factor: Deltaprime -> gamma + gammafac = (my_rs * Re_dp_diagonal(ising)) / tau_r + +c----------------------------------------------------------------------- +c populate sl_in for this surface from params() globals. +c Q_e/Q_i use spline-derived omega_e/omega_i normalised by +c tauk (= Fitzpatrick S^(1/3) * tau_H, from params()). +c----------------------------------------------------------------------- + sl_in%qval_arr(ising) = resm(ising) + sl_in%q_arr(ising) = my_qval + sl_in%lu_arr(ising) = lu +c Use the same -tauk prefactor for both so opposite-signed +c spline omega_e, omega_i (electron vs ion diamagnetic +c frequencies) produce opposite-signed Q_e, Q_i. Matches +c params.f convention: Q_e = -Qconv*omega_e, Q_i = -Qconv*omega_i. +c Using +tauk on Q_i cancels the physical sign and drives +c Q_e - Q_i -> 0 for ni=ne, Ti=Te (iota_e blows up). + sl_in%Q_e_arr(ising) = -tauk * omega_e + sl_in%Q_i_arr(ising) = -tauk * omega_i + sl_in%c_beta_arr(ising) = c_beta + sl_in%d_beta_arr(ising) = d_beta +c D_norm uses iota_e from the Q_e, Q_i just set above + sl_in%D_fac_arr(ising) = D_fac + sl_in%D_norm_arr(ising) = slayer_d_norm(D_fac, + $ sl_in%Q_e_arr(ising),sl_in%Q_i_arr(ising)) + sl_in%tau_arr(ising) = tau + sl_in%omegas_arr(ising) = omega + sl_in%psi_n_arr(ising) = respsi + sl_in%gammafac_arr(ising)= gammafac + sl_in%Re_dp_arr(ising) = Re_dp_diagonal(ising) + sl_in%Im_dp_arr(ising) = Im_dp_diagonal(ising) + sl_in%d_crit_arr(ising) = dc_tmp + sl_in%P_perp_arr(ising) = P_perp + sl_in%P_tor_arr(ising) = P_tor + sl_in%Qconv_arr(ising) = tauk +c the shear, Cole ds and unperturbed ExB Q as params() formed +c them, for the torque-balance (br_th_flag) path + sl_in%shear_arr(ising) = s_fitz + sl_in%ds_arr(ising) = ds + sl_in%Qexb_arr(ising) = REAL(Q,r8) + ENDDO + +c convert the full Deltaprime matrix to the r_s reference length + DO jsing = 1, msing + DO ising = 1, msing + sl_in%dp_matrix(ising,jsing) = + $ sl_in%k_ref_arr(ising) + $ **(0.5_r8 + sl_in%alpha_arr(ising)) + $ * sl_in%dp_matrix(ising,jsing) + $ * sl_in%k_ref_arr(jsing) + $ **(sl_in%alpha_arr(jsing) - 0.5_r8) + END DO + END DO +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + + END SUBROUTINE build_inputs + + END MODULE layerinputs_mod \ No newline at end of file diff --git a/slayer/params.f b/slayer/params.f index 2457fef3d..140bcf093 100644 --- a/slayer/params.f +++ b/slayer/params.f @@ -1,76 +1,337 @@ +c======================================================================= +c MODULE params_mod +c +c Computes the normalised layer-physics parameters needed by the +c SLAYER dispersion solver from dimensional equilibrium and +c kinetic-profile inputs. +c +c Subprograms contained: +c 1. params -- derive all normalised quantities (Q, ds, c_beta, +c D_norm, P_perp, P_tor, lu, delta_crit, ...) for a +c single rational surface and store them in sglobal_mod. +c 2. slayer_d_norm -- D_norm from D_fac and the diamagnetic +c frequencies. +c======================================================================= MODULE params_mod - USE sglobal_mod + USE sglobal_mod ! SLAYER global scalars, types, constants IMPLICIT NONE CONTAINS c----------------------------------------------------------------------- -c calculate parameters. +c subprogram 1. params. +c Compute all derived layer-physics parameters for a single +c rational surface from dimensional equilibrium / kinetic inputs. +c Results are written to module-level variables in sglobal_mod +c (tau, tau_r, tauk, lu, Q, Q_e, Q_i, ds, c_beta, d_beta, +c D_norm, P_perp, P_tor, delta_n, dc_tmp, eta, visc, rho_s, ...). +c +c Note: `pr` (magnetic Prandtl number) is read from sglobal_mod +c and must be set by the caller (e.g. pr = inpr from namelist). +c If pr == 0, tau_v and visc are set to zero (no viscosity). c----------------------------------------------------------------------- - SUBROUTINE params(n_e,t_e,t_i,omega, + SUBROUTINE params(n_e,t_e,t_i,omega,chis,dr_val,dgeo_val, $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) - REAL(r8), INTENT(IN) :: n_e,t_e,t_i,omega, - $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff +c --- arguments + REAL(r8), INTENT(IN) :: n_e ! electron density [m^-3] + REAL(r8), INTENT(IN) :: t_e ! electron temperature [eV] + REAL(r8), INTENT(IN) :: t_i ! ion temperature [eV] + REAL(r8), INTENT(IN) :: omega ! toroidal rotation [rad/s] + REAL(r8), INTENT(IN) :: dr_val ! radial width dr at surface + REAL(r8), INTENT(IN) :: dgeo_val ! unused; kept for the interface + REAL(r8), INTENT(IN) :: l_n ! density gradient length + REAL(r8), INTENT(IN) :: l_t ! temperature gradient length + REAL(r8), INTENT(IN) :: qval ! safety factor + REAL(r8), INTENT(IN) :: sval ! magnetic shear + REAL(r8), INTENT(IN) :: bt ! toroidal field [T] + REAL(r8), INTENT(IN) :: rs ! minor radius [m] + REAL(r8), INTENT(IN) :: R0 ! major radius [m] + REAL(r8), INTENT(IN) :: mu_i ! ion mass ratio to proton + REAL(r8), INTENT(IN) :: zeff ! effective charge + REAL(r8), DIMENSION(3), INTENT(IN) :: chis + ! (1) chi_perp [m^2/s] + ! (2) chi_tor [m^2/s] + ! (3) kappa [m^2/s] + LOGICAL, INTENT(IN) :: params_check ! .TRUE. = print diagnostics - LOGICAL, INTENT(IN) :: params_check +c --- local variables: basic plasma + REAL(r8) :: rho ! mass density [kg/m^3] + REAL(r8) :: b_l ! characteristic magnetic field [T] + REAL(r8) :: v_a ! Alfvén velocity [m/s] + REAL(r8) :: lbeta ! local beta-related quantity +c --- local variables: electron collision time (Braginskii) + REAL(r8) :: tau_ee_num ! numerator of tau_ee formula + REAL(r8) :: tau_ee_denom ! denominator of tau_ee formula + REAL(r8) :: tau_ee ! electron-electron collision time [s] +c --- local variables: parallel conductivity (Spitzer-Härm) + REAL(r8) :: sigma_par_1 ! neoclassical correction factor + REAL(r8) :: sigma_par_2 ! classical conductivity [1/(Ohm*m)] + REAL(r8) :: sigma_par ! parallel conductivity [1/(Ohm*m)] +c --- local variables: timescales + REAL(r8) :: tau_i ! ion collision time [s] + REAL(r8) :: tau_h ! Alfvén transit time [s] + REAL(r8) :: tau_v ! viscous time [s] + REAL(r8) :: tau_tor ! toroidal diffusion time [s] + REAL(r8) :: tau_perp ! perp. diffusion time [s] +c --- local variables: delta_crit iteration + REAL(r8) :: vte ! thermal electron speed [m/s] + REAL(r8) :: chi_par_smfp ! chi_par in short mfp limit + REAL(r8) :: chi_par_lmfp ! chi_par in long mfp limit + REAL(r8) :: chi_par ! effective parallel thermal cond. + REAL(r8) :: Wd ! magnetic island width proxy + REAL(r8) :: Wd_new ! updated Wd for convergence check + INTEGER :: wit ! iteration counter +c --- local variables: unused intermediates + REAL(r8) :: K_val, Csq ! kappa/eta and composite quantity - REAL(r8) :: rho,b_l,v_a,Qconv, - $ lbeta,tau_i,tau_h,tau_r,tau_v +c----------------------------------------------------------------------- +c stop on inputs that would give NaN or Inf below. +c----------------------------------------------------------------------- + IF (.NOT. (n_e > 0.0_r8 .AND. t_e > 0.0_r8 .AND. + $ t_i > 0.0_r8)) THEN + WRITE(*,*) 'ERROR: params needs positive n_e, t_e and t_i,' + WRITE(*,*) ' got n_e =',n_e,' [m^-3], t_e =',t_e, + $ ' [eV], t_i =',t_i,' [eV]' + WRITE(*,*) ' at q =',qval,'. Check the kinetic profiles.' + STOP 'params: non-positive density or temperature' + END IF + IF (sval == 0.0_r8) THEN + WRITE(*,*) 'ERROR: params needs a nonzero magnetic shear,' + WRITE(*,*) ' got sval = 0 at q =',qval,'. The Alfven time' + WRITE(*,*) ' and the layer widths are undefined.' + STOP 'params: zero magnetic shear' + END IF +c----------------------------------------------------------------------- +c Coulomb logarithm, basic plasma quantities, and Spitzer +c resistivity. +c----------------------------------------------------------------------- + lnLamb = 24 + 3.0_r8*LOG(10.0_r8) - 0.5_r8*LOG(n_e) + LOG(t_e) + + tau = t_i / t_e ! T_i / T_e + tau_i = 6.6e17_r8*mu_i**0.5_r8*(t_i/1e3_r8)**1.5_r8 + $ / (n_e*lnLamb) ! ion collision time [s] + eta = 1.65e-9_r8*lnLamb / (t_e/1e3_r8)**1.5_r8 ! Spitzer resistivity (Wesson) + rho = (mu_i*m_p)*n_e ! mass density [kg/m^3] + +c----------------------------------------------------------------------- +c Electron-electron collision time (Braginskii) and Spitzer-Härm +c parallel conductivity. +c----------------------------------------------------------------------- + tau_ee_num = 6.0_r8*SQRT(2.0_r8)*(pi**1.5_r8) + $ *(eps0**2.0_r8)*(m_e**0.5_r8)*(t_e**1.5_r8) + tau_ee_denom = lnLamb*(chag**2.5_r8)*n_e + tau_ee = tau_ee_num / tau_ee_denom + + sigma_par_1 = ( SQRT(2.0_r8) + 13.0_r8*(Zeff/4.0_r8) ) + $ / (Zeff*(SQRT(2.0_r8) + Zeff)) + sigma_par_2 = (n_e * (chag**2.0_r8) * tau_ee) / m_e + sigma_par = sigma_par_1 * sigma_par_2 + +c----------------------------------------------------------------------- +c Characteristic field, Alfvén speed, length scales, and +c fundamental timescales. +c Note: mr, nr (rational m/n) and nn (toroidal mode number) are +c module-level variables set before calling params(). +c----------------------------------------------------------------------- + b_l = (nr/mr)*rs*sval*bt/R0 ! characteristic B [T] + v_a = b_l / (mu0*rho)**0.5_r8 ! Alfvén velocity [m/s] + rho_s = 1.02e-4_r8*(mu_i*t_e)**0.5_r8 / bt ! ion Larmor at T_e [m] + d_i = ( (mu_i*m_p) / (n_e*(chag**2)*mu0) )**0.5_r8 + ! ion skin depth [m] + +c The layer timescales and widths below depend on the magnitude of +c the shear, not its sign: ABS(sval) keeps them positive on a +c reverse-shear surface. + tau_h = R0*(mu0*rho)**0.5_r8 / (nn*ABS(sval)*bt) ! Alfvén time [s] + tau_r = mu0*(rs**2.0_r8)*sigma_par ! resistive time [s] (Fitzpatrick) + IF (pr > 0.0_r8) THEN + tau_v = tau_r / pr ! viscous time [s] + visc = rho*rs**2.0_r8 / tau_v + ELSE + tau_v = 0.0_r8 + visc = 0.0_r8 + END IF + +c Lundquist number + lu = tau_r / tau_h + +c----------------------------------------------------------------------- +c Diamagnetic frequencies and normalised Q parameters. +c Qconv converts dimensional frequencies to the normalised Q +c used in the SLAYER dispersion relation (Cole scaling). +c----------------------------------------------------------------------- + Qconv = lu**(1.0_r8/3.0_r8) * tau_h ! frequency normalisation (Cole) + tauk = Qconv ! stored in sglobal_mod +c the layer sees the mode at the Doppler frequency n*omega_E, as in +c the coupled determinant's shift (slayer.f) and gpec_slayer + Q = Qconv * REAL(nn,r8) * omega ! normalised ExB frequency + +c Diamagnetic frequencies require gradient lengths; skip when the +c caller does not provide them (l_n <= 0 or l_t <= 0). Callers that +c compute omega_e/omega_i by other means must set Q_e and Q_i +c themselves afterwards. +c omega_* = k_theta*v_* carries the toroidal mode number through +c m = n*q: the bare gradient-length expression is multiplied by +c n, as in build_inputs (layerinputs.f). + IF (l_n > 0.0_r8 .AND. l_t > 0.0_r8) THEN + omega_e = -REAL(nn,r8)*t_e/(bt*R0) + $ *(1.0_r8/l_n + 1.0_r8/l_t)*qval + omega_i = REAL(nn,r8)*t_i/(bt*R0) + $ *(1.0_r8/l_n + 1.0_r8/l_t)*qval + Q_e = -Qconv * omega_e + Q_i = -Qconv * omega_i + END IF - ! mu_i: ion mass ratio to proton - tau= t_i/t_e ! ratio of ion to electron temperature - tau_i = 6.6e17*mu_i**0.5*(t_i/1e3)**1.5/(n_e*lnLamb) ! ion colls. - eta= 1.65e-9*lnLamb/(t_e/1e3)**1.5 ! spitzer resistivity (wesson) - rho=(mu_i*m_p)*n_e ! mass density +c normalised ion Larmor radius (critical stability parameter) + ds = lu**(1.0_r8/3.0_r8) * rho_s / rs - b_l=(nr/mr)*rs*sval*bt/R0 ! characteristic magnetic field - v_a=b_l/(mu0*rho)**0.5 ! alfven velocity - rho_s=1.02e-4*(mu_i*t_e)**0.5/bt ! ion Lamour by elec. Temp. +c----------------------------------------------------------------------- +c Plasma beta and Prandtl-number-like transport ratios. +c----------------------------------------------------------------------- + lbeta = (5.0_r8/3.0_r8)*mu0*n_e*chag*(t_e+t_i) / bt**2.0_r8 + c_beta = (lbeta / (1.0_r8+lbeta))**0.5_r8 - tau_h=R0*(mu0*rho)**0.5/(nn*sval*bt) ! alfven time across surface - tau_r=mu0*rs**2.0/eta ! resistive time scale - tau_v=tau_r/pr !rho*rs**2.0/visc ! viscous time scale +c kappa-based P_perp (JKP's definition, to implement later if desired) +c K_val = chis(3) / eta +c Csq = c_beta**2.0 + (1.0 - c_beta**2.0)*K_val - ! this one must be anomalous. calculated back from pr. - visc= rho*rs**2.0/tau_v +c IF (ABS(Csq) > 0.0) THEN +c P_perp = Csq +c ELSE +c tau_perp = (rs**2.0) / chis(1) +c END IF - lu=tau_r/tau_h ! Lundquist number -! pr=tau_r/tau_v ! Prandtl number. Only place needs viscosity. +c effective perpendicular and toroidal Prandtl numbers + tau_perp = (rs**2.0_r8) / chis(1) + P_perp = tau_r / tau_perp ! perp magnetic Prandtl number - omega_e=-t_e/(bt*R0)*(1.0/l_n+1.0/l_t)*qval ! elec. diamag - omega_i=t_i/(bt*R0)*(1.0/l_n+1.0/l_t)*qval ! ion diamag + tau_tor = (rs**2.0_r8) / chis(2) + P_tor = tau_r / tau_tor ! toroidal magnetic Prandtl number - ! now calculate the main 7 normalized parameters. +c----------------------------------------------------------------------- +c Normalised beta-related width and Delta norm factor. +c----------------------------------------------------------------------- +c d_beta uses Fitzpatrick's definition (tau' form). +c D = S^(1/3) * iota_e^(1/2) * d_beta/r_s with +c iota_e = Q_e/(Q_e - Q_i), the electron share of the diamagnetic +c frequency (see slayer_d_norm). When this call does not set Q_e +c and Q_i, the equal-gradient-ratio value iota_e = T_e/(T_e + T_i) +c = 1/(1 + tau) is used provisionally; callers that set Q_e and Q_i +c afterwards must reset D_norm = slayer_d_norm(D_fac,Q_e,Q_i). + d_beta = c_beta * d_i + D_fac = (d_beta/rs) * lu**(1.0_r8/3.0_r8) + IF (l_n > 0.0_r8 .AND. l_t > 0.0_r8) THEN + D_norm = slayer_d_norm(D_fac,Q_e,Q_i) + ELSE + D_norm = D_fac * SQRT(1.0_r8/(1.0_r8+tau)) + END IF - Qconv=lu**(1.0/3.0)*tau_h ! conversion to Qs based on Cole + delta_n = lu**(1.0_r8/3.0_r8) / rs ! normalisation for Delta values - ! note Q depends on Qconv even if omega is fixed. - Q=Qconv*omega - Q_e=-Qconv*omega_e - Q_i=-Qconv*omega_i +c----------------------------------------------------------------------- +c Critical Deltaprime (dc_tmp) via iterative chi_parallel calculation. +c The island-width Wd is iterated to converge the +c short-mfp / long-mfp interpolation for chi_parallel. +c dc_type (from sglobal_mod) selects the formula: +c 'rfitzp' -- R. Fitzpatrick +c default -- dc_tmp = 0 +c----------------------------------------------------------------------- + IF (ABS(dr_val) > 0.0_r8) THEN + + IF (.NOT. (chis(1) > 0.0_r8)) THEN + WRITE(*,*) 'ERROR: the critical Delta-prime needs' + WRITE(*,*) ' chi_perp > 0, got chi_perp =',chis(1), + $ ' [m^2/s] at q =',qval,'.' + WRITE(*,*) ' Set chi_p_prof in slayer.in.' + STOP 'params: chi_perp must be positive' + END IF + vte = SQRT((2.0_r8*(t_e*chag)) / m_e) + chi_par_smfp = (1.581_r8*tau_ee*(vte**2.0_r8)) + $ / (1.0_r8 + 0.2535_r8*Zeff) - ! This is the most critical parameter - ds=lu**(1.0/3.0)*rho_s/rs ! conversion based on Cole. + Wd = 0.1_r8 ! initial guess + DO wit = 1, 100 + chi_par_lmfp = (2.0_r8*R0*vte) + $ / (SQRT(pi)*nr*ABS(sval)*Wd) + chi_par = (chi_par_smfp*chi_par_lmfp) + $ / (chi_par_smfp + chi_par_lmfp) + Wd_new = SQRT(8.0_r8)*((chis(1)/chi_par)**0.25_r8) + $ * (1.0_r8/SQRT((rs/R0)*ABS(sval)*nr)) + IF (ABS(Wd_new - Wd) / MAX(ABS(Wd), 1.0e-30_r8) + $ < 1.0e-10_r8) THEN + Wd = Wd_new + EXIT + END IF + Wd = Wd_new + END DO + IF (wit > 100) THEN + WRITE(*,*) 'params: Wd iteration failed to converge' + STOP 'params: Wd iteration did not converge' + END IF - lbeta=(5.0/3.0)*mu0*n_e*chag*(t_e+t_i)/bt**2.0 - c_beta=(lbeta/(1.0+lbeta))**0.5 + SELECT CASE(dc_type) + CASE('rfitzp') + dc_tmp = -(SQRT(2.0_r8)*(pi**(1.5_r8))*dr_val) / Wd + CASE default + dc_tmp = 0.0_r8 + END SELECT - delta_n=lu**(1.0/3.0)/rs ! norm factor for delta primes + ELSE + dc_tmp = 0.0_r8 + END IF - ! quick diagnostics. +c----------------------------------------------------------------------- +c optional diagnostics (guarded by params_check flag). +c----------------------------------------------------------------------- IF (params_check) THEN - WRITE(*,*)"eta=",eta - WRITE(*,*)"S=",lu - WRITE(*,*)"Q=",Q - WRITE(*,*)"Q_e=",Q_e - WRITE(*,*)"Q_i=",Q_i - WRITE(*,*)"ds=",ds - WRITE(*,*)"c_beta=",c_beta + WRITE(*,*) 'eta = ', eta + WRITE(*,*) 'S = ', lu + WRITE(*,*) 'Q = ', Q + WRITE(*,*) 'Q_e = ', Q_e + WRITE(*,*) 'Q_i = ', Q_i + WRITE(*,*) 'ds = ', ds + WRITE(*,*) 'c_beta = ', c_beta ENDIF - +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- RETURN END SUBROUTINE params +c----------------------------------------------------------------------- +c subprogram 2. slayer_d_norm. +c Normalised ion-sound-radius parameter of the Fitzpatrick layer +c equations, D = S^(1/3) * iota_e^(1/2) * d_beta/r_s, from +c D_fac = S^(1/3)*d_beta/r_s and iota_e = Q_e/(Q_e - Q_i). +c iota_e must be positive, which holds when the electron and ion +c diamagnetic frequencies have opposite signs. +c----------------------------------------------------------------------- + FUNCTION slayer_d_norm(D_fac_in,Q_e_in,Q_i_in) + + REAL(r8), INTENT(IN) :: D_fac_in ! S^(1/3)*d_beta/r_s + REAL(r8), INTENT(IN) :: Q_e_in ! normalised electron diamag. + REAL(r8), INTENT(IN) :: Q_i_in ! normalised ion diamag. + REAL(r8) :: slayer_d_norm + REAL(r8) :: iota ! Q_e/(Q_e - Q_i) + + IF (Q_e_in == Q_i_in) THEN + WRITE(*,*) 'ERROR: Q_e = Q_i =',Q_e_in, + $ ', so iota_e = Q_e/(Q_e - Q_i) is undefined.' + STOP 'slayer_d_norm: Q_e equals Q_i' + END IF + iota = Q_e_in / (Q_e_in - Q_i_in) + IF (.NOT. (iota > 0.0_r8)) THEN + WRITE(*,*) 'ERROR: iota_e = Q_e/(Q_e - Q_i) =',iota, + $ ' is not positive (Q_e =',Q_e_in,', Q_i =',Q_i_in,').' + WRITE(*,*) ' D_norm = S^(1/3)*sqrt(iota_e)*d_beta/r_s is' + WRITE(*,*) ' undefined. The electron and ion diamagnetic' + WRITE(*,*) ' frequencies are expected to have opposite signs.' + STOP 'slayer_d_norm: iota_e is not positive' + END IF + slayer_d_norm = D_fac_in * SQRT(iota) + + RETURN + END FUNCTION slayer_d_norm - END MODULE params_mod + END MODULE params_mod \ No newline at end of file diff --git a/slayer/sglobal.f b/slayer/sglobal.f index d2a96e36b..1abd62098 100644 --- a/slayer/sglobal.f +++ b/slayer/sglobal.f @@ -1,25 +1,226 @@ +c======================================================================= +c MODULE sglobal_mod +c +c Global shared state for the SLAYER layer-physics package. +c Contains: +c - physical and mathematical constants +c - module-level scalar variables written by params() and +c consumed by the dispersion solvers (delta, slayer, gslayer) +c - AMR scanner storage (hash-based v1, cell-based v2) +c - derived types: slayer_inputs_type, slayer_outputs_type, +c deltas_outputs_type, result_type, amr_cell_type +c +c TODO: `pr` and `pe` (magnetic Prandtl number and its +c electron analogue) are declared but never explicitly set +c within this module or in params(). They must be initialised +c before params() is called (tau_v = tau_r / pr). Consider +c adding a dedicated init routine or making them INTENT(IN) +c arguments of params(). +c +c======================================================================= MODULE sglobal_mod - USE local_mod, ONLY: r8 + + USE local_mod, ONLY: r8 ! double-precision kind parameter IMPLICIT NONE - INTEGER :: mm,nn - INTEGER :: in_unit,out_unit,out2_unit,out3_unit, - $ bin_unit,bin_2d_unit,input_unit - REAL(r8) :: mr,nr - REAL(r8) :: Q_e,Q_i,pr,pe,c_beta,ds,tau - REAL(r8) :: eta,visc,rho_s,lu,omega_e,omega_i, - $ delta_n,layfac - COMPLEX(r8) :: Q +c----------------------------------------------------------------------- +c I/O unit numbers and trace control. +c----------------------------------------------------------------------- + INTEGER :: in_unit ! standard input unit + INTEGER :: out_unit ! primary output unit + INTEGER :: out2_unit ! secondary output unit + INTEGER :: out3_unit ! tertiary output unit + INTEGER :: bin_unit ! binary output unit + INTEGER :: bin_2d_unit ! 2-D binary output unit + INTEGER :: input_unit ! namelist / input unit + INTEGER :: n_trace ! trace / debug verbosity level + +c----------------------------------------------------------------------- +c Mode numbers (integer and real representations). +c----------------------------------------------------------------------- + INTEGER :: mm ! poloidal mode number (integer) + INTEGER :: nn ! toroidal mode number (integer) + REAL(r8) :: mr ! poloidal mode number (real copy) + REAL(r8) :: nr ! toroidal mode number (real copy) + CHARACTER(2) :: sn_str ! toroidal n as string + CHARACTER(2) :: sm_str ! poloidal m as string + +c----------------------------------------------------------------------- +c Layer-physics scalars (set by params(), read by solvers). +c These are the normalised parameters that enter the SLAYER +c dispersion relation. +c----------------------------------------------------------------------- +c --- temperature / collisionality + REAL(r8) :: tau ! T_i / T_e + REAL(r8) :: eta ! Spitzer resistivity [Ohm*m] + REAL(r8) :: visc ! anomalous viscosity [m^2/s] + REAL(r8) :: lnLamb = 17.0_r8 ! Coulomb log, set by params() +c --- length scales + REAL(r8) :: rho_s ! ion Larmor radius at T_e [m] + REAL(r8) :: d_i ! ion skin depth [m] + REAL(r8) :: d_beta ! beta-weighted ion scale d_beta = c_beta * d_i +c --- timescales + REAL(r8) :: tau_r ! resistive diffusion time [s] + REAL(r8) :: tauk ! Q-conversion factor (= Qconv) +c --- Lundquist and Prandtl numbers + REAL(r8) :: lu ! Lundquist number S = tau_r / tau_h + REAL(r8) :: pr ! magnetic Prandtl number (TODO: see header) + REAL(r8) :: pe ! electron Prandtl number (TODO: see header) + REAL(r8) :: P_perp ! perpendicular magnetic Prandtl number + REAL(r8) :: P_tor ! toroidal magnetic Prandtl number +c --- normalised layer parameters + REAL(r8) :: ds ! normalised ion Larmor radius + REAL(r8) :: c_beta ! compressional beta parameter + REAL(r8) :: D_norm ! normalised beta-weighted ion scale + REAL(r8) :: D_fac ! D_norm / sqrt(iota_e) + REAL(r8) :: delta_n ! Delta normalisation factor + REAL(r8) :: Qconv ! frequency normalisation (Cole) +c --- diamagnetic and rotation frequencies + REAL(r8) :: omega_e ! electron diamagnetic frequency [rad/s] + REAL(r8) :: omega_i ! ion diamagnetic frequency [rad/s] + COMPLEX(r8) :: Q ! normalised complex rotation frequency + REAL(r8) :: Q_e ! normalised electron diamagnetic Q + REAL(r8) :: Q_i ! normalised ion diamagnetic Q +c --- stability / Delta_crit + REAL(r8) :: dc_tmp ! computed Delta_crit + COMPLEX(r8) :: delta_eff ! (Deltaprime - Delta_crit)/S^(1/3) + CHARACTER(20) :: dc_type ! dc formula: 'rfitzp','none' +c --- solver workspace / results + COMPLEX(r8) :: g_tmp ! temporary complex growth rate + REAL(r8) :: gamma_fac ! growth-rate conversion factor +c --- miscellaneous + REAL(r8) :: iota_e ! electron iota + REAL(r8) :: layfac ! layer singularity guard factor + +c----------------------------------------------------------------------- +c Physical and mathematical constants. +c----------------------------------------------------------------------- + REAL(r8), PARAMETER :: pi = 3.1415926535897932385d0 + REAL(r8), PARAMETER :: mu0 = 4.0d-7 * pi ! vacuum permeability [H/m] + REAL(r8), PARAMETER :: m_e = 9.1094d-31 ! electron mass [kg] + REAL(r8), PARAMETER :: m_p = 1.6726d-27 ! proton mass [kg] + REAL(r8), PARAMETER :: chag = 1.6021917d-19 ! elementary charge [C] + REAL(r8), PARAMETER :: kval = 1.3807d-23 ! Boltzmann constant [J/K] + REAL(r8), PARAMETER :: eps0 = 8.8542d-12 ! vacuum permittivity [F/m] + COMPLEX(r8), PARAMETER :: ifac = (0.0d0, 1.0d0) ! imaginary unit + +c----------------------------------------------------------------------- +c AMR scanner storage -- hash-based deduplication (v1). +c----------------------------------------------------------------------- + INTEGER, PARAMETER :: MAX_PTS = 1000000 ! max unique eval points + INTEGER, PARAMETER :: HASH_SZ = 1000003 ! hash table size (prime) + REAL(r8), PARAMETER :: HASH_SCALE = 1.0d5 ! Re/Im quantisation scale + INTEGER, ALLOCATABLE :: hash_head(:) ! bucket heads (HASH_SZ) + INTEGER, ALLOCATABLE :: hash_next(:) ! chain pointers (MAX_PTS) + +c----------------------------------------------------------------------- +c AMR scanner storage -- cell-based refinement (v2). +c----------------------------------------------------------------------- + INTEGER, PARAMETER :: MAX_CELLS = 1000000 ! max AMR cells + + TYPE :: amr_cell_type + COMPLEX(r8) :: Q(4) ! corner Q-values (BL, BR, TL, TR) + COMPLEX(r8) :: D(4) ! corner dispersion values + LOGICAL :: needs_refine ! flagged for subdivision + END TYPE amr_cell_type + + TYPE(amr_cell_type), ALLOCATABLE :: amr_cells(:) + INTEGER :: n_amr_cells ! current number of active cells + +c----------------------------------------------------------------------- +c AMR output arrays (shared by v1 and v2). +c----------------------------------------------------------------------- + COMPLEX(r8), ALLOCATABLE :: Q_store(:) ! unique Q-points + COMPLEX(r8), ALLOCATABLE :: D_store(:) ! corresponding D-values + INTEGER :: n_pts ! number of stored points + +c----------------------------------------------------------------------- +c Derived types: solver I/O and scan results. +c----------------------------------------------------------------------- + +c result_type -- torque-scan output bucket + TYPE result_type + REAL(r8), ALLOCATABLE :: inQs(:) ! Re(Q) scan values + REAL(r8), ALLOCATABLE :: iinQs(:) ! Im(Q) scan values + REAL(r8), ALLOCATABLE :: Re_deltas(:) ! Re(Delta) results + REAL(r8), ALLOCATABLE :: Im_deltas(:) ! Im(Delta) results + INTEGER :: count ! number of entries + END TYPE result_type + +c slayer_inputs_type -- per-surface input arrays for SLAYER + TYPE slayer_inputs_type + INTEGER, ALLOCATABLE :: qval_arr(:) ! resonant poloidal m (labels) + REAL(r8), ALLOCATABLE :: q_arr(:) ! safety factor m/n + REAL(r8), ALLOCATABLE :: chi_p_arr(:) ! chi_perp [m^2/s] + REAL(r8), ALLOCATABLE :: chi_t_arr(:) ! chi_tor [m^2/s] + REAL(r8), ALLOCATABLE :: kappa_arr(:) ! kappa (thermal cond.) + REAL(r8), ALLOCATABLE :: psi_n_arr(:) ! normalised psi + REAL(r8), ALLOCATABLE :: lu_arr(:) ! Lundquist number + REAL(r8), ALLOCATABLE :: Qconv_arr(:) ! Q-conversion factor + REAL(r8), ALLOCATABLE :: Q_e_arr(:) ! normalised Q_e + REAL(r8), ALLOCATABLE :: Q_i_arr(:) ! normalised Q_i + REAL(r8), ALLOCATABLE :: c_beta_arr(:) ! compressional beta + REAL(r8), ALLOCATABLE :: d_beta_arr(:) ! beta-related width + REAL(r8), ALLOCATABLE :: D_norm_arr(:) ! normalised D + REAL(r8), ALLOCATABLE :: D_fac_arr(:) ! D / sqrt(iota_e) + REAL(r8), ALLOCATABLE :: tau_arr(:) ! T_i / T_e + REAL(r8), ALLOCATABLE :: P_perp_arr(:) ! perp Prandtl number + REAL(r8), ALLOCATABLE :: P_tor_arr(:) ! toroidal Prandtl + REAL(r8), ALLOCATABLE :: omegas_arr(:) ! ExB rotation per unit n [rad/s] + REAL(r8), ALLOCATABLE :: Qexb_arr(:) ! normalised ExB Q + REAL(r8), ALLOCATABLE :: shear_arr(:) ! shear for params() + REAL(r8), ALLOCATABLE :: ds_arr(:) ! Cole normalised rho_s + REAL(r8), ALLOCATABLE :: q_shift_arr(:) ! ExB Doppler shift of Q + REAL(r8), ALLOCATABLE :: omegas_e_arr(:) ! omega_e [rad/s] + REAL(r8), ALLOCATABLE :: omegas_i_arr(:) ! omega_i [rad/s] + REAL(r8), ALLOCATABLE :: gammafac_arr(:) ! gamma conversion +c Deltaprime used in the layer matching is referenced to the +c minor radius r_s (unit x = (r - r_s)/r_s), like every layer +c quantity. The *_psi arrays keep the outer-region value as +c STRIDE writes it, referenced to unit normalised poloidal flux. + REAL(r8), ALLOCATABLE :: Re_dp_arr(:) ! Re(Deltaprime), r_s ref. + REAL(r8), ALLOCATABLE :: Im_dp_arr(:) ! Im(Deltaprime), r_s ref. + REAL(r8), ALLOCATABLE :: Re_dp_psi_arr(:) ! Re(Deltaprime), psi_N ref. + REAL(r8), ALLOCATABLE :: Im_dp_psi_arr(:) ! Im(Deltaprime), psi_N ref. + REAL(r8), ALLOCATABLE :: k_ref_arr(:) ! K = r_s*dpsi_N/dr + REAL(r8), ALLOCATABLE :: alpha_arr(:) ! Mercier exponent sqrt(-D_I) + REAL(r8), ALLOCATABLE :: d_crit_arr(:) ! Delta_crit + COMPLEX(r8), ALLOCATABLE :: dp_matrix(:,:) ! full Deltaprime matrix, r_s ref. + END TYPE slayer_inputs_type - REAL(r8), PARAMETER :: pi=3.1415926535897932385, mu0=4e-7*pi, - $ m_e=9.1094e-31,m_p=1.6726e-27,chag=1.6022e-19, - $ kval=1.3807e-23,lnLamb=17.0 +c slayer_outputs_type -- per-surface solver results + TYPE slayer_outputs_type + COMPLEX(r8), ALLOCATABLE :: dels_db_arr(:) ! Delta from d_beta + COMPLEX(r8), ALLOCATABLE :: gamma_sol_arr(:) ! solved growth rate + COMPLEX(r8), ALLOCATABLE :: gamma_est_arr(:) ! estimated growth rate + REAL(r8), ALLOCATABLE :: br_th_arr(:) ! Br threshold + END TYPE slayer_outputs_type - ! lnLamb will be updated. +c deltas_outputs_type -- scan output (Q vs Delta) + TYPE deltas_outputs_type + REAL(r8), ALLOCATABLE :: inQs(:) ! Re(Q) values + REAL(r8), ALLOCATABLE :: iinQs(:) ! Im(Q) values + REAL(r8), ALLOCATABLE :: real_deltas(:) ! Re(Delta) + REAL(r8), ALLOCATABLE :: imag_deltas(:) ! Im(Delta) + END TYPE deltas_outputs_type - COMPLEX(r8), PARAMETER :: ifac=(0,1) +c roots_outputs_type -- zeros and poles of the dispersion function +c located in one AMR scan. Q values are stored in the scan plane, +c q = omega_hat - i*gamma_hat (so gamma = -AIMAG(q)/tauk and +c omega = REAL(q)/tauk), sorted by descending growth rate. + INTEGER, PARAMETER :: MAX_ROOTS = 32 ! max stored per scan + TYPE roots_outputs_type + INTEGER :: n_root = 0 ! zeros found + INTEGER :: n_pole = 0 ! poles found + COMPLEX(r8) :: q_root(MAX_ROOTS) ! zeros (scan plane) + COMPLEX(r8) :: q_pole(MAX_ROOTS) ! poles (scan plane) + LOGICAL :: polished(MAX_ROOTS) ! secant converged? + REAL(r8) :: tauk = 0.0_r8 ! Q -> [1/s] factor + END TYPE roots_outputs_type - CHARACTER(2) :: sn +c per-scan roots, filled by extract_roots (growthrates_mod) and +c written by slayer_netcdf_out + TYPE(roots_outputs_type), ALLOCATABLE :: roots_out(:) END MODULE sglobal_mod diff --git a/slayer/slayer.f b/slayer/slayer.f index 090af9d51..9d199a3a3 100644 --- a/slayer/slayer.f +++ b/slayer/slayer.f @@ -1,11 +1,10 @@ c----------------------------------------------------------------------- -c Slab LAYER based on linear drift MHD -c SLAYER: main program -c----------------------------------------------------------------------- -c----------------------------------------------------------------------- -c code organization. -c----------------------------------------------------------------------- -c slayer. +c SLAYER: Slab LAYER linear drift-MHD code. +c Main driver program. +c +c Computes tearing-mode layer quantities (inner layer Delta, +c growth rates, torque balance, field thresholds) from +c slab-geometry drift-MHD matching via a Riccati integration method. c----------------------------------------------------------------------- c----------------------------------------------------------------------- c declarations. @@ -13,125 +12,340 @@ PROGRAM slayer USE sglobal_mod - USE params_mod - USE delta_mod, ONLY: riccati,riccati_out, - $ parflow_flag,PeOhmOnly_flag + USE delta_mod, ONLY: riccati,riccati_f,riccati_del_s, + $ riccati_out,parflow_flag,PeOhmOnly_flag + USE gslayer_mod + USE growthrates_mod + USE layerinputs_mod + USE ieee_arithmetic, ONLY: ieee_value, ieee_quiet_nan IMPLICIT NONE - - CHARACTER(128) :: infile - INTEGER :: i,j,k,inum,jnum,knum,inn - INTEGER, DIMENSION(1) :: index - - LOGICAL :: params_flag,QPscan_flag,QPescan_flag,QPscan2_flag, - $ QDscan2_flag,Qbscan_flag,Qscan_flag, - $ onscan_flag,otscan_flag,ntscan_flag,nbtscan_flag, - $ Pe_flag,verbose,ascii_flag,bin_flag,netcdf_flag, - $ bal_flag,stability_flag,riccatiscan_flag,input_flag, - $ params_check - - REAL(r8) :: n_e,t_e,t_i,omega,omega0, - $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff - REAL(r8) :: inQ,inQ_e,inQ_i,inpr,inpe,inc_beta,inds,intau,inlu - REAL(r8) :: psi0,jxb,Q0,Q_sol,br_th - COMPLEX(r8) :: delta,delta_n_p - - REAL(r8) :: inQ_min,inQ_max,j_min,j_max,jpower,k_min,k_max,kpower, - $ Qratio - - INTEGER, DIMENSION(:), ALLOCATABLE :: mms,nns - - REAL(r8), DIMENSION(8) :: inpr_prof - - REAL(r8), DIMENSION(:), ALLOCATABLE :: inQs,iinQs,jxbl,bal, - $ prs,n_es,t_es,t_is,omegas,l_ns,l_ts,qvals,svals, - $ bts,rss,R0s,mu_is,zeffs,Q_soll,br_thl - REAL(r8), DIMENSION(:,:), ALLOCATABLE :: - $ js,ks,psis,jxbs,Q_sols,br_ths +c----------------------------------------------------------------------- +c local scalars — loop indices and counters. +c----------------------------------------------------------------------- + INTEGER :: i,j,k ! general loop indices + INTEGER :: inn ! number of input-file surfaces + INTEGER :: count ! generic counter + INTEGER :: nbad ! namelist inputs left unset +c----------------------------------------------------------------------- +c local scalars — numerical resolution. +c----------------------------------------------------------------------- + INTEGER :: inum ! resolution for 1-D scans + INTEGER :: jnum ! resolution for 2-D scan axis 1 + INTEGER :: knum ! resolution for 2-D scan axis 2 + INTEGER :: Q_num ! resolution for stab. scan Re(Q) + INTEGER :: msing_max ! max number of singular surfaces + INTEGER :: n_k ! number of rational surfaces +c----------------------------------------------------------------------- +c local scalars — MAXLOC result holder. +c----------------------------------------------------------------------- + INTEGER, DIMENSION(1) :: iloc ! result of MAXLOC +c----------------------------------------------------------------------- +c control flags — workflow. +c----------------------------------------------------------------------- + LOGICAL :: params_flag ! compute params from kinetic data + LOGICAL :: input_flag ! read multi-surface input file + LOGICAL :: read_eq ! read equilibrium files + LOGICAL :: verbose ! enable progress messages + LOGICAL :: params_check ! print diagnostic output in params +c----------------------------------------------------------------------- +c control flags — physics modes. +c----------------------------------------------------------------------- + LOGICAL :: est_gamma_flag ! estimate growth rate + LOGICAL :: match_gamma_flag ! asymptotically matched gamma + LOGICAL :: Pperp_Ptor_flag ! use Fitzpatrick layer model + LOGICAL :: coupling_flag ! coupled rational surfaces + LOGICAL :: rotation_flag ! ExB Doppler shift when coupled + LOGICAL :: br_th_flag ! Br threshold test scan + LOGICAL :: bal_flag ! torque balance scan + LOGICAL :: stability_flag ! complex-Q delta scan + LOGICAL :: Pe_flag ! include electron pressure +c----------------------------------------------------------------------- +c control flags — parameter-space scans. +c----------------------------------------------------------------------- + LOGICAL :: QPscan_flag ! (Q,P) scan + LOGICAL :: QPescan_flag ! (Q,Pe) scan + LOGICAL :: QPscan2_flag ! (Q,P) scan variant 2 + LOGICAL :: QDscan2_flag ! (Q,D) scan variant 2 + LOGICAL :: Qbscan_flag ! (Q,beta) scan + LOGICAL :: Qscan_flag ! 1-D Q scan + LOGICAL :: onscan_flag ! (omega,n) scan + LOGICAL :: otscan_flag ! (omega,T) scan + LOGICAL :: ntscan_flag ! (n,T) scan + LOGICAL :: nbtscan_flag ! (n,Bt) scan + LOGICAL :: riccatiscan_flag ! Riccati-variable scan + LOGICAL :: stabscan_flag ! stability scan (single surface) + LOGICAL :: coupled_stabscan_flag ! stability scan (coupled) + LOGICAL :: amr_flag ! adaptive mesh refinement scan +c----------------------------------------------------------------------- +c control flags — output format. +c----------------------------------------------------------------------- + LOGICAL :: ascii_flag ! write ASCII output files + LOGICAL :: bin_flag ! write binary output files + LOGICAL :: netcdf_flag ! write NetCDF output files +c----------------------------------------------------------------------- +c local scalars — physical input quantities. +c----------------------------------------------------------------------- + REAL(r8) :: n_e ! electron density [m^-3] + REAL(r8) :: t_e ! electron temperature [eV] + REAL(r8) :: t_i ! ion temperature [eV] + REAL(r8) :: omega ! toroidal rotation [rad/s] + REAL(r8) :: omega0 ! unused + REAL(r8) :: l_n ! density gradient scale length + REAL(r8) :: l_t ! temperature gradient scale length + REAL(r8) :: qval ! safety factor at surface + REAL(r8) :: sval ! magnetic shear at surface + REAL(r8) :: bt ! toroidal field [T] + REAL(r8) :: rs ! minor radius of surface [m] + REAL(r8) :: R0 ! major radius [m] + REAL(r8) :: mu_i ! ion mass number + REAL(r8) :: zeff ! effective charge + REAL(r8) :: dr_val ! radial derivative parameter + REAL(r8) :: dgeo_val ! unused; still read so old inputs parse + REAL(r8) :: scan_width ! half-width of complex-Q scan + REAL(r8) :: scan_center ! Re(Q) centre of complex-Q scan + REAL(r8) :: q_lo,q_hi,q_br ! expected coupled-branch Re(Q) +c----------------------------------------------------------------------- +c local scalars — normalized layer parameters (namelist overrides). +c----------------------------------------------------------------------- + REAL(r8) :: inQ ! normalized ExB rotation freq. + REAL(r8) :: inQ_e ! normalized electron diamagnetic + REAL(r8) :: inQ_i ! normalized ion diamagnetic + REAL(r8) :: inQ_e_ovr(8) ! per-surface Q_e overrides (read_eq) + REAL(r8) :: inQ_i_ovr(8) ! per-surface Q_i overrides (read_eq) + REAL(r8) :: inpr ! Prandtl number + REAL(r8), DIMENSION(8) :: inpr_prof ! read for OMFIT; unused + REAL(r8) :: inpe ! normalized electron pressure + REAL(r8) :: inc_beta ! normalized beta + REAL(r8) :: inds ! normalized D (magnetic diffusion) + REAL(r8) :: intau ! normalized tau = T_i/T_e + REAL(r8) :: inlu ! Lundquist number (params_flag=f) +c----------------------------------------------------------------------- +c local scalars — derived / scratch quantities. +c----------------------------------------------------------------------- + REAL(r8) :: psi0 ! reconnected flux (a.u.) + REAL(r8) :: jxb ! j x B torque (a.u.) + REAL(r8) :: Q0 ! unperturbed rotation frequency + REAL(r8) :: Q_sol ! solved rotation frequency + REAL(r8) :: br_th ! radial field threshold + REAL(r8) :: Qratio ! Q_e/Q ratio for scan2 variants +c----------------------------------------------------------------------- +c local scalars — scan grid helpers. +c----------------------------------------------------------------------- + REAL(r8) :: inQ_min,inQ_max ! rotation scan bounds + REAL(r8) :: j_min,j_max,jpower ! 2-D scan axis 1 + REAL(r8) :: k_min,k_max,kpower ! 2-D scan axis 2 + REAL(r8) :: ing_step ! growth-rate grid step + REAL(r8) :: ing_coarse ! Re(gamma) grid value + REAL(r8) :: iing_coarse ! Im(gamma) grid value +c----------------------------------------------------------------------- +c local scalars — complex quantities. +c----------------------------------------------------------------------- + COMPLEX(r8) :: delta ! layer Delta (tearing index) + COMPLEX(r8) :: delta_n_p ! Deltaprime scale factor + COMPLEX(r8) :: dels_db ! delta_s / d_beta + COMPLEX(r8) :: del_s ! delta_s + COMPLEX(r8) :: ingamma ! initial gamma guess (namelist) + COMPLEX(r8) :: delta_prime ! external Deltaprime (namelist) +c----------------------------------------------------------------------- +c local arrays — transport profile coefficients. +c----------------------------------------------------------------------- + REAL(r8) :: chis(3) ! chi_perp, chi_tor, kappa + REAL(r8), DIMENSION(8) :: chi_p_prof ! chi_perp radial profile + REAL(r8), DIMENSION(8) :: chi_t_prof ! chi_tor radial profile + REAL(r8), DIMENSION(8) :: kappa_prof ! kappa radial profile +c----------------------------------------------------------------------- +c local arrays — multi-surface input-file storage. +c----------------------------------------------------------------------- + INTEGER, DIMENSION(:), ALLOCATABLE :: mms,nns + REAL(r8), DIMENSION(:), ALLOCATABLE :: prs,n_es,t_es,t_is, + $ omegas,l_ns,l_ts,svals,qvals,bts,rss,R0s,mu_is,zeffs, + $ Q_soll,br_thl,pes +c----------------------------------------------------------------------- +c local arrays — scan workspace. +c----------------------------------------------------------------------- + REAL(r8), DIMENSION(:), ALLOCATABLE :: inQs,iinQs + REAL(r8), DIMENSION(:), ALLOCATABLE :: jxbl,bal + REAL(r8), DIMENSION(:,:), ALLOCATABLE :: js,ks,psis,jxbs, + $ Q_sols,br_ths REAL(r8), DIMENSION(:,:,:), ALLOCATABLE :: Q_solss,br_thss - COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: deltal + COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: deltal,outer_deltas COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: deltas +c----------------------------------------------------------------------- +c AMR (adaptive mesh refinement) growth-rate scan. +c----------------------------------------------------------------------- + INTEGER :: AMR_passes ! number of AMR refinement levels + INTEGER :: m_AMR ! effective number of surfaces +c----------------------------------------------------------------------- +c structured input/output types (defined in sglobal_mod). +c----------------------------------------------------------------------- + TYPE(slayer_inputs_type) :: sl_in + TYPE(slayer_outputs_type) :: sl_out + TYPE(deltas_outputs_type), ALLOCATABLE :: all_deltas_out(:) +c----------------------------------------------------------------------- +c file-path strings. +c----------------------------------------------------------------------- + CHARACTER(512) :: infile ! multi-surface input file path + CHARACTER(512) :: ncfile ! NetCDF equilibrium file path - NAMELIST/slayer_input/params_flag,input_flag,infile, - $ mm,nn,n_e,t_e,t_i,omega,l_n,l_t, - $ qval,sval,bt,rs,R0,zeff,mu_i,inQ,inQ_e,inQ_i, - $ inpr,inpr_prof,inpe,inc_beta,inds,intau,inlu,Q0,delta_n_p - NAMELIST/slayer_control/inum,jnum,knum,QPscan_flag,QPscan2_flag, - $ QPescan_flag,QDscan2_flag,Qbscan_flag,Qscan_flag, - $ onscan_flag,otscan_flag,ntscan_flag,nbtscan_flag, - $ layfac,Qratio,parflow_flag,peohmonly_flag,Pe_flag +c----------------------------------------------------------------------- +c namelist groups. +c----------------------------------------------------------------------- + NAMELIST/slayer_input/input_flag,infile, + $ ncfile,params_flag,mm,nn,n_e,t_e,t_i,sval,bt,rs,R0,omega, + $ l_t,l_n,qval,mu_i,zeff,dr_val,dgeo_val,chi_p_prof, + $ chi_t_prof,kappa_prof,inpr,inpe,inQ,inQ_e,inQ_i,inc_beta, + $ inds,intau,inlu,inpr_prof,Q0,delta_prime,delta_n_p,ingamma, + $ inQ_e_ovr,inQ_i_ovr + NAMELIST/slayer_control/inum,jnum,knum,Q_num,scan_width, + $ scan_center,rotation_flag, + $ AMR_passes,msing_max,dc_type,read_eq,Pperp_Ptor_flag, + $ coupling_flag,QPscan_flag,Qscan_flag,QPescan_flag, + $ Qbscan_flag,onscan_flag,otscan_flag,ntscan_flag, + $ nbtscan_flag,parflow_flag,peohmonly_flag,Pe_flag,layfac, + $ QPscan2_flag,QDscan2_flag,Qratio NAMELIST/slayer_output/verbose,ascii_flag,bin_flag,netcdf_flag, - $ stability_flag + $ est_gamma_flag,match_gamma_flag,stability_flag, + $ stabscan_flag,coupled_stabscan_flag,amr_flag,br_th_flag, + $ bal_flag + ! bal_flag is documented in slayer_diagnose, as on develop; it + ! is also accepted in slayer_output so either input layout reads NAMELIST/slayer_diagnose/riccati_out,riccatiscan_flag, - $ params_check,bal_flag + $ params_check,bal_flag c----------------------------------------------------------------------- c set initial values. +c defaults are overridden by namelist reads below. c----------------------------------------------------------------------- - mm=2 - nn=1 - mr = real(mm,4) - nr = real(nn,4) - n_e=1e19 - t_e=1e3 - t_i=1e3 - omega=1e4 - l_n=2e-1 - l_t=2e-1 - qval=2.0 - sval=2.0 - bt=1.0 - rs=0.5 - R0=1.0 - mu_i=2.0 - zeff=2.0 - inQ=20.0 ! Q=23.0 for DIII-D example. - inQ_e=2.0 ! Q_e=2.0 for DIII-D example. - inQ_i=-2.6 ! Q_i=-2.6 for DIII-D example. - inpr=1.0 ! - inpr_prof=1.0 - inpe=0.0 - inc_beta=0.7 ! c_beta=0.7 for DIII-D example. - inds=6.0 ! 6.0 for DIII-D example. - intau=1.0 ! 1.0 for DIII-D example. - inlu=1e8 - Q0=4.0 - delta_n_p=(1e-2,1e-2) - inum=400 ! resolution to find error field thresholds. - jnum=500 ! resolution for 2d scan along with Q,omega. - knum=100 ! resolution for 2d scan alont with the other. - in_unit=1 - out_unit=2 - out2_unit=3 - out3_unit=4 - bin_unit=5 - bin_2d_unit=6 - input_unit=7 - QPscan_flag=.FALSE. ! scan (Q,P) space for delta and torque. - QPescan_flag=.FALSE. ! scan (Q,Pe) space for delta and torque. - Qbscan_flag=.FALSE. ! scan (Q,beta) space for delta and torque. - onscan_flag=.FALSE. ! scan (omega,n) space for error fields. - otscan_flag=.FALSE. ! scan (omega,t) space for error fields. - ntscan_flag=.FALSE. ! scan (n,te) space for error fields. - nbtscan_flag=.FALSE. ! scan (n,bt) space for error fields. - layfac=0.02 - Qratio=0.5 - parflow_flag=.FALSE. - PeOhmOnly_flag=.TRUE. - Pe_flag=.FALSE. - params_flag=.TRUE. - input_flag=.FALSE. - infile="input_params.dat" - verbose=.TRUE. - ascii_flag=.TRUE. - bin_flag=.TRUE. - netcdf_flag=.FALSE. - riccati_out=.FALSE. - riccatiscan_flag=.FALSE. - params_check=.FALSE. - bal_flag=.FALSE. - stability_flag=.FALSE. + + ! mode numbers (sglobal_mod: INTEGER mm,nn; REAL mr,nr) + mm = 0 + nn = 0 + mr = 0.0 + nr = 0.0 + + ! kinetic / equilibrium inputs + n_e = 0.0 + t_e = 0.0 + t_i = 0.0 + omega = 0.0 + l_n = 0.0 + l_t = 0.0 + qval = 0.0 + sval = 0.0 + bt = 0.0 + rs = 0.0 + R0 = 0.0 + mu_i = 0.0 + zeff = 0.0 + dr_val = 0.0 + dgeo_val= 0.0 + + ! normalized layer-parameter overrides + inQ = 0.0 + inQ_e = 0.0 + inQ_i = 0.0 + inQ_e_ovr = 0.0 + inQ_i_ovr = 0.0 + inpr = 0.0 + inpr_prof = 1.0_r8 + inpe = 0.0 + inc_beta = 0.0 + inds = 0.0 + intau = 0.0 + inlu = 1.0e8_r8 + Q0 = 0.0 + + ! transport profile coefficients + chi_p_prof = 0.0 + chi_t_prof = 0.0 + kappa_prof = 0.0 + chis = 0.0 + + ! complex namelist inputs + delta_prime = (0.0,0.0) + delta_n_p = (0.0,0.0) + ingamma = (0.0,0.0) + + ! global module scalars (sglobal_mod) + gamma_fac = 0.0 + dc_type = "rfitzp" + + ! scan resolution defaults + inum = 400 ! 1-D resolution (error-field threshold scans) + jnum = 500 ! 2-D scan axis-1 resolution + knum = 100 ! 2-D scan axis-2 resolution + Q_num = 100 ! stability scan Re(Q) resolution + scan_width = 2.0 ! half-width for complex-Q scans + ! Re(Q) center of the scans. Left at HUGE, it is chosen by the + ! code: zero, or the middle of the Doppler-shifted branches in + ! a coupled scan with rotation. + scan_center = HUGE(scan_center) + AMR_passes = 4 ! AMR refinement levels + m_AMR = 0 ! no AMR scan stored yet + msing_max = 2 ! max singular surfaces to process + + ! I/O unit numbers (sglobal_mod) + in_unit = 1 + out_unit = 2 + out2_unit = 3 + out3_unit = 4 + bin_unit = 5 + bin_2d_unit= 6 + input_unit = 7 + + ! workflow flags + read_eq = .FALSE. + est_gamma_flag = .FALSE. + match_gamma_flag = .FALSE. + Pperp_Ptor_flag = .FALSE. + coupling_flag = .FALSE. + rotation_flag = .TRUE. + params_flag = .TRUE. + input_flag = .FALSE. + + ! parameter-space scan flags + QPscan_flag = .FALSE. + QPescan_flag = .FALSE. + QPscan2_flag = .FALSE. + QDscan2_flag = .FALSE. + Qbscan_flag = .FALSE. + onscan_flag = .FALSE. + otscan_flag = .FALSE. + ntscan_flag = .FALSE. + nbtscan_flag = .FALSE. + + ! physics / model flags + layfac = 0.02 + Qratio = 0.5 + parflow_flag = .FALSE. + PeOhmOnly_flag = .TRUE. + Pe_flag = .FALSE. + + ! file paths + infile = "" + ncfile = "" + + ! output control + verbose = .TRUE. + ascii_flag = .TRUE. + bin_flag = .TRUE. + netcdf_flag = .FALSE. + + ! diagnostic flags + riccati_out = .FALSE. + riccatiscan_flag = .FALSE. + params_check = .FALSE. + + ! remaining physics-mode flags + bal_flag = .FALSE. + stability_flag = .FALSE. + stabscan_flag = .FALSE. + coupled_stabscan_flag= .FALSE. + amr_flag = .FALSE. + br_th_flag = .FALSE. c----------------------------------------------------------------------- c read slayer.in. +c four namelist groups: input, control, output, diagnose. c----------------------------------------------------------------------- IF(verbose) WRITE(*,*)"" IF(verbose) WRITE(*,*)"SLAYER START" @@ -143,18 +357,128 @@ PROGRAM slayer READ(in_unit,NML=slayer_diagnose) CLOSE(UNIT=in_unit) + ! Critical-Delta model: fail on a misspelt name rather than + ! silently running without the stabilising offset. + SELECT CASE (TRIM(dc_type)) + CASE ('rfitzp','none','') + CASE ('toroidal') + WRITE(*,*) "ERROR: dc_type='toroidal' is not available." + WRITE(*,*) " Its geometric factor was not correct in the" + WRITE(*,*) " r_s reference length and has been removed." + WRITE(*,*) " Use dc_type='rfitzp'." + STOP 'slayer: dc_type toroidal removed' + CASE ('lar') + WRITE(*,*) "ERROR: dc_type='lar' is not available." + WRITE(*,*) " It was the rfitzp formula without its r_s" + WRITE(*,*) " normalization (rfitzp/r_s) and has been" + WRITE(*,*) " removed. Use dc_type='rfitzp' or 'none'." + STOP 'slayer: dc_type lar removed' + CASE DEFAULT + WRITE(*,*) "ERROR: unknown dc_type='",TRIM(dc_type),"'." + WRITE(*,*) " Use 'rfitzp' or 'none'." + STOP 'slayer: unknown dc_type' + END SELECT + +c----------------------------------------------------------------------- +c check that the inputs this run needs were set. Since v1.6 the +c namelist has no silent defaults: zero means "not set", or for +c inQ, inQ_e and inQ_i "compute from the profiles". A v1.5 input +c that omitted a variable is reported with the old default. +c----------------------------------------------------------------------- + nbad = 0 + ! dimensional inputs: needed when params() runs on the namelist + ! (params_flag does so even with read_eq) + IF (params_flag .OR. (.NOT. read_eq .AND. (est_gamma_flag .OR. + $ match_gamma_flag .OR. br_th_flag .OR. onscan_flag .OR. + $ otscan_flag .OR. ntscan_flag .OR. nbtscan_flag))) THEN + IF (mm <= 0) CALL slayer_unset('mm',REAL(mm,r8),'2',nbad) + IF (nn <= 0) CALL slayer_unset('nn',REAL(nn,r8),'1',nbad) + IF (.NOT. n_e > 0.0_r8) + $ CALL slayer_unset('n_e',n_e,'1e19 [m^-3]',nbad) + IF (.NOT. t_e > 0.0_r8) + $ CALL slayer_unset('t_e',t_e,'1e3 [eV]',nbad) + IF (.NOT. t_i > 0.0_r8) + $ CALL slayer_unset('t_i',t_i,'1e3 [eV]',nbad) + IF (sval == 0.0_r8) CALL slayer_unset('sval',sval,'2.0',nbad) + IF (.NOT. bt > 0.0_r8) CALL slayer_unset('bt',bt,'1.0',nbad) + IF (.NOT. rs > 0.0_r8) CALL slayer_unset('rs',rs,'0.5',nbad) + IF (.NOT. R0 > 0.0_r8) CALL slayer_unset('R0',R0,'1.0',nbad) + IF (.NOT. mu_i > 0.0_r8) + $ CALL slayer_unset('mu_i',mu_i,'2.0',nbad) + IF (.NOT. zeff > 0.0_r8) + $ CALL slayer_unset('zeff',zeff,'2.0',nbad) + ! the gradient lengths and q give omega_e, omega_i unless a + ! growth-rate run overrides both Q_e and Q_i + IF (.NOT. ((est_gamma_flag .OR. match_gamma_flag) .AND. + $ inQ_e /= 0.0_r8 .AND. inQ_i /= 0.0_r8)) THEN + IF (.NOT. l_n > 0.0_r8) + $ CALL slayer_unset('l_n',l_n,'0.2 [m]',nbad) + IF (.NOT. l_t > 0.0_r8) + $ CALL slayer_unset('l_t',l_t,'0.2 [m]',nbad) + IF (.NOT. qval > 0.0_r8) + $ CALL slayer_unset('qval',qval,'2.0',nbad) + END IF + END IF + ! Cole-model inputs: needed by the single-surface riccati() call + ! below, which runs unless a per-surface mode is on, and by the + ! legacy scans, which run unless a mode that stops first is on. + ! inpr is a namelist input in all of them. + IF (.NOT. (est_gamma_flag .OR. match_gamma_flag .OR. + $ br_th_flag .OR. input_flag) .OR. + $ (.NOT. (match_gamma_flag .OR. br_th_flag) .AND. + $ (bal_flag .OR. stability_flag .OR. riccatiscan_flag .OR. + $ QPscan_flag .OR. QPescan_flag .OR. QPscan2_flag .OR. + $ QDscan2_flag .OR. Qscan_flag .OR. onscan_flag .OR. + $ otscan_flag .OR. ntscan_flag .OR. nbtscan_flag))) THEN + IF (.NOT. inpr > 0.0_r8) + $ CALL slayer_unset('inpr',inpr,'1.0',nbad) + IF (.NOT. params_flag) THEN + IF (.NOT. inds > 0.0_r8) + $ CALL slayer_unset('inds',inds,'6.0',nbad) + IF (.NOT. intau > 0.0_r8) + $ CALL slayer_unset('intau',intau,'1.0',nbad) + END IF + END IF + IF (nbad > 0) THEN + WRITE(*,*) ' Since v1.6 SLAYER sets none of these silently;' + WRITE(*,*) ' set them in slayer.in. (With read_eq=t and' + WRITE(*,*) ' params_flag=f the dimensional inputs come from' + WRITE(*,*) ' STRIDE and the kinetic profile file instead.)' + WRITE(*,*) ' Other v1.5 defaults, for reference: inQ=20,' + WRITE(*,*) ' inQ_e=2, inQ_i=-2.6, inc_beta=0.7, Q0=4,' + WRITE(*,*) ' delta_n_p=(1e-2,1e-2). Zero is a valid choice' + WRITE(*,*) ' for those (no rotation or drift).' + STOP 'slayer: required inputs unset' + END IF + + ! Real copies of the mode numbers, read by params() in every + ! namelist-driven branch below. build_inputs resets them from + ! the STRIDE run when read_eq is set. + mr = mm + nr = nn + + ! Build toroidal-mode-number string for output filenames. IF (nn<10) THEN - WRITE(UNIT=sn,FMT='(I1)') nn - sn=ADJUSTL(sn) + WRITE(UNIT=sn_str,FMT='(I1)') nn + sn_str=ADJUSTL(sn_str) ELSE - WRITE(UNIT=sn,FMT='(I2)') nn + WRITE(UNIT=sn_str,FMT='(I2)') nn ENDIF c----------------------------------------------------------------------- -c calculate parameters as needed. +c compute normalized layer parameters from kinetic inputs. +c params() (params_mod) converts dimensional plasma profiles into +c the normalized quantities (Q, Q_e, Q_i, c_beta, ds, tau, lu) +c used by the Riccati solver. c----------------------------------------------------------------------- IF (params_flag) THEN - CALL params(n_e,t_e,t_i,omega, + nr = nn ! ensure toroidal mode number is set for Wd iteration + mr = mm ! ensure poloidal mode number is set + chis(1) = chi_p_prof(1) ! chi_perp + chis(2) = chi_t_prof(1) ! chi_tor + chis(3) = kappa_prof(1) ! kappa + CALL params(n_e,t_e,t_i,omega,chis,dr_val,dgeo_val, $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) + ! Copy module-level results into local working variables. inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -163,28 +487,43 @@ PROGRAM slayer intau=tau Q0=Q ELSE - lu=inlu + lu=inlu ! manual Lundquist number when params not computed ENDIF c----------------------------------------------------------------------- -c calculate basic delta, torque, balance, error fields. +c baseline single-surface delta, reconnected flux, & torque. +c skipped when the matched-gamma path is active (it computes +c its own delta internally). c----------------------------------------------------------------------- - delta=riccati(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe) - psi0=1.0/ABS(delta+delta_n_p) ! a.u. - jxb=-AIMAG(1.0/(delta+delta_n_p)) ! a.u. - WRITE(*,*)"delta=",delta - WRITE(*,*)"psi0=",psi0 - WRITE(*,*)"jxb=",jxb +c Also skipped for the growth-rate estimate, the Br threshold and +c the multi-surface input file, which compute their own deltas per +c surface and are usually run with these namelist inputs unset. + IF (.NOT. (match_gamma_flag .OR. est_gamma_flag .OR. + $ br_th_flag .OR. input_flag)) THEN + delta=riccati(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe) + psi0=1.0/ABS(delta+delta_n_p) ! reconnected flux [a.u.] + jxb=-AIMAG(1.0/(delta+delta_n_p)) ! j x B torque [a.u.] + IF (verbose) THEN + WRITE(*,*)"delta=",delta + WRITE(*,*)"psi0=",psi0 + WRITE(*,*)"jxb=",jxb + ENDIF + ENDIF c----------------------------------------------------------------------- -c calculate parameters as needed. +c multi-surface input-file mode. +c reads an external file of (m,n) surfaces with per-surface kinetic +c profiles, computes delta and the error-field threshold for each. c----------------------------------------------------------------------- IF (input_flag) THEN OPEN(UNIT=input_unit,FILE=infile,STATUS="old") READ(input_unit,*)inn - ALLOCATE(mms(inn),nns(inn),prs(inn), - $ n_es(inn),t_es(inn),t_is(inn),omegas(inn), - $ l_ns(inn),l_ts(inn),qvals(inn),svals(inn), - $ bts(inn),rss(inn),R0s(inn),mu_is(inn),zeffs(inn), - $ Q_soll(inn),br_thl(inn)) + ALLOCATE(mms(0:inn-1),nns(0:inn-1),prs(0:inn-1), + $ n_es(0:inn-1),t_es(0:inn-1),t_is(0:inn-1), + $ omegas(0:inn-1), + $ l_ns(0:inn-1),l_ts(0:inn-1),qvals(0:inn-1), + $ svals(0:inn-1), + $ bts(0:inn-1),rss(0:inn-1),R0s(0:inn-1), + $ mu_is(0:inn-1),zeffs(0:inn-1), + $ Q_soll(0:inn-1),br_thl(0:inn-1)) DO k=0,inn-1 READ(input_unit,'(2(1x,I2),14(1x,e12.4))') $ mms(k),nns(k),prs(k), @@ -195,13 +534,13 @@ PROGRAM slayer CLOSE(input_unit) DO k=0,inn-1 - WRITE(*,*)k + WRITE(*,*)k ! surface index mr=REAL(mms(k)) nr=REAL(nns(k)) inpr=prs(k) - CALL params(n_es(k),t_es(k),t_is(k),omegas(k), - $ l_ns(k),l_ts(k),qvals(k),svals(k),bts(k),rss(k),R0s(k), - $ mu_is(k),zeffs(k),params_check) + CALL params(n_es(k),t_es(k),t_is(k),omegas(k),chis,dr_val, + $ dgeo_val,l_ns(k),l_ts(k),qvals(k),svals(k),bts(k), + $ rss(k),R0s(k),mu_is(k),zeffs(k),params_check) inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -234,16 +573,17 @@ PROGRAM slayer bal(i)=2.0*inpr*(Q0-inQs(i))/jxbl(i) ENDDO - index=MAXLOC(bal) - Q_soll(k)=inQs(index(1)) + iloc=MAXLOC(bal) + iloc(1)=iloc(1)-1 ! MAXLOC counts from 1, inQs from 0 + Q_soll(k)=inQs(iloc(1)) br_thl(k)=sqrt(MAXVAL(bal)/lu*(svals(k)**2.0/2.0))*1e4 - WRITE(*,*)"Q_sol=",Q_soll(k) - WRITE(*,*)"br_th=",br_thl(k) + IF (verbose) WRITE(*,*)"Q_sol=",Q_soll(k) + IF (verbose) WRITE(*,*)"br_th=",br_thl(k) DEALLOCATE(inQs,deltal,jxbl,bal) ENDDO OPEN(UNIT=out_unit,FILE="slayer_input_bal_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,(2a17))') "Q_sol","br_th" DO k=0,inn-1 @@ -255,6 +595,785 @@ PROGRAM slayer $ bts,rss,R0s,mu_is,zeffs,Q_soll,br_thl,mms,nns) ENDIF c----------------------------------------------------------------------- +c estimate growth rate via resistive-layer thickness. +c Uses riccati_del_s to get delta_s/d_beta, then scales by +c d_beta to obtain the layer thickness delta_s and the estimated +c gamma. Inputs may come from equilibrium files (read_eq) or +c namelist. Subroutines: build_inputs (layerinputs_mod), +c allocate_inputs, allocate_outputs (growthrates_mod). +c----------------------------------------------------------------------- + IF (est_gamma_flag) THEN + WRITE(*,*)"------------------------------------------" + WRITE(*,*)">>> Estimating growth rate" + + IF (read_eq) THEN + ! Read equilibrium files for multi-surface inputs. + ! build_inputs (layerinputs_mod) reads STRIDE NetCDF data. + sl_in%chi_p_arr = chi_p_prof + sl_in%chi_t_arr = chi_t_prof + sl_in%kappa_arr = kappa_prof + + CALL build_inputs(infile,ncfile,sl_in,zeff,mu_i) + + n_k = SIZE(sl_in%qval_arr) + CALL allocate_outputs(n_k,sl_out) + +c Override Q_e/Q_i with namelist values if nonzero. +c inQ_e overrides all surfaces; per-surface arrays take +c precedence (surface k gets inQ_e_ovr(k) if nonzero). + IF (ABS(inQ_e) > 0.0) THEN + WRITE(*,*) 'Overriding all Q_e with inQ_e=',inQ_e + sl_in%Q_e_arr(:) = inQ_e + END IF + IF (ABS(inQ_i) > 0.0) THEN + WRITE(*,*) 'Overriding all Q_i with inQ_i=',inQ_i + sl_in%Q_i_arr(:) = inQ_i + END IF + DO k = 1, MIN(n_k, SIZE(inQ_e_ovr)) + IF (ABS(inQ_e_ovr(k)) > 0.0) THEN + WRITE(*,*) ' Q_e override surface',k,':', + $ sl_in%Q_e_arr(k),' ->',inQ_e_ovr(k) + sl_in%Q_e_arr(k) = inQ_e_ovr(k) + END IF + END DO + DO k = 1, MIN(n_k, SIZE(inQ_i_ovr)) + IF (ABS(inQ_i_ovr(k)) > 0.0) THEN + sl_in%Q_i_arr(k) = inQ_i_ovr(k) + END IF + END DO +c D_norm carries sqrt(iota_e): keep it consistent with the +c overridden Q_e, Q_i. + DO k = 1, n_k + sl_in%D_norm_arr(k) = slayer_d_norm( + $ sl_in%D_fac_arr(k),sl_in%Q_e_arr(k), + $ sl_in%Q_i_arr(k)) + END DO + + ELSE + ! Single-surface mode: build inputs from namelist. + n_k = 1 + mr = mm + nr = nn + + chis(1) = chi_p_prof(1) ! chi_perp + chis(2) = chi_t_prof(1) ! chi_tor + chis(3) = kappa_prof(1) ! kappa (thermal cond.) + + CALL params(n_e,t_e,t_i,omega,chis,dr_val,dgeo_val, + $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) + + ! Override computed parameters with nonzero namelist values. + IF (ABS(inQ) > 0.0) THEN + Q = inQ ! NAMELIST + END IF + IF (ABS(inQ_e) > 0.0) THEN + Q_e = inQ_e ! NAMELIST + END IF + IF (ABS(inQ_i) > 0.0) THEN + Q_i = inQ_i ! NAMELIST + END IF + IF (inpr > 0.0) THEN + pr = inpr ! NAMELIST + END IF + IF (intau > 0.0) THEN + tau = intau ! NAMELIST + END IF + IF (inds > 0.0) THEN + D_norm = inds ! NAMELIST + ELSE + D_norm = slayer_d_norm(D_fac,Q_e,Q_i) + END IF + ! Deltaprime -> gamma conversion, formed as in build_inputs + ! from the r_s-referenced Deltaprime and the resistive time + gamma_fac = rs*REAL(delta_prime)/tau_r + + CALL allocate_inputs(n_k,sl_in) ! growthrates_mod + CALL allocate_outputs(n_k,sl_out) ! growthrates_mod + + sl_in%qval_arr = (/ mm /) + sl_in%q_arr = (/ qval /) + sl_in%omegas_arr = (/ omega /) + sl_in%Qexb_arr = (/ REAL(Q,r8) /) + sl_in%shear_arr = (/ sval /) + sl_in%ds_arr = (/ ds /) + sl_in%Q_e_arr = (/ Q_e /) + sl_in%Q_i_arr = (/ Q_i /) + sl_in%psi_n_arr = (/ 0.0 /) + sl_in%Re_dp_arr = (/ REAL(delta_prime) /) + sl_in%Im_dp_arr = (/ AIMAG(delta_prime) /) + sl_in%Re_dp_psi_arr = (/ REAL(delta_prime) /) + sl_in%Im_dp_psi_arr = (/ AIMAG(delta_prime) /) + sl_in%d_crit_arr = (/ dc_tmp /) + sl_in%P_perp_arr = (/ P_perp /) + sl_in%P_tor_arr = (/ P_tor /) + sl_in%tau_arr = (/ tau /) + sl_in%D_norm_arr = (/ D_norm /) + sl_in%D_fac_arr = (/ D_fac /) + sl_in%d_beta_arr = (/ d_beta /) + sl_in%gammafac_arr = (/ gamma_fac /) + sl_in%c_beta_arr = (/ c_beta /) + sl_in%lu_arr = (/ lu /) + sl_in%Qconv_arr = (/ tauk /) + END IF + + ! Loop over rational surfaces to estimate growth rates. + DO k=1,n_k + WRITE(*,*) + WRITE(*,'(A,I0,A)') 'Calculating growth rate '// + $ 'estimate on m = ', + $ sl_in%qval_arr(k),' rational surface' + + D_norm = sl_in%D_norm_arr(k) + tau = sl_in%tau_arr(k) + P_tor = sl_in%P_tor_arr(k) + ! First arg is Q_e (electron diamagnetic freq), not Q + ! (ExB freq). This is intentional per riccati_del_s API. + dels_db=riccati_del_s(sl_in%Q_e_arr(k), + $ sl_in%Q_i_arr(k),sl_in%P_perp_arr(k), + $ 5.0*sl_in%D_norm_arr(k)) + + del_s = dels_db * sl_in%d_beta_arr(k) + + sl_out%gamma_est_arr(k) = sl_in%gammafac_arr(k)/del_s + sl_out%dels_db_arr(k) = dels_db + WRITE(*,*) + WRITE(*,'(A,F0.3,A)')'Growth rate estimate = ', + $ REAL(sl_out%gamma_est_arr(k)),' [Hz]' + + ENDDO + + IF (.NOT. (match_gamma_flag)) THEN + ! no matched solve in this run + sl_out%gamma_sol_arr = CMPLX( + $ ieee_value(1.0_r8, ieee_quiet_nan), + $ ieee_value(1.0_r8, ieee_quiet_nan), KIND=r8) + CALL output_gamma(est_gamma_flag,m_AMR,sl_in,sl_out, + $ all_deltas_out) + END IF + ENDIF +c----------------------------------------------------------------------- +c asymptotically matched growth rate. +c Matches the inner-layer Delta to the outer-region Delta' to +c find the self-consistent complex growth rate. Supports both +c single-surface and coupled multi-surface (AMR) modes. +c Subroutines: dispersion_AMR_v2, dispersion_det (growthrates_mod), +c riccati_f (delta_mod). +c----------------------------------------------------------------------- + IF (match_gamma_flag) THEN + WRITE(*,*)"------------------------------------------" + WRITE(*,*)">>> Calculating asymptotically matched growth rate" + + IF (read_eq) THEN + + IF (.NOT. est_gamma_flag) THEN + sl_in%chi_p_arr = chi_p_prof + sl_in%chi_t_arr = chi_t_prof + sl_in%kappa_arr = kappa_prof + + CALL build_inputs(infile,ncfile,sl_in,zeff,mu_i) + + n_k = SIZE(sl_in%qval_arr) + CALL allocate_outputs(n_k,sl_out) + END IF +c Apply Q_e/Q_i overrides (same logic as est_gamma path) + DO k = 1, MIN(n_k, SIZE(inQ_e_ovr)) + IF (ABS(inQ_e_ovr(k)) > 0.0) THEN + sl_in%Q_e_arr(k) = inQ_e_ovr(k) + END IF + END DO + DO k = 1, MIN(n_k, SIZE(inQ_i_ovr)) + IF (ABS(inQ_i_ovr(k)) > 0.0) THEN + sl_in%Q_i_arr(k) = inQ_i_ovr(k) + END IF + END DO +c D_norm carries sqrt(iota_e): keep it consistent with the +c overridden Q_e, Q_i. + DO k = 1, n_k + sl_in%D_norm_arr(k) = slayer_d_norm( + $ sl_in%D_fac_arr(k),sl_in%Q_e_arr(k), + $ sl_in%Q_i_arr(k)) + END DO + ELSE + n_k = 1 + + IF (.NOT. est_gamma_flag) THEN + + chis(1) = chi_p_prof(1) ! chi_perp + chis(2) = chi_t_prof(1) ! chi_tor + chis(3) = kappa_prof(1) ! kappa (thermal cond.) + + WRITE(*,*)"chis(1) (chi_perp): ",chis(1) + WRITE(*,*)"chis(2) (chi_tor): ",chis(2) + WRITE(*,*)"chis(3) (kappa): ",chis(3) + + ! Use namelist kinetic inputs instead of equilibrium files + CALL params(n_e,t_e,t_i,omega,chis,dr_val,dgeo_val, + $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) + + ! Override desired normalized parameters + IF (ABS(inQ) > 0.0) THEN + Q = inQ ! NAMELIST + END IF + IF (ABS(inQ_e) > 0.0) THEN + Q_e = inQ_e ! NAMELIST + END IF + IF (ABS(inQ_i) > 0.0) THEN + Q_i = inQ_i ! NAMELIST + END IF + IF (inpr > 0.0) THEN + P_perp = inpr ! NAMELIST + END IF + IF (intau > 0.0) THEN + tau = intau ! NAMELIST + END IF + IF (inds > 0.0) THEN + D_norm = inds ! NAMELIST + ELSE + D_norm = slayer_d_norm(D_fac,Q_e,Q_i) + END IF + ! Deltaprime -> gamma conversion, formed as in build_inputs + ! from the r_s-referenced Deltaprime and the resistive time + gamma_fac = rs*REAL(delta_prime)/tau_r + + IF (.NOT. est_gamma_flag) THEN + CALL allocate_inputs(n_k,sl_in) + CALL allocate_outputs(n_k,sl_out) + END IF + + sl_in%qval_arr = (/ mm /) + sl_in%q_arr = (/ qval /) + sl_in%omegas_arr = (/ omega /) + sl_in%Qexb_arr = (/ REAL(Q,r8) /) + sl_in%shear_arr = (/ sval /) + sl_in%ds_arr = (/ ds /) + sl_in%Q_e_arr = (/ Q_e /) + sl_in%Q_i_arr = (/ Q_i /) + sl_in%psi_n_arr = (/ 0.0 /) + sl_in%Re_dp_arr = (/ REAL(delta_prime) /) + sl_in%Im_dp_arr = (/ AIMAG(delta_prime) /) + sl_in%Re_dp_psi_arr = (/ REAL(delta_prime) /) + sl_in%Im_dp_psi_arr = (/ AIMAG(delta_prime) /) + sl_in%d_crit_arr = (/ dc_tmp /) + sl_in%P_perp_arr = (/ P_perp /) + sl_in%P_tor_arr = (/ P_tor /) + sl_in%tau_arr = (/ tau /) + sl_in%D_norm_arr = (/ D_norm /) + sl_in%D_fac_arr = (/ D_fac /) + sl_in%d_beta_arr = (/ d_beta /) + sl_in%gammafac_arr = (/ gamma_fac /) + sl_in%c_beta_arr = (/ c_beta /) + sl_in%lu_arr = (/ lu /) + sl_in%Qconv_arr = (/ tauk /) + + END IF + END IF + +c----------------------------------------------------------------------- +c sanity checks on the surface count and the coupling request. +c----------------------------------------------------------------------- + IF (msing_max < 1) THEN + WRITE(*,*) 'WARNING: msing_max < 1; using msing_max = 1' + msing_max = 1 + END IF + IF (msing_max > n_k) THEN + WRITE(*,'(A,I0,A,I0,A)') ' NOTE: msing_max = ',msing_max, + $ ' exceeds the ',n_k,' rational surfaces available;' + $ //' using all of them.' + msing_max = n_k + END IF + IF (coupling_flag .AND. .NOT. read_eq) THEN + WRITE(*,*) 'ERROR: coupling_flag=t needs the STRIDE' + WRITE(*,*) ' Deltaprime matrix; set read_eq=t, infile' + WRITE(*,*) ' and ncfile, or set coupling_flag=f.' + STOP 'slayer: coupling_flag requires read_eq' + END IF + +c----------------------------------------------------------------------- +c ExB Doppler shift and scan-window centre. +c Each uncoupled layer is solved in its own ExB frame, where a +c rotation only relabels the frequency, so no shift is applied. +c In the coupled determinant the surfaces share one lab-frame +c eigenvalue and surface k sees it Doppler-shifted by +c n*Omega_E,k: q_shift = -tau_k*n*Omega_E,k. In the scanned Q +c (normalised by tau_1) branch k then sits near +c tau_1*(n*Omega_E,k - Q_e,k/tau_k), taking the electron-drift +c frequency as the estimate of its own-frame root. +c----------------------------------------------------------------------- + sl_in%q_shift_arr = 0.0_r8 + IF (coupling_flag .AND. rotation_flag .AND. read_eq) THEN + q_lo = HUGE(q_lo) + q_hi = -HUGE(q_hi) + DO k = 1, MIN(n_k,msing_max) + sl_in%q_shift_arr(k) = -sl_in%Qconv_arr(k) + $ * REAL(nn,r8) * sl_in%omegas_arr(k) + q_br = sl_in%Qconv_arr(1)*( REAL(nn,r8) + $ * sl_in%omegas_arr(k) + $ - sl_in%Q_e_arr(k)/sl_in%Qconv_arr(k) ) + q_lo = MIN(q_lo,q_br) + q_hi = MAX(q_hi,q_br) + WRITE(*,'(A,I0,A,ES11.4,A,ES11.4)') + $ ' ExB rotation, surface ',k,': n*Omega_E [rad/s] =', + $ REAL(nn,r8)*sl_in%omegas_arr(k), + $ ', branch expected near Re(Q) =',q_br + END DO + IF (scan_center >= HUGE(scan_center)) + $ scan_center = 0.5_r8*(q_lo + q_hi) + IF (q_lo < scan_center - scan_width .OR. + $ q_hi > scan_center + scan_width) THEN + WRITE(*,'(A)') ' WARNING: the Doppler-shifted branches' + $ //' do not all fit in the scan window' + WRITE(*,'(A,2ES11.3,A,2ES11.3,A)') + $ ' Re(Q) = [',scan_center-scan_width, + $ scan_center+scan_width,'], branches span [', + $ q_lo,q_hi,'].' + WRITE(*,'(A)') ' Increase scan_width or set ' + $ //'scan_center, or a coupled root may be missed.' + END IF + END IF + IF (scan_center >= HUGE(scan_center)) scan_center = 0.0_r8 + WRITE(*,'(A,ES11.4,A,ES11.4)') + $ ' Scan window: Re(Q) center =',scan_center, + $ ', half-width =',scan_width + +c----------------------------------------------------------------------- +c allocate output arrays for AMR delta storage. +c----------------------------------------------------------------------- + IF (AMR_flag .AND. .NOT. coupling_flag) THEN + ALLOCATE(all_deltas_out(n_k)) + ALLOCATE(roots_out(n_k)) + ELSEIF (AMR_flag .AND. coupling_flag) THEN + ALLOCATE(all_deltas_out(1)) + ALLOCATE(roots_out(1)) + END IF + +c growth rates stay NaN on surfaces where no root is extracted + sl_out%gamma_sol_arr = CMPLX( + $ ieee_value(1.0_r8, ieee_quiet_nan), + $ ieee_value(1.0_r8, ieee_quiet_nan), KIND=r8) + + IF (AMR_flag) THEN + IF (coupling_flag) THEN + m_AMR = 1 + ELSE + m_AMR = MIN(n_k,msing_max) + END IF + END IF + + WRITE(*,*) "Resonant m of the surfaces: ",sl_in%qval_arr +c----------------------------------------------------------------------- +c loop over rational surfaces to find matched growth rates. +c----------------------------------------------------------------------- + DO k=1,MIN(n_k,msing_max) + WRITE(*,*) + WRITE(*,'(A,I0,A)') 'Calculating growth rate on m = ', + $ sl_in%qval_arr(k),' rational surface:' + + ! Load per-surface parameters into module-level scalars. + Q_e = sl_in%Q_e_arr(k) + Q_i = sl_in%Q_i_arr(k) + P_perp = sl_in%P_perp_arr(k) + P_tor = sl_in%P_tor_arr(k) + tau = sl_in%tau_arr(k) + D_norm = sl_in%D_norm_arr(k) + c_beta = sl_in%c_beta_arr(k) + tauk = sl_in%Qconv_arr(k) + iota_e = Q_e / (Q_e - Q_i) + + WRITE(*,*)"Q_e: ",Q_e + WRITE(*,*)"Q_i: ",Q_i + WRITE(*,*)"P_perp: ",P_perp + WRITE(*,*)"P_tor: ",P_tor + WRITE(*,*)"tau: ",tau + WRITE(*,*)"D_norm: ",D_norm + WRITE(*,*)"tauk: ",tauk + WRITE(*,*)"iota_e: ",iota_e + WRITE(*,*)"Delta_prime (r_s ref.): ",sl_in%Re_dp_arr(k), + $ sl_in%Im_dp_arr(k) + WRITE(*,*)"Delta_crit: ",sl_in%d_crit_arr(k) + + ! (Deltaprime - Delta_crit)/S^(1/3), with the complex + ! Deltaprime diagonal so that the uncoupled matching is + ! the diagonal limit of the coupled determinant. + delta_eff = (CMPLX(sl_in%Re_dp_arr(k),sl_in%Im_dp_arr(k), + $ KIND=r8) - sl_in%d_crit_arr(k)) + $ /(sl_in%lu_arr(k)**(1.0_r8/3.0_r8)) + pe = 0.0 + +c----------------------------------------------------------------------- +c uncoupled AMR scan (one surface at a time). +c dispersion_AMR_v2 (growthrates_mod) populates Q_store, D_store. +c----------------------------------------------------------------------- + IF (AMR_flag .AND. .NOT. coupling_flag) THEN + + WRITE(*,'(A,I0,A)') 'Calling uncoupled AMR scan on m = ', + $ sl_in%qval_arr(k),' rational surface:' + ! The tearing root usually sits near the electron + ! drift frequency, Re(Q) = -Q_e. + IF (ABS(-Q_e - scan_center) > scan_width) THEN + WRITE(*,'(A,ES10.3,A)') + $ ' WARNING: Re(Q) = -Q_e =',-Q_e, + $ ' is outside the scan window.' + WRITE(*,'(A)') ' Increase scan_width or set ' + $ //'scan_center, or the root may be missed.' + END IF + CALL dispersion_AMR_v2(n_k,sl_in,msing_max,scan_width, + $ Q_num,AMR_passes,coupling_flag,scan_center) + + ! Zeros of the dispersion function in this surface's + ! own ExB frame; the most unstable one is reported as + ! gamma_sol = gamma [1/s] + i*omega [rad/s]. + CALL extract_roots(n_k,sl_in,msing_max,coupling_flag, + $ sl_in%Qconv_arr(k),roots_out(k)) + IF (roots_out(k)%n_root > 0) THEN + sl_out%gamma_sol_arr(k) = CMPLX( + $ -AIMAG(roots_out(k)%q_root(1)), + $ REAL(roots_out(k)%q_root(1)), KIND=r8) + $ / sl_in%Qconv_arr(k) + END IF + + ! Re-allocate output arrays for this surface. + IF (ALLOCATED(all_deltas_out(k)%inQs)) + $ DEALLOCATE(all_deltas_out(k)%inQs) + IF (ALLOCATED(all_deltas_out(k)%iinQs)) + $ DEALLOCATE(all_deltas_out(k)%iinQs) + IF (ALLOCATED(all_deltas_out(k)%real_deltas)) + $ DEALLOCATE(all_deltas_out(k)%real_deltas) + IF (ALLOCATED(all_deltas_out(k)%imag_deltas)) + $ DEALLOCATE(all_deltas_out(k)%imag_deltas) + + ALLOCATE(all_deltas_out(k)%inQs(n_pts), + $ all_deltas_out(k)%iinQs(n_pts)) + ALLOCATE(all_deltas_out(k)%real_deltas(n_pts), + $ all_deltas_out(k)%imag_deltas(n_pts)) + + ! Flatten unique AMR points into 1-D output arrays. + DO i = 1, n_pts + all_deltas_out(k)%inQs(i) = REAL(Q_store(i)) + all_deltas_out(k)%iinQs(i) = -AIMAG(Q_store(i)) + + all_deltas_out(k)%real_deltas(i) = REAL(D_store(i)) + all_deltas_out(k)%imag_deltas(i) = AIMAG(D_store(i)) + END DO + + ! Clean up temporary AMR memory. + DEALLOCATE(Q_store, D_store) + + WRITE(*,'(A,I2,A,I7,A,2ES14.6)') + $ ' Surface', k, ': n_pts=', + $ SIZE(all_deltas_out(k)%real_deltas), + $ ' out_chksum=', + $ SUM(all_deltas_out(k)%real_deltas), + $ SUM(all_deltas_out(k)%imag_deltas) + + END IF + +c----------------------------------------------------------------------- +c single-surface stability scan on [Re(Q), Im(Q)] grid. +c Uses riccati_f() (new Fitzpatrick TJ-like formalism) or +c. riccati() (orig. SLAYER/Waelbroeck). +c----------------------------------------------------------------------- + IF ((stabscan_flag)) THEN + WRITE(*,*)"------------------------------------------" + WRITE(*,'(A,F0.1)')' >>> Running [Re(Q),'// + $ 'Im(Q)] scan with Q width = ', + $ scan_width + + ing_step = (2.0*scan_width) / (Q_num - 1) + count = 0 + + ALLOCATE(inQs(1:(Q_num+1)),iinQs(1:Q_num)) + ALLOCATE(deltas(1:(Q_num+1),1:Q_num)) + + DO i = 1, (Q_num+1) + DO j = 1, Q_num + ing_coarse = -scan_width + (i - 1) * ing_step ! added 0.5* + iing_coarse = -scan_width + (j - 1) * ing_step + ! Evaluate riccati function + g_tmp = CMPLX(ing_coarse,iing_coarse) + IF (Pperp_Ptor_flag) THEN + delta=riccati_f() + ELSE + ! the Cole riccati() takes ds, not D_norm + delta=riccati(iing_coarse,Q_e,Q_i,P_perp, + $ c_beta,sl_in%ds_arr(k),tau,pe, + $ iinQ=ing_coarse) + END IF + inQs(i) = ing_coarse + IF (Pperp_Ptor_flag) THEN + iinQs(j) = iing_coarse + ELSE + iinQs(j) = -iing_coarse + END IF + deltas(i,j) = delta + ENDDO + ENDDO + + + IF (k<10) THEN + WRITE(UNIT=sm_str,FMT='(I1)') sl_in%qval_arr(k) + sm_str=ADJUSTL(sm_str) + ELSE + WRITE(UNIT=sm_str,FMT='(I2)') sl_in%qval_arr(k) + ENDIF + + OPEN(UNIT=out_unit,FILE="slayer_stability_n"// + $ TRIM(sn_str)//"m"//TRIM(sm_str)//".out", + $ STATUS="UNKNOWN") + WRITE(out_unit,'(1x,4(a17))') "RE(Q)", + $ "IM(Q)","RE(delta)","IM(delta)" + DO i=1,Q_num+1 + DO j=1,Q_num + WRITE(out_unit,'(1x,4(es17.8e3))') + $ inQs(i),iinQs(j), + $ REAL(deltas(i,j)),AIMAG(deltas(i,j)) + ENDDO + ENDDO + CLOSE(out_unit) + + DEALLOCATE(inQs,iinQs,deltas) + ENDIF + ENDDO + + IF (.NOT. (est_gamma_flag)) THEN + sl_in%d_beta_arr = 0.0_r8 + sl_out%dels_db_arr = CMPLX(0.0_r8, 0.0_r8, KIND=r8) + END IF + +c----------------------------------------------------------------------- +c coupled AMR scan (all surfaces simultaneously). +c dispersion_AMR_v2 (growthrates_mod) with coupling_flag = .TRUE. +c----------------------------------------------------------------------- + IF (AMR_flag .AND. coupling_flag) THEN + + CALL dispersion_AMR_v2(n_k,sl_in,msing_max,scan_width, + $ Q_num,AMR_passes,coupling_flag,scan_center) + + ! Zeros of the coupled determinant. Q is normalized by + ! the first surface's tau_k. The coupled eigenvalue is + ! global, so it is reported through roots_out only and + ! the per-surface gamma_sol_arr is left as NaN. + CALL extract_roots(n_k,sl_in,msing_max,coupling_flag, + $ sl_in%Qconv_arr(1),roots_out(1)) + + ! Re-allocate output arrays. + IF (ALLOCATED(all_deltas_out(1)%inQs)) + $ DEALLOCATE(all_deltas_out(1)%inQs) + IF (ALLOCATED(all_deltas_out(1)%iinQs)) + $ DEALLOCATE(all_deltas_out(1)%iinQs) + IF (ALLOCATED(all_deltas_out(1)%real_deltas)) + $ DEALLOCATE(all_deltas_out(1)%real_deltas) + IF (ALLOCATED(all_deltas_out(1)%imag_deltas)) + $ DEALLOCATE(all_deltas_out(1)%imag_deltas) + + ALLOCATE(all_deltas_out(1)%inQs(n_pts), + $ all_deltas_out(1)%iinQs(n_pts)) + ALLOCATE(all_deltas_out(1)%real_deltas(n_pts), + $ all_deltas_out(1)%imag_deltas(n_pts)) + + ! Flatten unique AMR points into 1-D output arrays. + DO i = 1, n_pts + all_deltas_out(1)%inQs(i) = REAL(Q_store(i)) + + ! scan plane q = omega_hat - i*gamma_hat: store the + ! growth-rate axis as +gamma_hat, as in the uncoupled + ! output above. + all_deltas_out(1)%iinQs(i) = -AIMAG(Q_store(i)) + + all_deltas_out(1)%real_deltas(i) = REAL(D_store(i)) + all_deltas_out(1)%imag_deltas(i) = AIMAG(D_store(i)) + END DO + + ! Clean up temporary AMR memory. + DEALLOCATE(Q_store, D_store) + + END IF + +c----------------------------------------------------------------------- +c coupled-surface stability scan on [Re(Q), Im(Q)] grid. +c Uses dispersion_det (growthrates_mod) for the full dispersion +c determinant including inter-surface coupling. +c----------------------------------------------------------------------- + IF (coupled_stabscan_flag) THEN + WRITE(*,*)"------------------------------------------" + WRITE(*,'(A,F0.1)')' >>> Running [Re(Q),'// + $ 'Im(Q)] determinant scan with radius = ', + $ scan_width + + ing_step = (2.0 * scan_width) / (Q_num - 1) + count = 0 + + !IF (.NOT. stabscan_flag) THEN + ALLOCATE(inQs(1:(Q_num+1)),iinQs(1:Q_num)) + ALLOCATE(deltas(1:(Q_num+1),1:Q_num)) + !END IF + + inQs=0.0; iinQs=0.0; deltas=(0.0,0.0) + + DO i = 1, (Q_num+1) + DO j = 1, Q_num + ! This scan passes g = gamma_hat + i*omega_hat + ! directly, so its imaginary axis is the frequency + ! and carries the scan-window center. + ing_coarse = -scan_width + (i - 1) * ing_step + iing_coarse = scan_center - scan_width + $ + (j - 1) * ing_step + inQs(i) = ing_coarse + iinQs(j) = iing_coarse + + ! Evaluate determinant + g_tmp = CMPLX(ing_coarse,iing_coarse) + deltas(i,j)=dispersion_det(g_tmp,n_k,sl_in,msing_max) + ENDDO + ENDDO + + OPEN(UNIT=out_unit,FILE="slayer_determinants_n"// + $ TRIM(sn_str)//".out", STATUS="UNKNOWN") + WRITE(out_unit,'(1x,4(a17))') "RE(Q)", + $ "IM(Q)","RE(det)","IM(det)" + DO i=1,Q_num+1 + DO j=1,Q_num + WRITE(out_unit,'(1x,4(es17.8e3))') + $ inQs(i),iinQs(j), + $ REAL(deltas(i,j)),AIMAG(deltas(i,j)) + ENDDO + ENDDO + CLOSE(out_unit) + + DEALLOCATE(inQs,iinQs,deltas) + END IF + + CALL output_gamma(est_gamma_flag,m_AMR,sl_in,sl_out, + $ all_deltas_out) + stop + ENDIF ! match_gamma_flag +c----------------------------------------------------------------------- +c Br threshold test scan (analytic). +c For testing & verification only. Scans rotation to find the +c critical radial-field threshold from a simple torque balance. +c----------------------------------------------------------------------- + IF (br_th_flag) THEN + WRITE(*,*)"------------------------------------------" + WRITE(*,*)">>> Computing Br threshold" + + IF (read_eq) THEN + sl_in%chi_p_arr = chi_p_prof + sl_in%chi_t_arr = chi_t_prof + sl_in%kappa_arr = kappa_prof + CALL build_inputs(infile,ncfile,sl_in,zeff,mu_i) + n_k = SIZE(sl_in%qval_arr) + CALL allocate_outputs(n_k,sl_out) + ELSE + n_k = 1 + + chis(1) = chi_p_prof(1) + chis(2) = chi_t_prof(1) + chis(3) = kappa_prof(1) + + CALL params(n_e,t_e,t_i,omega,chis,dr_val,dgeo_val, + $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) + + inQ=Q + inQ_e=Q_e + inQ_i=Q_i + inc_beta=c_beta + inds=ds + intau=tau + + CALL allocate_inputs(n_k,sl_in) + CALL allocate_outputs(n_k,sl_out) + + sl_in%qval_arr = (/ mm /) + sl_in%q_arr = (/ qval /) + sl_in%omegas_arr = (/ omega /) + sl_in%Qexb_arr = (/ REAL(Q,r8) /) + sl_in%shear_arr = (/ sval /) + sl_in%ds_arr = (/ ds /) + sl_in%Q_e_arr = (/ Q_e /) + sl_in%Q_i_arr = (/ Q_i /) + sl_in%psi_n_arr = (/ 0.0 /) + sl_in%Re_dp_arr = (/ 0.0 /) + sl_in%Im_dp_arr = (/ 0.0 /) + sl_in%Re_dp_psi_arr = (/ 0.0 /) + sl_in%Im_dp_psi_arr = (/ 0.0 /) + sl_in%d_crit_arr = (/ 0.0 /) + sl_in%P_perp_arr = (/ P_perp /) + sl_in%P_tor_arr = (/ P_tor /) + sl_in%tau_arr = (/ tau /) + sl_in%D_norm_arr = (/ D_norm /) + sl_in%D_fac_arr = (/ D_fac /) + sl_in%d_beta_arr = (/ d_beta /) + sl_in%gammafac_arr = (/ gamma_fac /) + sl_in%c_beta_arr = (/ c_beta /) + sl_in%lu_arr = (/ lu /) + sl_in%Qconv_arr = (/ tauk /) + END IF +c----------------------------------------------------------------------- +c loop over rational surfaces to compute Br threshold. +c----------------------------------------------------------------------- + delta_n_p = 1e-2 + inum = 200 + inQ_max = 10.0 + inQ_min = -10.0 + + DO k=1,n_k + WRITE(*,*) + WRITE(*,'(A,I0,A)') 'Computing Br threshold on m = ', + $ sl_in%qval_arr(k),' rational surface' + + Q_e = sl_in%Q_e_arr(k) + Q_i = sl_in%Q_i_arr(k) + inQ_e = Q_e + inQ_i = Q_i + inpr = sl_in%P_perp_arr(k) + inc_beta = sl_in%c_beta_arr(k) + ! riccati() takes the Cole ds and the torque balance the + ! unperturbed ExB Q, as at the other call sites above + inds = sl_in%ds_arr(k) + intau = sl_in%tau_arr(k) + Q0 = sl_in%Qexb_arr(k) + + ALLOCATE(inQs(0:inum),deltal(0:inum), + $ jxbl(0:inum),bal(0:inum)) + DO i=0,inum + inQs(i)=inQ_min+(REAL(i)/inum)*(inQ_max-inQ_min) + deltal(i)=riccati(inQs(i),inQ_e,inQ_i, + $ inpr,inc_beta,inds,intau,inpe) + jxbl(i)=-AIMAG(1.0/(deltal(i)+delta_n_p)) + bal(i)=2.0*inpr*(Q0-inQs(i))/jxbl(i) + ENDDO + + ! as in gpec_slayer: ignore non-finite entries (0/0, or a + ! near-zero jxb) and floor at zero with a warning when the + ! balance has no positive maximum, rather than a NaN + iloc=MAXLOC(bal,MASK=(bal==bal .AND. ABS(bal)= 0 .AND. bal(MAX(iloc(1),0))>0.0_r8) THEN + br_th=SQRT(bal(iloc(1))/sl_in%lu_arr(k) + $ *(sl_in%shear_arr(k)**2.0/2.0)) + ELSE + br_th=0.0_r8 + WRITE(*,'(1x,a,i0,a)') "!! WARNING: torque balance "// + $ "has no positive maximum at m=",sl_in%qval_arr(k), + $ "; br_th set to 0" + ENDIF + + sl_out%br_th_arr(k) = br_th + sl_out%gamma_sol_arr(k) = 0.0 + sl_out%gamma_est_arr(k) = 0.0 + sl_out%dels_db_arr(k) = 0.0 + + WRITE(*,'(A,ES12.4)') ' br_th = ', br_th + DEALLOCATE(inQs,deltal,jxbl,bal) + ENDDO +c----------------------------------------------------------------------- +c write output. +c----------------------------------------------------------------------- + CALL output_gamma(est_gamma_flag,m_AMR,sl_in,sl_out, + $ all_deltas_out) + STOP + ENDIF ! br_th_flag +c----------------------------------------------------------------------- c find solutions based on simple torque balance. c----------------------------------------------------------------------- IF (bal_flag)THEN @@ -288,7 +1407,7 @@ PROGRAM slayer ! write components of torque balance IF(ascii_flag)THEN OPEN(UNIT=out_unit,FILE="slayer_bal_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,5(a17))') "inQ","RE(delta)", $ "IM(delta)","jxb","bal" @@ -299,9 +1418,18 @@ PROGRAM slayer CLOSE(out_unit) ENDIF - index=MAXLOC(bal) - Q_sol=inQs(index(1)) - br_th=sqrt(MAXVAL(bal)/lu*(sval**2.0/2.0))*1e4 + ! as in gpec_slayer: ignore non-finite entries and floor at + ! zero with a warning when the balance has no positive maximum + iloc=MAXLOC(bal,MASK=(bal==bal .AND. ABS(bal)= 0 .AND. bal(MAX(iloc(1),0))>0.0_r8) THEN + br_th=sqrt(bal(iloc(1))/lu*(sval**2.0/2.0))*1e4 + ELSE + br_th=0.0_r8 + WRITE(*,*)"!! WARNING: torque balance has no positive "// + $ "maximum; br_th set to 0" + ENDIF WRITE(*,*)"Q_sol=",Q_sol WRITE(*,*)"br_th=",br_th DEALLOCATE(inQs,deltal,jxbl,bal) @@ -327,7 +1455,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_stability_n"// - $ TRIM(sn)//".out", STATUS="UNKNOWN") + $ TRIM(sn_str)//".out", STATUS="UNKNOWN") WRITE(out_unit,'(1x,4(a17))') "RE(Q)", $ "IM(Q)","RE(delta)","IM(delta)" DO i=0,inum @@ -364,7 +1492,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_riccatiscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,5(a17))') "x","yphs","yamp", $ "RE(delta)","IM(delta)" DO j=0,jnum @@ -406,7 +1534,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_QPescan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Q","Pe","RE(delta)", $ "IM(delta)","psi","jxb" DO j=0,jnum @@ -421,7 +1549,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE='slayer_QPescan_n' - $ //TRIM(sn)//'.bin', + $ //TRIM(sn_str)//'.bin', $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -461,7 +1589,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_QPscan"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Q","Pr","RE(delta)", $ "IM(delta)","psi","jxb" DO j=0,jnum @@ -476,7 +1604,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE="slayer_QPscan_" - $ //TRIM(sn)//".bin", + $ //TRIM(sn_str)//".bin", $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -510,7 +1638,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_Qscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Q","Pr","RE(delta)", $ "IM(delta)","psi","jxb" DO j=0,jnum @@ -552,7 +1680,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_QPscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Q","Pr","RE(delta)", $ "IM(delta)","psi","jxb" DO j=0,jnum @@ -567,7 +1695,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE="slayer_QPscan_n"// - $ TRIM(sn)//".bin", + $ TRIM(sn_str)//".bin", $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -620,7 +1748,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE='slayer_QDscan_n'// - $ TRIM(sn)//'.bin', + $ TRIM(sn_str)//'.bin', $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -651,8 +1779,9 @@ PROGRAM slayer DO k=0,knum ks(j,k)=k_min+(k_max-k_min)*(REAL(k)/knum) - CALL params(n_e*ks(j,k),t_e,t_i,omega*js(j,k), - $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) + CALL params(n_e*ks(j,k),t_e,t_i,omega*js(j,k),chis, + $ dr_val,dgeo_val,l_n,l_t,qval,sval,bt,rs, + $ R0,mu_i,zeff,params_check) inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -683,8 +1812,9 @@ PROGRAM slayer jxb=-AIMAG(1.0/(delta+delta_n_p)) bal(i)=2.0*inpr*(Q0-inQs(i))/jxb ENDDO - index=MAXLOC(bal) - Q_sols(j,k)=inQs(index(1)) + iloc=MAXLOC(bal) + iloc(1)=iloc(1)-1 ! MAXLOC counts from 1, inQs from 0 + Q_sols(j,k)=inQs(iloc(1)) br_ths(j,k)=sqrt(MAXVAL(bal)/lu)*1e4 WRITE(*,*)"br_ths=",br_ths(j,k) ENDDO @@ -692,7 +1822,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_onscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Omega","Density", $ "Omega_i","Omega_e","Omega_sol","Field_Threshold" DO j=0,jnum @@ -709,7 +1839,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE='slayer_onscan_n'// - $ TRIM(sn)//'.bin', + $ TRIM(sn_str)//'.bin', $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -739,8 +1869,8 @@ PROGRAM slayer ks(j,k)=k_min+(k_max-k_min)*(REAL(k)/knum) CALL params(n_e,t_e*ks(j,k),t_i*ks(j,k), - $ omega*js(j,k),l_n,l_t,qval,sval,bt, - $ rs,R0,mu_i,zeff,params_check) + $ omega*js(j,k),chis,dr_val,dgeo_val,l_n,l_t,qval, + $ sval,bt,rs,R0,mu_i,zeff,params_check) inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -769,8 +1899,9 @@ PROGRAM slayer jxb=-AIMAG(1.0/(delta+delta_n_p)) bal(i)=2.0*inpr*(Q0-inQs(i))/jxb ENDDO - index=MAXLOC(bal) - Q_sols(j,k)=inQs(index(1)) + iloc=MAXLOC(bal) + iloc(1)=iloc(1)-1 ! MAXLOC counts from 1, inQs from 0 + Q_sols(j,k)=inQs(iloc(1)) br_ths(j,k)=sqrt(MAXVAL(bal)/lu)*1e4 WRITE(*,*)"t_e=",t_e*ks(j,k),"br_ths=",br_ths(j,k) ENDDO @@ -778,7 +1909,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_otscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Omega","Temperature", $ "Omega_i","Omega_e","Omega_sol","Field_Threshold" DO j=0,jnum @@ -795,7 +1926,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE='slayer_otscan_n'// - $ TRIM(sn)//'.bin', + $ TRIM(sn_str)//'.bin', $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -825,7 +1956,8 @@ PROGRAM slayer ks(j,k)=k_min+(k_max-k_min)*(REAL(k)/knum) CALL params(n_e*ks(j,k),t_e*js(j,k),t_i*js(j,k),omega, - $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) + $ chis,dr_val,dgeo_val,l_n,l_t,qval,sval,bt,rs,R0, + $ mu_i,zeff,params_check) inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -854,8 +1986,9 @@ PROGRAM slayer jxb=-AIMAG(1.0/(delta+delta_n_p)) bal(i)=2.0*inpr*(Q0-inQs(i))/jxb ENDDO - index=MAXLOC(bal) - Q_sols(j,k)=inQs(index(1)) + iloc=MAXLOC(bal) + iloc(1)=iloc(1)-1 ! MAXLOC counts from 1, inQs from 0 + Q_sols(j,k)=inQs(iloc(1)) br_ths(j,k)=sqrt(MAXVAL(bal)/lu)*1e4 WRITE(*,*)"br_ths=",br_ths(j,k) ENDDO @@ -863,7 +1996,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_ntscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Temperature","Density", $ "Omega_i","Omega_e","Omega_sol","Field_Threshold" DO j=0,jnum @@ -879,7 +2012,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE='slayer_ntscan_n'// - $ TRIM(sn)//'.bin', + $ TRIM(sn_str)//'.bin', $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -909,9 +2042,9 @@ PROGRAM slayer ks(j,k)=k_min+(k_max-k_min)*(REAL(k)/knum) - CALL params(n_e*ks(j,k),t_e,t_i,omega, - $ l_n,l_t,qval,sval,bt*js(j,k),rs,R0,mu_i,zeff, - $ params_check) + CALL params(n_e*ks(j,k),t_e,t_i,omega,chis,dr_val, + $ dgeo_val,l_n,l_t,qval,sval,bt*js(j,k),rs,R0,mu_i, + $ zeff,params_check) inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -941,8 +2074,9 @@ PROGRAM slayer jxb=-AIMAG(1.0/(delta+delta_n_p)) bal(i)=2.0*inpr*(Q0-inQs(i))/jxb ENDDO - index=MAXLOC(bal) - Q_sols(j,k)=inQs(index(1)) + iloc=MAXLOC(bal) + iloc(1)=iloc(1)-1 ! MAXLOC counts from 1, inQs from 0 + Q_sols(j,k)=inQs(iloc(1)) br_ths(j,k)=sqrt(MAXVAL(bal)/lu)*1e4 WRITE(*,*)"br_ths=",br_ths(j,k) ENDDO @@ -950,7 +2084,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_nbtscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,4(a17))') "Bt","Density", $ "Omega_sol","Field_Threshold" DO j=0,jnum @@ -966,7 +2100,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE='slayer_nbtscan_n'// - $ TRIM(sn)//'.bin', + $ TRIM(sn_str)//'.bin', $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -978,5 +2112,27 @@ PROGRAM slayer DEALLOCATE(inQs,bal,js,ks,Q_sols,br_ths) ENDIF + CONTAINS +c----------------------------------------------------------------------- +c report a namelist input that this run needs but was left unset, +c with the default v1.5 used silently in its place. +c----------------------------------------------------------------------- + SUBROUTINE slayer_unset(name,val,v15,nbad) + + CHARACTER(*), INTENT(IN) :: name ! namelist variable + REAL(r8), INTENT(IN) :: val ! value read + CHARACTER(*), INTENT(IN) :: v15 ! v1.5 default + INTEGER, INTENT(INOUT) :: nbad ! running count + + IF (nbad == 0) THEN + WRITE(*,*) 'ERROR: slayer.in leaves required inputs unset:' + END IF + WRITE(*,'(4x,a,a,es10.3,a,a,a)') name,' = ',val, + $ ' (v1.5 default: ',v15,')' + nbad = nbad+1 + + RETURN + END SUBROUTINE slayer_unset + END PROGRAM slayer diff --git a/slayer/slayer_netcdf.f b/slayer/slayer_netcdf.f new file mode 100644 index 000000000..e2ce42a98 --- /dev/null +++ b/slayer/slayer_netcdf.f @@ -0,0 +1,562 @@ +c======================================================================= +c file slayer_netcdf.f +c Writes SLAYER solver results to a NetCDF output file. +c +c The per-surface scalar inputs (Lundquist number, Q-normalisation, +c Prandtl numbers, …) and the solver outputs (growth rates, Delta +c values) are stored in a single NetCDF-3/64-bit-offset file named +c slayer_output_n.nc. +c +c Ragged AMR scan data (variable number of evaluation points per +c surface) are zero-padded into rectangular arrays before writing. +c +c (Resolved: FLAG 1 -- early RETURN for msing == 0. +c FLAG 2 -- buffers now use fill_val; _FillValue attributes added. +c FLAG 3 -- Q_id removed (unused). +c FLAG 4 -- c_b_id removed (unused). +c FLAG 5 -- stale unused-variable comment removed. +c FLAG 6 -- version from INCLUDE "version.inc". +c FLAG 7 -- local sn_local replaces global sn_str.) +c======================================================================= +c----------------------------------------------------------------------- +c code organisation. +c----------------------------------------------------------------------- +c 0. slayer_netcdf_mod -- module declarations +c 1. sl_check -- NetCDF status checker +c 2. slayer_netcdf_out -- main output routine +c----------------------------------------------------------------------- +c +c----------------------------------------------------------------------- +c subprogram 0. slayer_netcdf_mod. +c Module wrapper — imports sglobal_mod (shared types and globals) +c and the NetCDF Fortran-90 API. +c----------------------------------------------------------------------- + MODULE slayer_netcdf_mod + + USE sglobal_mod + USE netcdf + USE ieee_arithmetic, ONLY: ieee_value, ieee_quiet_nan + + IMPLICIT NONE + + CONTAINS +c +c----------------------------------------------------------------------- +c subprogram 1. sl_check. +c Assert that a NetCDF operation succeeded; abort with a message +c if it did not. +c----------------------------------------------------------------------- + SUBROUTINE sl_check(stat) +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + INTEGER, INTENT(IN) :: stat ! return code from any nf90_* call +c----------------------------------------------------------------------- +c check status and abort on error. +c----------------------------------------------------------------------- + IF (stat /= nf90_noerr) THEN + PRINT *, TRIM(nf90_strerror(stat)) + STOP "ERROR: failed to write/read netcdf file" + ENDIF + + RETURN + END SUBROUTINE sl_check +c +c----------------------------------------------------------------------- +c subprogram 2. slayer_netcdf_out. +c Write per-surface SLAYER inputs and solver outputs to the +c NetCDF file slayer_output_n.nc . +c +c Arguments: +c msing -- number of rational surfaces +c m_AMR -- number of AMR-scanned surfaces +c est_gamma_flag -- .TRUE. to include estimated growth rates +c sl_in -- slayer_inputs_type (per-surface inputs) +c sl_out -- slayer_outputs_type (solver results) +c all_deltas_out -- array(m_AMR) of deltas_outputs_type (AMR +c scan results, potentially ragged) +c----------------------------------------------------------------------- + SUBROUTINE slayer_netcdf_out(msing, m_AMR, est_gamma_flag, + $ sl_in, sl_out, all_deltas_out) +c----------------------------------------------------------------------- +c declarations -- subroutine arguments. +c----------------------------------------------------------------------- + INTEGER, INTENT(IN) :: msing ! number of rational surfaces + INTEGER, INTENT(IN) :: m_AMR ! number of AMR surfaces + LOGICAL, INTENT(IN) :: est_gamma_flag ! include estimated gammas? + + TYPE(slayer_inputs_type), INTENT(IN) :: sl_in + TYPE(slayer_outputs_type), INTENT(IN) :: sl_out + TYPE(deltas_outputs_type), INTENT(IN) :: all_deltas_out(m_AMR) + +c----------------------------------------------------------------------- +c declarations -- NetCDF file and dimension IDs. +c----------------------------------------------------------------------- + INTEGER :: ncid ! NetCDF file ID + INTEGER :: qsing_dim ! dim: rational surfaces (msing) + INTEGER :: nAMR_dim ! dim: AMR surfaces (m_AMR) + INTEGER :: i_dim ! dim: Re/Im component (2) + INTEGER :: dim_pts_id ! dim: max AMR eval points + +c----------------------------------------------------------------------- +c declarations -- NetCDF variable IDs. +c Each *_id holds the handle returned by nf90_def_var and is +c later passed to the matching nf90_put_var call. +c----------------------------------------------------------------------- + INTEGER :: qsing_id ! "r" — surface index + INTEGER :: qr_id ! "q_rational" — safety factor + INTEGER :: omegas_id ! "omegas" — rotation freq + INTEGER :: qc_id ! "tau_k" — Q-conversion + INTEGER :: Q_e_id ! "Q_e" — norm Q_e + INTEGER :: Q_i_id ! "Q_i" — norm Q_i + INTEGER :: S_id ! "S" — Lundquist + INTEGER :: pr_id ! "psi_n_rational" — norm psi + INTEGER :: p_perp_id ! "P_perp" + INTEGER :: p_tor_id ! "P_tor" + INTEGER :: Dnorm_id ! "D" — normalised D + INTEGER :: dpp_id ! "Delta_prime_rational" (complex) + INTEGER :: dprs_id ! "Delta_prime_rs_rational" (complex) + INTEGER :: kref_id ! "k_ref" + INTEGER :: qsh_id ! "Q_shift" + INTEGER :: alpha_id ! "alpha_mercier" + INTEGER :: dc_id ! "Delta_crit_rational" + INTEGER :: dels_db_id ! "delta_s_d_b" (complex) + INTEGER :: d_b_id ! "d_beta" + INTEGER :: gs_id ! "growth rate" (complex) + INTEGER :: ge_id ! "est. growth rate" (complex) + INTEGER :: br_th_id ! "br_th" — Br threshold + +c AMR variable IDs + INTEGER :: var_q_id ! "Q_AMR" — scan Q-points + INTEGER :: var_d_id ! "Deltas_AMR" — scan Delta values + INTEGER :: var_npts_id ! "n_amr_pts" — points per surface + +c root variable IDs (zeros and poles of the dispersion function) + INTEGER :: root_dim ! dim: stored roots per scan + INTEGER :: nroot_id ! "n_roots" + INTEGER :: npole_id ! "n_poles" + INTEGER :: qroot_id ! "Q_roots" + INTEGER :: qpole_id ! "Q_poles" + INTEGER :: groot_id ! "gamma_roots" + INTEGER :: oroot_id ! "omega_roots" + INTEGER :: proot_id ! "root_polished" + LOGICAL :: have_roots ! roots_out holds this run's roots + LOGICAL :: have_amr ! AMR scan data are present + REAL(r8), ALLOCATABLE :: buffer_r(:,:,:) ! (root, surf, Re/Im) + REAL(r8), ALLOCATABLE :: buffer_p(:,:,:) ! (root, surf, Re/Im) + REAL(r8), ALLOCATABLE :: buffer_g(:,:) ! gamma (root, surf) + REAL(r8), ALLOCATABLE :: buffer_o(:,:) ! omega (root, surf) + INTEGER, ALLOCATABLE :: buffer_f(:,:) ! polished (root, surf) + +c----------------------------------------------------------------------- +c declarations -- AMR rectangular-buffer workspace. +c The ragged per-surface scan data are padded into fixed-size +c rectangular arrays before writing to NetCDF. +c----------------------------------------------------------------------- + INTEGER :: max_pts_all ! max points across surfaces + INTEGER :: s ! surface loop index + INTEGER :: n_curr ! points on current surface + REAL(r8), ALLOCATABLE :: buffer_q(:,:,:) ! (pts, surf, Re/Im) + REAL(r8), ALLOCATABLE :: buffer_d(:,:,:) ! (pts, surf, Re/Im) + INTEGER, ALLOCATABLE :: n_pts_arr(:) ! points per surface + REAL(r8) :: fill_val ! NaN padding for ragged arrays + +c----------------------------------------------------------------------- +c declarations -- miscellaneous locals. +c----------------------------------------------------------------------- + CHARACTER(64) :: ncfile ! output file name + CHARACTER(2) :: sn_local ! local n-string for filename + LOGICAL, PARAMETER :: debug_flag = .FALSE. ! verbose trace + INCLUDE "version.inc" + +c----------------------------------------------------------------------- +c build the output filename from the toroidal mode number. +c----------------------------------------------------------------------- + IF (debug_flag) PRINT *, "Called slayer_netcdf_out" + + IF (nn < 10) THEN + WRITE(UNIT=sn_local, FMT='(I1)') nn + sn_local = ADJUSTL(sn_local) + ELSE + WRITE(UNIT=sn_local, FMT='(I2)') nn + ENDIF + ncfile = "slayer_output_n"//TRIM(sn_local)//".nc" + IF (debug_flag) PRINT *, ncfile + +c----------------------------------------------------------------------- +c create the NetCDF file (clobber any existing file). +c----------------------------------------------------------------------- + IF (debug_flag) PRINT *, " - Creating netcdf file" + CALL sl_check( nf90_create(ncfile, + $ cmode=OR(NF90_CLOBBER, NF90_64BIT_OFFSET), ncid=ncid) ) + +c----------------------------------------------------------------------- +c reform ragged AMR Delta outputs into rectangular arrays. +c +c Each surface may have a different number of AMR scan points. +c We find the maximum, allocate rectangular buffers of that size, +c and copy in the per-surface data. Unused trailing slots are +c filled with fill_val (-9.99E33). +c----------------------------------------------------------------------- + +c step 1: find the maximum AMR grid size across all surfaces. + max_pts_all = 0 + have_amr = .FALSE. + IF (m_AMR > 0) THEN + have_amr = ALLOCATED(all_deltas_out(1)%inQs) + END IF + IF (have_amr) THEN + DO s = 1, m_AMR + max_pts_all = MAX(max_pts_all, + $ SIZE(all_deltas_out(s)%inQs)) + END DO + END IF + +c step 2: allocate rectangular buffers (pts × surfaces × Re/Im). + ALLOCATE(buffer_q(max_pts_all, m_AMR, 2)) + ALLOCATE(buffer_d(max_pts_all, m_AMR, 2)) + ALLOCATE(n_pts_arr(m_AMR)) + + fill_val = ieee_value(1.0d0, ieee_quiet_nan) + buffer_q = fill_val ! NaN padding for ragged arrays + buffer_d = fill_val + n_pts_arr = 0 + +c step 3: flatten the ragged data into the buffers. + IF (have_amr) THEN + DO s = 1, m_AMR + n_curr = SIZE(all_deltas_out(s)%inQs) + n_pts_arr(s) = n_curr + +c Q-coordinate (Re and Im parts) + buffer_q(1:n_curr, s, 1) = + $ all_deltas_out(s)%inQs(1:n_curr) + buffer_q(1:n_curr, s, 2) = + $ all_deltas_out(s)%iinQs(1:n_curr) + +c Delta result (Re and Im parts) + buffer_d(1:n_curr, s, 1) = + $ all_deltas_out(s)%real_deltas(1:n_curr) + buffer_d(1:n_curr, s, 2) = + $ all_deltas_out(s)%imag_deltas(1:n_curr) + END DO + END IF + +c step 4: roots and poles, in the same (omega_hat, gamma_hat) +c convention as Q_AMR. + have_roots = .FALSE. + IF (ALLOCATED(roots_out) .AND. have_amr) + $ have_roots = SIZE(roots_out) >= m_AMR + IF (have_roots) THEN + ALLOCATE(buffer_r(MAX_ROOTS, m_AMR, 2)) + ALLOCATE(buffer_p(MAX_ROOTS, m_AMR, 2)) + ALLOCATE(buffer_g(MAX_ROOTS, m_AMR)) + ALLOCATE(buffer_o(MAX_ROOTS, m_AMR)) + ALLOCATE(buffer_f(MAX_ROOTS, m_AMR)) + buffer_f = 0 + DO s = 1, m_AMR + buffer_r(:, s, 1) = REAL(roots_out(s)%q_root) + buffer_r(:, s, 2) = -AIMAG(roots_out(s)%q_root) + buffer_p(:, s, 1) = REAL(roots_out(s)%q_pole) + buffer_p(:, s, 2) = -AIMAG(roots_out(s)%q_pole) + buffer_g(:, s) = fill_val + buffer_o(:, s) = fill_val + IF (roots_out(s)%tauk > 0.0_r8) THEN + buffer_g(:, s) = buffer_r(:, s, 2)/roots_out(s)%tauk + buffer_o(:, s) = buffer_r(:, s, 1)/roots_out(s)%tauk + END IF + WHERE (roots_out(s)%polished) buffer_f(:, s) = 1 + END DO + END IF + +c----------------------------------------------------------------------- +c define global file attributes. +c----------------------------------------------------------------------- + IF (debug_flag) PRINT *, " - Defining netcdf globals" + CALL sl_check( nf90_put_att(ncid, nf90_global, + $ "title", "SLAYER outputs") ) + CALL sl_check( nf90_put_att(ncid, nf90_global, + $ "version", version) ) + +c----------------------------------------------------------------------- +c define dimensions and per-surface NetCDF variables. +c----------------------------------------------------------------------- + IF (debug_flag) PRINT *, " - Defining dimensions in netcdf" + WRITE(*,*) ">>> Writing results to NetCDF output file" + + IF (msing == 0) THEN + WRITE(*,*) "WARNING: msing == 0, skipping NetCDF output" + DEALLOCATE(buffer_q, buffer_d, n_pts_arr) + CALL sl_check( nf90_close(ncid) ) + RETURN + END IF + +c -- core dimensions -- + CALL sl_check( nf90_def_dim(ncid, "r", msing, qsing_dim) ) + IF (have_amr) THEN + CALL sl_check( nf90_def_dim(ncid, "r_AMR", m_AMR, + $ nAMR_dim) ) + END IF + CALL sl_check( nf90_def_dim(ncid, "i", 2, i_dim) ) + +c -- scalar per-surface variables -- + CALL sl_check( nf90_def_var(ncid, "r", nf90_int, + $ qsing_dim, qsing_id) ) + CALL sl_check( nf90_put_att(ncid, qsing_id, "long_name", + $ "Resonant poloidal mode number m of each surface") ) + CALL sl_check( nf90_def_var(ncid, "q_rational", + $ nf90_double, qsing_dim, qr_id) ) + CALL sl_check( nf90_def_var(ncid, "omegas", nf90_double, + $ qsing_dim, omegas_id) ) + CALL sl_check( nf90_put_att(ncid, omegas_id, "long_name", + $ "ExB rotation frequency per unit n [rad/s]") ) + CALL sl_check( nf90_def_var(ncid, "Q_shift", nf90_double, + $ qsing_dim, qsh_id) ) + CALL sl_check( nf90_put_att(ncid, qsh_id, "long_name", + $ "ExB Doppler shift -tau_k*n*Omega_E applied to Q in "// + $ "the coupled determinant (zero when not applied)") ) + CALL sl_check( nf90_def_var(ncid, "tau_k", nf90_double, + $ qsing_dim, qc_id) ) + CALL sl_check( nf90_def_var(ncid, "Q_e", nf90_double, + $ qsing_dim, Q_e_id) ) + CALL sl_check( nf90_def_var(ncid, "Q_i", nf90_double, + $ qsing_dim, Q_i_id) ) + CALL sl_check( nf90_def_var(ncid, "S", nf90_double, + $ qsing_dim, S_id) ) + CALL sl_check( nf90_def_var(ncid, "psi_n_rational", + $ nf90_double, qsing_dim, pr_id) ) + CALL sl_check( nf90_def_var(ncid, "P_perp", nf90_double, + $ qsing_dim, p_perp_id) ) + CALL sl_check( nf90_def_var(ncid, "P_tor", nf90_double, + $ qsing_dim, p_tor_id) ) + +c----------------------------------------------------------------------- +c define additional variables (D, Delta', growth rates). +c----------------------------------------------------------------------- + CALL sl_check( nf90_def_var(ncid, "D", nf90_double, + $ qsing_dim, Dnorm_id) ) + CALL sl_check( nf90_def_var(ncid, "Delta_prime_rational", + $ nf90_double, (/qsing_dim, i_dim/), dpp_id) ) + CALL sl_check( nf90_put_att(ncid, dpp_id, "long_name", + $ "Outer-region Delta_prime diagonal, referenced to unit "// + $ "normalized poloidal flux (as written by STRIDE)") ) + CALL sl_check( nf90_def_var(ncid, "Delta_prime_rs_rational", + $ nf90_double, (/qsing_dim, i_dim/), dprs_id) ) + CALL sl_check( nf90_put_att(ncid, dprs_id, "long_name", + $ "Delta_prime diagonal referenced to the minor radius "// + $ "r_s, K^(2*alpha)*Delta_prime; the value used in the "// + $ "layer matching") ) + CALL sl_check( nf90_def_var(ncid, "k_ref", nf90_double, + $ qsing_dim, kref_id) ) + CALL sl_check( nf90_put_att(ncid, kref_id, "long_name", + $ "Reference-length factor K = r_s*dpsi_N/dr") ) + CALL sl_check( nf90_def_var(ncid, "alpha_mercier", nf90_double, + $ qsing_dim, alpha_id) ) + CALL sl_check( nf90_put_att(ncid, alpha_id, "long_name", + $ "Mercier exponent sqrt(-D_I)") ) + CALL sl_check( nf90_def_var(ncid, "Delta_crit_rational", + $ nf90_double, qsing_dim, dc_id) ) + + IF (est_gamma_flag) THEN + CALL sl_check( nf90_def_var(ncid, "delta_s_d_b", + $ nf90_double, (/qsing_dim, i_dim/), dels_db_id) ) + CALL sl_check( nf90_def_var(ncid, "d_beta", nf90_double, + $ qsing_dim, d_b_id) ) + CALL sl_check( nf90_def_var(ncid, "est. growth rate", + $ nf90_double, (/qsing_dim, i_dim/), ge_id) ) + END IF + + CALL sl_check( nf90_def_var(ncid, "growth rate", + $ nf90_double, (/qsing_dim, i_dim/), gs_id) ) + CALL sl_check( nf90_put_att(ncid, gs_id, "long_name", + $ "Most unstable matched root on each surface: gamma [1/s] "// + $ "(i=1) and omega [rad/s] (i=2), in that surface's ExB "// + $ "frame; NaN where no root was extracted") ) + + IF (ALLOCATED(sl_out%br_th_arr)) THEN + CALL sl_check( nf90_def_var(ncid, "br_th", + $ nf90_double, qsing_dim, br_th_id) ) + END IF + +c----------------------------------------------------------------------- +c define AMR scan dimensions and variables. +c----------------------------------------------------------------------- + IF (have_amr) THEN + CALL sl_check( nf90_def_dim(ncid, "amr_pts", + $ max_pts_all, dim_pts_id) ) + CALL sl_check( nf90_def_var(ncid, "n_amr_pts", NF90_INT, + $ (/nAMR_dim/), var_npts_id) ) + CALL sl_check( nf90_def_var(ncid, "Q_AMR", NF90_DOUBLE, + $ (/dim_pts_id, nAMR_dim, i_dim/), var_q_id) ) + CALL sl_check( nf90_put_att(ncid, var_q_id, "long_name", + $ "AMR scan points, (omega, gamma)*tau_k") ) + CALL sl_check( nf90_def_var(ncid, "Deltas_AMR", NF90_DOUBLE, + $ (/dim_pts_id, nAMR_dim, i_dim/), var_d_id) ) + CALL sl_check( nf90_put_att(ncid, var_d_id, "long_name", + $ "Dispersion function at the AMR scan points") ) + +c roots and poles of the dispersion function from each AMR scan. + IF (have_roots) THEN + CALL sl_check( nf90_def_dim(ncid, "root", MAX_ROOTS, + $ root_dim) ) + CALL sl_check( nf90_def_var(ncid, "n_roots", NF90_INT, + $ (/nAMR_dim/), nroot_id) ) + CALL sl_check( nf90_def_var(ncid, "n_poles", NF90_INT, + $ (/nAMR_dim/), npole_id) ) + CALL sl_check( nf90_def_var(ncid, "Q_roots", NF90_DOUBLE, + $ (/root_dim, nAMR_dim, i_dim/), qroot_id) ) + CALL sl_check( nf90_put_att(ncid, qroot_id, "long_name", + $ "Zeros of the dispersion function, (omega, gamma)*"// + $ "tau_k, sorted by descending gamma") ) + CALL sl_check( nf90_def_var(ncid, "Q_poles", NF90_DOUBLE, + $ (/root_dim, nAMR_dim, i_dim/), qpole_id) ) + CALL sl_check( nf90_put_att(ncid, qpole_id, "long_name", + $ "Poles of the dispersion function, (omega, gamma)*"// + $ "tau_k, AMR cell centers") ) + CALL sl_check( nf90_def_var(ncid, "gamma_roots", + $ NF90_DOUBLE, (/root_dim, nAMR_dim/), groot_id) ) + CALL sl_check( nf90_put_att(ncid, groot_id, "units", "1/s") ) + CALL sl_check( nf90_def_var(ncid, "omega_roots", + $ NF90_DOUBLE, (/root_dim, nAMR_dim/), oroot_id) ) + CALL sl_check( nf90_put_att(ncid, oroot_id, "units", + $ "rad/s") ) + CALL sl_check( nf90_def_var(ncid, "root_polished", NF90_INT, + $ (/root_dim, nAMR_dim/), proot_id) ) + CALL sl_check( nf90_put_att(ncid, proot_id, "long_name", + $ "1 if the secant refinement of the root converged") ) + CALL sl_check( nf90_put_att(ncid, qroot_id, + $ "_FillValue", fill_val) ) + CALL sl_check( nf90_put_att(ncid, qpole_id, + $ "_FillValue", fill_val) ) + CALL sl_check( nf90_put_att(ncid, groot_id, + $ "_FillValue", fill_val) ) + CALL sl_check( nf90_put_att(ncid, oroot_id, + $ "_FillValue", fill_val) ) + END IF + +c set _FillValue attribute on AMR arrays for proper padding. + CALL sl_check( nf90_put_att(ncid, var_q_id, + $ "_FillValue", fill_val) ) + CALL sl_check( nf90_put_att(ncid, var_d_id, + $ "_FillValue", fill_val) ) + END IF + +c----------------------------------------------------------------------- +c end NetCDF define mode. +c----------------------------------------------------------------------- + CALL sl_check( nf90_enddef(ncid) ) + +c----------------------------------------------------------------------- +c write per-surface scalar variables. +c----------------------------------------------------------------------- + CALL sl_check( nf90_put_var(ncid, qsing_id, + $ sl_in%qval_arr) ) + CALL sl_check( nf90_put_var(ncid, qr_id, + $ sl_in%q_arr) ) + CALL sl_check( nf90_put_var(ncid, pr_id, + $ sl_in%psi_n_arr) ) + CALL sl_check( nf90_put_var(ncid, omegas_id, + $ sl_in%omegas_arr)) + CALL sl_check( nf90_put_var(ncid, S_id, + $ sl_in%lu_arr) ) + CALL sl_check( nf90_put_var(ncid, qsh_id, + $ sl_in%q_shift_arr) ) + CALL sl_check( nf90_put_var(ncid, qc_id, + $ sl_in%Qconv_arr) ) + CALL sl_check( nf90_put_var(ncid, Q_e_id, + $ sl_in%Q_e_arr) ) + CALL sl_check( nf90_put_var(ncid, Q_i_id, + $ sl_in%Q_i_arr) ) + CALL sl_check( nf90_put_var(ncid, p_perp_id, + $ sl_in%P_perp_arr)) + CALL sl_check( nf90_put_var(ncid, p_tor_id, + $ sl_in%P_tor_arr) ) + CALL sl_check( nf90_put_var(ncid, Dnorm_id, + $ sl_in%D_norm_arr)) + +c----------------------------------------------------------------------- +c write complex Delta' as a RESHAPE'd (msing, 2) real array. +c----------------------------------------------------------------------- + CALL sl_check( nf90_put_var(ncid, dpp_id, + $ RESHAPE( (/sl_in%Re_dp_psi_arr, sl_in%Im_dp_psi_arr/), + $ (/msing, 2/) )) ) + CALL sl_check( nf90_put_var(ncid, dprs_id, + $ RESHAPE( (/sl_in%Re_dp_arr, sl_in%Im_dp_arr/), + $ (/msing, 2/) )) ) + CALL sl_check( nf90_put_var(ncid, kref_id, + $ sl_in%k_ref_arr) ) + CALL sl_check( nf90_put_var(ncid, alpha_id, + $ sl_in%alpha_arr) ) + CALL sl_check( nf90_put_var(ncid, dc_id, + $ sl_in%d_crit_arr) ) + +c----------------------------------------------------------------------- +c write estimated growth-rate outputs (only when requested). +c----------------------------------------------------------------------- + IF (est_gamma_flag) THEN + CALL sl_check( nf90_put_var(ncid, dels_db_id, + $ RESHAPE( (/REAL(sl_out%dels_db_arr), + $ AIMAG(sl_out%dels_db_arr)/), + $ (/msing, 2/) )) ) + + CALL sl_check( nf90_put_var(ncid, d_b_id, + $ sl_in%d_beta_arr) ) + + CALL sl_check( nf90_put_var(ncid, ge_id, + $ RESHAPE( (/REAL(sl_out%gamma_est_arr), + $ AIMAG(sl_out%gamma_est_arr)/), + $ (/msing, 2/) )) ) + END IF + +c----------------------------------------------------------------------- +c write solved growth rate (always present). +c----------------------------------------------------------------------- + CALL sl_check( nf90_put_var(ncid, gs_id, + $ RESHAPE( (/REAL(sl_out%gamma_sol_arr), + $ AIMAG(sl_out%gamma_sol_arr)/), + $ (/msing, 2/) )) ) + + IF (ALLOCATED(sl_out%br_th_arr)) THEN + CALL sl_check( nf90_put_var(ncid, br_th_id, + $ sl_out%br_th_arr) ) + END IF + +c----------------------------------------------------------------------- +c write AMR scan arrays. +c----------------------------------------------------------------------- + IF (have_amr) THEN + CALL sl_check( nf90_put_var(ncid, var_npts_id, n_pts_arr) ) + CALL sl_check( nf90_put_var(ncid, var_q_id, buffer_q) ) + CALL sl_check( nf90_put_var(ncid, var_d_id, buffer_d) ) + END IF + + IF (have_roots) THEN + CALL sl_check( nf90_put_var(ncid, nroot_id, + $ roots_out(1:m_AMR)%n_root) ) + CALL sl_check( nf90_put_var(ncid, npole_id, + $ roots_out(1:m_AMR)%n_pole) ) + CALL sl_check( nf90_put_var(ncid, qroot_id, buffer_r) ) + CALL sl_check( nf90_put_var(ncid, qpole_id, buffer_p) ) + CALL sl_check( nf90_put_var(ncid, groot_id, buffer_g) ) + CALL sl_check( nf90_put_var(ncid, oroot_id, buffer_o) ) + CALL sl_check( nf90_put_var(ncid, proot_id, buffer_f) ) + DEALLOCATE(buffer_r, buffer_p, buffer_g, buffer_o, buffer_f) + END IF + +c----------------------------------------------------------------------- +c deallocate AMR rectangular buffers. +c----------------------------------------------------------------------- + DEALLOCATE(buffer_q, buffer_d, n_pts_arr) + +c----------------------------------------------------------------------- +c close the NetCDF file. +c----------------------------------------------------------------------- + IF (debug_flag) PRINT *, " - Closing netcdf file" + CALL sl_check( nf90_close(ncid) ) + +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE slayer_netcdf_out + END MODULE slayer_netcdf_mod \ No newline at end of file diff --git a/stride/free.f b/stride/free.f index 6e4ddcab0..e0142a675 100755 --- a/stride/free.f +++ b/stride/free.f @@ -168,7 +168,6 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,op_netcdf_out) INTEGER :: ipert,jpert,isol,info,lwork INTEGER, DIMENSION(mpert) :: m INTEGER, DIMENSION(1) :: imax - REAL(r8) :: v1 REAL(r8), DIMENSION(mpert) :: ep,ev,et REAL(r8), DIMENSION(3*mpert-1) :: rwork COMPLEX(r8) :: phase,norm @@ -176,6 +175,10 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,op_netcdf_out) COMPLEX(r8), DIMENSION(mpert,mpert) :: wt,wpt,wvt COMPLEX(r8), DIMENSION(mpert,mpert) :: nmat,smat CHARACTER(24), DIMENSION(mpert) :: message + REAL(r8), DIMENSION(sq%mx+1) :: ln_q + REAL(r8), DIMENSION(msing) :: shr + TYPE(spline_type) :: shr_spl + REAL(r8) :: v1,respsi c----------------------------------------------------------------------- c write formats. c----------------------------------------------------------------------- @@ -295,11 +298,26 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,op_netcdf_out) WRITE(out_unit,90)(isol,ep(isol),ev(isol),isol=1,mpert) WRITE(out_unit,80) c----------------------------------------------------------------------- +c magnetic shear s = psi_N*d(ln q)/d(psi_N) at the rational surfaces +c (written to the netcdf file for SLAYER). +c----------------------------------------------------------------------- + CALL spline_alloc(shr_spl,SIZE(sq%xs(:))-1,1) + shr_spl%xs=sq%xs(:) + ln_q=LOG(sq%fs(:,4)) + shr_spl%fs(:,1)=ln_q + CALL spline_fit(shr_spl,"extrap") + DO ising=1,msing + respsi=sing(ising)%psifac + CALL spline_eval(shr_spl,respsi,1) + shr(ising) = respsi * shr_spl%f1(1) + ENDDO + CALL spline_dealloc(shr_spl) +c----------------------------------------------------------------------- c optionally write netcdf file. c----------------------------------------------------------------------- IF(present(op_netcdf_out))THEN IF(op_netcdf_out) CALL stride_netcdf_out(wp,wv,wt,ep,ev,et, - $ delta_prime_mat,plasma1,vacuum1,total1) + $ delta_prime_mat,plasma1,vacuum1,total1,shr) ENDIF c----------------------------------------------------------------------- c deallocate diff --git a/stride/stride.F b/stride/stride.F index 585058ce7..10cd23f45 100644 --- a/stride/stride.F +++ b/stride/stride.F @@ -349,7 +349,8 @@ PROGRAM stride total1=0_r8 ALLOCATE(mx0(mpert,mpert),vx0(mpert)) CALL stride_netcdf_out(mx0,mx0,mx0,vx0,vx0,vx0, - $ delta_prime_mat,plasma1,vacuum1,total1) + $ delta_prime_mat,plasma1,vacuum1,total1, + $ SPREAD(0._r8,1,msing)) CALL stride_dealloc ENDIF IF(mat_flag .OR. ode_flag)DEALLOCATE(amat,bmat,cmat,ipiva,jmat) diff --git a/stride/stride_netcdf.f b/stride/stride_netcdf.f index 08a7995d2..939ebf075 100644 --- a/stride/stride_netcdf.f +++ b/stride/stride_netcdf.f @@ -48,23 +48,30 @@ END SUBROUTINE check c----------------------------------------------------------------------- c declarations. c----------------------------------------------------------------------- - SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp,pl1,va1,to1) + SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp, + $ pl1,va1,to1,shr) REAL(r8), DIMENSION(mpert), INTENT(IN) :: epi,evi,eti COMPLEX(r8), DIMENSION(mpert,mpert), INTENT(IN) :: wp,wv,wt COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: dp REAL(r8), INTENT(IN) :: pl1,va1,to1 + REAL(r8), DIMENSION(msing), INTENT(IN) :: shr INTEGER :: i, ncid, $ i_dim, m_dim, mo_dim, p_dim, i_id, m_id, mo_id, p_id, $ f_id, q_id, dv_id, mu_id, di_id, dr_id, ca_id, $ wp_id, wpv_id, wv_id, wvv_id, wt_id, wtv_id, $ r_dim, rp_dim, l_dim, lp_dim, r_id, rp_id, l_id, lp_id, - $ pr_id, qr_id, dp_id, ap_id, bp_id, gp_id, dpp_id + $ pr_id, qr_id, dp_id, ap_id, bp_id, gp_id, dpp_id, + $ shear_id,resm_id,drr_id,dir_id COMPLEX(r8), DIMENSION(mpert) :: ep,ev,et CHARACTER(2) :: sn CHARACTER(64) :: ncfile + REAL(r8) :: respsi + REAL(r8), DIMENSION(msing) :: dr_rationals,di_rationals + INTEGER, DIMENSION(msing) :: resm + INTEGER :: ising,jsing COMPLEX(r8), DIMENSION(msing,msing) :: ap,bp,gammap,deltap @@ -85,6 +92,15 @@ SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp,pl1,va1,to1) ep = CMPLX(epi, 0.0) ev = CMPLX(evi, 0.0) et = CMPLX(eti, 0.0) + + ! evaluate the Mercier (D_I) and resistive interchange (D_R) + ! indices on rational surfaces + DO i=1,msing + respsi = sing(i)%psifac + CALL spline_eval(locstab,respsi,0) + di_rationals(i)=locstab%f(1)/respsi + dr_rationals(i)=locstab%f(2)/respsi + END DO c----------------------------------------------------------------------- c open files c----------------------------------------------------------------------- @@ -109,6 +125,7 @@ SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp,pl1,va1,to1) CALL check( nf90_put_att(ncid,nf90_global,'mhigh', mhigh)) CALL check( nf90_put_att(ncid,nf90_global,'mpert', mpert)) CALL check( nf90_put_att(ncid,nf90_global,'mband', mband)) + CALL check( nf90_put_att(ncid,nf90_global,'msing', msing)) CALL check( nf90_put_att(ncid,nf90_global,'psilow', psilow)) CALL check( nf90_put_att(ncid,nf90_global,'amean', amean)) CALL check( nf90_put_att(ncid,nf90_global,'rmean', rmean)) @@ -158,6 +175,9 @@ SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp,pl1,va1,to1) CALL check( nf90_def_dim(ncid, "psi_n", sq%mx+1, p_dim) ) CALL check( nf90_def_var(ncid, "psi_n", nf90_double, p_dim, p_id)) IF(msing>0)THEN + DO i=1,msing + resm(i)=NINT(sing(i)%q*nn) + ENDDO CALL check( nf90_def_dim(ncid,"lr_index",2*msing,l_dim) ) CALL check( nf90_def_var(ncid,"lr_index",nf90_int,l_dim,l_id)) CALL check( nf90_put_att(ncid,l_dim,"long_name", @@ -183,6 +203,18 @@ SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp,pl1,va1,to1) $ r_dim,qr_id) ) CALL check( nf90_put_att(ncid,qr_id,"long_name", $ "Safety Factor at Rational Surfaces") ) + CALL check( nf90_def_var(ncid, "shear", nf90_double, r_dim, + $ shear_id) ) + CALL check( nf90_def_var(ncid, "resm", nf90_int, r_dim, + $ resm_id) ) + CALL check( nf90_def_var(ncid, "dr_rational", nf90_double, + $ r_dim, drr_id) ) + CALL check( nf90_put_att(ncid,drr_id,"long_name", + $ "Resistive Interchange Index D_R at Rational Surfaces") ) + CALL check( nf90_def_var(ncid, "di_rational", nf90_double, + $ r_dim, dir_id) ) + CALL check( nf90_put_att(ncid,dir_id,"long_name", + $ "Mercier Index D_I at Rational Surfaces") ) ENDIF ! define variables IF(debug_flag) PRINT *," - Defining variables in netcdf" @@ -246,6 +278,10 @@ SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp,pl1,va1,to1) $ i=1,msing)/)) ) CALL check( nf90_put_var(ncid,qr_id, (/(sing(i)%q, $ i=1,msing)/)) ) + CALL check( nf90_put_var(ncid,shear_id, shr) ) + CALL check( nf90_put_var(ncid,resm_id, resm) ) + CALL check( nf90_put_var(ncid,drr_id, dr_rationals)) + CALL check( nf90_put_var(ncid,dir_id, di_rationals)) ENDIF IF(debug_flag) PRINT *," - Putting profile variables in netcdf"