Fix culvert discharge chatter, make drain volume transfer mass-exact, limit withdrawal per step

[keesnederhoff]

Fix culvert discharge chatter, make drain volume transfer mass-exact, limit withdrawal per step

Tim reported wiggles in the culvert discharges when running the refreshed testbed after the Galibier release. In the Newark drainage model, the discharge on volume-limited culverts (drains 108 and 15) oscillates full-amplitude (0 to -16 m3/s) across the whole event, while drains that are not volume-limited (e.g. 10, 100) stay smooth. Bisected across the release archive: v2.2.0 is smooth, v2.3.0 (#211, moving gates) swapped the order of the volume clamp and the relaxation in the shared drainage tail, so the on/off clamp output is applied directly and fed back as the relaxation memory.

Changes

1. sfincs_discharges.f90: clamp first, relax last, so the applied discharge is a smooth low-pass again (restores the v2.0-v2.2 order).
2. sfincs_continuity.f90: drainage structures are applied as pairs. The donor cell is floored at zero volume (subgrid) or at the bed (regular), and the receiving cell gets exactly the volume that was actually removed. v2.2 delivered the full discharge to the outfall even when the intake had run dry.
3. sfincs_discharges.f90: the volume limit is now a fraction of the donor cell volume per step, new keyword drainage_volfrac (default 0.5, range 0.05-1, 1 recovers the old hard limit). With the hard limit q <= V/dt a volume-limited culvert emptied its cell in one step and the limited discharge switched between zero and full demand; with a fraction the cell volume decays geometrically and the discharge converges smoothly to the inflow.
4. sfincs_continuity.f90: plain point sources are applied and floored (previously silently dropped when the cell volume was negative).
5. New end-of-run log line: Drain/source withdrawal: X m3, not removed (dry cell): Y m3 (Z%), the volume drains nominally withdrew and the part the dry donor never gave up.

Newark testbed, 8 h event (nTV = total variation / 2 range, about 1 for smooth, about 40 for chatter; reversals = direction changes > 5 % of range; jump = largest step change / range)

variant	drain 108 nTV / rev. / jump	drain 15	drain 100	water created	not delivered
main v2.4.1	41.3 / 187 / -	39.2 / 130 / -	4.2 / 0 / -	4,570 m3 (0.16 %)	-
clamp-then-relax only (first push of this PR)	0.87 / 0 / 0.01	0.70 / 0 / 0.01	3.9 / 0 / 0.02	66,392 m3 (2.20 %)	-
relax unclamped, clamp last (review suggestion)	1.09 / 0, spikes	1.18 / 0, spikes	6.98 / 11	75,541 m3 (2.50 %)	-
pair delivery, drainage_volfrac = 1	0.87 / 0 / 0.01	1.20 / 0 / 0.02	4.2 / 0 / 0.03	0	3.56 %
this PR: pair delivery, drainage_volfrac = 0.5	0.89 / 0 / 0.01	0.87 / 0 / 0.01	4.2 / 0 / 0.02	0	0.31 %
pair delivery, drainage_volfrac = 0.25	0.90 / 0 / 0.01	0.90 / 0 / 0.01	4.0 / 0 / 0.02	0	0.002 %

drainage_volfrac = 0.5 also removes the sharp dip and the deep late peak that = 1 still shows on drain 15. Flood map: zsmax at 0.5 is within 8 mm of the value at 1 in every cell; at 0.25 the limiter starts to throttle (11 cells change by more than 5 cm, up to 0.31 m), which is why 0.5 is the default.

At the SFINCS default structure_relax = 10 (Newark sets 600) the volume-limited drains stay fixed (drain 108 nTV 33 on main, 1.8 here; drain 15 15 on main, 1.3 here; undelivered 3.7 % at drainage_volfrac = 1, 0.07 % at 0.5). The small control drains 100 and 10 chatter at relax 10 on main and on this branch alike (drain 100 nTV 66 vs 64, about 570 reversals in both), so that is a separate, pre-existing weakness of the 10 s default against a 1.8 s time step and is not addressed here.

Behaviour change to be aware of. Because the outfall no longer receives phantom water, zsmax drops by up to 1.1 m in 173 cells around outfalls relative to the first push of this PR, and relative to v2.2 which had the same phantom water. Regular-grid gate case 01_Implementation/controlled_gates_group/controlled_gates: water levels at the observation points agree to 1 mm with the first push, and the gate intake never dries; the only discharge difference is a phase shift of the +-0.5 m3/s flicker both versions show while the gate sits at near-zero head.

[Leynse]

Nice fix Kees!
2 questions;

1. Did you think about implications for mass conservation?
2. And possible impact on other point sources as disfile etc, that are also affected by this?

Copilot review - One real concern:

1. Reordering relax after the clamp breaks the volume limit (and mass conservation)

sfincs_discharges.f90:626-645 — the clamp now applies to the target discharge, but the value actually written to qtsrc is a blend with the previous (possibly much larger) discharge. So the applied withdrawal can exceed the water available in nmin. The new max(..., 0.0) floor in continuity absorbs that on the donor side, but qtsrc(jout) = qq still delivers the full discharge to the receiving cell → net water creation whenever the drain runs a cell dry. Under the old order this couldn't happen.

Suggested alternative that fixes chatter and keeps the limiter authoritative: the chatter comes from relaxing toward the clamped previous discharge (qq0 = -qtsrc(jin) at :440/:540 feeds the clamped value back into the Bates update). Keep an unclamped qq_prev per structure, relax against that, then clamp last:

qq = w*qq_target + (1.0 - w)*qq_prev(idrn)   ! smooth, unclamped
qq_prev(idrn) = qq
... volume clamp ...                          ! authoritative, applied last
qtsrc(jin) = -qq ; qtsrc(jout) = qq

2. Scope of the continuity change is wider than drainage structures

sfincs_continuity.f90:319 — the loop runs over nsrcdrn, i.e. all point sources, not just drains. Previously a positive source into a cell with z_volume < 0 was silently dropped; now it is applied. That is arguably the correct fix, but it's a behaviour change for plain src points that isn't mentioned in the PR title/description — worth calling out in the release notes. Also, momentum:697, infiltration:649 and meteo:1321 all explicitly test for negative z_volume, so silently flooring it at source cells hides mass errors that were previously visible.

[keesnederhoff]

Thanks Tim, both points were worth checking.

1. Mass conservation: you are right, and the suggested alternative does not fix it.

variant	drain 108 nTV / reversals	drain 15	drain 100 (control)	water created
main v2.4.1 (relax, then clamp)	41.3 / 187	39.2 / 130	4.2 / 0	4,570 m3 (0.16 %)
this PR as first pushed (clamp, then relax)	0.87 / 0	0.70 / 0	3.9 / 0	66,392 m3 (2.20 %)
Copilot suggestion (relax unclamped, clamp last)	1.09 / 0, spikes	1.18 / 0, spikes	6.98 / 11	75,541 m3 (2.50 %)
pair delivery, drainage_volfrac = 1	0.87 / 0	1.20 / 0	4.2 / 0	0
pair delivery, drainage_volfrac = 0.5 (now in the PR)	0.89 / 0	0.87 / 0	4.2 / 0	0

nTV is total variation over twice the range (about 1 for a smooth hydrograph, about 40 for full-amplitude chatter); reversals counts direction changes larger than 5 % of the range.

• The PR as first pushed creates 2.2 % of the drained volume:
• The "clamp last" design still creates 2.5 %.
• What does work, and is now pushed, is two things. First, continuity treats each drain as a pair: the source cell is floored, and the receiver gets exactly the volume that was actually removed, so the pair conserves by construction. Second, the volume limit is a fraction of the source cell volume per step (drainage_volfrac, default 0.5). The old hard limit let a culvert empty its cell in one step, so the available volume flipped between zero and full
• On your default-relaxation hypothesis: Newark sets structure_relax = 600, so I also ran everything at the default 10 s. The volume-limited drains stay fixed there.

2. Scope of the continuity change. Correct, the floor also applies to plain src points: a source into a cell with negative volume used to be dropped and is now applied and floored. I kept that, because dropping a prescribed inflow was itself a silent mass error, but it is now in the PR description, and the new log line makes it visible in every run.