Take the collision profile radius from a tabulated label at the physical position (#613)

Collisional runs on geometric REFCOORDS chartmaps (`rho_convention =
unknown`, for example map2disc charts that reach the wall) cannot use
radially varying `&profiles` today, because the chart radius carries no
flux meaning and `init_profile_radius` stops
(src/profile_radius.f90:114). The profiles are, however, often known as
a function of a flux label from a separate equilibrium or field-line
model.

This adds the `&config` key `profile_label_file`. It names a netCDF
table `label(phi, Z, R)` over one field period (uniform `R`, `Z` in m or
cm, `phi` at `(i-1)*2*pi/(nfp*nphi)`, global attribute `nfp`). With
`integ_coords = 5`, `profile_s` evaluates the reference coordinates at
the marker state to obtain `(R, phi, Z)` and interpolates the label
trilinearly, periodic in phi. A position outside the table stops the run
instead of extrapolating. Without the file, the existing error and all
other radius mappings are unchanged. The format is documented in
`src/profile_label.f90` and `docs/model/collisions.rst`, and the input
reference is regenerated.

Tests (synthetic data only):
- `test_profile_radius_label` builds an analytic label (bilinear in R, Z
plus a periodic triangle wave in phi with kinks on grid nodes, which
trilinear interpolation reproduces exactly) on the existing test
chartmap. It compares `profile_s` and the collision coefficients with
the closed form at the independently computed physical position to
1e-12, and it also checks the metre-to-cm conversion and phi
periodicity.
- `test_profile_radius_label_outside` checks that a table which does not
cover the marker stops the run.

Native verification (sCluster, Release, libneo 924ba77, CMake/CTest),
head bcfec2c:
- job 3356096 (node20): `ctest -R
'profile_radius|profiles|input_reference|input_validation|config_reference|collision'`
21/21 pass. `ctest -L unit` 76/77, where the only failure was a 120 s
timeout of `test_splined_kernel_bench` on a loaded node. Rerun alone in
job 3356160, it passes in 42 s, and this PR does not touch that code.
- The same commit on top of #608-#611 (job 3356095, node10):
`golden_record_collision` passes and 37/37 focused native tests pass.

Chris&AI

Author: krystophny

docs/input-reference.json

@@ -4,7 +4,7 @@
   "namelists": {
     "config": {
       "source": "src/params.f90",
-      "source_line": 240,
+      "source_line": 243,
       "variables": [
         {
           "name": "notrace_passing",
@@ -1312,6 +1312,19 @@
           "source": "src/params.f90",
           "source_line": 198
         },
+        {
+          "name": "profile_label_file",
+          "type": "character",
+          "kind": "",
+          "length": "1000",
+          "shape": [],
+          "extent": [],
+          "default": "''",
+          "default_value": "''",
+          "documentation": "Optional netCDF profile label s(R, phi, Z) for the collision profiles with integ_coords = REFCOORDS (format: src/profile_label.f90).",
+          "source": "src/params.f90",
+          "source_line": 201
+        },
         {
           "name": "wall_input",
           "type": "character",
@@ -1323,7 +1336,7 @@
           "default_value": "''",
           "documentation": "",
           "source": "src/params.f90",
-          "source_line": 199
+          "source_line": 202
         },
         {
           "name": "wall_units",
@@ -1336,7 +1349,7 @@
           "default_value": "'m'",
           "documentation": "",
           "source": "src/params.f90",
-          "source_line": 200,
+          "source_line": 203,
           "validation": {
             "choices": [
               "m",
@@ -1355,7 +1368,7 @@
           "default_value": "'auto'",
           "documentation": "",
           "source": "src/params.f90",
-          "source_line": 204,
+          "source_line": 207,
           "validation": {
             "choices": [
               "auto",
@@ -1376,7 +1389,7 @@
           "default_value": "-1000",
           "documentation": "Sentinel: -1000 means user did not set it",
           "source": "src/params.f90",
-          "source_line": 206
+          "source_line": 209
         },
         {
           "name": "crossing_level",
@@ -1389,7 +1402,7 @@
           "default_value": "3",
           "documentation": "SPECTRE interface treatment: only full orbit (3) is supported.",
           "source": "src/params.f90",
-          "source_line": 208,
+          "source_line": 211,
           "validation": {
             "choices": [
               3
@@ -1407,7 +1420,7 @@
           "default_value": "1.5_dp",
           "documentation": "Full orbit near SPECTRE interfaces: entry distance in radial gyroradii and the relative tolerance of RK45 on the sharp per-volume field. 1.5 keeps the switch's p_zeta invariant as at 3 (entry 1 breaks it) at about half the FO time (concept/spectre_loss_study/results/fo-residence).",
           "source": "src/params.f90",
-          "source_line": 213,
+          "source_line": 216,
           "validation": {
             "positive": true
           }
@@ -1423,7 +1436,7 @@
           "default_value": "1.0e-9_dp",
           "documentation": "",
           "source": "src/params.f90",
-          "source_line": 214,
+          "source_line": 217,
           "validation": {
             "positive": true
           }
@@ -1443,7 +1456,7 @@
           "default_value": "(/32, 64, 32/)",
           "documentation": "Per-volume FO spline points: rho, theta, zeta per field period.",
           "source": "src/params.f90",
-          "source_line": 216
+          "source_line": 219
         },
         {
           "name": "spectre_sbeg_is_toroidal_flux",
@@ -1456,7 +1469,7 @@
           "default_value": ".false.",
           "documentation": "",
           "source": "src/params.f90",
-          "source_line": 217
+          "source_line": 220
         },
         {
           "name": "spectre_ncon_r",
@@ -1469,7 +1482,7 @@
           "default_value": "48",
           "documentation": "SPECTRE per-volume Meiss construction grid (n_r, n_th, n_phi). Each volume allocates rank-3 arrays plus batch splines of this size, so the total peak memory scales as Mvol*n_r*n_th*n_phi. spectre_ncon_phi = -1 means automatic: 32 for fields with toroidal harmonics and a minimal phi grid for axisymmetric fields (all n = 0), the dominant memory saver for tokamak cases. A positive value forces that phi count verbatim, bypassing the axisymmetric clamp.",
           "source": "src/params.f90",
-          "source_line": 226,
+          "source_line": 229,
           "validation": {
             "min": 1
           }
@@ -1485,7 +1498,7 @@
           "default_value": "48",
           "documentation": "",
           "source": "src/params.f90",
-          "source_line": 227,
+          "source_line": 230,
           "validation": {
             "min": 1
           }
@@ -1501,7 +1514,7 @@
           "default_value": "-1",
           "documentation": "",
           "source": "src/params.f90",
-          "source_line": 228
+          "source_line": 231
         },
         {
           "name": "spectre_ncon_order",
@@ -1514,7 +1527,7 @@
           "default_value": "5",
           "documentation": "",
           "source": "src/params.f90",
-          "source_line": 229,
+          "source_line": 232,
           "validation": {
             "min": 1
           }
@@ -1530,7 +1543,7 @@
           "default_value": "1000000",
           "documentation": "",
           "source": "src/params.f90",
-          "source_line": 230,
+          "source_line": 233,
           "validation": {
             "min": 1
           }
@@ -1546,7 +1559,7 @@
           "default_value": "1.0e-2_dp",
           "documentation": "",
           "source": "src/params.f90",
-          "source_line": 231,
+          "source_line": 234,
           "validation": {
             "positive": true
           }
@@ -1646,7 +1659,7 @@
           "default_value": ".true.",
           "documentation": "Removed &config keys. They are still read so that old input files run; a value other than the historical default prints a warning.",
           "source": "src/params.f90",
-          "source_line": 235
+          "source_line": 238
         },
         {
           "name": "loopskip",
@@ -1659,7 +1672,7 @@
           "default_value": "0",
           "documentation": "",
           "source": "src/params.f90",
-          "source_line": 236
+          "source_line": 239
         },
         {
           "name": "bmod_ref",
@@ -1672,7 +1685,7 @@
           "default_value": "5.0d4",
           "documentation": "",
           "source": "src/params.f90",
-          "source_line": 237
+          "source_line": 240
         },
         {
           "name": "macrostep_log_ratio",
@@ -1685,7 +1698,7 @@
           "default_value": "1.0d3",
           "documentation": "",
           "source": "src/params.f90",
-          "source_line": 238
+          "source_line": 241
         }
       ]
     },

docs/input-reference.rst

@@ -25,7 +25,7 @@ and ``DOC/config.md`` for additional usage, units and option interactions.
 &config
 -------
 
-Source: ``src/params.f90:240``.
+Source: ``src/params.f90:243``.
 
 .. list-table::
    :header-rows: 1
@@ -565,6 +565,12 @@ Source: ``src/params.f90:240``.
      - ``scalar``
      - ``''``
      -
+   * - ``profile_label_file``
+     - ``character(len=1000)``
+     - ``default``
+     - ``scalar``
+     - ``''``
+     - Optional netCDF profile label s(R, phi, Z) for the collision profiles with integ_coords = REFCOORDS (format: src/profile_label.f90).
    * - ``wall_input``
      - ``character(len=1000)``
      - ``default``

docs/model/collisions.rst

@@ -98,7 +98,18 @@ radius; for 6 the enclosed toroidal flux of the SPECTRE chart
 radius gives :math:`s = \rho^2`, and ``vmec_extended``/``vmec_to_wall``
 charts give :math:`s = (\rho/\rho_\mathrm{lcfs})^2`. Radii without
 toroidal-flux meaning (geometric or poloidal chartmaps) are accepted only
-for radially uniform ``&profiles``.
+for radially uniform ``&profiles``, unless ``&config profile_label_file``
+supplies the profile radius as a table :math:`s(R,\varphi,Z)`. With
+``integ_coords = 5`` the label is then interpolated trilinearly at the
+marker's physical position, which the reference coordinates give, and
+periodically in :math:`\varphi` modulo :math:`2\pi/N_\mathrm{fp}`. The
+netCDF file has dimensions ``R``, ``Z`` and ``phi``, uniform coordinate
+variables ``R(R)`` and ``Z(Z)`` with ``units`` ``"m"`` or ``"cm"``,
+``phi(phi)`` at :math:`(i-1)\,2\pi/(N_\mathrm{fp} n_\varphi)`, the
+variable ``label(phi, Z, R)`` and the integer global attribute ``nfp``.
+The table has to cover the whole tracing domain, and a position outside
+it stops the run. How the label is constructed, for example a flux label
+of a separate equilibrium, is the user's model choice.
 
 Discrete step
 -------------

examples/simple_full.in

@@ -66,6 +66,8 @@ field_input = ''                 ! field file: VMEC wout (.nc), Boozer chartmap
 !                                  SPECTRE (.h5) or coils file; '' = use netcdffile
 coord_input = ''                 ! reference coordinate file: VMEC wout or chartmap (.nc);
 !                                  '' = netcdffile, then field_input
+profile_label_file = ''          ! netCDF label s(R,phi,Z) as &profiles radius for
+!                                  integ_coords = 5 (src/profile_label.f90); '' = chart radius
 netcdffile = ''                  ! VMEC wout file; deprecated fallback for field_input and
 !                                  coord_input; '' = 'wout.nc' where a VMEC file is required
 radial_electric_potential_slope = 0d0 ! Boozer integration: Phi(s) = slope*s added to H; 0 = off

src/CMakeLists.txt

@@ -55,6 +55,7 @@
     parse_ants.f90
     profiles.f90
     collis_alphas.f90
+    profile_label.f90
     profile_radius.f90
     netcdf_orbit_output.f90
     netcdf_results_output.f90

src/params.f90

@@ -196,6 +196,9 @@ module params
 
     character(1000) :: field_input = ''
     character(1000) :: coord_input = ''
+    !> Optional netCDF profile label s(R, phi, Z) for the collision profiles
+    !> with integ_coords = REFCOORDS (format: src/profile_label.f90).
+    character(1000) :: profile_label_file = ''
     character(1000) :: wall_input = ''
     character(16) :: wall_units = 'm'
     !> Radial threshold for querying an external STL wall. A negative value
@@ -261,7 +264,7 @@ module params
         am1, am2, Z1, Z2, &
         densi1, densi2, tempi1, tempi2, tempe, &
         batch_size, ran_seed, particle_index_offset, reuse_batch, &
-        field_input, coord_input, &
+        field_input, coord_input, profile_label_file, &
         wall_input, wall_units, chart_boundary_kind, integ_coords, crossing_level, &
         spectre_fo_entry_larmor, spectre_fo_tolerance, spectre_fo_grid, &
         spectre_sbeg_is_toroidal_flux, spectre_ncon_r, spectre_ncon_th, &