|
| 1 | +#!/usr/bin/env python3 |
| 2 | +"""Write KIN6D's ``nested_input`` namelist for a fixed-boundary VMEC input. |
| 3 | +
|
| 4 | +KIN6D's ``nested_equilibrium`` prescribes iota(s) and p(s) as monomials in the |
| 5 | +normalized toroidal flux s, with a positive toroidal flux. VMEC flips theta |
| 6 | +internally when the boundary runs clockwise in (R, Z) and reverses iota with it; |
| 7 | +the same flip is applied here, so iota keeps its handedness relative to the |
| 8 | +boundary. |
| 9 | +
|
| 10 | +usage: convert_vmec_to_kin6d.py INPUT NAMELIST |
| 11 | +""" |
| 12 | +from __future__ import annotations |
| 13 | + |
| 14 | +import argparse |
| 15 | +from pathlib import Path |
| 16 | + |
| 17 | +import numpy as np |
| 18 | + |
| 19 | +from convert_vmec_to_gvec import _dense_vmec_namelist, _mark_unsupported |
| 20 | + |
| 21 | +# KIN6D's monomial profiles hold at most 33 coefficients |
| 22 | +_MAX_COEFFICIENTS = 33 |
| 23 | + |
| 24 | + |
| 25 | +def _monomials(nml: dict, kind: str, prefix: str) -> np.ndarray: |
| 26 | + kind_type = str(nml.get(f"{kind}_type", "power_series")).strip().lower() |
| 27 | + if kind_type != "power_series": |
| 28 | + raise NotImplementedError(f"KIN6D takes power-series profiles, not {kind.upper()}_TYPE = '{kind_type}'") |
| 29 | + values = np.atleast_1d(np.asarray(nml.get(prefix, [0.0]), dtype=float)) |
| 30 | + values = np.trim_zeros(np.nan_to_num(values), "b") |
| 31 | + if len(values) > _MAX_COEFFICIENTS: |
| 32 | + raise NotImplementedError(f"KIN6D holds at most {_MAX_COEFFICIENTS} {prefix.upper()} coefficients") |
| 33 | + return values if len(values) else np.zeros(1) |
| 34 | + |
| 35 | + |
| 36 | +def _boundary(nml: dict) -> dict[str, dict[tuple[int, int], float]]: |
| 37 | + """Boundary coefficients keyed (n, m) as in VMEC, oriented counterclockwise in (R, Z).""" |
| 38 | + ntor = int(nml.get("ntor", 0)) |
| 39 | + out = {} |
| 40 | + for field in ("rbc", "rbs", "zbc", "zbs"): |
| 41 | + dense = nml.get(field) |
| 42 | + out[field] = {} |
| 43 | + if not isinstance(dense, list): |
| 44 | + continue |
| 45 | + for m, row in enumerate(dense): |
| 46 | + for j, value in enumerate(row or []): |
| 47 | + if value: |
| 48 | + out[field][(j - ntor, m)] = float(value) |
| 49 | + theta = np.linspace(0.0, 2.0 * np.pi, 512, endpoint=False) |
| 50 | + |
| 51 | + def at_zeta0(field, derivative=False): |
| 52 | + total = np.zeros_like(theta) |
| 53 | + for (_, m), c in out[field].items(): |
| 54 | + if field in ("rbc", "zbc"): |
| 55 | + total += -c * m * np.sin(m * theta) if derivative else c * np.cos(m * theta) |
| 56 | + else: |
| 57 | + total += c * m * np.cos(m * theta) if derivative else c * np.sin(m * theta) |
| 58 | + return total |
| 59 | + |
| 60 | + area = np.sum((at_zeta0("rbc") + at_zeta0("rbs")) * (at_zeta0("zbc", True) + at_zeta0("zbs", True))) |
| 61 | + out["flipped"] = area < 0 |
| 62 | + if out["flipped"]: |
| 63 | + # theta -> -theta: (n, m) -> (-n, m) for m > 0; sine coefficients change sign |
| 64 | + for field in ("rbc", "rbs", "zbc", "zbs"): |
| 65 | + sign = -1.0 if field in ("rbs", "zbs") else 1.0 |
| 66 | + out[field] = {(-n if m else n, m): (sign * c if m else c) for (n, m), c in out[field].items()} |
| 67 | + return out |
| 68 | + |
| 69 | + |
| 70 | +def convert(source: Path, destination: Path) -> None: |
| 71 | + try: |
| 72 | + nml = _dense_vmec_namelist(source) |
| 73 | + except ValueError as exc: |
| 74 | + _mark_unsupported(destination, str(exc)) |
| 75 | + if int(nml.get("ncurr", 0)) == 1: |
| 76 | + _mark_unsupported(destination, "KIN6D prescribes iota; current-constrained inputs (NCURR = 1) are not supported") |
| 77 | + for key, neutral in (("bloat", 1.0), ("spres_ped", 1.0)): |
| 78 | + if abs(float(nml.get(key, neutral))) != neutral: |
| 79 | + _mark_unsupported(destination, f"KIN6D does not take {key.upper()} = {nml[key]}") |
| 80 | + try: |
| 81 | + scale = float(nml.get("pres_scale", 1.0)) |
| 82 | + pressure = scale * _monomials(nml, "pmass", "am") if scale else np.zeros(1) |
| 83 | + iota = _monomials(nml, "piota", "ai") |
| 84 | + except NotImplementedError as exc: |
| 85 | + _mark_unsupported(destination, str(exc)) |
| 86 | + if bool(nml.get("lrfp", False)): |
| 87 | + _mark_unsupported(destination, "KIN6D prescribes iota, not the safety factor of LRFP = T") |
| 88 | + boundary = _boundary(nml) |
| 89 | + if boundary["flipped"]: |
| 90 | + iota = -iota |
| 91 | + phiedge = abs(float(nml.get("phiedge", 1.0))) |
| 92 | + mpol = max([m for f in ("rbc", "rbs", "zbc", "zbs") for (_, m) in boundary[f]] + [1]) |
| 93 | + ntor = max([abs(n) for f in ("rbc", "rbs", "zbc", "zbs") for (n, _) in boundary[f]] + [0]) |
| 94 | + lines = [ |
| 95 | + "&nested_input", |
| 96 | + f" nfp = {int(nml.get('nfp', 1))}, mpol = {mpol}, ntor = {ntor}, phiedge = {phiedge!r},", |
| 97 | + ] |
| 98 | + for field in ("rbc", "rbs", "zbc", "zbs"): |
| 99 | + lines += [f" {field}({n},{m}) = {float(c)!r}," for (n, m), c in sorted(boundary[field].items())] |
| 100 | + lines += [f" iota({k}) = {float(c)!r}," for k, c in enumerate(iota) if c] |
| 101 | + lines += [f" pressure({k}) = {float(c)!r}," for k, c in enumerate(pressure) if c] |
| 102 | + # VMEC's boundary angle is not a field-line angle, so the boundary field angle is free |
| 103 | + lines += [" boundary_nodes = 16, degree = 3, angle_mpol = 4, tolerance = 1e-8, max_iterations = 400,", "/"] |
| 104 | + destination.write_text("\n".join(lines) + "\n") |
| 105 | + |
| 106 | + |
| 107 | +def main() -> int: |
| 108 | + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) |
| 109 | + parser.add_argument("input", type=Path) |
| 110 | + parser.add_argument("namelist", type=Path) |
| 111 | + args = parser.parse_args() |
| 112 | + convert(args.input, args.namelist) |
| 113 | + return 0 |
| 114 | + |
| 115 | + |
| 116 | +if __name__ == "__main__": |
| 117 | + raise SystemExit(main()) |
0 commit comments