diff --git a/docs/examples.rst b/docs/examples.rst index e9f79fa2..f83b9830 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 bc4e65dc..04b3fbb6 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:: @@ -39,4 +39,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_kinetic_example/README b/docs/examples/DIIID_kinetic_example/README index beb77909..26653414 100644 --- a/docs/examples/DIIID_kinetic_example/README +++ b/docs/examples/DIIID_kinetic_example/README @@ -30,4 +30,4 @@ Finally, we run PENTRC for the full neoclassical torque calculation:: PENTRC START => ... PENTRC STOP => Normal termination. -The full general-aspect-ratio (``fgar``) and trapped-only (``tgar``) torque results can be compared to assess the trapped-particle dominance. The ``data_dir`` in ``pentrc.in`` points to the precomputed special function matrix used by the reduced large-aspect-ratio (``rlar``) method, which serves as a cross-check. +The full general-aspect-ratio (``fgar``) and trapped-only (``tgar``) torque results can be compared to assess the trapped-particle dominance. The ``data_dir`` in ``pentrc.in`` points to the precomputed special function matrix used by the reduced large-aspect-ratio (``rlar``) method, which can be enabled with ``rlar_flag=t`` as a cross-check. diff --git a/docs/index.rst b/docs/index.rst index aa1df6a0..09be708f 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 f7b48839..43813b04 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 4629dade..491cd44d 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 a0d80b2a..22d35140 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 00000000..5da0559f --- /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/sign_conventions.rst b/docs/sign_conventions.rst index 40dfccb4..25290474 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 59e5c59c..7265bca2 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 038d6918..5c376751 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 9f1c7e3c..57b50304 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 069ae193..7eb32a60 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -3047,6 +3047,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)+ diff --git a/input/dcon.in b/input/dcon.in index 8db58813..57ec14a1 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 86f3b7f3..6189d78c 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/install/README b/install/README index 992eaeb0..3b3dd094 100644 --- a/install/README +++ b/install/README @@ -45,8 +45,8 @@ Personal copies of GPEC should be cloned from the `Github repository &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