Equilibrium/KineticForces - BUGFIX! - Port the Fortran 1.6.0 kinetic profile, bounce-integral and torque-method corrections

docs/src/conventions.md

@@ -379,7 +379,7 @@ Here:
 
 - ``\omega_b`` and ``\omega_D`` are the bounce and magnetic-precession frequency **coefficients at
   thermal energy** (``x = 1``): ``\omega_b \propto v_\mathrm{th}`` and
-  ``\omega_D \propto v_\mathrm{th}^2`` (`wbhat`, `wdhat` in `Torque.jl`). Their velocity-space
+  ``\omega_D \propto v_\mathrm{th}^2`` (`wbbar`, `wdbar` times `bhat`, `dhat` in `BounceAveraging.jl`). Their velocity-space
   dependence enters explicitly as ``\sqrt{x}`` (bounce ``\propto v``) and ``x`` (precession
   ``\propto E``).
 - ``\ell_\mathrm{eff}`` is the effective bounce harmonic (`PitchIntegration.jl`):

examples/DIIID-like_SLAYER_example/gpec.toml

@@ -18,7 +18,7 @@ psihigh = 0.9995               # Upper limit of normalized flux coordinate
 mpsi = 0                       # Number of radial grid points (0 = auto-compute from psi_accuracy)
 psi_accuracy = 0.001           # Target absolute error in q for auto-mpsi
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 
 # Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning]

examples/DIIID-like_error_field_example/gpec.toml

@@ -22,7 +22,7 @@ psihigh = 0.995                # Upper limit of normalized poloidal flux (captur
 mpsi = 0                       # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 impurity_charge = 6                     # Impurity charge (carbon)
 impurity_mass = 12                      # Impurity mass [proton masses]

examples/DIIID-like_gal_resistive_example/gpec.toml

@@ -17,7 +17,7 @@ psihigh = 0.993                # Upper limit of normalized flux coordinate (0.99
 mpsi = 128                     # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 
 # Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning]

examples/DIIID-like_gal_resistive_pe_example/gpec.toml

@@ -18,7 +18,7 @@ psihigh = 0.993                # Upper limit of normalized flux coordinate (0.99
 mpsi = 128                     # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 
 # Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning]

examples/DIIID-like_ideal_example/gpec.toml

@@ -17,7 +17,7 @@ psihigh = 0.995                # Upper limit of normalized poloidal flux (captur
 mpsi = 0                       # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 impurity_charge = 6                     # Impurity charge (carbon)
 impurity_mass = 12                      # Impurity mass [proton masses]

examples/DIIID-like_ideal_multi_n_example/gpec.toml

@@ -17,7 +17,7 @@ psilow = 1e-4                  # Lower limit of normalized poloidal flux
 psihigh = 0.9                  # Upper limit of normalized poloidal flux
 mpsi = 128                     # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 
 # Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning]

examples/DIIID-like_riccati_deltaprime_example/gpec.toml

@@ -17,7 +17,7 @@ psihigh = 0.995                # Upper limit of normalized poloidal flux (captur
 mpsi = 0                       # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Override for on-axis safety factor (0 = use input value)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 
 # Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning]

examples/DIIID-like_truncation_example/gpec.toml

@@ -20,7 +20,7 @@ psihigh = 0.9999               # Upper limit of normalized poloidal flux (past t
 mpsi = 0                       # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace each flux surface
 
 [ForceFreeStates]

examples/Solovev_ideal_example/gpec.toml

@@ -13,7 +13,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 [Wall]

examples/Solovev_ideal_example_3D/gpec.toml

@@ -12,7 +12,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 # Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning]

examples/Solovev_ideal_example_multi_n/gpec.toml

@@ -13,7 +13,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 [Wall]

examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml

@@ -13,7 +13,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 [Wall]

examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml

@@ -13,7 +13,7 @@ psilow = 1e-4                           # Lower limit of normalized poloidal flu
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
 mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                               # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7                # Relative tolerance of the ODE solves that trace each flux surface
 
 [Wall]

examples/a10_kinetic_example/gpec.toml

@@ -19,7 +19,7 @@ psilow = 1e-3                  # Lower limit of normalized poloidal flux
 psihigh = 0.99                 # Upper limit of normalized poloidal flux
 mpsi = 16                      # Number of radial grid points (low value for fast iteration)
 mtheta = 256                   # Number of poloidal grid points
-newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value; q is positive, so negative values are rejected)
 flux_surface_rtol = 1e-7       # Relative tolerance of the ODE solves that trace each flux surface
 impurity_charge = 6            # Impurity charge
 impurity_mass = 12             # Impurity mass

src/Deck.jl

@@ -49,6 +49,8 @@ const _MOVED_KEYS = Dict{String,Dict{String,String}}(
         "impurity_charge" => "moved to [Equilibrium]",
         "impurity_mass" => "moved to [Equilibrium]",
         "nn" => "removed; KineticForces takes the toroidal mode number from the run",
+        "rlar_flag" => "removed; the large-aspect-ratio methods are not ported, use fgar_flag or tgar_flag",
+        "clar_flag" => "removed; the large-aspect-ratio methods are not ported, use fgar_flag or tgar_flag",
         "kinetic_file" => _KINETIC_FILE_HOME),
     "SLAYER" => Dict{String,String}(
         "enabled" => "a [Tearing] section with inner_model = \"SLAYER\" runs the tearing analysis",

src/Equilibrium/Axisymmetric/Direct.jl

@@ -609,17 +609,13 @@ Direct equilibrium driver: locates the axis, traces each flux surface with `fiel
     if q0 <= 0.0
         @warn "q0 extrapolation to axis gives q0 = $(@sprintf("%.3f", q0)) ≤ 0 — likely a spline artifact from psilow being too large; check psilow or use newq0 to override."
     end
-    # newq0 = -1 flips the extrapolated q0; resolved locally so the config keeps the user's request.
     newq0 = equil_params.newq0
-    if newq0 == -1
-        newq0 = -q0
-    end
     if newq0 != 0.0
         @info "Revising q-profile for newq0 = $(@sprintf("%.3f", newq0))"
         f0 = profiles.F_spline.y[1] - profiles.F_deriv(psi_nodes[1]; hint=Ref(1)) * psi_nodes[1]
         f0fac = f0^2 * ((newq0 / q0)^2 - 1.0)
         for i in 1:(mpsi+1)
-            ffac = sqrt(1.0 + f0fac / profiles.F_spline.y[i]^2) * sign(newq0)
+            ffac = sqrt(1.0 + f0fac / profiles.F_spline.y[i]^2)
             sq_nodes[i, 1] *= ffac
             sq_nodes[i, 4] *= ffac
             rzphi_nodes[i, :, 3] .*= ffac

src/Equilibrium/Axisymmetric/Inverse.jl

@@ -283,16 +283,13 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No
     f1_sq_lo = sq_deriv(sq_xs[1])
     f1_sq_hi = sq_deriv(sq_xs[end])
     q0 = f_sq[1, 4] - f1_sq_lo[4] * sq_xs[1]
-    if newq0 == -1
-        newq0 = -q0
-    end
 
     if newq0 != 0
-        f0 = f_sq[1, 2] - f1_sq_lo[2] * sq_xs[1]
+        f0 = f_sq[1, 1] - f1_sq_lo[1] * sq_xs[1]
         f0fac = f0^2 * ((newq0 / q0)^2 - 1)
         q0 = newq0
         for ipsi in 0:mpsi
-            ffac = sqrt(1 + f0fac / f_sq[ipsi+1, 1]^2) * sign(newq0)
+            ffac = sqrt(1 + f0fac / f_sq[ipsi+1, 1]^2)
             sq_fs[ipsi+1, 1] *= ffac
             sq_fs[ipsi+1, 4] *= ffac
             rzphi_fs[ipsi+1, :, 3] *= ffac

src/Equilibrium/Config.jl

@@ -53,7 +53,7 @@ specified in the input.
     two-pass auto grid (knot count scales as τ^(-1/3))
   - `mtheta::Int` - Number of poloidal grid points
   - `newq0::Float64` - Target on-axis safety factor q(0); the q and F profiles are rescaled to
-    meet it (0 = use input value, -1 = use the axis extrapolation with its sign flipped)
+    meet it (0 = use input value); q is positive, so negative values are rejected
   - `flux_surface_rtol::Float64` - Relative tolerance of the ODE solves that trace each flux surface; their absolute tolerance
     is `equil_abstol(flux_surface_rtol)`, i.e. `flux_surface_rtol` capped at `EQUIL_ABSTOL_MAX`
   - `imas_cocos::Int` - Expected COCOS convention of an IMAS `dd.equilibrium` input (11 = IMAS standard, 2 = GPEC internal)
@@ -174,6 +174,7 @@ specified in the input.
             error("Cannot recognize jac_type = $(jac_type)")
         end
         (isfinite(flux_surface_rtol) && flux_surface_rtol > 0) || error("flux_surface_rtol = $flux_surface_rtol must be a positive, finite relative tolerance")
+        newq0 >= 0 || error("newq0 = $newq0 is not supported: q must be positive. Use newq0 = 0 (keep the input q) or newq0 > 0.")
         if psihigh > 1.0
             @warn "psihigh = $psihigh exceeds 1.0 (separatrix); clamping to 1.0"
         end