From 11a3e937a176e7c10d7d67a3eb1ab56ca0317acc Mon Sep 17 00:00:00 2001 From: Sergio Date: Wed, 22 Jul 2026 13:25:06 +0100 Subject: [PATCH] Add resonance flag diagnostics --- ldsfl/flowfield.py | 15 ++- ldsfl/flowfield_periodic.py | 15 ++- ldsfl/main.py | 2 + ldsfl/resonance.py | 139 ++++++++++++++++++++++++++++ tests/test_resonance_diagnostics.py | 115 +++++++++++++++++++++++ 5 files changed, 282 insertions(+), 4 deletions(-) create mode 100644 ldsfl/resonance.py create mode 100644 tests/test_resonance_diagnostics.py diff --git a/ldsfl/flowfield.py b/ldsfl/flowfield.py index d3e7b6b..06717c4 100644 --- a/ldsfl/flowfield.py +++ b/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( *, diff --git a/ldsfl/flowfield_periodic.py b/ldsfl/flowfield_periodic.py index 23c4cc5..285f123 100644 --- a/ldsfl/flowfield_periodic.py +++ b/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) diff --git a/ldsfl/main.py b/ldsfl/main.py index a69a1f2..8a61756 100644 --- a/ldsfl/main.py +++ b/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, } diff --git a/ldsfl/resonance.py b/ldsfl/resonance.py new file mode 100644 index 0000000..5046d54 --- /dev/null +++ b/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), + } diff --git a/tests/test_resonance_diagnostics.py b/tests/test_resonance_diagnostics.py new file mode 100644 index 0000000..6bd0a1d --- /dev/null +++ b/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