Fix solver validation, history semantics, and Numba parity

.github/workflows/tests.yml

@@ -9,6 +9,23 @@ permissions:
   contents: read
 
 jobs:
+  numba-equivalence:
+    name: Numba backend equivalence / Python 3.12
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout repository
+        uses: actions/checkout@v4
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+      - name: Install development and Numba dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install -e ".[dev,numba]"
+      - name: Run backend guardrails
+        run: python -m pytest tests/test_flowfield_backend_guardrails.py -q
+
   smoke-tests:
     name: Tests / Python ${{ matrix.python-version }}
     runs-on: ubuntu-latest

CHANGELOG.md

@@ -8,6 +8,12 @@ This project follows a lightweight changelog style inspired by [Keep a Changelog
 
 ### Fixed
 
+- Corrected the optional Numba `SL=1` sub-resonant routing so its lambda-2 and lambda-3 terms use the same downstream `SEMIANA2` kernels as the NumPy reference; dedicated CI now exercises Numba parity.
+- Centralized scientific parameter and core run-control validation, including exact integer checks for `Id`, `flagbed`, `Mdat`, save intervals, worker counts, and other discrete controls.
+- Rejected dimensional output unless finite positive physical length and velocity scales are available; the CLI now requires both scales explicitly.
+- Made `--no-plots` suppress every PNG, while sinuosity CSV rows are appended incrementally and exactly once.
+- Added `state_step` to variable histories and included the final post-step state, aligning variable and sinuosity histories from step 0 through the final completed step.
+- Made missing final numerical outputs fatal and recorded finalization state and errors in `run_status.json`.
 - Included freshwater density when converting bed shear stress in pascals to Shields stress.
 - Preserved previous simulations by reserving a separate output directory for every run, including continuation segments and colliding historical labels.
 - Preserved dimensional output scaling during GUI continuation and aligned initial overlays with output units.
@@ -23,6 +29,7 @@ This project follows a lightweight changelog style inspired by [Keep a Changelog
 
 ### Changed
 
+- Variable histories now contain `N + 1` state rows for an `N`-iteration run; `jt` remains as the legacy counter and `state_step` is the completed-iteration count.
 - Documentation and release hygiene are now tracked explicitly through contributor and changelog files.
 
 ## [0.6.5] - documented release

USER_MANUAL.md

@@ -16,6 +16,12 @@ an unknown ID is rejected before the batch starts.
 
 Use `gui_ldsfl.py` when you want help building input files, checking conversions, scaling and filtering geometry, selecting stop criteria, and viewing the final planform.
 
+All entry points share the same scientific validation. `Parameter.csv` must
+contain `Id`, `Beta`, `ds`, `Thetha`, `flagbed`, `r`, and `Mdat`. Numeric values
+must be finite; positive quantities must be greater than zero; `Id`, `flagbed`,
+and `Mdat` must be exact integers. `flagbed` must be 1 or 2, IDs must be unique,
+and fractional values such as `Mdat=6.5` are rejected rather than truncated.
+
 ## 3. Core notation
 
 In this project:
@@ -60,6 +66,11 @@ mode. At least one criterion must have a positive limit, or sinuosity stability
 stopping must be enabled. Negative and non-finite stopping limits are rejected
 by the solver.
 
+These rules also apply to direct `run_case()` calls. Save intervals, iteration
+labels, worker counts, smoothing and resampling settings, stability settings,
+flow options, and output units are checked before a run creates its output
+directory. For example, `--nprint 0` fails immediately with a validation error.
+
 ## 6. Bank erodibility and resonance
 
 The advanced GUI includes **Bank erodibility / erosion rate**. The historical
@@ -84,6 +95,7 @@ LDSFL-Meander writes case outputs under `Output/<id_files>/`, including:
 - `plot/`
 - `files/` (saved run variables)
 - `run_config.json` (solver options, initial parameters, and input SHA-256 hashes)
+- `run_status.json` (simulation completion, output completeness, stop reason, and output errors)
 - `run_manifest.json`
 - `gui_final_overlay.png`
 
@@ -98,6 +110,38 @@ Dimensional length outputs use the physical half-width `B_0` independently of
 how the input geometry was prepared. Already dimensionless coordinates and
 continuation geometry therefore still produce metres when requested.
 
+Dimensional output is accepted only when both `B_0` and a finite positive
+physical reference velocity `U0` are known. GUI input modes that cannot derive
+`U0` must be supplemented with velocity/friction information or use
+dimensionless output. CLI and direct API users must explicitly provide
+`output_length_scale` and `output_velocity_scale`; no scale of 1 is silently
+treated as physical.
+
+With plotting disabled, including through `--no-plots`, no PNG is created. The
+sinuosity CSV is still required and is appended in blocks without rewriting
+the complete history. Its `step` values run from 0 through the final completed
+iteration exactly once. The optional sinuosity PNG is generated only at final
+completion when plots are enabled.
+
+The `var_*.csv` files use two counters:
+
+- `state_step` is the number of completed migration iterations represented by the row;
+- `jt` is the retained legacy/internal loop counter and equals `state_step + 1`.
+
+The initial row has `state_step=0`, `jt=1`, `dt=0`, and `dt_cum=0`. The `dt`
+field is the adaptive timestep used to reach the represented state, so a run
+that completes `N` iterations records states 0 through `N`, including the final
+post-step geometry and cumulative time.
+
+The final `xyu` snapshot, variable-history flush, and sinuosity CSV are required
+scientific outputs. If any cannot be written, the solver attempts the other
+final outputs, records the errors in `run_status.json`, and raises `RuntimeError`
+instead of returning a successful result. Plot failures remain non-fatal.
+
+The optional Numba backend now follows the NumPy reference for both
+sub-resonant and super-resonant `SL=1` flow evaluation, and CI checks this
+equivalence when installing the `numba` extra.
+
 ## 8. Full manual
 
 See `docs/LDSFL_Meander_user_manual.pdf` for the full LaTeX manual.

docs/timestep_notes.md

@@ -21,6 +21,23 @@ contains:
 | `step` | Solver iteration index. |
 | `sinuo` | Dimensionless sinuosity at that iteration. |
 
+Variable-history files also contain `state_step` and the legacy `jt` counter.
+`state_step` is the completed migration-iteration count and uses the same
+0-through-N convention as `sinuosity_history.step`. `jt` is retained for
+backward compatibility and equals `state_step + 1`.
+
+Each variable-history row describes a state rather than an iteration about to
+run. The initial state is:
+
+```text
+state_step = 0, jt = 1, dt = 0, dt_cum = 0
+```
+
+After the first migration update, the next row is `state_step=1`, `jt=2`.
+Here `dt` is the adaptive timestep used to reach that represented state and
+`dt_cum` includes it. A run completing N iterations writes N+1 state rows and
+the final row's `dt_cum` equals the value returned in the run result.
+
 The relationship between model iterations and dimensional physical time depends on the nondimensionalisation and physical scaling used for the study. Do not interpret `step` as seconds, days, years, or any dimensional unit unless the corresponding scaling has been defined.
 
 ## Adaptive computational timestep

examples/reproducible_case1_short/expected_tree.txt

@@ -5,7 +5,7 @@ expected_output/1_9_0005_03_2_05/files/var_1_9_0005_03_2_05_000000_2.csv
 expected_output/1_9_0005_03_2_05/files/var_1_9_0005_03_2_05_000000_4.csv
 expected_output/1_9_0005_03_2_05/files/var_1_9_0005_03_2_05_000000_6.csv
 expected_output/1_9_0005_03_2_05/plot/
-expected_output/1_9_0005_03_2_05/plot/sinuosity_history_1_9_0005_03_2_05.png
+expected_output/1_9_0005_03_2_05/run_status.json
 expected_output/1_9_0005_03_2_05/xy_cut/
 expected_output/1_9_0005_03_2_05/xyu/
 expected_output/1_9_0005_03_2_05/xyu/xyu_1_9_0005_03_2_05_000000_2.csv

gui_smoke_test.py

@@ -5,10 +5,10 @@
 from pathlib import Path
 
 from ldsfl.gui_utils import (
-    DimensionlessInputs,
     DimensionalInputs,
-    GuiCaseConfig,
+    DimensionlessInputs,
     GeometrySettings,
+    GuiCaseConfig,
     RunControls,
     config_from_dict,
     config_to_dict,
@@ -69,6 +69,7 @@ def main() -> None:
                 d50=0.05,
                 mobility_mode='direct_shields',
                 theta0=0.3,
+                velocity=1.2,
                 flagbed=2,
                 rpic_0=0.5,
                 Mdat=6,
@@ -155,6 +156,7 @@ def main() -> None:
                 d50=0.05,
                 mobility_mode='direct_shields',
                 theta0=0.3,
+                velocity=1.2,
                 flagbed=2,
                 rpic_0=0.5,
                 Mdat=6,
@@ -176,6 +178,7 @@ def main() -> None:
                 d50=0.05,
                 mobility_mode='direct_shields',
                 theta0=0.3,
+                velocity=1.2,
                 flagbed=2,
                 rpic_0=0.5,
                 Mdat=6,

ldsfl/evolution.py

@@ -26,6 +26,8 @@ def dxdy2(ERT: float, U: np.ndarray, x: np.ndarray, y: np.ndarray,
         raise FloatingPointError(
             f"Non-finite migration speed encountered in dxdy2: ERT={ERT!r} at jt={jt}"
         )
+    if ERT <= 0.0:
+        raise ValueError(f"ERT must be > 0 in dxdy2; got {ERT!r} at jt={jt}")
     if not np.all(np.isfinite(U)):
         raise FloatingPointError(
             f"Non-finite migration speed encountered in dxdy2: U contains non-finite values at jt={jt}"
@@ -38,6 +40,8 @@ def dxdy2(ERT: float, U: np.ndarray, x: np.ndarray, y: np.ndarray,
         raise FloatingPointError(
             f"Non-finite timestep stability coefficient encountered in dxdy2: cstab={cstab!r} at jt={jt}"
         )
+    if cstab <= 0.0:
+        raise ValueError(f"cstab must be > 0 in dxdy2; got {cstab!r} at jt={jt}")
 
     dxdtl = np.empty_like(U)
     dydtl = np.empty_like(U)

ldsfl/flowfield_numba.py

@@ -16,7 +16,6 @@
 
 from __future__ import annotations
 
-
 import numpy as np
 
 # --- Compatibility shim for a known `numba` ↔ `coverage` API mismatch ---
@@ -28,12 +27,12 @@
     import coverage  # type: ignore
 
     if not hasattr(coverage, "types"):
-        class _Types:  # noqa: WPS431 (nested class)
+        class _Types:
             pass
         coverage.types = _Types()  # type: ignore[attr-defined]
 
     if not hasattr(coverage.types, "Tracer"):
-        class _Tracer:  # noqa: WPS431
+        class _Tracer:
             pass
         coverage.types.Tracer = _Tracer  # type: ignore[attr-defined]
 
@@ -52,7 +51,7 @@ class _Tracer:  # noqa: WPS431
     pass
 
 try:
-    from numba import njit
+    from numba import njit, prange
 except Exception as e:  # pragma: no cover
     # Provide a clearer error when imported without numba.
     raise ImportError(
@@ -263,10 +262,6 @@ def _jtoll_for_lam(lam: complex, s_pad: np.ndarray, toll: float, *, upwind: bool
 # Public entry points
 # -----------------------------
 
-
-
-from numba import prange
-
 _MODE_FN_CACHE = {}
 _MODE_ADDER_CACHE = {}
 
@@ -492,7 +487,8 @@ def jtoll_sl1(real_lam_abs: float) -> int:
             j3toll = jtoll_sl1(abs(lm3r))
             j4toll = jtoll_sl1(abs(lm4r))
 
-            # um1: -Am*(g10*sum1 + g20*sum2 + g30*sum3) using SEMIANA1
+            # Sub-resonant reference algebra: lambda1 is upstream/SEMIANA1;
+            # lambda2, lambda3 and lambda4 are downstream/SEMIANA2.
             if use_cached:
                 # lm1 (allow_pos_k False)
                 kmax1 = max(1, min(N, j1toll - 1))
@@ -502,20 +498,20 @@ def jtoll_sl1(real_lam_abs: float) -> int:
                 # lm2/lm3 in the reference use allow_pos_k True
                 kmax2 = max(1, min(N, j2toll - 1))
                 cA2, cB2, lmds2, exp2, kmax2 = _cached_tables_for_lam_np(lm2, c_pad, deltas, kmax2, allow_pos_k=True)
-                fill_upstr_cached_real(out, cA2, cB2, lmds2, exp2, kmax2, j2toll, N, complex(mode_scale * (-Aj), 0.0) * g20jm)
+                fill_dwstr_cached_real(out, cA2, cB2, lmds2, exp2, kmax2, j2toll, N, complex(mode_scale * Aj, 0.0) * g20jm)
 
                 kmax3 = max(1, min(N, j3toll - 1))
                 cA3, cB3, lmds3, exp3, kmax3 = _cached_tables_for_lam_np(lm3, c_pad, deltas, kmax3, allow_pos_k=True)
-                fill_upstr_cached_real(out, cA3, cB3, lmds3, exp3, kmax3, j3toll, N, complex(mode_scale * (-Aj), 0.0) * g30jm)
+                fill_dwstr_cached_real(out, cA3, cB3, lmds3, exp3, kmax3, j3toll, N, complex(mode_scale * Aj, 0.0) * g30jm)
 
                 # um2: +Am*g40*sum4 using SEMIANA2
                 kmax4 = max(1, min(N, j4toll - 1))
                 cA4, cB4, lmds4, exp4, kmax4 = _cached_tables_for_lam_np(lm4, c_pad, deltas, kmax4, allow_pos_k=True)
                 fill_dwstr_cached_real(out, cA4, cB4, lmds4, exp4, kmax4, j4toll, N, complex(mode_scale * Aj, 0.0) * g40jm)
             else:
                 fill_upstr_direct_real(out, c_pad, float(deltas), lm1, j1toll, N, complex(mode_scale * (-Aj), 0.0) * g10jm)
-                fill_upstr_direct_real(out, c_pad, float(deltas), lm2, j2toll, N, complex(mode_scale * (-Aj), 0.0) * g20jm)
-                fill_upstr_direct_real(out, c_pad, float(deltas), lm3, j3toll, N, complex(mode_scale * (-Aj), 0.0) * g30jm)
+                fill_dwstr_direct_real(out, c_pad, float(deltas), lm2, j2toll, N, complex(mode_scale * Aj, 0.0) * g20jm)
+                fill_dwstr_direct_real(out, c_pad, float(deltas), lm3, j3toll, N, complex(mode_scale * Aj, 0.0) * g30jm)
                 fill_dwstr_direct_real(out, c_pad, float(deltas), lm4, j4toll, N, complex(mode_scale * Aj, 0.0) * g40jm)
 
         elif lm2r > 0.0: