Grid Transformations and Profile Flow

PORTALS optimizes plasma profiles by working in normalized-gradient space (\(a/L_T\), \(a/L_n\)) on a coarse radial grid of ~5 control points. Transport codes (TGLF, CGYRO, NEO) need full profiles on a fine grid. This page documents every grid transformation, interpolation, and derivative computation that happens along the way — and the subtle numerical issues that arise from them.

Grid types

MITIM uses several radial grids throughout the workflow. Understanding which grid is active at each stage is essential for debugging discrepancies.

Grid	Typical size	Description
Original	~30–120 pts	The user-provided `input.gacode` file. Profiles live on a \(\rho\) (square root of normalized toroidal flux) coordinate.
Enhanced fine	~100+ pts	After `improve_resolution_profiles()`. Uniform fill between control points, plus extra points at \(\pm 10^{-3}\) and \(\pm 2 \times 10^{-3}\) around each control point.
Coarse powerstate	~5 pts	The PORTALS prediction radii (e.g., \(\rho = 0.2, 0.4, 0.55, 0.7, 0.85\)). These are the free-parameter locations for the optimization. A zero at the axis is prepended internally.
Fine targets (optional)	variable	A higher-resolution grid used only inside `calculateTargets()` for more accurate volume integration of power sources (fusion, radiation, exchange).
Scalar extraction	1 pt per surface	What TGLF/CGYRO/NEO receives: local values at each flux surface, obtained by cubic-spline interpolation from the enhanced fine grid.

Workflow overview

The following diagram traces the input.gacode data from the user’s file all the way to what TGLF/CGYRO receives. Grid changes and derivative computations are highlighted.

USER's input.gacode  (~30-120 pts, in rho)
|
+---> input.gacode_original  (untouched backup)
|
v
profiles.correct(recalculate_ptot=True)
|
===== powerstate.__init__() ==========================================
|
|  improve_resolution_profiles()
|     ~40 pts --> ~100+ pts
|     Cubic-spline interpolation of all quantities to new grid
|     Adds extra points at +/-1e-3, +/-2e-3 of each control point
|     ** Modifies profiles IN-PLACE **
|
|  gacode_to_powerstate()
|     ~100+ pts --> ~5 coarse pts
|     Cubic-spline interpolation of profiles to powerstate grid
|     DERIVATIVE: parameterizers.piecewise_linear()
|       Computes a/LT via derivation_into_Lx() on the ~100+ pt grid
|       Samples at coarse control points
|       Creates profile_constructor functions (fine/coarse/middle)
|
===== PORTALSinit: First roundtrip ===================================
|
|  update_var(): a/LT (coarse) --> integration --> T (coarse)
|     ** Profile now differs from original (roundtrip error) **
|
|  from_powerstate(): a/LT (coarse) --> profile_constructors_fine
|     --> T on ~100+ pt grid
|     Writes: input.gacode (parameterized), input.gacode_modified
|
===== Each PORTALS evaluation ========================================
|
|  1. update_var(): insert new a/LT, integrate to T on coarse grid
|  2. calculateProfileFunctions(): gyro-Bohm norms, collisionality, etc.
|  3. calculateTargets(): power balance (optionally on fine targets grid)
|  4. _produce_profiles():
|       profile_constructors_fine: a/LT (coarse) --> T on ~100+ pt grid
|       derive_quantities(): recompute aLTe, aLne, etc. on fine grid
|       Write: input.gacode (~100+ pts)
|
|  5. Transport code input preparation:
|       to_tglf() / to_cgyro() / to_neo():
|         Cubic-spline interpolation in r/a space to exact flux surface
|         --> RLTS, RLNS, DLNTDR, DLNNDR (scalar values per surface)
|
======================================================================

Key operations

Derivative computation

The normalized gradient scale length is defined as:

\[\frac{a}{L_X} = -\frac{a}{X} \frac{dX}{dr}\]

MITIM computes this via derivation_into_Lx(r, p), which uses a 3-point Lagrange interpolating polynomial derivative of \(-\ln(p)\). This is the same scheme as GACODE’s bound_deriv in expro_util.f90, ensuring that MITIM and GACODE (TGLF, CGYRO, NEO) compute identical gradient values from the same input.gacode profiles.

The Lagrange polynomial derivative evaluates at each point \(r_i\) using its two nearest neighbors \((r_{i-1}, r_i, r_{i+1})\) for interior points, and a one-sided 3-point stencil at the boundaries.

An alternative central-difference variant is available as derivation_into_Lx_central(). It computes:

\[z_i = -\frac{\ln p_{i+1} - \ln p_{i-1}}{r_{i+1} - r_{i-1}} \quad \text{(interior points)}\]

This is structurally more consistent with integration_Lx’s trapezoidal log-space formula and reduces the fine-grid roundtrip error by ~30%. However, it uses a different stencil than GACODE’s bound_deriv, which introduces small discrepancies (~0.01%) when transport codes recompute gradients from the same input.gacode. For this reason, the Lagrange method is the default.

Integration (gradient to profile)

integration_Lx(x, z, f_bound) reconstructs a profile from its normalized gradient, adapted from TGYRO’s tgyro_profile_functions.f90:

\[\frac{T_i}{T_{i+1}} = \exp\!\Big[\tfrac{1}{2}(z_i + z_{i+1})(r_{i+1} - r_i)\Big]\]

where \(z = a/L_T\) and \(f_\text{bound}\) is the boundary condition at the last grid point. This trapezoidal rule in log-space is exact for piecewise-linear \(z(r)\), which is the representation used after coarse-grid sampling and interpolation.

Interpolation coordinate: \(r/a\) vs \(\rho\)

When extracting local quantities at a specific flux surface (e.g., for TGLF input), MITIM interpolates profiles using a cubic spline. The choice of which radial coordinate to use for the spline matters:

GACODE’s approach: expro_locsim in GACODE uses natural cubic splines in \(r/a\) (rmin) space (see cub_spline1 in expro_locsim.f90).
MITIM’s approach: The to_tglf(), to_cgyro(), to_neo(), and to_gx() functions now also interpolate in \(r/a\) space, matching GACODE’s convention. Input \(\rho\) values are first converted to \(r/a\) before interpolation.

Why this matters: The mapping \(\rho \to r/a\) is nonlinear (it depends on the magnetic equilibrium geometry). For a point that falls between grid nodes, a cubic spline in \(\rho\) and a cubic spline in \(r/a\) will give slightly different interpolated values. Near the plasma edge where the mapping is most nonlinear, the discrepancy can reach ~0.1–0.2%.

By using the same interpolation coordinate as GACODE, MITIM ensures that the gradient values written to input.tglf / input.cgyro match what GACODE’s internal routines would compute from input.gacode via PROFILE_MODEL=2.

Profile reconstruction and the trailing edge

When PORTALS reconstructs a full profile from coarse gradients (via profile_constructors_fine), the following steps are applied:

Piecewise-linear gradient interpolation: The gradient values at the coarse control points are linearly interpolated to the full fine grid.
smoothAroundCoarsing (controlled by transport.flatten_gradients_at_control_points in the PORTALS namelist): The gradient between control points is piecewise-linear, which creates a slope change (kink in gradient space) at each control point. If a code recomputes gradients from the reconstructed profile using a higher-order stencil (e.g., TGYRO’s 3-point Lagrange), it will “see” the slope change through neighboring points and return a value that differs from the intended control point gradient.

When enabled (default: true), the gradient at \(\pm 1\) neighboring fine-grid points around each control point is flattened to match the control point value. This creates a tiny plateau so that any derivative stencil returns the correct value regardless of its order. The extra points at \(\pm 10^{-3}\) added by improve_resolution_profiles() exist specifically to provide these close neighbors; they are skipped when smoothing is disabled.

When disabled (false), the piecewise-linear gradient is left as-is, the \(\pm 10^{-3}\) extra grid points are not added, and the fine grid is more uniform. This may be preferable when the PORTALS evaluation does not rely on codes that recompute gradients from profiles.
Integration with direct BC: The gradient profile is integrated via integration_Lx with the boundary condition applied directly at the last control point. No extra points or rescaling are needed.
Hermite trailing edge blending: Beyond the last control point, the reconstructed profile must transition back to the original (unmodified) profile. This is done with a C1-continuous Hermite smoothstep blend (\(w(t) = 3t^2 - 2t^3\)) over a small region (a few grid points). The reconstructed profile is linearly extrapolated from the junction using its last gradient, then blended with the original profile. This eliminates the derivative discontinuity (kink) that previously caused ~0.1% gradient discrepancies when transport codes recomputed gradients from the reconstructed profiles.

The roundtrip problem

Converting a profile to gradients and back does not reproduce the original profile:

\[T \xrightarrow{\text{derivative}} \frac{a}{L_T} \xrightarrow{\text{coarsen to 5 pts}} \xrightarrow{\text{interpolate}} \xrightarrow{\text{integrate}} T' \neq T\]

The error has two sources:

Derivative–integration mismatch: No point-value derivative is the exact mathematical inverse of integration_Lx. The integration uses interval averages \(\tfrac{1}{2}(z_i + z_{i+1})\), which cannot be uniquely inverted to point values \(z_i\).
Coarsening: Sampling the gradient at ~5 points and piecewise-linearly interpolating loses the original gradient shape between control points.

Current policy: PORTALS allows exactly one roundtrip pass at initialization (in PORTALSinit). This “bakes in” the roundtrip error but ensures the starting point is self-consistent with the parameterization. Subsequent evaluations operate purely in gradient space and do not introduce additional roundtrip degradation — the profile is always reconstructed from the current gradients, never re-differentiated.

GACODE consistency

GACODE’s internal derivative (bound_deriv in expro_util.f90) and MITIM’s derivation_into_Lx are both 2nd-order accurate but use different stencils. When GACODE (e.g., CGYRO with PROFILE_MODEL=2) reads the same input.gacode file and recomputes gradients, small differences can arise:

Source	Typical magnitude	Cause
Derivative stencil	< 0.01%	Both methods are 2nd-order on the same grid; differences are \(O(h^2)\) where \(h\) is the fine grid spacing.
Interpolation coordinate	~0.1–0.2% at edge	Fixed by interpolating in \(r/a\) instead of \(\rho\) (now implemented).
Text format truncation	~0.001%	`input.gacode` text format has ~7 significant digits. Derivatives of truncated values differ from derivatives of full-precision values.
Trailing edge blending	negligible	The Hermite smoothstep blend ensures C1 continuity at the trailing edge junction.

Known caveats and future work

Grid overdimensioning: The \(\pm 10^{-3}\) extra points added by improve_resolution_profiles() exist to support smoothAroundCoarsing. They are automatically skipped when transport.flatten_gradients_at_control_points: false is set in the PORTALS namelist.
Multi-pass idempotency: Achieving perfect stability (pass 2 = pass 1) through the derivative function alone is not possible because integration_Lx uses interval averages that cannot be uniquely inverted to point values. True idempotency would require changes at the parameterization level (e.g., self-consistent initialization via iterative gradient correction).