Non-hydrostatic solver rewrite, velocity-form momentum scheme, and moving-bed source by maartenvanormondt · Pull Request #352 · Deltares/SFINCS

maartenvanormondt · 2026-07-13T09:47:19Z

Non-hydrostatic solver rewrite, velocity-form momentum scheme, and moving-bed source

SUMMARY

This PR introduces a rewritten, GPU-ready non-hydrostatic pressure solver, a new
velocity-form momentum scheme, and a moving-seafloor source term for
landslide/earthquake tsunami generation. Together these make SFINCS suitable for
dispersive tsunami and short-wave modelling at significantly improved performance
and robustness. The default momentum scheme remains the flux form
(momentum_scheme = bates), so existing models are unaffected; the new scheme is
opt-in via momentum_scheme = velocity.

VELOCITY-FORM MOMENTUM SCHEME (sfincs_momentum_velocity.f90, new)

An alternative to the flux-form (Bates et al., 2010) momentum update, selected
with momentum_scheme = velocity. Velocity uv is the primary state; forcing is
assembled directly as acceleration, with the continuity flux derived as q = u*h.
Upwind surface/bed selection at the face avoids the average-bed depth overshoot
on steep slopes and improves run-up. Advection uses the momentum-conserving
formulation of Yamazaki, Kowalik & Cheung (2009) with second-order upwind
surface reconstruction, which is important for accurate bore and tsunami
propagation. Manning friction uses a fast cube-root evaluation (no pow calls,
fully vectorizable). Dual OpenMP/OpenACC directives throughout. Note for
reviewers: the advection limiter (advlim clamp and mask_adv gate) is disabled in
this scheme, so the full advective tendency is applied; this is intentional for
tsunami applications.

NON-HYDROSTATIC SOLVER REWRITE (sfincs_nonhydrostatic.f90)

The pressure projection (single-layer Keller-box) is rebuilt around an assembled
5-point stencil with a Chronopoulos-Gear single-reduction conjugate gradient,
replacing the previous BiCGSTAB-ILU approach. Symmetric Jacobi scaling is baked
into the operator (unit diagonal, no preconditioner work per iteration), a
ghost-element trick gives a branch-free coalescing matvec, and each CG iteration
runs as two fused kernels with a single host synchronisation. The solver runs
entirely device-resident on GPU and is equally effective with OpenMP on CPU.
Wave breaking is handled with either a steepness/HFA criterion (Smit et al.,
2013) or the NEOWAVE Froude criterion, with hysteresis, neighbour spreading,
hole filling and an optional reformation time scale. Robustness controls include
an open-boundary fade-in, global and localized 2dx pressure filters, bed-slope
caps, a depth-based pressure limiter and a tunable Keller-box dispersion factor.
Works on regular and quadtree meshes (refinement transitions stay hydrostatic)
and on projected and geographic grids.

MOVING-BED SOURCE (nonh_movingbed)

Adds d(zb)/dt to the bottom kinematic condition, fed through the new BMI/XMI
variable dzbext and the update_bed_level() entry point. This enables
landslide-generated tsunami simulation from an external bed-motion time series
or a coupled model.

NEW / CHANGED KEYWORDS

momentum_scheme (default bates): bates (flux form) or velocity
nonhydrostatic (default false): enable the NH solver (renamed from nonh)
nonh_movingbed (default false): moving-seafloor source via dzbext
nonh_tol (0.01), nonh_itermax (100): CG convergence
nonh_fnudge (1.0), nonh_tstop (-999): correction nudge / stop time
nonh_brsteep (0.4), nonh_brfr (0.5), nonh_brsmooth (0), nonh_treform (1.0):
breaking onset, smoothing and reformation
nonh_filter (0.5), nonh_smoothbnd (0.5), nonh_smoothdep (0), nonh_pmax (0):
2dx filters and pressure limiter
nonh_fadein (0), nonh_dzbmax (0.1), nonh_slsmooth (0), nonh_disp (1.0):
boundary fade-in, slope caps, dispersion factor

COMPATIBILITY NOTES

The keyword nonh has been renamed to nonhydrostatic: old input files using
nonh = 1 will silently no longer activate the solver. Consider a legacy alias
or startup warning.
Hydrostatic runs with the default (bates) scheme are unaffected.

TESTING

Validated against dispersive wave benchmarks (dispersion relation flat to
Airy up to kd of about 2.5 with nonh_disp = 1.0).
CPU (OpenMP) / GPU (OpenACC) consistency checked; CG iteration statistics
reported in the timing summary.
(add specific validation cases/figures here, e.g. solitary wave run-up,
landslide benchmark)

Generated with Claude Code (https://claude.com/claude-code)

commit 37b7969 Author: Maarten van Ormondt <maarten.vanormondt@deltares.nl> Date: Thu Jun 18 12:48:12 2026 +0200 Remove commented legacy upwind velocity form Remove an old, commented-out velocity-average single-sided upwind implementation from sfincs_momentum_velocity.f90. This cleans up dead reference code left next to the current blended upwind calculation (using vp/vn), making the file clearer; there is no functional change to the algorithm. commit b6d689d Author: Maarten van Ormondt <maarten.vanormondt@deltares.nl> Date: Thu Jun 18 12:34:38 2026 +0200 Remove legacy advection scheme options Remove handling and log messages for several advection scheme keywords (upw_div, sd03, mca, upw_div_upw, upw_flux, mca_sg). Keep default advection_scheme = 1 (upw1) and the 'original' option; unrecognized values still produce a warning. Cleans up deprecated/alternate scheme branches and related comments. commit d275452 Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Thu Jun 18 11:50:58 2026 +0200 Remove duplicate sfincs_data.f90 entry Remove duplicated <File> entry for ..\src\sfincs_data.f90 from source/sfincs_lib/sfincs_lib.vfproj to avoid duplicate inclusion in the project and potential build warnings. commit 0899dcc Author: Maarten van Ormondt <maarten.vanormondt@deltares.nl> Date: Wed Apr 29 13:55:11 2026 +0200 Support external dzbext and dynamic zb output Add support for an externally-supplied delta bed-level array (dzbext) and optional time-varying bed level output. Changes include: - Introduce logical use_dzbext and allocatable target dzbext in sfincs_data, and deallocate dzbext on shutdown. - Expose dzbext through the BMI getters (shape/type/rank) and lazily allocate/initialize dzbext when the BMI flag is enabled. - Replace the old update_zbuv binding with a BMI-facing bmi_update_bed_level (bound as update_bed_level) that calls a new update_bed_level routine. - Implement update_bed_level to apply dzbext to cell-centre or subgrid bed arrays and rebuild uv-derived quantities (zbuvmx), with separate logic for subgrid and non-subgrid modes. - Add support in NetCDF output to define zb either with or without a time dimension when store_dynamic_bed_level is enabled, and adjust writes so static writes are skipped when the variable is time-varying; per-step writer and map outputs now use subgrid_z_zmin for subgrid runs when writing dynamic bed level. Also includes minor whitespace/cleanup changes. These changes enable external drivers (via BMI) to modify bed level increments and optionally record time-varying bed topography in output files. (cherry picked from commit a87a934) commit bd513ea Author: Maarten van Ormondt <maarten.vanormondt@deltares.nl> Date: Thu Jun 18 11:29:32 2026 +0200 Normalize .vfproj File tags, add momentum file Replace multi-line empty <File>...</File> entries with self-closing <File .../> tags across source/sfincs_lib/sfincs_lib.vfproj (notably in snapwave group and several src/*.f90 entries) and add a new project entry for src/sfincs_momentum_velocity.f90. This is a project file cleanup (no source code changes) to keep the .vfproj consistent and reduce noise in diffs. commit 9832035 Merge: adc9089 30b5c0d Author: Maarten van Ormondt <maarten.vanormondt@deltares.nl> Date: Thu Jun 18 11:27:49 2026 +0200 Merge branch 'main' into sd03 commit adc9089 Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Thu Jun 18 08:54:50 2026 +0200 Add nonh breaking smoothing & flux-based advection Introduce configurable non-hydrostatic breaking/activity controls and a momentum-conserving flux-based cross-advection scheme. Changes include: - sfincs_input.f90: adjust several nonh defaults (nonh_fnudge, nonh_tol, nonh_filter, nonh_brsteep, nonh_brfr) and add new input parameters nonh_brsmooth, nonh_slsmooth, nonh_treform with comments. - sfincs_momentum_velocity.f90: add flux-based two-sided upwind cross-advection (Yamazaki eq.22) using fluxes (q0) and new vp/vn speeds; copy q->q0 each step and include q0 in OpenACC present lists; keep old velocity-average form commented for reference. - sfincs_nonhydrostatic.f90: add declarations and allocation for nonh_brsmooth, nonh_slsmooth, nonh_treform and new nh_bract array; implement optional frozen-bed-slope smoothing (nonh_slsmooth), gradual breaking reformation (nonh_treform), and optional smoothing of breaking ACTIVITY over neighbours (nonh_brsmooth). Use nh_bract (smoothed/activity) to scale pressure operator rows/columns and pnh contribution while keeping nh_brfac as the hysteresis STATE. Overall this commit smooths breaking/bed-slope transitions, enables gradual breaking recovery, and replaces the previous single-sided velocity cross-term with a momentum-conserving flux-based approach. commit b3844f1 Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Thu Jun 11 11:06:15 2026 +0200 Switch to assembled non-hydrostatic solver Refactor non-hydrostatic solver to an assembled-stencil, GPU-ready conjugate-gradient implementation and related API/initialization changes. Renamed input/config variables from nh_* to nonh_* and moved/removed related data (mask_nonh now in sfincs_nonhydrostatic). Removed legacy sfincs_nonhydrostatic_v2 and bicgstab_solver_ilu references from project/Makefile. Initialization now builds per-row/face assembled operator, frozen bed slopes use subgrid cell-mean bed for subgrid models, and the previous initialize_zb_effective helper was removed in favor of computing the cell-mean bed at domain init. Updated calls in sfincs_lib and sfincs_input to use the new init/compute routines and parameter names; memory layouts and device (OpenACC/OpenMP) residency prepared for per-step kernels. commit 17ae05b Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Thu Jun 11 10:24:14 2026 +0200 Add nonhydrostatic v2 (assembled-stencil CG) Introduce a GPU-ready assembled-stencil non-hydrostatic solver (sfincs_nonhydrostatic_v2.f90) and runtime selection between solver versions. Updates: add new source file and include it in the Visual Studio project and Makefile, add nh_version input to choose v1 (matrix-free CG) or v2 (assembled-stencil CG), and call the appropriate init/compute routines. Implement effective per-cell bed handling for subgrid models (zb_effective): new flag, initialization routine (initialize_zb_effective) and logic to recompute zb = zs - z_volume/area when needed; update continuity and momentum-velocity code paths to use the effective bed consistently. In the non-hydrostatic v1 module, precompute and freeze per-face capped bed slopes (nh_slbed) from the file bed to avoid feeding subgrid slope jitter into the bottom kinematic condition, and replace inline slope calculations with nh_slbed usage. Overall changes improve solver performance (reduced per-iteration memory traffic, Jacobi scaling, branch-free assembly) and add correct handling for subgrid bathymetry when using the velocity-form momentum scheme or non-hydrostatic solver. commit 950b51d Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Wed Jun 10 09:33:16 2026 +0200 Per-row inverse grid spacing for NH Introduce per-row inverse grid spacing arrays nh_dxr and nh_dyr and allocate them. Compute nh_dxr/nh_dyr per cell (handling crsgeo lat-dependent dx via dxm and quadtree refinement levels) and replace hard-coded dxrinv(1)/dyrinv(1) uses with these per-row values in bed-slope, kinematic (w_b), breaking divergence and time-derivative calculations. Update nh_cf face spacing to use dxminv for crsgeo or per-refinement dxrinv otherwise. Add explanatory comments about spacing choices to ensure correct operators across refinement and geographic grids. commit 95baa5c Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Wed Jun 10 09:00:36 2026 +0200 Remove adv_wiggle and conveyance options Remove the deprecated adv_wiggle parameter and the conveyance option/logic from the input and data modules. Deleted adv_wiggle variable and its read from sfincs_input.f90, removed convstr and conveyance_scheme declarations and the conveyance-selection block, and cleaned up related comments. This simplifies advection/conveyance handling by eliminating the NEOWAVE/upwind conveyance option and the TVD-style adv_wiggle limiter. commit 8864995 Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Wed Jun 10 07:38:47 2026 +0200 Update sfincs_domain.f90 commit 537ed14 Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Wed Jun 10 07:29:20 2026 +0200 Use upwind bed (zbup) instead of avg bed Replace zbavg with zbup and compute face wetness and flow depth from the upwind cell bed rather than the average of neighboring beds. zbup is chosen based on the sign of the previous-step velocity (uv0) with a deterministic tie-breaker when velocity is near zero. Update comments and the OpenMP private variable list accordingly. This avoids depth overshoot on steep downslopes while preserving run-up behavior. commit 66a3e9d Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Wed Jun 10 07:16:03 2026 +0200 Refactor momentum advection and subgrid handling Simplify and clean up momentum routines to improve subgrid support and reduce legacy advection complexity. Changes include: - Remove unused variables and legacy NEOWAVE temporaries. - Unify computation of zsu (u-point surface) and replace zwet/zbavg logic with a boolean iok wet/dry test; simplify hu calculation to use zbuvmx for conveyance. - Rework advection handling: remove/condense old advection branches and unify advection limiting to a single advlim clamp. - In velocity solver, introduce zbnm/zbnmu so Mader upwind reconstruction uses an effective bed for subgrid (or real bed for non-subgrid), ensuring consistent FLU computation. - Add wiggle-suppression limiter (mdrv) applied when subgrid and wiggle_suppression are enabled. - Miscellaneous formatting and minor performance-oriented reorderings. These changes aim to make the solver more robust with subgrid tables, simplify the advection logic, and reduce numerical wiggle artifacts. commit 0a64124 Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Tue Jun 9 23:22:13 2026 +0200 Make velocity momentum scheme default and refactor Switch the default momentum scheme from Bates (flux) to velocity form and update input parsing/logging accordingly. Reverse scheme enum semantics and branch to call compute_fluxes_velocity for the velocity scheme. Implement extensive refactor of the velocity-form momentum routine: treat pressure/viscosity/Coriolis/atm/wave as accelerations (remove multiplying hu), rename qfr->ufr, replace hu**(7/3) logic with a fast cube-root + hu^(4/3) factor (hu43) for Manning friction, and reuse the cube-root when wetting cells. Improve NEOWAVE advection stencil handling at refinement boundaries, adjust advective speed computations, apply advective acceleration directly to frc, and tighten wind/wave stress handling to use acceleration (stress/depth) with huvmin protection. Update OpenMP/ACC pragmas and private lists, remove the old power7over3 lookup function, and add small cleanups (variable renames, comments, and warnings/log messages). commit 9aa632a Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Tue Jun 9 07:48:46 2026 +0200 Switch momentum code to velocity form Remove the obsolete flux2nd input and variable and migrate momentum/velocity routines to the velocity form (uv) instead of carrying fluxes (q). Update input reading, OpenMP/ACC present/private lists, and loop copies to use uv/uv0. Rework friction and continuity computations: compute friction proxy qfr from velocities (hu*|u| or hu*sqrt(u^2+v^2)) and set q = uv*hu for continuity. Clean up many now-unused flux-related variables and simplify NEOWAVE-related comments and work variables. commit b87e1f3 Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Tue Jun 9 07:38:37 2026 +0200 Remove theta smoothing and qsm variable Remove legacy theta-smoothing logic and the associated qsm variable from sfincs_momentum_velocity.f90. qsm was deleted from the declarations and the OpenMP private list, and the entire theta-smoothing block (marked as not recommended) was removed along with a redundant commented adv line. Also simplified a VELOCITY-form comment. These changes clean up dead/unused code and do not alter the core velocity update logic. commit 58dac5e Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Mon Jun 8 18:11:35 2026 +0200 Add conveyance modes and subgrid-safe advection Introduce a conveyance_scheme option and a subgrid-safe upwind-flux advection mode. Changes include: - Add integer conveyance_scheme to sfincs_data and new input keyword 'conveyance' (convstr) parsed in sfincs_input.f90; default 'max' (Bates), 'upw' selects NEOWAVE-style upwind-surface + average-bed. - Add advection scheme 4 (advstr 'upw_flux' or 'mca_sg'): a subgrid-safe, momentum-conservative upwind-flux that upwinds the stored face fluxes (q0) instead of reconstructing depth zs - zb. - Update domain handling to allocate and populate zb from quadtree when reading subgrid files; avoid double allocation in nonhydrostatic branch. - Modify momentum routines to support conveyance modes: - New vars zbavg and zwet; decide surface at u-point (zsu) either as max(zs) (default) or upwind(zs) for conveyance=upw. - Determine wet/dry threshold zwet based on subgrid, average-bed (upwind conveyance) or max-bed (default). - Compute hu using average-bed depth for conveyance=upw, else use max-bed as before. - Implement advection_scheme==4 in momentum to upwind cell-centered mass flux from stored face fluxes and use a scheme-2 cross term (keeps momentum conservation and is subgrid-safe). - In velocity-form momentum (sfincs_momentum_velocity): - Extend OpenACC data list with cell_area, cell_area_m2, z_flags_iref and use z_volume/cell_area to compute cell-averaged depths when subgrid is active (avoids using zb there). - For non-subgrid, preserve upwind zeta reconstruction and max(zs - zb,0) behavior. - Adjust advection computation to use dnm/dnmu from either subgrid volumes or traditional zs-zb and update dnm+dnmu handling accordingly. - Simplify flux computation q(ip) to use uv*hu (commented-out the flux2nd branch). These changes enable NEOWAVE-style conveyance handling and a subgrid-compatible advection option to improve bore capture and correct behavior with subgrid tables. commit 2ba1783 Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Mon Jun 8 14:53:08 2026 +0200 Add NEOWAVE 2nd-order continuity flux option Introduce a new logical input 'flux2nd' (default .true.) and data flag to select Yamazaki/NEOWAVE continuity flux: true = 2nd-order central flux (eq.15), false = 1st-order upwind (eq.13). Changes: declare flux2nd in sfincs_data, read it in sfincs_input, and implement NEOWAVE upwind zeta logic, average bed (zbavg) and central flux depth (hufl) in sfincs_momentum_velocity. The patch updates the OpenMP private list, uses zbavg for flow-depth checks and hu computation, and computes q either as uv*hu (subgrid or upwind) or uv*0.5*(D_nm+D_nmu) for the 2nd-order option. Comments added to explain behavior and rationale (reduce O(dx) diffusion for propagating waves). commit 488eb6c Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Fri Jun 5 12:27:22 2026 +0200 Add NEOWAVE velocity-form momentum option Introduce an alternative velocity-form momentum solver (Yamazaki/NEOWAVE) and related options. - Add new module source/src/sfincs_momentum_velocity.f90 implementing compute_fluxes_velocity (velocity-form update, NEOWAVE-style advection and MCA option). - Register new source file in project and build (sfincs_lib.vfproj, Makefile.am). - Add momentum_scheme and nh_brfr parameters to sfincs_data and parse new 'scheme' (momstr) and optional 'nh_brfr' inputs in sfincs_input.f90; choose momentum solver at runtime. - Preserve initial velocities in sfincs_domain (don't zero uv) so velocity-form scheme gets correct startup state. - Wire up conditional dispatch in sfincs_lib.f90 to call compute_fluxes_velocity when momentum_scheme==1. - Extend advection options and momentum solver (sfincs_momentum.f90) with support for advection_scheme==3 (MCA / NEOWAVE-style) and refactor/comments; add Froude-based nonhydrostatic breaking criterion (nh_brfr) in sfincs_nonhydrostatic.f90. These changes enable using a velocity-based momentum formulation and NEOWAVE-compatible advection/breaking behavior while keeping the original Bates flux-form as default. commit aa31cda Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Thu Jun 4 17:50:05 2026 +0200 Nonhydrostatic: parallel gather CG & diagnostics Rework the non-hydrostatic solver to a race-free, gather-form, OpenMP-friendly implementation and add CG diagnostics. Introduces per-row incident-face map (nh_cellface) and per-face scratch nh_hf, plus nh_bb for the RHS and cg_w for the single-reduction (Chronopoulos–Gear) CG. apply_A_nh2 was converted to a two-pass gather operator (per-face flux then per-row sum). The CG was rewritten to run inside one parallel region with fused reductions to reduce synchronization, and iteration counters (nh_iter_total, nh_solve_count, nh_iter_max) were added and printed in the timing summary. Numerous loops were parallelized with OpenMP directives and neighbor-scatter patterns were converted to gather-style operations (filtering/smoothing, RHS assembly, diag computation) for thread safety and performance. commit b99a0a8 Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Thu Jun 4 12:44:31 2026 +0200 Simplify non-hydrostatic solver, remove sponges Remove legacy open-boundary sponge, dry-front, run-up taper, sparse matrix assembly and BiCGSTAB-ILU plumbing from the non-hydrostatic module and input. Many related variables/arrays and input parameters (nh_sponge_width, nh_sponge_coef, nh_relax, nh_dryfront, nh_runuph0/1, nh_brwidth, nh_brdilate, sponge ramps/diag, sparse matrix indices, etc.) were deleted and their allocation/usage stripped out. The code and comments now describe a matrix-free, SPD pressure-projection solved with a Jacobi-preconditioned conjugate-gradient (compute_nonhydrostatic2); boundary handling and breaking logic were simplified accordingly. Also changed default nh_disp to 1.0 and updated explanatory comments. These changes simplify the solver, remove non-parallelizable legacy pieces, and streamline the non-hydrostatic operator. commit 8b8ab37 Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Thu Jun 4 11:10:28 2026 +0200 Unify non-hydrostatic solver to matrix-free CG Remove legacy non-hydrostatic solver selection and related parameters, and make the matrix-free SPD conjugate-gradient solver the sole non-hydrostatic code path. Changes: dropped nonh_solver and removed input reads for nh_reformsteep, nh_subiter and nh_omega; simplified call site to always invoke compute_nonhydrostatic2; deleted the old compute_nonhydrostatic implementation; and cleaned up comments and allocations to remove solver-specific qualifiers. compute_nonhydrostatic2 now uses nh_disp for the Keller-box factor, uses 0.25*nh_brsteep as the reform threshold (XBeach default), and always runs the Jacobi-preconditioned CG solve (no separate Jacobi-only path). This simplifies control flow and removes the legacy BiCGSTAB-ILU and explicit pressure-relaxation branches. commit 14f059c Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Wed Jun 3 23:31:58 2026 +0200 Add adv_wiggle and enhanced nonhydrostatic Introduce a selective 2dx TVD-style adv_wiggle limiter and extend non-hydrostatic handling and solvers. Changes include: - Data & input: add adv_wiggle and new nonhydrostatic parameters (nh_brwidth, nh_disp, nh_pmax, nh_subiter, nh_brdilate, nh_omega) and input readers/comments for the expanded solver options. - Momentum: implement adv_wiggle limiter (selective numerical diffusion ~|q|*lap(u)) applied to dqxudx/dqyudy to damp grid-scale zigzags while preserving accuracy elsewhere. - Library: dispatch compute_nonhydrostatic2 for nonh_solver values 2..4. - Nonhydrostatic module: replace integer nh_breaking with real nh_brfac (smooth 1->0 breaking factor), add kbfac tunable vertical factor (nh_disp), refine breaking detection by smoothing dzdt, hysteretic state machine, hole-filling and optional dilation of breaking region, and keep face coupling active while enforcing pnh=0 as an internal Dirichlet constraint. Add solver-4 (weighted-Jacobi relaxation) as a reduction-free alternative, warm-start logic, and a depth-based pnh limiter (nh_pmax). Also adjust vertical source scaling and a few waveform/geometry fixes (hb/zb usage and hu calc). The goal is improved stability and physical behaviour for breaking fronts, reduced 2dx noise, tunable dispersion, and an alternative parallel-friendly pressure solver. commit 5c47134 Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Tue Jun 2 02:28:07 2026 +0200 Add matrix-free nonhydrostatic CG solver and options Introduce a new matrix-free, SPD conjugate-gradient non-hydrostatic solver (nonh_solver=2) and a suite of related options and data structures. Added many new input keywords (nonh_solver, nh_sponge_width, nh_sponge_coef, nh_relax, nh_filter, nh_dryfront, nh_runuph0/1, nh_dzbmax, nh_fadein, nh_brsteep, nh_reformsteep, nh_smoothbnd, nh_smoothdep) and the corresponding module variables and work arrays. Implemented compute_nonhydrostatic2 (Jacobi-preconditioned CG) with per-face coefficients, open-boundary sponge & fade-in ramps, hydrostatic-breaking switch, wet/dry-front handling, spatial filtering and localized smoothing. Kept legacy compute_nonhydrostatic and call it by default; select the new solver when nonh_solver==2. Also updated domain neighbor checks to allow open-boundary neighbors (kcs 2/3) and handle true 1-D grids, fixed boundary pnh handling to avoid out-of-bounds reads, added a new advection scheme keyword (upw_div) with alias sd03 and adjusted advection limiting behavior for the momentum-conservative scheme, and made a few small code cleanups/format fixes. commit 1244a89 Merge: 2948fc9 142cc4b Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Fri May 29 16:26:55 2026 +0200 Merge branch 'main' of https://github.com/Deltares/SFINCS commit 2948fc9 Merge: 196a287 83ddc38 Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Sun Apr 26 11:24:24 2026 +0200 Merge branch 'main' of https://github.com/Deltares/SFINCS commit 196a287 Merge: 595f914 fde0a93 Author: Maarten van Ormondt <82946105+maartenvanormondt@users.noreply.github.com> Date: Thu Feb 12 14:56:49 2026 +0100 Merge branch 'main' of https://github.com/Deltares/SFINCS commit 595f914 Merge: 0e88d61 82c8c4c Author: Maarten van Ormondt <maarten.vanormondt@deltares.nl> Date: Wed Jan 14 07:51:56 2026 +0100 Merge branch 'main' of https://github.com/Deltares/SFINCS commit 0e88d61 Merge: 18d1e91 2076235 Author: Maarten van Ormondt <maarten.vanormondt@deltares.nl> Date: Tue Dec 9 07:57:33 2025 +0100 Merge branch 'main' of https://github.com/Deltares/SFINCS commit 18d1e91 Author: Maarten van Ormondt <maarten.vanormondt@deltares.nl> Date: Fri Aug 1 12:49:48 2025 -0400 added build.log to gitignore

Introduce an optional moving-seafloor source for the nonhydrostatic solver: - Add a new logical input flag nonh_movingbed (read in sfincs_input.f90) and a module variable nonh_movingbed in sfincs_nonhydrostatic.f90. - When enabled, the code adds d(zb)/dt = dzbext/dt to the bottom kinematic wb (and subtracts it from ws) so a moving bed radiates a depth-filtered (Kajiura-like) surface response instead of being stamped onto zs. The update is applied in a separate OpenMP-parallel loop and is guarded by use_dzbext, dt>0, and a depth threshold so the original step remains byte-for-byte identical when the feature is off. - Adjust subgrid level bookkeeping in sfincs_domain.f90: translate subgrid_z_dep by dzbext for each level (clamped at -20 m) so volume->water-level tables remain consistent with bed shifts and partially-wet cells track the moving bed. - Add the debug Intel OpenMP binary libiomp5md.dll under source/sfincs/x64/Debug. These changes enable correctly modeling moving seafloor effects on nonhydrostatic surface dynamics while preserving existing behavior when the feature is disabled.

Split from tsunami-branch commit a345d07, where this change was bundled with the SnapWave 1-D/radstress work: the mask_adv gate and the advlim clamp on the advective acceleration are commented out so the full advection tendency is applied. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Copilot

Pull request overview

This PR upgrades SFINCS’ dispersive/short-wave modeling capabilities by adding a GPU-ready rewritten non-hydrostatic pressure projection solver, an opt-in velocity-form momentum scheme, and infrastructure for moving-bed forcing via BMI/XMI (dzbext) and time-varying bed output.

Changes:

Replaces the previous non-hydrostatic linear solver with an assembled-stencil CG-based projection designed for CPU/GPU execution, and adds moving-bed source support.
Adds a new velocity-form momentum update path and wires it into the time stepping via momentum_scheme.
Extends I/O and coupling hooks (NetCDF map output and BMI) to support non-hydrostatic pressure output and time-varying bed level / external bed motion.

Reviewed changes

Copilot reviewed 14 out of 15 changed files in this pull request and generated 7 comments.

Show a summary per file

File	Description
source/third_party_open/netcdf/netcdf-fortran-4.6.1/netcdff_c.vcxproj	Updates Visual Studio toolset selection for the bundled NetCDF Fortran project.
source/src/sfincs.f90	Initializes new BMI flag `use_dzbext`.
source/src/sfincs_output.f90	Ensures `pnh` is copied to host for map output in OpenACC builds.
source/src/sfincs_nonhydrostatic.f90	Major rewrite of the non-hydrostatic solver (assembled stencil + CG), breaking controls, GPU residency, and moving-bed forcing.
source/src/sfincs_ncoutput.F90	Adds/adjusts dynamic bed-level (`zb`) output shape and writing for regular/quadtree maps.
source/src/sfincs_momentum_velocity.f90	Introduces the new velocity-form momentum scheme implementation.
source/src/sfincs_lib.f90	Selects between momentum schemes; updates non-hydrostatic stop-time variable usage; adds CG iteration reporting.
source/src/sfincs_input.f90	Adds `momentum_scheme` parsing; renames/expands non-hydrostatic input parameters.
source/src/sfincs_domain.f90	Moves `mask_nonh` ownership to the non-hydrostatic module; adds bed update entry point and preserves initial velocities.
source/src/sfincs_data.f90	Adds `use_dzbext`, `dzbext`, `momentum_scheme`, and `zb_effective`; removes legacy `nh_*` parameters from this module.
source/src/sfincs_continuity.f90	Updates effective bed handling for subgrid runs when using velocity scheme and/or non-hydrostatics.
source/src/sfincs_bmi.f90	Exposes `dzbext` and adds a BMI entry point to apply bed updates.
source/src/Makefile.am	Adds new/updated sources for build system integration.
source/sfincs_lib/sfincs_lib.vfproj	Updates Visual Fortran project file entries to include new sources and remove old ones.
.gitignore	Ignores `source/build.log`.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>

Updates the inline comment in `sfincs_input.f90` describing momentum scheme options. This is a documentation-only wording change with no functional impact.

Temporarily disable writing of non-hydrostatic parameters (nonh, nonh_fnudge, nonh_tstop, nonh_tol, nonh_itermax) to NetCDF output. Also updates variable names from nh_ prefix to nonh_ prefix.

The merge of main (netcdf output redo #326) into qsrc lost pieces of the drainage output: his_type/map_type components, the his variable definitions, the river/urban-drainage his writes, the map-side cumulative urban drainage depth, the bucket-capacity map output, and left a duplicated early-return + implicit none in ncoutput_his_finalize and a duplicated toml-f filter in sfincs_lib.vfproj. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Introduce a local characteristic speed cap (`umax = sqrt(g*hu) + |u|`) for u- and v-direction advective velocities in `sfincs_momentum_velocity.f90`, and add `umax` to the OpenMP private list. This prevents `q/D` from exceeding the CFL-resolvable bound when subgrid cell-mean depth is small versus conveyance depth (e.g., steep valleys), improving stability. With speeds now bounded, the explicit advection acceleration limiter is disabled and the full computed advection term is applied.

mdrv was missing from the OpenMP private list of the main momentum loop: a correctness race (threads could read each other's wiggle damping driver, giving non-reproducible results), and severe false sharing - every thread rewriting the same shared scalar per wet uv point made the whole kernel ~1.7x slower than the bates scheme, independent of the advection setting. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

maartenvanormondt and others added 4 commits July 13, 2026 11:17

Update sfincs_input.f90

5f833ef

maartenvanormondt requested review from Leynse and Copilot July 13, 2026 09:47

Copilot started reviewing on behalf of maartenvanormondt July 13, 2026 09:47 View session

Copilot AI reviewed Jul 13, 2026

View reviewed changes

maartenvanormondt and others added 7 commits July 13, 2026 11:57

Potential fix for pull request finding

e147b11

Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>

Potential fix for pull request finding

f3c878c

Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>

Potential fix for pull request finding

49f47c7

Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>

Potential fix for pull request finding

13e1f0c

Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>

Clarify momentum scheme comment text

04e1e8e

Updates the inline comment in `sfincs_input.f90` describing momentum scheme options. This is a documentation-only wording change with no functional impact.

Merge branch 'main' into nonh_velocity_scheme

41b90f4

Comment out non-hydrostatic NetCDF attributes

03a2674

Temporarily disable writing of non-hydrostatic parameters (nonh, nonh_fnudge, nonh_tstop, nonh_tol, nonh_itermax) to NetCDF output. Also updates variable names from nh_ prefix to nonh_ prefix.

Leynse added this to the v2.4.1 milestone Jul 15, 2026

maartenvanormondt and others added 3 commits July 15, 2026 09:54

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Non-hydrostatic solver rewrite, velocity-form momentum scheme, and moving-bed source - #352

Non-hydrostatic solver rewrite, velocity-form momentum scheme, and moving-bed source#352
maartenvanormondt wants to merge 14 commits into
mainfrom
nonh_velocity_scheme

maartenvanormondt commented Jul 13, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Uh oh!

Conversation

maartenvanormondt commented Jul 13, 2026

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants