diff --git a/POTATO/CMakeLists.txt b/POTATO/CMakeLists.txt index 492540e7..d19436c4 100644 --- a/POTATO/CMakeLists.txt +++ b/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 diff --git a/POTATO/README.md b/POTATO/README.md index 365fe7ac..f5b8584e 100644 --- a/POTATO/README.md +++ b/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 diff --git a/POTATO/SRC/CMakeLists.txt b/POTATO/SRC/CMakeLists.txt index 40381191..82ae5d78 100644 --- a/POTATO/SRC/CMakeLists.txt +++ b/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 diff --git a/POTATO/SRC/resonance_mode_bounds_mod.f90 b/POTATO/SRC/resonance_mode_bounds_mod.f90 new file mode 100644 index 00000000..635d0840 --- /dev/null +++ b/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 diff --git a/POTATO/SRC/resonant_int.f90 b/POTATO/SRC/resonant_int.f90 index 8375a69f..ada36fa6 100644 --- a/POTATO/SRC/resonant_int.f90 +++ b/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)) diff --git a/POTATO/SRC/test_signed_toroidal_bound.f90 b/POTATO/SRC/test_signed_toroidal_bound.f90 new file mode 100644 index 00000000..4ad5c7fa --- /dev/null +++ b/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 diff --git a/POTATO/test/bmod_converter/test_bmod_converter.py b/POTATO/test/bmod_converter/test_bmod_converter.py new file mode 100644 index 00000000..9ebef96c --- /dev/null +++ b/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()) diff --git a/POTATO/test/profile_converter/test_profile_converter.py b/POTATO/test/profile_converter/test_profile_converter.py new file mode 100644 index 00000000..2f5567a4 --- /dev/null +++ b/POTATO/test/profile_converter/test_profile_converter.py @@ -0,0 +1,111 @@ +#!/usr/bin/env python3 +"""Regression for explicit NEO-RT-to-POTATO profile and frequency signs.""" + +from __future__ import annotations + +import importlib.util +import sys +import tempfile +from pathlib import Path + +import numpy as np + + +TOOL = Path(__file__).parents[2] / "tools" / "neo_rt_profiles_to_potato.py" +SPEC = importlib.util.spec_from_file_location("neo_rt_profiles_to_potato", 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) + profile = root / "profile.in" + plasma = root / "plasma.in" + geometry = root / "components.npz" + profile.write_text( + "0.0 0.02\n" + "0.5 0.03\n" + "1.0 0.04\n" + ) + plasma_grid = np.linspace(0.0, 1.0, 239) + plasma.write_text( + "% N am1 am2 Z1 Z2\n" + f"{plasma_grid.size} 2.0 3.0 1.0 1.0\n" + "% s n1 n2 T1 T2 Te\n" + + "".join( + f"{s:.16e} {(1.0-0.4*s)*1.0e14:.16e} 0.0 " + f"{(1.0-0.4*s)*1.0e4:.16e} 1.0 " + f"{(1.2-0.4*s)*1.0e4:.16e}\n" + for s in plasma_grid + ) + ) + s = np.linspace(0.0, 1.0, 11) + np.savez(geometry, s_sqrt_poloidal=np.sqrt(s), s_toroidal=s) + + metadata = {} + for sign in (-1, 0, 1): + output = root / f"profile_{sign:+d}.in" + metadata[sign] = MODULE.convert( + profile, plasma, geometry, output, r0_cm=100.0, + psi_span_tm2=2.0, relation_sign=sign, degree=3, + ) + if metadata[sign]["selected_ion_index_one_based"] != 1: + raise AssertionError("physical ion was not selected by charge and density") + if metadata[sign]["plasma_grid_count"] != 239: + raise AssertionError("non-default NEO-RT grid count was not preserved") + + plus = metadata[1]["potential_statv_edge"] + minus = metadata[-1]["potential_statv_edge"] + zero = metadata[0]["potential_statv_edge"] + if not plus > 0.0 or not np.isclose(minus, -plus) or zero != 0.0: + raise AssertionError("potential-gradient relation signs were not preserved") + zero_coefficients = np.asarray( + MODULE.numeric_rows(root / "profile_+0.in")[3] + ) + if not np.array_equal(zero_coefficients, np.zeros(4)): + raise AssertionError("zero-frequency control produced a nonzero potential") + + profile_with_vth = root / "profile_with_vth.in" + temperatures = (1.0-0.4*plasma_grid)*1.0e4 + vth = np.sqrt( + 2.0*temperatures*MODULE.EV_CGS/(2.0*MODULE.AMU_CGS) + ) + mach = np.interp(plasma_grid, [0.0, 0.5, 1.0], [0.02, 0.03, 0.04]) + profile_with_vth.write_text( + "".join( + f"{s:.16e} {m:.16e} {speed:.16e}\n" + for s, m, speed in zip(plasma_grid, mach, vth, strict=True) + ) + ) + with_vth = MODULE.convert( + profile_with_vth, plasma, geometry, root / "profile_with_vth.out", + r0_cm=100.0, psi_span_tm2=2.0, relation_sign=1, degree=3, + ) + if with_vth["optional_profile_vth_max_relative_error"] > 1.0e-5: + raise AssertionError("equivalent optional vth column did not close") + if not np.isclose( + with_vth["potential_statv_edge"], + metadata[1]["potential_statv_edge"], + rtol=1.0e-12, + ): + raise AssertionError("two- and three-column profile schemas differ") + + truncated = root / "plasma_truncated.in" + truncated.write_text("\n".join(plasma.read_text().splitlines()[:-1]) + "\n") + try: + MODULE.convert( + profile, truncated, geometry, root / "invalid.in", + r0_cm=100.0, psi_span_tm2=2.0, relation_sign=1, degree=3, + ) + except ValueError as error: + if "declares 239 rows, found 238" not in str(error): + raise + else: + raise AssertionError("truncated NEO-RT plasma grid was accepted") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/POTATO/tools/boozer_npz_to_bmod_n.py b/POTATO/tools/boozer_npz_to_bmod_n.py new file mode 100755 index 00000000..6750d7c1 --- /dev/null +++ b/POTATO/tools/boozer_npz_to_bmod_n.py @@ -0,0 +1,329 @@ +#!/usr/bin/env python3 +"""Map a signed Boozer ``Delta|B|`` spectrum to POTATO's R-Z input grid. + +The converter deliberately does not smooth or fit the perturbation. It uses +the Boozer geometry already stored in a NEO-RT chartmap, reconstructs the +complex single-n amplitude on those surfaces, and linearly interpolates that +amplitude in the poloidal R-Z plane. Values outside the mapped surface are +zero, and the rectangular grid includes an explicit margin so POTATO does not +silently clamp ordinary orbit excursions to a nonzero boundary value. +""" + +from __future__ import annotations + +import argparse +import hashlib +import json +import struct +from pathlib import Path + +import h5py +import numpy as np +from scipy.interpolate import CubicSpline, LinearNDInterpolator, RegularGridInterpolator + + +def _sha256(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as stream: + for chunk in iter(lambda: stream.read(1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def _record(payload: bytes) -> bytes: + size = len(payload) + return struct.pack("=i", size) + payload + struct.pack("=i", size) + + +def _text_attribute(value: object) -> str: + """Return a scalar HDF5 text attribute as a normal Python string.""" + array = np.asarray(value) + if array.ndim != 0: + raise ValueError("chartmap text attributes must be scalar") + scalar = array.item() + if isinstance(scalar, bytes): + return scalar.decode("utf-8") + return str(scalar) + + +def write_bmod_n(path: Path, rad_cm: np.ndarray, zet_cm: np.ndarray, + amplitude_tesla: np.ndarray) -> None: + """Write POTATO's three-record, Fortran-sequential ``bmod_n.dat``.""" + if amplitude_tesla.shape != (rad_cm.size, zet_cm.size): + raise ValueError("amplitude grid does not match the R-Z axes") + amplitude_gauss = np.asarray(amplitude_tesla, dtype=np.complex128) * 1.0e4 + with path.open("wb") as stream: + stream.write(_record(struct.pack("=ii", rad_cm.size, zet_cm.size))) + stream.write(_record( + np.asarray(rad_cm, dtype=np.float64).tobytes() + + np.asarray(zet_cm, dtype=np.float64).tobytes() + )) + stream.write(_record( + amplitude_gauss.real.tobytes(order="F") + + amplitude_gauss.imag.tobytes(order="F") + )) + + +def read_bmod_n(path: Path) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Read a POTATO perturbation grid; primarily used for round-trip checks.""" + data = path.read_bytes() + offset = 0 + + def record() -> bytes: + nonlocal offset + (size,) = struct.unpack_from("=i", data, offset) + offset += 4 + payload = data[offset:offset + size] + offset += size + (trailer,) = struct.unpack_from("=i", data, offset) + offset += 4 + if trailer != size: + raise ValueError("Fortran record markers differ") + return payload + + nrad, nzet = struct.unpack("=ii", record()) + axes = np.frombuffer(record(), dtype=np.float64) + values = np.frombuffer(record(), dtype=np.float64) + if offset != len(data) or axes.size != nrad + nzet or values.size != 2*nrad*nzet: + raise ValueError("invalid bmod_n.dat record sizes") + rad = axes[:nrad].copy() + zet = axes[nrad:].copy() + count = nrad*nzet + real = values[:count].reshape((nrad, nzet), order="F") + imag = values[count:].reshape((nrad, nzet), order="F") + return rad, zet, (real + 1j*imag)/1.0e4 + + +def convert(chartmap: Path, components: Path, output: Path, *, component: str, + n_tor: int, s_max: float, nrad: int, nzet: int, + ntheta: int, margin_fraction: float) -> dict: + if n_tor == 0: + raise ValueError("the perturbation toroidal mode must be nonzero") + if not 0.0 < s_max <= 1.0: + raise ValueError("s_max must be in (0, 1]") + if nrad < 8 or nzet < 8: + raise ValueError("the R-Z grid must have at least 8 points per direction") + if ntheta < 2: + raise ValueError("ntheta must be at least two") + if margin_fraction <= 0.0: + raise ValueError("margin_fraction must be positive") + + with h5py.File(chartmap, "r") as handle: + s_geometry = np.asarray(handle["s"], dtype=float) + theta = np.asarray(handle["theta"], dtype=float) + zeta = np.asarray(handle["zeta"], dtype=float) + x = np.asarray(handle["x"], dtype=float) + y = np.asarray(handle["y"], dtype=float) + z = np.asarray(handle["z"], dtype=float) + if "zeta_convention" not in handle.attrs: + raise ValueError("chartmap lacks the required zeta_convention attribute") + zeta_convention = _text_attribute(handle.attrs["zeta_convention"]) + if zeta_convention != "boozer": + raise ValueError( + f"chartmap zeta_convention must be 'boozer', got {zeta_convention!r}" + ) + expected = (zeta.size, theta.size, s_geometry.size) + if x.shape != expected or y.shape != expected or z.shape != expected: + raise ValueError(f"chartmap x/y/z must have shape {expected}") + if zeta.size < 1: + raise ValueError("chartmap zeta axis must not be empty") + if (np.any(np.diff(s_geometry) <= 0.0) + or np.any(np.diff(theta) <= 0.0) + or np.any(np.diff(zeta) <= 0.0)): + raise ValueError("chartmap s, theta, and zeta axes must increase strictly") + + with np.load(components) as data: + s_spectrum = np.asarray(data["boozer_s"], dtype=float) + modes = np.asarray(data["boozer_m"], dtype=int) + key = f"boozer_{component}" + if key not in data: + raise KeyError(f"{components} has no {key}") + coefficients = np.asarray(data[key], dtype=np.complex128) + if coefficients.shape != (s_spectrum.size, modes.size): + raise ValueError("Boozer coefficient matrix has inconsistent dimensions") + if np.any(np.diff(s_spectrum) <= 0.0): + raise ValueError("Boozer spectrum surfaces must increase strictly") + + selected = s_geometry <= s_max + if np.count_nonzero(selected) < 3: + raise ValueError("s_max leaves fewer than three chartmap surfaces") + s_map = s_geometry[selected] + r_chart = np.hypot(x[0][:, selected], y[0][:, selected]).T + z_chart = z[0][:, selected].T + theta_map = np.linspace(0.0, 2.0*np.pi, ntheta, endpoint=False) + theta_closed = np.append(theta, 2.0*np.pi) + r_map = CubicSpline( + theta_closed, np.column_stack((r_chart, r_chart[:, 0])), axis=1, + bc_type="periodic", + )(theta_map) + z_map = CubicSpline( + theta_closed, np.column_stack((z_chart, z_chart[:, 0])), axis=1, + bc_type="periodic", + )(theta_map) + + # The spectrum amplitude multiplies exp(i*n*phi_B), whereas POTATO's R-Z + # amplitude multiplies exp(i*n*phi_geom). The chartmap is sampled on the + # first constant-Boozer-zeta plane, whose cylindrical points generally do + # not all have phi_geom=zeta[0]. Transform the coefficient independently + # of the signed mode/helicity choice: + # + # A_RZ = A_B * exp(i*n*(phi_B - phi_geom)). + # + # Interpolate the unit phasor rather than a wrapped angle, then explicitly + # restore unit modulus so the coordinate map cannot alter the amplitude. + phi_geom = np.arctan2(y[0][:, selected], x[0][:, selected]).T + toroidal_shift_phasor = np.exp(1j*n_tor*(zeta[0] - phi_geom)) + toroidal_shift_map = CubicSpline( + theta_closed, + np.column_stack((toroidal_shift_phasor, toroidal_shift_phasor[:, 0])), + axis=1, + bc_type="periodic", + )(theta_map) + shift_modulus = np.abs(toroidal_shift_map) + if np.any(shift_modulus <= np.finfo(float).tiny): + raise ValueError("toroidal-angle phase interpolation produced zero modulus") + toroidal_shift_map /= shift_modulus + + # Piecewise-linear radial interpolation is exact at the input surfaces and + # introduces no radial smoothing. The unresolved axis interval is set to + # zero, consistent with regular nonaxisymmetric harmonics there. + coeff_map = np.zeros((s_map.size, modes.size), dtype=np.complex128) + in_spectrum = s_map >= s_spectrum[0] + for index in range(modes.size): + coeff_map[in_spectrum, index] = ( + np.interp(s_map[in_spectrum], s_spectrum, coefficients[:, index].real) + + 1j*np.interp(s_map[in_spectrum], s_spectrum, coefficients[:, index].imag) + ) + amplitude_map = ( + coeff_map @ np.exp(1j*np.outer(modes, theta_map)) + ) * toroidal_shift_map + + points = np.column_stack((r_map.ravel(), z_map.ravel())) + values = amplitude_map.ravel() + r_span = float(np.ptp(r_map)) + z_span = float(np.ptp(z_map)) + rad_cm = np.linspace( + float(np.min(r_map) - margin_fraction*r_span), + float(np.max(r_map) + margin_fraction*r_span), nrad, + ) + zet_cm = np.linspace( + float(np.min(z_map) - margin_fraction*z_span), + float(np.max(z_map) + margin_fraction*z_span), nzet, + ) + rr, zz = np.meshgrid(rad_cm, zet_cm, indexing="ij") + interpolator = LinearNDInterpolator(points, values, fill_value=0.0j) + amplitude_grid = np.asarray(interpolator(rr, zz), dtype=np.complex128) + if not np.all(np.isfinite(amplitude_grid)): + raise ValueError("non-finite values produced by R-Z interpolation") + + output.parent.mkdir(parents=True, exist_ok=True) + write_bmod_n(output, rad_cm, zet_cm, amplitude_grid) + round_rad, round_zet, round_values = read_bmod_n(output) + roundtrip_atol = np.finfo(float).eps*max(float(np.max(np.abs(amplitude_grid))), 1.0) + if not (np.array_equal(round_rad, rad_cm) and np.array_equal(round_zet, zet_cm) + and np.allclose(round_values, amplitude_grid, rtol=2.0e-16, + atol=roundtrip_atol)): + raise RuntimeError("bmod_n.dat round-trip check failed") + + # Quantify only gridding error at mapped input points. The no-filter + # radial and Fourier reconstruction happens before this interpolation. + grid_interp = RegularGridInterpolator( + (rad_cm, zet_cm), amplitude_grid, bounds_error=False, fill_value=0.0j + ) + mapped_back = grid_interp(points).reshape(amplitude_map.shape) + denominator = max(float(np.linalg.norm(amplitude_map)), np.finfo(float).tiny) + gridding_relative_l2 = float(np.linalg.norm(mapped_back-amplitude_map)/denominator) + gridding_absolute_l2_by_surface = np.linalg.norm( + mapped_back-amplitude_map, axis=1 + ) + surface_reference_l2 = np.linalg.norm(amplitude_map, axis=1) + reference_floor = np.finfo(float).eps*max( + float(np.max(surface_reference_l2)), 1.0 + ) + gridding_relative_l2_by_surface = [ + None if reference <= reference_floor else float(error/reference) + for error, reference in zip( + gridding_absolute_l2_by_surface, surface_reference_l2, strict=True + ) + ] + metadata = { + "format": "POTATO bmod_n.dat, Fortran sequential, complex amplitude in gauss", + "inputs": { + "chartmap": {"path": str(chartmap), "sha256": _sha256(chartmap)}, + "components": {"path": str(components), "sha256": _sha256(components)}, + }, + "output": {"path": str(output), "sha256": _sha256(output)}, + "component": component, + "signed_toroidal_mode": n_tor, + "source_fourier_convention": "Delta|B|=Re[A_B*exp(i*n*phi_B)]", + "output_fourier_convention": "Delta|B|=Re[A_RZ(R,Z)*exp(i*n*phi_geom)]", + "toroidal_angle_transform": "A_RZ=A_B*exp(i*n*(phi_B-phi_geom))", + "chartmap_zeta_convention": zeta_convention, + "chartmap_zeta_slice": float(zeta[0]), + "toroidal_shift_radians_max_abs": float(np.max(np.abs(np.angle( + np.exp(1j*(zeta[0] - phi_geom)) + )))), + "toroidal_phase_correction_radians_max_abs": float(np.max(np.abs( + np.angle(toroidal_shift_map) + ))), + "radial_coordinate": "s_tor", + "s_map_min": float(s_map[0]), + "s_map_max": float(s_map[-1]), + "spectrum_s_min": float(s_spectrum[0]), + "spectrum_s_max": float(s_spectrum[-1]), + "rad_cm": [float(rad_cm[0]), float(rad_cm[-1]), nrad], + "zet_cm": [float(zet_cm[0]), float(zet_cm[-1]), nzet], + "ntheta_reconstruction": ntheta, + "chartmap_ntheta": int(theta.size), + "margin_fraction": margin_fraction, + "outside_mapped_surface": "zero", + "radial_interpolation": "piecewise linear; no smoothing or fit", + "rz_interpolation": "piecewise linear Delaunay; no smoothing or fit", + "gridding_relative_l2": gridding_relative_l2, + "gridding_relative_l2_by_surface": { + "s_tor": s_map.tolist(), + "absolute_l2_tesla": gridding_absolute_l2_by_surface.tolist(), + "reference_l2_tesla": surface_reference_l2.tolist(), + "relative_l2": gridding_relative_l2_by_surface, + "relative_l2_null_policy": ( + "null when reference_l2_tesla <= eps*max(max_reference,1 T)" + ), + }, + "grid_zero_fraction": float(np.count_nonzero(amplitude_grid == 0.0j)/amplitude_grid.size), + "amplitude_tesla_max": float(np.max(np.abs(amplitude_grid))), + } + metadata_path = output.with_suffix(output.suffix + ".json") + metadata_path.write_text(json.dumps(metadata, indent=2, sort_keys=True) + "\n") + return metadata + + +def parser() -> argparse.ArgumentParser: + result = argparse.ArgumentParser(description=__doc__) + result.add_argument("chartmap", type=Path) + result.add_argument("components", type=Path) + result.add_argument("output", type=Path) + result.add_argument("--component", default="total") + result.add_argument("--n-tor", type=int, required=True, + help="signed native toroidal harmonic used by POTATO") + result.add_argument("--s-max", type=float, default=0.704, + help="largest mapped s_tor surface (default: 0.704)") + result.add_argument("--nrad", type=int, default=801) + result.add_argument("--nzet", type=int, default=801) + result.add_argument("--ntheta", type=int, default=1024) + result.add_argument("--margin-fraction", type=float, default=0.05) + return result + + +def main() -> None: + args = parser().parse_args() + metadata = convert( + args.chartmap, args.components, args.output, + component=args.component, n_tor=args.n_tor, s_max=args.s_max, + nrad=args.nrad, nzet=args.nzet, ntheta=args.ntheta, + margin_fraction=args.margin_fraction, + ) + print(json.dumps(metadata, indent=2, sort_keys=True)) + + +if __name__ == "__main__": + main() diff --git a/POTATO/tools/neo_rt_profiles_to_potato.py b/POTATO/tools/neo_rt_profiles_to_potato.py new file mode 100755 index 00000000..f50bdfe3 --- /dev/null +++ b/POTATO/tools/neo_rt_profiles_to_potato.py @@ -0,0 +1,184 @@ +#!/usr/bin/env python3 +"""Convert NEO-RT s_tor profiles to POTATO polynomials in explicit s_pol.""" + +from __future__ import annotations + +import argparse +import hashlib +import json +from pathlib import Path + +import numpy as np +from numpy.polynomial import Polynomial + +C_CGS = 2.99792458e10 +EV_CGS = 1.602176e-12 +AMU_CGS = 1.660538e-24 + + +def sha256(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as stream: + for chunk in iter(lambda: stream.read(1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def numeric_rows(path: Path) -> list[list[float]]: + rows = [] + for line in path.read_text().splitlines(): + stripped = line.strip() + if not stripped or stripped.startswith(("#", "%")): + continue + try: + rows.append([float(value) for value in stripped.split()]) + except ValueError: + continue + if not rows: + raise ValueError(f"{path} contains no numeric rows") + return rows + + +def fit_descending(x: np.ndarray, y: np.ndarray, degree: int) -> tuple[np.ndarray, float]: + polynomial = Polynomial.fit(x, y, degree, domain=[0.0, 1.0]).convert() + coefficients = np.pad(polynomial.coef, (0, degree + 1-polynomial.coef.size)) + scale = max(float(np.linalg.norm(y)), np.finfo(float).tiny) + return coefficients[::-1], float(np.linalg.norm(polynomial(x)-y)/scale) + + +def convert(profile: Path, plasma: Path, geometry: Path, output: Path, *, + r0_cm: float, psi_span_tm2: float, relation_sign: int, + degree: int = 9) -> dict: + if r0_cm <= 0.0 or psi_span_tm2 <= 0.0: + raise ValueError("r0_cm and psi_span_tm2 must be positive") + if relation_sign not in (-1, 0, 1): + raise ValueError("relation_sign must be -1, 0, or +1") + + rotation = np.asarray(numeric_rows(profile), dtype=float) + plasma_rows = numeric_rows(plasma) + if rotation.shape[1] < 2 or len(plasma_rows) < 2: + raise ValueError("truncated NEO-RT profile input") + species = np.asarray(plasma_rows[0], dtype=float) + if species.size < 5: + raise ValueError("expected the NEO-RT two-ion header with leading grid count") + grid_count = int(round(species[0])) + if grid_count <= 1 or not np.isclose(species[0], grid_count): + raise ValueError(f"invalid NEO-RT plasma grid count {species[0]}") + masses = species[1:3] + charges = species[3:5] + data = np.asarray(plasma_rows[1:], dtype=float) + if data.shape[0] != grid_count: + raise ValueError( + f"NEO-RT plasma header declares {grid_count} rows, found {data.shape[0]}" + ) + if data.shape[1] < 6: + raise ValueError("expected s, two densities, two ion temperatures, and Te") + candidates = [index for index, charge in enumerate(charges) + if charge > 0.0 and np.any(data[:, 1+index] > 0.0)] + if len(candidates) != 1: + raise ValueError(f"cannot select one physical ion from charges {charges}") + ion = candidates[0] + + with np.load(geometry) as mapping: + s_pol_map = np.asarray(mapping["s_sqrt_poloidal"], dtype=float)**2 + s_tor_map = np.asarray(mapping["s_toroidal"], dtype=float) + if np.any(np.diff(s_pol_map) <= 0.0) or np.any(np.diff(s_tor_map) <= 0.0): + raise ValueError("geometry flux map must increase strictly") + + s_pol = np.linspace(0.0, 1.0, 1001) + s_tor = np.interp(s_pol, s_pol_map, s_tor_map) + density = np.interp(s_tor, data[:, 0], data[:, 1+ion]) + temperature_ev = np.interp(s_tor, data[:, 0], data[:, 3+ion]) + mach = np.interp(s_tor, rotation[:, 0], rotation[:, 1]) + vth_cm_s = np.sqrt( + 2.0*temperature_ev*EV_CGS/(masses[ion]*AMU_CGS) + ) + supplied_vth_relative_error = None + if rotation.shape[1] >= 3: + supplied_vth = np.interp(s_tor, rotation[:, 0], rotation[:, 2]) + supplied_vth_relative_error = float( + np.max(np.abs(supplied_vth/vth_cm_s-1.0)) + ) + if supplied_vth_relative_error > 1.0e-5: + raise ValueError( + "profile vth column disagrees with charge-selected plasma " + f"temperature and mass by {supplied_vth_relative_error:.3e}" + ) + omega_e_s = mach*vth_cm_s/r0_cm + + # Phi is statvolt and psi_pol is gauss cm^2. The sign is explicit because + # NEO-RT and POTATO encode their poloidal/toroidal orientation differently. + dphi_ds_pol = relation_sign*omega_e_s*(psi_span_tm2*1.0e8)/C_CGS + potential_statv = np.zeros_like(s_pol) + potential_statv[1:] = np.cumsum( + 0.5*(dphi_ds_pol[1:]+dphi_ds_pol[:-1])*np.diff(s_pol) + ) + + density_coef, density_error = fit_descending(s_pol, density, degree) + temperature_coef, temperature_error = fit_descending(s_pol, temperature_ev, degree) + potential_coef, potential_error = fit_descending(s_pol, potential_statv, degree) + output.parent.mkdir(parents=True, exist_ok=True) + with output.open("w") as stream: + stream.write("% NEO-RT profiles mapped from s_tor to s_pol\n") + stream.write("% density, dummy, ion temperature, potential; descending powers\n") + for values in (density_coef, np.zeros(degree+1), + temperature_coef, potential_coef): + stream.write(" ".join(f"{value:.16e}" for value in values) + "\n") + + metadata = { + "inputs": { + "profile": {"path": str(profile), "sha256": sha256(profile)}, + "plasma": {"path": str(plasma), "sha256": sha256(plasma)}, + "geometry": {"path": str(geometry), "sha256": sha256(geometry)}, + }, + "output": {"path": str(output), "sha256": sha256(output)}, + "input_abscissa": "s_tor", + "output_abscissa": "s_pol=rho_pol^2", + "plasma_grid_count": grid_count, + "selected_ion_index_one_based": ion + 1, + "selected_ion_mass_amu": float(masses[ion]), + "selected_ion_charge": float(charges[ion]), + "thermal_speed_source": ( + "sqrt(2*Ti_eV*eV_cgs/(mass_amu*amu_cgs)) from plasma.in" + ), + "optional_profile_vth_max_relative_error": supplied_vth_relative_error, + "r0_cm": r0_cm, + "psi_pol_span_tm2": psi_span_tm2, + "potential_relation": "dPhi/dpsi_pol = relation_sign*Omega_E/c", + "potential_relation_sign": relation_sign, + "potential_units": "statvolt", + "omega_e_s_range": [float(np.min(omega_e_s)), float(np.max(omega_e_s))], + "potential_statv_edge": float(potential_statv[-1]), + "polynomial_degree": degree, + "relative_l2_fit": { + "density": density_error, + "temperature": temperature_error, + "potential": potential_error, + }, + } + output.with_suffix(output.suffix + ".json").write_text( + json.dumps(metadata, indent=2, sort_keys=True) + "\n" + ) + return metadata + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("profile", type=Path) + parser.add_argument("plasma", type=Path) + parser.add_argument("geometry", type=Path) + parser.add_argument("output", type=Path) + parser.add_argument("--r0-cm", type=float, required=True) + parser.add_argument("--psi-span-tm2", type=float, required=True) + parser.add_argument("--relation-sign", type=int, choices=(-1, 0, 1), required=True) + parser.add_argument("--degree", type=int, default=9) + args = parser.parse_args() + print(json.dumps(convert( + args.profile, args.plasma, args.geometry, args.output, + r0_cm=args.r0_cm, psi_span_tm2=args.psi_span_tm2, + relation_sign=args.relation_sign, degree=args.degree, + ), indent=2, sort_keys=True)) + + +if __name__ == "__main__": + main()