Add resonance flag diagnostics

ldsfl/flowfield.py

@@ -507,9 +507,20 @@ def _mode_SL1(
 
 
 def _compute_flag(lamb2: np.ndarray) -> int:
-    idx = 5 if lamb2.size >= 6 else max(0, lamb2.size - 1)
-    return 1 if np.real(lamb2[idx]) < 0 else -1
+    """Resonance indicator: +1 sub-resonant, -1 super-resonant.
 
+    The sign of ``Re(lambda2)`` for the fundamental lateral mode distinguishes
+    the two regimes: it is negative below the resonant aspect ratio and
+    positive above it. This transition flips bend migration from
+    downstream- to upstream-dominated in the Zolezzi-Seminara model family.
+
+    The previous implementation indexed mode 6 (``idx = 5``) rather than mode
+    1. Mode 6 remains strongly negative across realistic aspect ratios, so the
+    reported flag was effectively constant. Mode 1 crosses zero at ``beta_R``.
+    """
+    if lamb2.size == 0:
+        return 1
+    return 1 if np.real(lamb2[0]) < 0 else -1
 
 def _run_modes(
     *,

ldsfl/flowfield_periodic.py

@@ -281,9 +281,20 @@ def _precompute_modes(
 
 
 def _compute_flag(lamb2: np.ndarray) -> int:
-    idx = 5 if lamb2.size >= 6 else max(0, lamb2.size - 1)
-    return 1 if np.real(lamb2[idx]) < 0 else -1
+    """Resonance indicator: +1 sub-resonant, -1 super-resonant.
 
+    The sign of ``Re(lambda2)`` for the fundamental lateral mode distinguishes
+    the two regimes: it is negative below the resonant aspect ratio and
+    positive above it. This transition flips bend migration from
+    downstream- to upstream-dominated in the Zolezzi-Seminara model family.
+
+    The previous implementation indexed mode 6 (``idx = 5``) rather than mode
+    1. Mode 6 remains strongly negative across realistic aspect ratios, so the
+    reported flag was effectively constant. Mode 1 crosses zero at ``beta_R``.
+    """
+    if lamb2.size == 0:
+        return 1
+    return 1 if np.real(lamb2[0]) < 0 else -1
 
 # --------------------------
 # Periodic BC core (dUZSBC1)

ldsfl/main.py

@@ -21,6 +21,7 @@
 from .outputs import ensure_dirs, plot_it, save_sinuosity_history, save_variables, save_xystcu
 from .profile import preprof_3
 from .resistance import resistance_function_flagbed
+from .resonance import resonance_report
 from .stability import sinuosity_equivalence_stability
 
 
@@ -880,6 +881,7 @@ def run_case(
         "output_length_scale": output_length_scale,
         "output_velocity_scale": output_velocity_scale,
         "sinuo_final": float(sinuo_hist[-1]),
+        "resonance": resonance_report(beta, theta0, ds, rpic_0, flagbed, Mdat),
         "sinuosity_stability": stability_info,
     }
 

ldsfl/resonance.py

@@ -0,0 +1,139 @@
+"""Resonance diagnostics for the linearised bend theory.
+
+The reduced model has a resonant aspect ratio ``beta_R`` at which the
+fundamental mode's streamwise decay rate changes sign. Below it the flow
+response is downstream-dominated (sub-resonant); above it, upstream-dominated
+(super-resonant). Which side a run sits on controls the qualitative behaviour
+of bend migration, so it is worth reporting rather than leaving implicit.
+
+``beta_R`` depends on the other physical inputs. For flagbed = 2, r = 0.5:
+
+    theta0 = 0.20            beta_R ~  8.97
+    theta0 = 0.30            beta_R ~  9.68
+    theta0 = 0.50            beta_R ~  9.50
+    ds     = 0.002           beta_R ~ 11.34
+    ds     = 0.005           beta_R ~  9.68
+    ds     = 0.020           beta_R ~  7.19
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from .flowfield import _precompute_modes
+from .resistance import resistance_function_flagbed
+
+#: Relative distance from beta_R within which a run is reported as near-resonant.
+NEAR_RESONANCE_BAND = 0.02
+
+
+def fundamental_decay_rate(
+    beta: float,
+    theta0: float,
+    ds: float,
+    rpic_0: float,
+    flagbed: int = 2,
+    Mdat: int = 6,
+) -> float:
+    """Return ``Re(lambda2)`` for the fundamental lateral mode.
+
+    Negative values are sub-resonant; positive values are super-resonant.
+    """
+    rpic, cf0, ct, cd, phit, phid, f0 = resistance_function_flagbed(
+        int(flagbed), float(theta0), float(ds), float(rpic_0)
+    )
+    result = _precompute_modes(
+        cf0,
+        ct,
+        cd,
+        phit,
+        phid,
+        float(beta),
+        rpic,
+        float(theta0),
+        f0,
+        int(Mdat),
+    )
+    lamb2 = result[2]
+    if len(lamb2) == 0:
+        return float("nan")
+    return float(np.real(lamb2[0]))
+
+
+def resonant_aspect_ratio(
+    theta0: float,
+    ds: float,
+    rpic_0: float,
+    flagbed: int = 2,
+    Mdat: int = 6,
+    bracket: tuple[float, float] = (2.0, 200.0),
+    tolerance: float = 1.0e-6,
+) -> float | None:
+    """Bisect for the ``beta`` at which the fundamental decay rate vanishes.
+
+    Returns ``None`` when the bracket does not contain a sign change.
+    """
+    lo, hi = float(bracket[0]), float(bracket[1])
+    args = (theta0, ds, rpic_0, flagbed, Mdat)
+    f_lo = fundamental_decay_rate(lo, *args)
+    f_hi = fundamental_decay_rate(hi, *args)
+    if not np.isfinite(f_lo) or not np.isfinite(f_hi) or f_lo * f_hi > 0.0:
+        return None
+
+    while (hi - lo) > tolerance * max(1.0, lo):
+        mid = 0.5 * (lo + hi)
+        f_mid = fundamental_decay_rate(mid, *args)
+        if not np.isfinite(f_mid):
+            return None
+        if f_mid * f_lo > 0.0:
+            lo = mid
+            f_lo = f_mid
+        else:
+            hi = mid
+
+    return 0.5 * (lo + hi)
+
+
+def resonance_report(
+    beta: float,
+    theta0: float,
+    ds: float,
+    rpic_0: float,
+    flagbed: int = 2,
+    Mdat: int = 6,
+) -> dict:
+    """Summarise where a run sits relative to resonance."""
+    beta = float(beta)
+    args = (theta0, ds, rpic_0, flagbed, Mdat)
+    decay = fundamental_decay_rate(beta, *args)
+    beta_r = resonant_aspect_ratio(*args)
+
+    if decay < 0.0:
+        state = "sub-resonant"
+    elif decay > 0.0:
+        state = "super-resonant"
+    else:
+        state = "resonant"
+
+    distance = None
+    if beta_r is not None and beta_r > 0.0:
+        distance = beta / beta_r - 1.0
+        if abs(distance) <= NEAR_RESONANCE_BAND:
+            state = "near-resonant"
+
+    if decay == 0.0:
+        influence_length = float("inf")
+    elif np.isfinite(decay):
+        influence_length = 1.0 / abs(decay)
+    else:
+        influence_length = float("nan")
+
+    return {
+        "state": state,
+        "flag": 1 if decay < 0.0 else -1,
+        "beta": beta,
+        "resonant_beta": beta_r,
+        "relative_distance_to_resonance": distance,
+        "fundamental_decay_rate": decay,
+        "influence_length_half_widths": float(influence_length),
+    }

tests/test_resonance_diagnostics.py

@@ -0,0 +1,115 @@
+"""The resonance flag must track the fundamental mode, not mode 6.
+
+``_compute_flag`` previously read ``lamb2[5]``. ``Re(lambda2)`` for mode 6 is
+roughly -5 to -13 across realistic aspect ratios, so the flag reported
+"sub-resonant" unconditionally. The fundamental mode is the one that changes
+sign at the resonant aspect ratio.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from ldsfl.flowfield import _compute_flag, _precompute_modes
+from ldsfl.flowfield_periodic import _compute_flag as _compute_periodic_flag
+from ldsfl.resistance import resistance_function_flagbed
+from ldsfl.resonance import (
+    fundamental_decay_rate,
+    resonance_report,
+    resonant_aspect_ratio,
+)
+
+THETA0 = 0.3
+DS = 0.005
+RPIC0 = 0.5
+
+
+def _lamb2(beta: float, mdat: int = 6):
+    rpic, cf0, ct, cd, phit, phid, f0 = resistance_function_flagbed(2, THETA0, DS, RPIC0)
+    return _precompute_modes(cf0, ct, cd, phit, phid, beta, rpic, THETA0, f0, mdat)[2]
+
+
+def test_free_and_periodic_flags_change_sign_across_the_resonant_aspect_ratio():
+    beta_r = resonant_aspect_ratio(THETA0, DS, RPIC0)
+    assert beta_r is not None
+
+    below = _lamb2(beta_r * 0.9)
+    above = _lamb2(beta_r * 1.1)
+
+    assert _compute_flag(below) == 1, "expected sub-resonant below beta_R"
+    assert _compute_flag(above) == -1, "expected super-resonant above beta_R"
+    assert _compute_periodic_flag(below) == 1, "periodic flag should use the same convention"
+    assert _compute_periodic_flag(above) == -1, "periodic flag should use the same convention"
+
+
+def test_flag_is_not_constant_over_the_realistic_aspect_ratio_range():
+    """The previous implementation returned +1 for every one of these."""
+    flags = {_compute_flag(_lamb2(beta)) for beta in (4.0, 6.0, 9.0, 12.0, 20.0, 40.0)}
+    assert flags == {1, -1}
+
+
+def test_mode_six_decay_rate_is_uninformative():
+    """Documents why indexing mode 6 could never work."""
+    for beta in (4.0, 9.0, 20.0, 40.0):
+        assert _lamb2(beta)[5].real < -1.0
+
+
+def test_resonant_aspect_ratio_matches_the_sign_change():
+    beta_r = resonant_aspect_ratio(THETA0, DS, RPIC0)
+    assert beta_r is not None
+    assert fundamental_decay_rate(beta_r * 0.99, THETA0, DS, RPIC0) < 0.0
+    assert fundamental_decay_rate(beta_r * 1.01, THETA0, DS, RPIC0) > 0.0
+
+
+def test_resonant_aspect_ratio_moves_with_grain_size():
+    coarse = resonant_aspect_ratio(THETA0, 0.02, RPIC0)
+    fine = resonant_aspect_ratio(THETA0, 0.002, RPIC0)
+    assert coarse is not None and fine is not None
+    assert fine > coarse, "beta_R should increase as relative grain size decreases"
+
+
+@pytest.mark.parametrize(
+    "beta, expected",
+    [
+        (6.0, "sub-resonant"),
+        (9.0, "sub-resonant"),
+        (12.0, "super-resonant"),
+        (20.0, "super-resonant"),
+    ],
+)
+def test_resonance_report_states(beta, expected):
+    assert resonance_report(beta, THETA0, DS, RPIC0)["state"] == expected
+
+
+def test_influence_length_grows_towards_resonance():
+    """A domain must be long compared with 1/|Re(lambda)|, which diverges at beta_R."""
+    beta_r = resonant_aspect_ratio(THETA0, DS, RPIC0)
+    assert beta_r is not None
+    far = resonance_report(6.0, THETA0, DS, RPIC0)["influence_length_half_widths"]
+    near = resonance_report(beta_r * 0.999, THETA0, DS, RPIC0)["influence_length_half_widths"]
+    assert np.isfinite(far)
+    assert near > 10.0 * far
+
+
+def test_resonance_report_is_included_in_run_summary(tmp_path):
+    from ldsfl.main import run_case
+
+    input_dir = tmp_path / "Input"
+    input_dir.mkdir()
+    (input_dir / "Parameter.csv").write_text(
+        "Id,Beta,ds,Thetha,flagbed,r,Mdat,flagbed=1 plane; flagbed=2 dunes\n"
+        "1,9.0,0.005,0.3,2,0.5,6,2\n",
+        encoding="utf-8",
+    )
+    (input_dir / "xy.csv").write_text(
+        "0,0\n1,0.01\n2,0.03\n3,0.06\n4,0.08\n5,0.09\n",
+        encoding="utf-8",
+    )
+
+    result = run_case(tmp_path, 1, max_steps=2, Nprint=2, do_plots=False)
+
+    report = result["resonance"]
+    assert report["state"] in {"sub-resonant", "near-resonant", "super-resonant", "resonant"}
+    assert report["flag"] in {1, -1}
+    assert report["resonant_beta"] is None or report["resonant_beta"] > 0.0