Fix/issue 1 junction surface area#4
Open
wiesnerfriedman wants to merge 193 commits intoSWMM-Project:mainfrom
Open
Fix/issue 1 junction surface area#4wiesnerfriedman wants to merge 193 commits intoSWMM-Project:mainfrom
wiesnerfriedman wants to merge 193 commits intoSWMM-Project:mainfrom
Conversation
Major routing fixes to align new engine with legacy SWMM behavior: - Add non-conduit flow callback inside DWSolver Picard iteration loop (pumps/orifices/weirs now computed per iteration, matching legacy) - Fix double-counting: updateNodeFlows skips non-conduits when callback active - Fix normal flow slope check (was reversed: h1>=h2 → y1<y2 matching legacy) - Remove EXTRAN slot width (legacy returns 0 for non-SLOT surcharge) - Add CrownCutoff depth capping for width computation in EXTRAN mode - Pre-build conduit_idx_ and nc_indices_ for O(conduits) inner loops Orifice/weir flow equations completely rewritten to match legacy: - Setting-adjusted effective opening (cOrif/cWeir recomputed dynamically) - BOTTOM vs SIDE orifice head computation with proper submergence - Villemonte submergence correction for weir-mode orifice flow - Weir setting raises effective crest height (legacy line 2257) - Orifice type (SIDE/BOTTOM), flap gate, close time now parsed - Weir type, flap gate, end contractions now parsed - Added param1 field to LinkData for structure sub-type storage Pump control fixes: - target_setting initialized from pump initial state (ON/OFF) - Pump setting = target_setting (matching legacy pump_getInflow) - Depth-based startup/shutoff hysteresis sets target_setting - Non-conduit callback only iterates non-conduit link indices Node initialization fix: - reset_state() now preserves init_depth (was zeroing all depths) - Initial volumes computed from init_depth after reset Report unit conversion fix: - All flow statistics (link max flow, node inflow, pump flow, outfall flow, flood rate) now converted from CFS to display units using Qcf factor Test constant renames: DEFAULT_HEAD_TOL, DEFAULT_MAX_TRIALS, MIN_TIMESTEP Result: Flow routing continuity error improved from -2983% to -6.1% (legacy: +1.0%). Remaining gap traced to 35 IRREGULAR cross-sections not yet supported in XSectBatch — blocks flow through trunk sewers.
Signed-off-by: cbuahin <caleb.buahin@gmail.com>
IRREGULAR cross-section support: - Per-link transect table lookup in XSectBatch (area/width/hrad) - New perlink_tabulated() kernel for shapes with individual tables - Transect name → index resolution deferred to PostParseResolver (handles XSECTIONS appearing before TRANSECTS in input files) - TransectData tables built from station-elevation data during init - ShapeGroup now holds per-link table pointers for IRREGULAR group - transect_tables stored in SimulationContext for lifetime management Performance optimizations: - Pre-allocated buf_d/buf_r in ShapeGroup (no allocation in hot path) - Removed thread_local vector allocations in DW width-capping - Changed OpenMP schedule from dynamic(64) to static for momentum Result: Flow routing continuity error +0.509% (legacy: +1.028%). Node 2883 (IRREGULAR conduit): 0.32/1.19 ft (legacy: 0.19/1.09 ft). Runtime: 475s (legacy: 302s). Remaining gap: WWTP pump activation.
CUSTOM cross-section support: - Shape curve name resolution deferred to PostParseResolver - buildCustomTables() computes area/width/hrad tables from depth-width shape curves (integrated area, computed perimeter) - CUSTOM tables stored as TransectData alongside IRREGULAR tables - XSectBatch dispatches CUSTOM shapes to perlink_tabulated kernel - SWR00362559_2 now matches legacy: a_full=95.37 (legacy 95.43) XSectBatch pre-allocated buffers: - ShapeGroup buf_d/buf_r pre-allocated in resize() for hot loop - computeAreas/HydRad/Widths use pre-allocated buffers instead of per-call std::vector allocation Note: CSO outfall flows still too high — orifice/weir equations at regulators need further investigation.
The DW momentum solver handles adverse slopes naturally through the St. Venant equations. Removing the reversal keeps conduit node ordering consistent with the input file and matches how legacy conduits like M04 (ARCH_BASCULE → REG_DWF_06) appear in reports. Note: -42% continuity error persists after CUSTOM xsect support. The mass balance violation traces to water being created in routing, likely from CUSTOM conduit geometry interaction with the DW solver. Requires investigation of the buildCustomTables() area/hrad tables.
Signed-off-by: cbuahin <caleb.buahin@gmail.com>
Critical fix: translateShape() used simple +1 offset which misaligned LinkData::XsectShape enums with XSectBatch::XSectShape enums beyond the first 8 shapes. IRREGULAR(18) mapped to 19 (SEMIELLIPTICAL) instead of 22, CUSTOM(19) to 20 (BASKETHANDLE) instead of 23, etc. Replaced with explicit switch mapping for all 22 shape types. Also simplified CUSTOM pre-computation in PostParseResolver to let buildCustomTables() handle all geometry (removed redundant inline area/perimeter integration that used wrong width scaling). Result: Continuity -42% → -26%, flooding 397 → 192 acre-ft. Runtime: 428s (down from 475s due to correct batch kernels).
Signed-off-by: cbuahin <caleb.buahin@gmail.com>
Signed-off-by: cbuahin <caleb.buahin@gmail.com>
Signed-off-by: cbuahin <caleb.buahin@gmail.com>
Signed-off-by: cbuahin <caleb.buahin@gmail.com>
Signed-off-by: cbuahin <caleb.buahin@gmail.com>
Signed-off-by: cbuahin <caleb.buahin@gmail.com>
Signed-off-by: cbuahin <caleb.buahin@gmail.com>
EXTRAN's dQ/dH surcharge formulation is discontinuous and not a smooth fixed-point map G(y). Anderson acceleration mixing can produce depths below the crown floor that setNodeDepth enforces, leading to non-physical results. This guard skips AA mixing when a node is surcharged under EXTRAN, falling back to the standard Picard iterate. AA remains active for SLOT and DYNAMIC_SLOT methods which produce smooth G(y) maps.
…uard Skip Anderson acceleration for surcharged nodes under EXTRAN
Only strip conduit/orifice half-areas at a STORAGE node when the storage curve provides area > MIN_SURFAREA. For degenerate storage nodes (zero or near-zero curve), keep the legacy pipe-half contribution so the Picard denominator stays bounded away from MIN_SURFAREA. Fixes: DynamicWave.cpp::updateNodeFlows (conduit path) Fixes: SWMMEngine.cpp non-conduit callback (orifice path) Tests: 5 new StorageHalfAreaGuard unit tests
… regimes
Anderson acceleration (AA) requires a smooth fixed-point operator G
between consecutive Picard iterates. Three surcharge formulations
violate this assumption:
EXTRAN: discontinuous dQ/dH formulation at crown
DYNAMIC_SLOT: per-iterate geometry rewrite (Sharior et al. 2023)
SLOT: C⁰ kink at the slot cutoff (~0.985·yFull)
The existing code only skipped AA for EXTRAN surcharged nodes (from
pr/aa-extran-guard). This extends the guard to all three regimes via
a per-node aa_skip_ flag vector, scatter-computed each Picard iteration
after geometry is known:
- EXTRAN: skip when xnode_.is_surcharged
- DPS: skip end-nodes of conduits with active slot area (As > 0)
- SLOT: skip end-nodes of conduits with depth_mid/yFull in [0.98, 1.02]
Free-surface nodes (95%+ of network) retain full AA speedup.
No physics is altered; only the numerical accelerator is gated.
Implementation:
- DynamicWave.hpp: add aa_skip_ (vector<uint8_t>), computeAASkipFlags()
declaration, aaSkipFlags() const accessor
- DynamicWave.cpp: allocate aa_skip_ in init(), implement
computeAASkipFlags() with O(n_conduits) scatter, call after
computeLinkGeometry() in execute(), guard AA block in
updateNodeDepths() with !aa_skip_[ui]
Tests: 6 new AASkipFlagTest cases in test_routing.cpp
- FreeSurfaceNoSkip, ExtranSurchargedSkips, DPSActiveSkipsEndNodes,
SlotNearKinkSkipsEndNodes, SlotFarFromKinkNoSkip, AADisabledNoFlags
All 61 routing tests pass.
…orage-conduit-halves fix(issue-2): conditionally zero storage-node conduit half-areas
…derson-acceleration-scope fix(issue-3): extend Anderson acceleration skip guard to DPS and SLOT…
…Node, wire option into DWSolver with unit conversion, add .inp-driven regression tests
Phase 1: Non-storage getSurfArea() returns 0.0 (physical area only)
Phase 2: Storage getSurfArea() purely geometric (no MIN_SURFAREA clamp)
Wire MIN_SURFAREA/HEAD_TOL options into DWSolver with UCF conversion
Phase 3: Add minimal_conduit.inp fixture and DWNodeContinuity tests
(DepthRisesWithForcedLateralInflow, MultipleStepsAccumulateDepth)
All 61 routing + 21 site-drainage tests pass.
wiesnerfriedman
force-pushed
the
fix/issue-1-junction-surface-area
branch
from
dac366e to
8f75098
Compare
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Node::getSurfAreareturnsMIN_SURFAREAas a constant for non-storage nodes, which gets added to conduit half-areas instead of acting as a floor. This inflates junction surface area, damping depth response. Fix moves the floor toDWSolver::setNodeDepthasmax(surfArea, MIN_SURFAREA).