Handle native signed modes in POTATO torque

POTATO/CMakeLists.txt

@@ -88,6 +88,16 @@ add_test(NAME test_find_all_roots_bracketed
    COMMAND test_find_all_roots_bracketed.x
 )
 
+add_executable(test_signed_toroidal_bound.x
+   SRC/test_signed_toroidal_bound.f90
+)
+target_link_libraries(test_signed_toroidal_bound.x
+   potato_base
+)
+add_test(NAME test_signed_toroidal_bound
+   COMMAND test_signed_toroidal_bound.x
+)
+
 add_executable(test_wall_loss.x
    SRC/test_wall_loss.f90
 )
@@ -153,8 +163,16 @@ if(Python3_Interpreter_FOUND)
         COMMAND test_neo2_reader.x "${_neo2_reader_fixture}")
     set_tests_properties(test_neo2_reader PROPERTIES
         DEPENDS generate_neo2_reader_oracle)
-endif()
-if(Python3_Interpreter_FOUND)
+    add_test(NAME potato_bmod_converter
+        COMMAND ${Python3_EXECUTABLE}
+            ${CMAKE_CURRENT_SOURCE_DIR}/test/bmod_converter/test_bmod_converter.py
+    )
+    set_tests_properties(potato_bmod_converter PROPERTIES TIMEOUT 60)
+    add_test(NAME potato_profile_converter
+        COMMAND ${Python3_EXECUTABLE}
+            ${CMAKE_CURRENT_SOURCE_DIR}/test/profile_converter/test_profile_converter.py
+    )
+    set_tests_properties(potato_profile_converter PROPERTIES TIMEOUT 60)
     add_test(NAME potato_invariant_handoff
         COMMAND ${Python3_EXECUTABLE}
             ${CMAKE_CURRENT_SOURCE_DIR}/test/invariant_handoff/test_invariant_handoff.py

POTATO/README.md

@@ -29,6 +29,51 @@ the LCFS into the scrape-off layer. Set it explicitly to `.false.` for a
 closed-flux-only run. Edge extension does not turn the convex wall into the
 LCFS or remove the outer-domain requirement.
 
+For a MARS field already converted to Boozer harmonics, generate POTATO's
+single-`n` cylindrical perturbation without the legacy converter's radial
+Gaussian filter:
+
+```bash
+python tools/boozer_npz_to_bmod_n.py chartmap.nc components.npz bmod_n.dat \
+    --component total --n-tor=+3 --s-max=0.95
+```
+
+Here `chartmap.nc` supplies the accepted Boozer-surface geometry and
+`components.npz` is the provenance product from `rmp_torque mars_to_boozer`.
+The signed `n` must also be used as `n_tor` in `potato.in`. The converter writes
+a JSON sidecar, performs no smoothing or fit, uses `s_tor` explicitly, writes
+zero outside the outer mapped surface, and adds a zero-valued rectangular
+margin so the POTATO spline is not normally evaluated at a clamped nonzero
+boundary. Production target orbits must remain inside that mapped surface.
+The displayed `n=+3` is the current TC24 MARS-to-CCW result:
+`phi_CCW=-phi_MARS` and `n_CCW=-RNTOR=+3`. It is not a universal MARS
+default; another field must supply its own signed laboratory-coordinate map.
+
+The matching profile converter also keeps the coordinate and electric-field
+conventions explicit:
+
+```bash
+python tools/neo_rt_profiles_to_potato.py profile.in plasma.in components.npz \
+    profile_poly.in --r0-cm=641.0903942 --psi-span-tm2=11.88279543 \
+    --relation-sign=1
+```
+
+It selects the physical ion using charge and nonzero density, maps
+`s_tor -> s_pol=rho_pol^2`, and records the polynomial residuals. The JSON
+sidecar also records the signed relation
+`dPhi/dpsi_pol = relation_sign*Omega_E/c`; opposite signs are separate
+convention diagnostics, never an unrecorded curve flip. For the displayed
+right-handed direct-EQDSK chart, the native field equations give
+`Omega_phi,E=c*dPhi/dpsi_pol`, so `relation_sign=+1` maps an already-CCW
+`Omega_E` profile. The current TC24 NEO profile has first been serialized as
+`Omega_E,CCW=-Omega_E,MARS`; the converter does not apply that MARS-to-CCW
+change a second time. The first value in the NEO `plasma.in` header is a
+validated radial-row count, not a required value of 50. The compiled NEO-RT
+profile schema needs only `s_tor,M_t`; the converter derives
+`v_th=sqrt(2*T_i/m_i)` from the charge-selected `plasma.in` species. If a
+legacy third `v_th` column is present, it is checked against that source rather
+than used as an independent temperature profile.
+
 ## Running with OpenMP
 
 The grid build and the per-mode root search run in parallel with OpenMP. Use one

POTATO/SRC/CMakeLists.txt

@@ -34,6 +34,7 @@ add_library(potato_base
   sorting.f90
   binsrc.f90
   bmod_pert.f90
+  resonance_mode_bounds_mod.f90
   resonant_int.f90
   eqmagprofs.f90
   box_counting.f90

POTATO/SRC/resonance_mode_bounds_mod.f90

@@ -0,0 +1,30 @@
+module resonance_mode_bounds_mod
+    implicit none
+
+contains
+
+    pure function resonant_delphi_bound(m_modes, n_modes) result(bound)
+        integer, intent(in) :: m_modes(:), n_modes(:)
+        double precision :: bound
+        double precision, parameter :: pi = 3.14159265358979d0
+
+        ! The signed toroidal mode remains in m*Omega_b+n*Omega_phi=0 and in the
+        ! Fourier phase.  Only this symmetric search extent is a magnitude.  Using
+        ! signed n here made the native n<0 case produce a negative extent and
+        ! silently discard every otherwise valid resonance contribution.
+        bound = 2.d0*pi*(maxval(abs(dble(m_modes))/abs(dble(n_modes))) &
+            + 1.d0/dble(minval(abs(n_modes))))
+    end function resonant_delphi_bound
+
+    pure logical function canonical_flux_outside_lcfs(psi_star, psi_axis, psi_edge)
+        double precision, intent(in) :: psi_star, psi_axis, psi_edge
+
+        ! Outside means beyond the edge in the axis-to-edge flux direction.  The
+        ! earlier psi_star < psi_edge test was valid only when psi decreased from
+        ! axis to edge and rejected every ITER resonance after field_eq selected
+        ! the opposite native GEQDSK gauge.
+        canonical_flux_outside_lcfs = &
+            (psi_star - psi_edge)*(psi_edge - psi_axis) > 0.d0
+    end function canonical_flux_outside_lcfs
+
+end module resonance_mode_bounds_mod

POTATO/SRC/resonant_int.f90

@@ -57,6 +57,7 @@ module resint_mod
 !$omp threadprivate(nperp_max,delint_mode,respoints_jp,respoints_all, &
 !$omp               respoints_all_tmp,respoint,resline_unit,resline_diag_unit, &
 !$omp               resline_unit_is_private,resline_diag_unit_is_private)
+
 end module resint_mod
 !
 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -170,7 +171,8 @@ subroutine integrate_class_resonances
     use logging_mod,        only : tee_message
     use interp_cache_mod,   only : interp_cache_reset
     use field_sub,          only : psif
-    use field_eq_mod,       only : psi_sep
+    use field_eq_mod,       only : psi_axis,psi_sep
+    use resonance_mode_bounds_mod, only : canonical_flux_outside_lcfs
     use, intrinsic :: ieee_arithmetic, only : ieee_is_finite
     !
     implicit none
@@ -286,7 +288,7 @@ subroutine integrate_class_resonances
             !
             dpsiastdx=dpsiast_dRst*delta_R*dxi_dx !$\difp{\psi^\ast}{x}$
             !
-            if(psiast_res.lt.psi_sep) then
+            if(canonical_flux_outside_lcfs(psiast_res,psi_axis,psi_sep)) then
                 ! SOL resonance (rho_pol > 1): the orbit leaves the field domain,
                 ! so pertham/find_bounce cannot close it and would only grind to
                 ! the integrator cap before returning zero.  It is also outside
@@ -484,6 +486,7 @@ subroutine resonant_torque
     use get_matrix_mod,    only : iclass,delphi_max
     use form_classes_doublecount_mod, only : nclasses
     use orbit_dim_mod,     only : numbasef
+    use resonance_mode_bounds_mod, only : resonant_delphi_bound
     use resint_mod,        only : nmodes,marr,narr,delint_mode,respoints_jp,respoints_all,nperp_max, &
         respoints_all_tmp,respoint,resline_unit,resline_diag_unit, &
         resline_unit_is_private,resline_diag_unit_is_private
@@ -534,7 +537,7 @@ subroutine resonant_torque
     ! Bound the class root search to the resonant range: |delphi_b| = 2*pi*|m|/n
     ! at a resonance, so nothing past max|m|/n can resonate.  One n-step margin
     ! keeps the extreme-m root safely inside the trimmed domain.
-    delphi_max=2.d0*pi*(maxval(abs(dble(marr))/dble(narr))+1.d0/dble(minval(narr)))
+    delphi_max=resonant_delphi_bound(marr,narr)
     write(msg, '(A,ES14.6)') &
         'class root search bounded to |delphi_b| <= ', delphi_max
     call tee_message(trim(msg))

POTATO/SRC/test_signed_toroidal_bound.f90

@@ -0,0 +1,27 @@
+program test_signed_toroidal_bound
+    use resonance_mode_bounds_mod, only: resonant_delphi_bound, &
+        canonical_flux_outside_lcfs
+    implicit none
+
+    integer, parameter :: m_modes(7) = [-3, -2, -1, 0, 1, 2, 3]
+    integer, parameter :: n_positive(7) = 3
+    integer, parameter :: n_negative(7) = -3
+    double precision :: positive_bound, negative_bound
+
+    positive_bound = resonant_delphi_bound(m_modes, n_positive)
+    negative_bound = resonant_delphi_bound(m_modes, n_negative)
+
+    if (positive_bound <= 0.d0) error stop "positive n produced nonpositive bound"
+    if (negative_bound <= 0.d0) error stop "negative n produced nonpositive bound"
+    if (positive_bound /= negative_bound) &
+        error stop "search bound depends on toroidal-mode sign"
+
+    if (canonical_flux_outside_lcfs(0.5d0, 0.d0, 1.d0)) &
+        error stop "inside point rejected for increasing flux"
+    if (.not. canonical_flux_outside_lcfs(1.1d0, 0.d0, 1.d0)) &
+        error stop "outside point accepted for increasing flux"
+    if (canonical_flux_outside_lcfs(-0.5d0, 0.d0, -1.d0)) &
+        error stop "inside point rejected for decreasing flux"
+    if (.not. canonical_flux_outside_lcfs(-1.1d0, 0.d0, -1.d0)) &
+        error stop "outside point accepted for decreasing flux"
+end program test_signed_toroidal_bound

POTATO/test/bmod_converter/test_bmod_converter.py

@@ -0,0 +1,121 @@
+#!/usr/bin/env python3
+"""Self-contained regression for the no-filter Boozer-to-POTATO converter."""
+
+from __future__ import annotations
+
+import importlib.util
+import json
+import sys
+import tempfile
+from pathlib import Path
+
+import h5py
+import numpy as np
+from scipy.interpolate import RegularGridInterpolator
+
+
+TOOL = Path(__file__).parents[2] / "tools" / "boozer_npz_to_bmod_n.py"
+SPEC = importlib.util.spec_from_file_location("boozer_npz_to_bmod_n", TOOL)
+MODULE = importlib.util.module_from_spec(SPEC)
+assert SPEC.loader is not None
+SPEC.loader.exec_module(MODULE)
+
+
+def main() -> int:
+    with tempfile.TemporaryDirectory() as directory:
+        root = Path(directory)
+        chartmap = root / "chartmap.nc"
+        components = root / "components.npz"
+        output = root / "bmod_n.dat"
+
+        s = np.linspace(0.01, 0.81, 9)
+        theta = np.linspace(0.0, 2.0*np.pi, 32, endpoint=False)
+        zeta = np.array([0.0, np.pi])
+        radius = 20.0*np.sqrt(s)
+        r = 160.0 + radius[:, None]*np.cos(theta)
+        z_plane = radius[:, None]*np.sin(theta)
+        toroidal_shift = 0.08*np.sin(theta)
+        x = np.empty((2, theta.size, s.size))
+        y = np.empty_like(x)
+        z = np.empty_like(x)
+        for index, angle in enumerate(zeta):
+            # The chart is at constant Boozer zeta, but geometric phi differs
+            # by a manufactured, poloidally varying toroidal shift.
+            phi_geom = angle-toroidal_shift
+            x[index] = (r*np.cos(phi_geom)).T
+            y[index] = (r*np.sin(phi_geom)).T
+            z[index] = z_plane.T
+        with h5py.File(chartmap, "w") as handle:
+            handle.attrs["zeta_convention"] = "boozer"
+            handle["s"] = s
+            handle["rho"] = np.sqrt(s)
+            handle["theta"] = theta
+            handle["zeta"] = zeta
+            handle["x"] = x
+            handle["y"] = y
+            handle["z"] = z
+
+        modes = np.array([-2, 2])
+        coefficients = np.column_stack((0.002*s, (0.01+0.003j)*s))
+        np.savez(
+            components, boozer_s=s, boozer_m=modes,
+            boozer_total=coefficients,
+        )
+        metadata = MODULE.convert(
+            chartmap, components, output, component="total", n_tor=-3,
+            s_max=0.81, nrad=257, nzet=257, ntheta=256,
+            margin_fraction=0.1,
+        )
+        rad, zet, values = MODULE.read_bmod_n(output)
+        if metadata["signed_toroidal_mode"] != -3:
+            raise AssertionError("signed toroidal harmonic was not preserved")
+        if json.loads(output.with_suffix(".dat.json").read_text()) != metadata:
+            raise AssertionError("metadata sidecar differs from returned provenance")
+        if values.shape != (257, 257) or not np.all(np.isfinite(values)):
+            raise AssertionError("invalid POTATO output grid")
+        surface_errors = metadata["gridding_relative_l2_by_surface"]
+        if not np.array_equal(surface_errors["s_tor"], s):
+            raise AssertionError("surface-error provenance lost the mapped radial grid")
+        if (
+            len(surface_errors["relative_l2"]) != s.size
+            or not np.all(np.isfinite(surface_errors["relative_l2"]))
+            or len(surface_errors["absolute_l2_tesla"]) != s.size
+            or len(surface_errors["reference_l2_tesla"]) != s.size
+        ):
+            raise AssertionError("invalid per-surface gridding-error provenance")
+
+        # Check the resolved interior surfaces after the R-Z gridding step.
+        sample = RegularGridInterpolator((rad, zet), values)
+        theta_test = np.linspace(0.0, 2.0*np.pi, 128, endpoint=False)
+        errors = []
+        references = []
+        for surface in (2, 4, 6):
+            points = np.column_stack((
+                160.0 + radius[surface]*np.cos(theta_test),
+                radius[surface]*np.sin(theta_test),
+            ))
+            expected = np.sum(
+                coefficients[surface, :, None]
+                * np.exp(1j*np.outer(modes, theta_test)), axis=0,
+            )
+            expected *= np.exp(1j*(-3)*0.08*np.sin(theta_test))
+            actual = sample(points)
+            errors.append(actual-expected)
+            references.append(expected)
+        relative_l2 = np.linalg.norm(np.concatenate(errors))/np.linalg.norm(
+            np.concatenate(references)
+        )
+        if relative_l2 > 2.0e-2:
+            raise AssertionError(f"R-Z reconstruction error is {relative_l2:.3e}")
+        if metadata["toroidal_angle_transform"] != (
+                "A_RZ=A_B*exp(i*n*(phi_B-phi_geom))"):
+            raise AssertionError("missing Boozer-to-geometric toroidal-angle provenance")
+        if not np.isclose(metadata["toroidal_shift_radians_max_abs"], 0.08):
+            raise AssertionError("manufactured toroidal shift was not recorded")
+        if values[0, 0] != 0.0j or values[-1, -1] != 0.0j:
+            raise AssertionError("outside-map grid margin must be zero")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())