Analytical tests batch 2

.gitignore

@@ -1,6 +1,10 @@
 
 .DS_Store
 
+# Working notes not intended for version control
+doc/FlashX_SWMM_crosswalk.md
+doc/Verification_Implementation_Plan.md
+
 # Eclipse Stuff
 .metadata/
 .settings/

src/engine/hydraulics/KinematicWave.cpp

@@ -82,8 +82,8 @@ int KWSolver::solveConduit(int idx, const XSectParams& xs,
     } else if (q_in_norm <= 0.0) {
         a_in_norm = 0.0;
     } else {
-        double s_needed = q_in_norm / beta1;
-        a_in_norm = xsect::getAofS(xs, s_needed * a_full) / a_full;
+        double s_needed = q_in_norm / beta1;  // dimensional section factor = Q_in/beta
+        a_in_norm = xsect::getAofS(xs, s_needed) / a_full;
     }
 
     // Finite-difference coefficients

src/engine/hydrology/LID.cpp

@@ -501,16 +501,29 @@ void LIDSolver::batchBarrelFlux(LIDGroupSoA& g, double rainfall, double dt) {
             new_depth -= drain * dt;
         }
 
+        // Exfiltration through storage floor (Euler is exact for constant-rate ODE)
+        double exfil = 0.0;
+        if (g.stor_void[ui] > 0.0 && new_depth > 0.0) {
+            double ksat_eff = getStorageExfil(g.stor_ksat[ui], g.stor_clog[ui],
+                                              g.wb_inflow[ui]);
+            if (ksat_eff > 0.0) {
+                double exfil_depth = (ksat_eff / g.stor_void[ui]) * dt;
+                exfil_depth = std::min(exfil_depth, new_depth);
+                new_depth -= exfil_depth;
+                exfil = exfil_depth * g.stor_void[ui];
+            }
+        }
+
         g.stor_depth[ui] = std::max(new_depth, 0.0);
         g.surface_runoff[ui] = overflow;
         g.drain_flow[ui] = drain;
         g.evap_loss[ui] = 0.0;
-        g.infil_loss[ui] = 0.0;
+        g.infil_loss[ui] = (dt > 0.0) ? exfil / dt : 0.0;
 
         // Water balance tracking
         g.wb_inflow[ui]     += unit_inflow * dt;
         g.wb_evap[ui]       += 0.0;
-        g.wb_infil[ui]      += 0.0;
+        g.wb_infil[ui]      += exfil;
         g.wb_surf_flow[ui]  += overflow * dt;
         g.wb_drain_flow[ui] += drain * dt;
         g.wb_final_vol[ui]   = g.stor_depth[ui];

tests/benchmarks/README.md

@@ -0,0 +1,160 @@
+# Benchmark Data Layout
+
+## Purpose
+
+On `swmm6_rel`, this directory is already the Google Benchmark performance subtree for `openswmm.engine`.
+
+It can also hold reference datasets and supporting metadata used for solver verification, but those assets should live in subdirectories so they do not interfere with the top-level benchmark executables and CMake files.
+
+The intended use case is to retain benchmark artifacts generated externally, including with internal Flash-X workflows, without importing Flash-X source code into this repository.
+
+## Principles
+
+1. Store generated reference data, not external solver source.
+2. Keep provenance next to the dataset.
+3. Prefer open, stable, text-based formats.
+4. Separate benchmark data from test code.
+5. Make each dataset reproducible from its metadata and generation notes.
+
+## Recommended Layout
+
+```text
+tests/
+  benchmarks/
+    CMakeLists.txt
+    bench_engine_vs_legacy.cpp
+    bench_timeseries_lookup.cpp
+    bench_hydraulics.cpp
+    README.md
+    provenance-template.md
+    generated/
+      <benchmark-name>/
+        README.md
+        provenance.yaml
+        reference.csv
+        reference.json
+        notes.md
+    manufactured/
+      <benchmark-name>/
+        README.md
+        definition.md
+        reference.csv
+    shared/
+      units.md
+      conventions.md
+```
+
+## Directory Roles
+
+### Top-level files
+
+The top level of `tests/benchmarks/` is reserved for Google Benchmark performance targets already used by this branch.
+
+Current examples:
+
+- `bench_engine_vs_legacy.cpp`
+- `bench_timeseries_lookup.cpp`
+- `bench_hydraulics.cpp`
+- `CMakeLists.txt`
+
+Do not place verification datasets directly alongside those files. Put verification assets in the subdirectories below.
+
+### `generated/`
+
+For benchmark outputs created by an external tool or solver run.
+
+Typical examples:
+
+- time histories,
+- reference hydrographs,
+- depth profiles,
+- infiltration trajectories,
+- storage-loss tables,
+- tabulated exact or semi-analytic data.
+
+Each benchmark directory should contain:
+
+- `README.md`: what the dataset represents and how tests consume it,
+- `provenance.yaml`: machine-readable generation metadata,
+- `reference.csv` or `reference.json`: the actual reference values,
+- `notes.md`: optional generation caveats, tolerances, or assumptions.
+
+### `manufactured/`
+
+For hand-defined or analytically derived verification cases where the reference solution is defined by formulas or compact tables rather than an external production run.
+
+Typical examples:
+
+- smooth manufactured dynamic-wave solution,
+- exact scalar ODE trajectories for integrator testing,
+- closed-form infiltration cases,
+- linear storage recession cases.
+
+### `shared/`
+
+For conventions reused by multiple datasets.
+
+Suggested contents:
+
+- unit conventions,
+- time origin conventions,
+- coordinate/sign conventions,
+- variable naming conventions,
+- acceptable interpolation rules between stored points and test query points.
+
+## File Format Guidance
+
+Prefer:
+
+- `CSV` for dense numeric tables,
+- `JSON` for structured reference datasets with metadata-rich records,
+- `YAML` for provenance and configuration,
+- `Markdown` for human-readable notes.
+
+Avoid:
+
+- opaque binary formats when a text format is practical,
+- external-tool-native formats that require the external tool to parse,
+- embedding provenance only in code comments.
+
+## Naming Guidance
+
+Benchmark names should describe both the physics and the scenario, for example:
+
+- `odesolve-exponential-decay`
+- `infil-greenampt-constant-rainfall`
+- `exfil-cylindrical-storage-greenampt`
+- `kinwave-step-inflow-rectangular-conduit`
+
+## Provenance Minimum
+
+Every externally generated benchmark should record:
+
+- generator name,
+- generator version or commit,
+- source problem/setup name,
+- input deck or parameter file identity,
+- date generated,
+- variables exported,
+- units,
+- spatial and temporal sampling,
+- any post-processing steps,
+- known limitations.
+
+For Flash-X-generated data, the metadata should say that the dataset was generated with Flash-X and cite the Flash-X publication, while keeping Flash-X source code outside this repository unless there is an explicit licensing decision to vendor code.
+
+## How Tests Should Consume Data
+
+Tests should:
+
+1. load a single benchmark dataset from `tests/benchmarks/generated/...` or `tests/benchmarks/manufactured/...`,
+2. run the SWMM code path under test,
+3. compare the computed result to the reference values,
+4. report max error, RMS or $L_2$, and conservation or mass-balance error where relevant,
+5. fail with benchmark-specific tolerances stored in code or metadata.
+
+The benchmark dataset should remain immutable once baselined. If regenerated, update provenance and explain why the baseline changed.
+
+## Branch-Specific Note
+
+On `swmm6_rel`, correctness tests belong under `tests/unit/engine/` and `tests/regression/`. This directory remains performance-first; the verification datasets stored here are inputs to those correctness tests, not replacements for them.

tests/benchmarks/manufactured/dynwave-gvf-backwater-m1/definition.md

@@ -0,0 +1,82 @@
+# Benchmark: dynwave-gvf-backwater-m1
+
+## Purpose
+
+Verify that the DW solver (St. Venant equations) converges to the analytically
+computed **M1 backwater profile** (GVF, subcritical, tailwater above normal depth)
+after sufficient time under steady inflow and a fixed tailwater boundary condition.
+
+## Physical setup
+
+A 1000 ft reach of rectangular open channel (b = 5 ft, y_full = 4 ft, S₀ = 0.001,
+n = 0.013) is discretised into 5 conduits of 200 ft each (6 nodes J0–J5).
+A constant lateral inflow Q = 10 cfs enters at the upstream junction J0.
+The downstream node J5 is a FIXED outfall at water-surface elevation y_d = 1.5 × y_n.
+
+## Parameters
+
+| Symbol    | Value     | Unit    | Notes                                      |
+|-----------|-----------|---------|--------------------------------------------|
+| b         | 5.0       | ft      | channel width (RECT_OPEN)                  |
+| y_full    | 4.0       | ft      | conduit full depth (prevents surcharge)    |
+| S₀        | 0.001     | ft/ft   | bed slope                                  |
+| n         | 0.013     | —       | Manning roughness                           |
+| L_conduit | 200.0     | ft      | per-conduit length (5 conduits total)      |
+| Q         | 10.0      | cfs     | steady lateral inflow at J0                |
+| dt        | 30.0      | s       | routing timestep                           |
+| N_steps   | 120       | —       | steps (T = 3600 s = 1 hr)                  |
+
+Derived quantities (PHI = 1.486, US customary):
+```
+beta      = PHI * sqrt(S₀) / n  ≈ 3.6163 ft^{1/3}/s
+y_n       = 0.781692 ft  (Manning normal depth at Q = 10 cfs)
+y_c       = 0.498963 ft  (critical depth at Q = 10 cfs)
+Fr_n      ≈ 0.510  (mild slope: y_n > y_c ✓)
+y_d       = 1.5 × y_n = 1.172539 ft  (fixed tailwater BC)
+```
+
+## Analytical reference: M1 GVF profile
+
+The GVF ODE (gradually varied flow, subcritical regime):
+
+    dy/dx = (S₀ - Sf(y)) / (1 - Fr²(y))
+
+where `Sf = (Q·n / (PHI·A·R^{2/3}))²` and `Fr² = Q² / (g·A²·(A/B))`.
+
+Integrated by RK4 from x = 1000 ft (y = y_d) upstream to x = 0 ft,
+sampling at each node:
+
+| Node | x (ft) | z_inv (ft) | y_GVF (ft) |
+|------|--------|------------|------------|
+| J0   | 0      | 1.000      | 0.792075   |
+| J1   | 200    | 0.800      | 0.809203   |
+| J2   | 400    | 0.600      | 0.848366   |
+| J3   | 600    | 0.400      | 0.921831   |
+| J4   | 800    | 0.200      | 1.032399   |
+| J5   | 1000   | 0.000      | 1.172539   |
+
+The M1 profile asymptotes toward y_n from above as x → 0. At J0, the depth
+is within 1.3% of normal depth (0.792075 vs 0.781692), reflecting the finite
+reach length.
+
+## Expected accuracy
+
+After T = 3600 s (≈ 18 wave-travel times), the DW solver should have converged
+to quasi-steady state. Test tolerance is **5% of y_n ≈ 0.039 ft** at J0–J4.
+J5 is not checked (it is the fixed boundary condition).
+
+The GVF reference was computed with g = 32.174 ft/s²; the DW solver uses
+g = 32.2 ft/s². The 0.08% difference in g is negligible versus the 5% tolerance.
+
+## Verification target
+
+`src/engine/hydraulics/DynamicWave.cpp` — `DWSolver::execute`
+
+## Consuming test
+
+`tests/unit/engine/test_routing.cpp` — `TEST(DWSolverGVF, BackwaterM1Benchmark)`
+
+## References
+
+- Henderson, F.M. (1966). *Open Channel Flow*. Macmillan. Chapter 5 (gradually varied flow, M1 backwater curve).
+- Manning, R. (1891). On the flow of water in open channels and pipes. *Transactions of the Institution of Civil Engineers of Ireland*, 20, 161–207.

tests/benchmarks/manufactured/dynwave-gvf-backwater-m1/provenance.yaml

@@ -0,0 +1,90 @@
+benchmark: dynwave-gvf-backwater-m1
+type: manufactured
+description: >
+  GVF M1 backwater profile on a 5-conduit rectangular open channel reach.
+  Under constant inflow Q = 10 cfs and fixed tailwater y_d = 1.5 y_n, the
+  St. Venant equations admit the M1 GVF profile as the unique subcritical
+  steady state.  Reference depths are computed by RK4 integration of the
+  GVF ODE from downstream to upstream.
+swmm_target: src/engine/hydraulics/DynamicWave.cpp (DWSolver::execute)
+
+problem:
+  name: gvf-m1-backwater
+  ode: "dy/dx = (S0 - Sf(y)) / (1 - Fr^2(y))"
+  ode_direction: "integrated from x=1000 (downstream) to x=0 (upstream)"
+  friction: "Manning: Sf = (Q*n / (PHI*A*R^{2/3}))^2"
+  froude:   "Fr^2 = Q^2 / (g*A^2*(A/B))"
+  integration_scheme: RK4 with step dx = -200 ft (upstream direction)
+  boundary_condition: "y(x=1000) = y_d = 1.5 * y_n"
+
+parameters:
+  b_ft: 5.0
+  y_full_ft: 4.0
+  S0: 0.001
+  n_mann: 0.013
+  L_conduit_ft: 200.0
+  Q_cfs: 10.0
+  PHI: 1.486
+  g_fts2: 32.174
+  unit_system: US
+  derived:
+    beta_ft13_per_s: 3.616272    # PHI * sqrt(S0) / n
+    y_n_ft: 0.781692              # normal depth at Q = 10 cfs
+    y_c_ft: 0.498963              # critical depth at Q = 10 cfs
+    Fr_n: 0.510                   # Froude at normal depth (mild slope: Fr_n < 1)
+    y_d_ft: 1.172539              # tailwater BC = 1.5 * y_n
+
+dataset:
+  file: reference.csv
+  columns:
+    node:      "junction label (J0–J5)"
+    x_ft:      "station from upstream end (ft)"
+    z_inv_ft:  "invert elevation (ft)"
+    y_gvf_ft:  "GVF steady-state depth (ft)"
+  n_points: 6
+  note: >
+    J5 is the fixed outfall boundary; its depth equals y_d by construction
+    and is not checked by the consuming test.
+
+tolerances:
+  y_max_abs_ft: 0.039            # 5% of y_n = 0.05 * 0.781692
+  nodes_checked: [J0, J1, J2, J3, J4]
+  basis: >
+    The DW solver uses a 5-conduit spatial discretisation (Δx = 200 ft) and
+    a θ-implicit Preissmann scheme.  After T = 3600 s the flow is within
+    numerical steady state.  Spatial discretisation error on such a coarse
+    grid is estimated at < 2% of y_n; 5% provides comfortable margin.
+
+attribution:
+  references:
+    gvf: >
+      Henderson, F.M. (1966). Open Channel Flow.  Macmillan.  Chapter 5
+      (Gradually Varied Flow — M1 backwater curve).
+    manning: >
+      Manning, R. (1891). On the flow of water in open channels and pipes.
+      Trans. Inst. Civil Eng. Ireland, 20, 161-207.
+  license_note: benchmark data are original to this project; no copied source code.
+
+provenance:
+  generated_by: >
+    Python script using fixed-step RK4 with dx = -200 ft to integrate
+    dy/dx = (S0 - Sf) / (1 - Fr^2) from x=1000 to x=0,
+    with boundary condition y(x=1000) = y_d = 1.5 * y_n.
+  date: "2026-05-14"
+  author: openswmm.engine development team
+  computation: |
+    # Normal depth: solve Q = beta * A(y) * R(y)^{2/3} for y
+    # beta = 1.486 * sqrt(0.001) / 0.013 = 3.616272
+    # A(y) = 5*y,  R(y) = 5*y/(2*y+5)
+    # y_n = 0.781692 ft  (bisection)
+    # y_d = 1.5 * y_n = 1.172539 ft
+    #
+    # GVF ODE integrated backwards (dx = -200 ft per step):
+    # dy/dx = f(y) = (0.001 - Sf(y)) / (1 - Fr^2(y))
+    # Sample at x = 1000, 800, 600, 400, 200, 0 ft
+    #
+    # Values rounded to 6 decimal places.
+  reproducible: true
+  notes: >
+    The DW solver uses g = 32.2 ft/s² (SWMM legacy); the reference used
+    g = 32.174 ft/s².  The 0.08% difference is negligible vs the 5% test tolerance.

tests/benchmarks/manufactured/dynwave-gvf-backwater-m1/reference.csv

@@ -0,0 +1,7 @@
+node,x_ft,z_inv_ft,y_gvf_ft
+J0,0,1.000,0.792075
+J1,200,0.800,0.809203
+J2,400,0.600,0.848366
+J3,600,0.400,0.921831
+J4,800,0.200,1.032399
+J5,1000,0.000,1.172539