Add Old River Control Structure Capability#113
Conversation
…VS Code from shell
…. conus_lakes test passing.
…handle the string IDs.
…raised. Add plots to give visual reassurance that the runs produced good results. Delete old outputs before tests run
…or of catch-all ciss creek
Reduce the gage crosswalk to exactly one routing link per gage. vfp_id is not unique in the link table because flowpaths longer than the discretization length are split into several routing links that all inherit the parent vfp_id, so the join was one-to-many and every gage was registered at every link of its flowpath, a median of 8 on the CONUS hydrofabric. That assimilated a single observation at many consecutive segments and split the cached execution plan at all of them. The gage's own segment_order is a different quantity from the link table's and matches only about a fifth of gages, so placement uses the flowpath outlet link, the same rule already applied to diversion gages. Unplaceable gages are now dropped with a warning instead of entering the map as a null key, and gages sharing a link are reported rather than silently collapsing. Reject a diversion configured without streamflow nudging. The kernel subtracts the observed diversion from the donor flowpath, but nothing adds it to the receiving river: that water arrives only because nudging imposes the observed discharge at the diversion gage on a headwater of the receiving system. With nudging off the subtraction still ran and the transferred water was destroyed. Give diversion_da a typed model instead of a free-form dict, so a misspelled key is rejected rather than silently ignored on a feature that moves water between river systems. Pass diversion_da from the BMI driver. It defaulted to an empty dict there, so the diversion was silently disabled on the path ngen drives. Build the historical-median fallback with one column per routing timestep. The columns were laid out on a fixed 5 minute grid while the kernel indexes them by routing timestep, so with dt=60 the fallback ran out about a fifth of the way through a run and with dt=900 the timestamps advanced three times too slowly. Substituted values are now logged with their share of the window, and a gage with no routing link is reported instead of raising KeyError. Warn when a donor segment is present in a compute job but its gage is not, which previously skipped the diversion with no signal. Hoist an empty-map test out of the kernel's segment loop so runs that divert nothing do not pay a dict membership check per segment per timestep.
…odies preprocess_data_assimilation returned early whenever the reservoir_da table or the waterbody table was empty, leaving self.gages empty. Reservoir DA needs lakes, but streamflow and diversion DA do not, so on any lake-free domain gage nudging was silently disabled even with a fully populated gages layer. On the Ohio benchmark subset, which has zero lakes, this resolved 0 gages instead of 318. Split the two concerns: the reservoir crosswalks are still skipped, with a warning when waterbodies exist but carry no reservoir_da records, and the streamflow and diversion crosswalk now runs on both paths. Rename the link table's segment_order to link_segment_order when joining. The gages layer carries a column of the same name holding a different quantity, so the merge produced suffixed columns and the outlet-link lookup raised KeyError.
The only diversion test was an integration case that was commented out entirely, so nothing in the suite exercised this feature. Add unit coverage for the parts that need no subset hydrofabric, forcing or timeslices, so the mechanism is guarded on every run rather than only where the generated data exists. test_diversion_da covers the job-local kernel map, including the two ways it can decline to divert: a donor outside the job, which guards the binary search against landing on an unrelated segment, and a donor whose gage is absent, which now warns instead of skipping silently. It also covers the climatological fallback: one column per routing timestep, calendar-month values, real observations never overwritten, substitution logged, and an unresolved gage reported rather than raising. test_diversion_config pins the three configuration shapes the report documents: nudging supplies observations, climatology alone supplies them in forecast mode, and both together. Neither source is not an error, it just means nothing is diverted, so it warns. Restore the Old River integration tests using fixtures rather than inline guards, and assert the behavior that defines the structure: turning the transfer on must lower the Mississippi peak below it and raise the Atchafalaya peak, and the climatology-only mode must still lower the Mississippi peak. Both are A/B comparisons against the no-diversion control, so they need no hardcoded bounds. Case inputs are still generated out of band by test.nhf.prep_tests, matching the other NHF cases. Relax the configuration guard for diversions. Requiring streamflow nudging would have rejected the historical-median-only run the report documents and the old_river case exercises: the climatological fill runs outside the nudging branch and writes the gage's row, so that mode does supply the transfer. Warn when neither source is configured instead of raising.
…path Place the diversion gage on the same routing link as the streamflow crosswalk. Reducing the streamflow crosswalk to one link per gage moved observations to the flowpath outlet, but the diversion site map still took the first joined sub-link, so on any multi-link flowpath the two disagreed, the kernel map could not find the gage's observations, and no flow was diverted. Both now use the same rule. Pass dt and nts into the BMI DataAssimilation. It was constructed with empty run parameters, so file-based nudging computed its resampling frequency from a null dt, and the climatological diversion fallback defaulted to dt=300 with nts=0 and built a single column that the kernel, indexing from timestep 1, never reads. Route the window-sliced observation frame in the BMI driver. The slice at the top of the loop was computed and then discarded in favour of the full frame, so every update restarted nudging and the donor subtraction from observation column zero. Fail with a clear error when nudging is enabled without a timeslice folder. The run died inside pathlib with a TypeError about NoneType, which named neither the setting nor the feature that needed it. Report how many routing links carry a gage. This count drives how finely the cached execution plan splits reaches, so a regression in it is a routing slowdown with no other visible cause; the new benchmark reads it back from the log. Add benchmark/bench_da.py, an A/B harness for the assimilation capabilities. Arms differ by exactly one configuration block and report wall time, memory and the gage link count. An arm whose feature cannot engage on the given domain, such as a diversion whose donor link is absent, is skipped with the reason rather than reporting a number that looks like a measurement. Stop the stale conftest from hiding the BMI tests, matching the reservoir DA branch.
… separate Rebuild the diversion's climatological frame at each forcing loop. With nudging off nothing else rebuilds it, so the frame stayed on the first loop's columns while the kernel restarts its timestep index at zero every loop. A run spanning several months therefore kept diverting the first month's climatology, which on the Old River case means diverting April's value through July, an over-subtraction from the Mississippi of up to 80 percent that nothing reported. Gate the reservoir DA shortcut on nudging rather than on the observation frame being non-empty. The climatological fill makes that frame non-empty even with nudging off, and it then holds only the diversion gage's rows, so USGS reservoir persistence derived its observations from a frame with no reservoirs in it and went silently dark. USACE and USBR were unaffected; they never take this shortcut. Add tests for the multi-loop case, including the reason the caller has to clear the frame: a fully populated frame has no gaps, so the fill alone cannot notice that the calendar advanced.
Validated end to end against real data: the 2011 spring flood, 2773 hourly NWM v3.0 retrospective forcing files and 11289 USGS timeslices over ten forcing windows. With observed diversion the Mississippi peak at Baton Rouge falls from 65794 to 44443 cms while the Atchafalaya at Simmesport rises from 5275 to 24314, which is the behavior the control structure exists to produce. The forecast-mode test compared the climatological run against the no-diversion run, but those differ in two ways: the latter also has streamflow nudging enabled. The nudged Mississippi is pulled down to its observed value while the unnudged one routes free on retrospective forcing that overestimates the peak, so the comparison reported the diverted run as higher and the diversion effect was invisible. Add a control with nudging off and no diversion, so the only difference is the diversion itself. Against it the climatological run lowers the Mississippi peak from 90942 to 85992 and raises the Atchafalaya from 5275 to 9235. Smaller than the observed case, as expected from monthly means standing in for a flood peak, and in the right direction.
…let nexus Only one observation can be assimilated per routing link, but the CONUS hydrofabric places 4096 USGS gages on 1729 shared virtual flowpaths. Which one survived was previously decided by dict insertion order, so re-exporting the geopackage with rows in a different order silently changed which gage drove assimilation at that link. Nothing in the gages layer separates co-located gages positionally: within every one of those 1729 groups their segment_order and dn_virtual_nex_id are identical, so they cannot be placed on distinct sub-links. Drainage area discriminates only 1069 groups and is null for half the colliding gages. Rank instead by an active gage over a discontinued one, then by the gage physically closest to the flowpath's downstream nexus, which is the one that best represents flow at the outlet link it is placed on, then by site number so the result cannot depend on row order at all. On CONUS that leaves 85 groups undecided after the first two and none after the third. Distance is computed only for colliding gages, so the geometry read covers roughly 1700 nexus points rather than the full 1.65 M, and a geometry failure degrades to status and site number with a warning. Re-sort after the per-flowpath outlet lookup. The groupby that selects the outlet link sorts by site number and discards the ranking, which without this made the deduplication keep the lowest site number and drop the active gages. The warning now names the sites that were dropped rather than only counting them.
The execution plan is a topological decomposition built once and reused for every forcing window, but its gage split points came from that window's observation frame. A first window that yielded no observations therefore produced a plan with no gage splits at all, and nothing rebuilds the plan when the feed recovers, so the whole run routed without splitting at gages and the in-kernel override could not land at the gaged nexus. Only a DEBUG line reported it. In normal operation this changes nothing. The observation frame is built by left-joining the gage crosswalk, so a gage with no data still gets an all-NaN row: measured on the Old River case the two sets are identical at 8 gages with 6 of the 8 rows entirely NaN, and CONUS splits at the same 21450 links either way. The sets diverge only when a window yields no observations whatsoever. Pass the network's gage set through to the plan and use it where present, falling back to the observation frame so callers that predate the argument are unaffected. The plan then depends on where gages are, which is static, rather than on which observations happened to arrive.
bench_da.py writes netCDF output and a timings json under benchmark/, neither of which belongs in the repository. The existing data/ rule does not cover a symlink named data, because a trailing slash matches directories only, and local benchmark datasets are commonly symlinked into place.
… supply The transfer conserves mass by construction, since the observed diversion is removed from the donor and the same value is imposed at the receiving system's headwater gage. The zero-flow clamp is the one exception: when the observed diversion exceeds the routed donor flow the donor is floored at zero and gives up less than the receiving river gains, and the difference is water the run created. Capping the transfer at the routed flow would fix it at the source but needs the donor's discharge inside the receiving reach's compute job, which crosses the parallel decomposition. Observed diversion has stayed below modelled Mississippi discharge across the overlapping record, so the condition is monitored rather than enforced. The previous check keyed on donor flow reaching zero, which reported a genuinely dry reach with nothing to divert as though water had been created, and gave a count of timesteps rather than anything about mass. It now fires only where a transfer was actually requested, and reports the volume involved as an upper bound, the donor's pre-clamp flow no longer being recoverable from the results. Document in the configuration that the full observed discharge is removed from the donor, which assumes the model routes no flow of its own down the diversion path. That holds while the diversion gage sits on a headwater flowpath, as at Old River; a gage with upstream contributing area would need the observed discharge minus the simulated flow already leaving the donor.
Applies the change from sclaw's commit 4b35b66, "Use DA at Vicksburg for historical median run", which gives CFG_HISTORICAL the same observation archive as CFG_NO_DIVERSION with streamflow nudging left on. The diversion gage is absent from that archive, so climatology still supplies it and the case measures what it is meant to, but every other gage is now assimilated and the run has the shape a real forecast has. Reapplied by hand rather than cherry-picked because this file was restructured to use fixtures in the meantime. This also makes CFG_NO_DIVERSION a matched control, since the only remaining difference between the two arms is the diversion itself, so the separate CFG_HISTORICAL_CONTROL case added earlier is redundant and is dropped along with its extra config and model run. Verified on the 2011 flood: against the CFG_NO_DIVERSION control the climatology driven transfer moves the Mississippi peak at Baton Rouge from 65794 to 60844 and the Atchafalaya at Simmesport from 5275 to 9235.
built_case cleared stale outputs but reused whatever config YAML was already on disk, and prep only writes that file when it is missing. Editing a Config in Python therefore left the previous settings live, and the run silently used settings the test no longer described. This is not a visible failure, it is a test that quietly measures the wrong thing. It cost a full diagnosis here: after CFG_HISTORICAL was changed to nudge the rest of the network, the case still ran with nudging off, so the forecast mode arm was compared against a nudged control, the Mississippi peak came out higher with the diversion enabled than without it, and the mass balance warning fired at every requested transfer. All of it was the old config. The config is cheap and deterministic to rebuild, so it is now written from the in-code Config on every run. The expensive inputs, domain and forcing and timeslices, are still reused.
4b35b66 to
e19b745
Compare
|
Issues found and fixed on this branch. Correctness
BMI path
Mass balance
Configuration and diagnostics
Performance
Tests and benchmarks
ValidationVerified end to end on the 2011 flood with real data: 2773 NWM retrospective forcing files and 11289 USGS timeslices. Peaks in cms at Baton Rouge on the Mississippi and Simmesport on the Atchafalaya:
Both arms move discharge in opposite directions at the two gages, which is the claim the report makes from this same case. CONUS gage placement goes from 10.70 links per gage to 1.00. Ohio streamflow DA goes from 0 to 318 gages. The execution plan change measures at plus 1.09 percent on CONUS, t(4) = 4.00. Full suite is unchanged against the base branch: same failure set, with 32 additional passing tests. |
The per-loop refresh reads reservoir_persistence_usgs directly off reservoir_da_parameters, one level above where it lives. The USACE branch a few lines down reads the same setting through reservoir_persistence_da, which is the correct nesting, so this branch could never fire. The giveaway is inside it: the reindex assigned to an undefined local rather than the instance frame, which would have raised NameError the first time it ran. It stayed hidden while the branch above refreshed the frame for most runs. Gating that branch on streamflow_nudging, which this PR does so the diversion's climatological fill cannot masquerade as gage observations, routes any run that assimilates reservoirs without gage nudging into this one instead, where _reservoir_usgs_df is never rebuilt. The kernel keeps recomputing observation offsets against the advancing t0, so from the second window onward type-2 reservoirs are handed observations dated before the window and quietly stop assimilating.
2 participants
Add Old River Control Structure support for NHF routing. Functionality is achieved through data assimilation alone (no physics-based model). By knowing divergence discharges a-priori, the current T-Route parallelization structure does not need to be modified.
Additions
test/nhf/old_river/run-troute-with-water-transfer.ipynbNew config file fields
The config file now accepts a new sub-dictionary for
data_assimilation_parameterswith keydiversion_da. Here is an example of the format.Here the
diversion_gage_crosswalkhas key value representing fp_id that water will be diverted from and value representing site_no for gage that tracks the divergence flowRemovals
None
Changes
None
Testing
All integration tests are passing. The test_conus_lakes.py test is passing. I have not run the regression check because 1) we expect different results and 2) the streamflow DA method that pre-existed in T-Route leads to substantial speed slowdowns.
Checklist