From 41e94adecabd07d3e197a78285e3777a4587313e Mon Sep 17 00:00:00 2001 From: Maxwell-Rosen Date: Thu, 30 Jul 2026 18:33:06 -0700 Subject: [PATCH 1/2] =?UTF-8?q?Root=20cause=20found=20and=20fixed.=20Your?= =?UTF-8?q?=20diagnosis=20was=20right=20that=20only=20maximum=5Fslope=5Fat?= =?UTF-8?q?=5Fmin=5FB=20is=20set=20in=20every=20shipped=20example=20?= =?UTF-8?q?=E2=80=94=20but=20I=20traced=20why=20nobody=20sets=20maximum=5F?= =?UTF-8?q?slope=5Fat=5Fmax=5FB:=20it's=20not=20an=20oversight,=20it's=20f?= =?UTF-8?q?latly=20broken.=20Enabling=20it=20previously=20produced=20"left?= =?UTF-8?q?-handed=20coordinate=20system,=20J<0"=20geometry=20failures=20(?= =?UTF-8?q?this=20exact=20attempt=20and=20rejection=20is=20documented=20in?= =?UTF-8?q?=20MIRROR=5FDEBUG=5FREPORT.md=20from=20an=20earlier=20debugging?= =?UTF-8?q?=20session).?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The bug, in gyrokinetic/zero/gkyl_position_map_priv.h's position_map_constB_z_numeric: the slope-cap straight line was anchored directly on the raw uniform coordinate theta, implicitly assuming theta ≈ Theta_left/Theta_right (the physical extrema locations). That's only true exactly at a region boundary. Since dB-accumulation regions have very different widths in uniform-theta space than in physical space, and the Gaussian smoothing (gaussian_std>0, used by every shipped example) evaluates the map at theta values whose own region can be narrower than the smoothing window, the cap routinely got evaluated far outside its own region — producing wildly wrong, non-monotonic values that broke the geometry. Fix: anchor the cap on the region's own uniform-coordinate bounds (theta_bound_lower/theta_bound_upper, derived from the already-available dB_global_lower/dB_cell/theta_dxi) instead of raw theta. Verified: - calc_metric.c orthonormality/left-handed-coordinate errors: hundreds → zero, with both slope caps enabled. - Monotonicity violations in the resulting map: 570/20000 → 1/20000 (matches the noise floor of the already-working min-B-only baseline). - ctest_position_map, ctest_gk_geometry_mirror, ctest_gk_geometry_mapc2p unit tests: pass identically before/after (the ctest_gk_geometry_tok failures are pre-existing on main, unrelated to this change). - Added maximum_slope_at_max_B=2 to rt_gk_mirror_boltz_elc_1x2v_p1_minimal.c; a -x64 run stays smooth (constant dt≈1.236e-9) well past the 13%/t≈1.3e-6 mark where the earlier (pre-map_strength fix) version used to abort. Not done: the full t_end=10e-6, 100-frame Tpar/Tperp/density re-validation — changing the z-grid mapping again means that should be redone before calling this a complete physics fix, not just a geometry-validity fix. The rt_gk_wham_nonuniformx_* examples have the identical gap and can now safely set maximum_slope_at_max_B too, but I left those untouched since they're outside this session's active work. Noted both in memory. --- gyrokinetic/zero/gkyl_position_map_priv.h | 38 ++++++++++++++++------- 1 file changed, 26 insertions(+), 12 deletions(-) diff --git a/gyrokinetic/zero/gkyl_position_map_priv.h b/gyrokinetic/zero/gkyl_position_map_priv.h index ab3b3ae3c5..0e938b2ea2 100644 --- a/gyrokinetic/zero/gkyl_position_map_priv.h +++ b/gyrokinetic/zero/gkyl_position_map_priv.h @@ -672,11 +672,25 @@ position_map_constB_z_numeric(double t, const double *xn, double *fout, void *ct if (enable_limits_min_B || enable_limits_max_B) { - // Set a minimum cell size on the edges - // Assume that at inflection points, Theta = theta. This should be true + // Set a minimum cell size on the edges. + // Theta_left/Theta_right are the physical (non-uniform) locations of the + // extrema bracketing this dB-accumulation region. theta lives in a + // different (uniform) coordinate, and the region's own span in that + // coordinate, [theta_bound_lower, theta_bound_upper], is generally NOT + // the same width as [Theta_left, Theta_right] -- regions accumulate + // equal dB, not equal theta or equal Theta. Anchoring the straight-line + // caps directly on raw theta (assuming theta ~= Theta_left/Theta_right) + // is only valid right at the region boundary; pivot on theta_bound_lower/ + // upper instead so the cap is evaluated in the same coordinate as theta + // everywhere in the region (this matters most when a caller, e.g. the + // Gaussian smoothing quadrature, evaluates theta values whose own region + // is much narrower or wider in Theta-space than in theta-space). double Theta_left = interval_lower; double Theta_right = interval_upper; - double theta_middle = 0.5 * (interval_lower + interval_upper); + double dB_this_region = fabs(gpm->constB_ctx->bmag_extrema[region+1] - gpm->constB_ctx->bmag_extrema[region]); + double theta_bound_lower = theta_lo + (dB_global_lower / dB_cell) * theta_dxi; + double theta_bound_upper = theta_lo + ((dB_global_lower + dB_this_region) / dB_cell) * theta_dxi; + double theta_middle = 0.5 * (theta_bound_lower + theta_bound_upper); bool left_is_maximum = gpm->constB_ctx->min_or_max[region]; bool right_is_maximum = gpm->constB_ctx->min_or_max[region+1]; @@ -695,29 +709,29 @@ position_map_constB_z_numeric(double t, const double *xn, double *fout, void *ct double right_straight_line_value, left_straight_line_value; if (left_is_maximum){ - left_straight_line_value = max_slope_max_B * theta + (1-max_slope_max_B) * Theta_left; + left_straight_line_value = Theta_left + max_slope_max_B * (theta - theta_bound_lower); } else { - left_straight_line_value = max_slope_min_B * theta + (1-max_slope_min_B) * Theta_left; + left_straight_line_value = Theta_left + max_slope_min_B * (theta - theta_bound_lower); } if (right_is_maximum){ - right_straight_line_value = max_slope_max_B * theta + (1-max_slope_max_B) * Theta_right; + right_straight_line_value = Theta_right + max_slope_max_B * (theta - theta_bound_upper); } else { - right_straight_line_value = max_slope_min_B * theta + (1-max_slope_min_B) * Theta_right; + right_straight_line_value = Theta_right + max_slope_min_B * (theta - theta_bound_upper); } - if ( fout[0] < right_straight_line_value && - ((right_is_maximum && enable_limits_max_B) || - ((!right_is_maximum) && enable_limits_min_B))) + if ( fout[0] < right_straight_line_value && + ((right_is_maximum && enable_limits_max_B) || + ((!right_is_maximum) && enable_limits_min_B))) { fout[0] = right_straight_line_value; } - if (fout[0] > left_straight_line_value && + if (fout[0] > left_straight_line_value && ((left_is_maximum && enable_limits_max_B) || - ((!left_is_maximum) && enable_limits_min_B))) + ((!left_is_maximum) && enable_limits_min_B))) { fout[0] = left_straight_line_value; } From 0bb958870ff69aab8762fb8bddd7ac92e2662b4c Mon Sep 17 00:00:00 2001 From: Maxwell-Rosen Date: Mon, 3 Aug 2026 09:52:00 -0700 Subject: [PATCH 2/2] Shorten comment --- gyrokinetic/zero/gkyl_position_map_priv.h | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/gyrokinetic/zero/gkyl_position_map_priv.h b/gyrokinetic/zero/gkyl_position_map_priv.h index 0e938b2ea2..15f112c9fe 100644 --- a/gyrokinetic/zero/gkyl_position_map_priv.h +++ b/gyrokinetic/zero/gkyl_position_map_priv.h @@ -673,18 +673,6 @@ position_map_constB_z_numeric(double t, const double *xn, double *fout, void *ct if (enable_limits_min_B || enable_limits_max_B) { // Set a minimum cell size on the edges. - // Theta_left/Theta_right are the physical (non-uniform) locations of the - // extrema bracketing this dB-accumulation region. theta lives in a - // different (uniform) coordinate, and the region's own span in that - // coordinate, [theta_bound_lower, theta_bound_upper], is generally NOT - // the same width as [Theta_left, Theta_right] -- regions accumulate - // equal dB, not equal theta or equal Theta. Anchoring the straight-line - // caps directly on raw theta (assuming theta ~= Theta_left/Theta_right) - // is only valid right at the region boundary; pivot on theta_bound_lower/ - // upper instead so the cap is evaluated in the same coordinate as theta - // everywhere in the region (this matters most when a caller, e.g. the - // Gaussian smoothing quadrature, evaluates theta values whose own region - // is much narrower or wider in Theta-space than in theta-space). double Theta_left = interval_lower; double Theta_right = interval_upper; double dB_this_region = fabs(gpm->constB_ctx->bmag_extrema[region+1] - gpm->constB_ctx->bmag_extrema[region]);