Skip to content

Nswing - A finite differences tsunami modeling supplement#9093

Open
joa-quim wants to merge 30 commits into
masterfrom
nswing
Open

Nswing - A finite differences tsunami modeling supplement#9093
joa-quim wants to merge 30 commits into
masterfrom
nswing

Conversation

@joa-quim

@joa-quim joa-quim commented Jul 7, 2026

Copy link
Copy Markdown
Member

NSWING is a non linear shallow water with nested grids tsunami modeling code. I have it for a long time compiled as a MEX file, which is used in Mirone. Now, whit help of Fable, Opus and Sonnet (credits were not enough to satisfy the hungry models) I converted it into a GMT supplement. In the transition steps Fable rewrote some parts using OpenMP and it now runs twice fast as it used to.

It's not easy to make tests for it, and the modeling result is stored in a 3D netcdf cube. GMT can be used to visualize it but a bit cumbersome. For the time being, Mirone is still the best tool I know to visualize data like this but that will, hopefully (so Claude credits help) change soon too.

The man page has a quick example of what we can get without the nesting grids and it's something like this.
Level_=_597

I'm pasting too the result of what is known as Benchmark 1 of the Catalina set of benchmarks that compares against the analytical solution of a wave inundating a beach with 10% slope. Results are excellent.

GMTplot

joa-quim and others added 30 commits June 25, 2026 18:28
- Exclude -R from the rebuilt legacy argv[] (already consumed by
  GMT_Parse_Common into GMT->common.R): fixes "Unknown option -R..."
- resamplegrid()/GMT_get_bcr_z() now skip land source nodes in the
  bicubic stencil, renormalizing over wet nodes only. Coastal terrain
  elevations stored in etad on dry nodes were bleeding into nearby
  wet-node interpolates as invented negative sea-surface values,
  triggered only via -J's (now -P's) nested-grid hold/resync path.
- Rename time-jump option -J to -P: GMT_Create_Options() unconditionally
  prefixes '+' to any -J argument starting with a digit (EPSG-code
  shorthand support in gmt_parse.c), which mangled a numeric
  -J<time_jump> into -J+<time_jump> before nswing's own parser ever
  saw it, making run_jump_time silently pick up time_jump's value.
- Remove the old -O<int>,<outfname> (unused maregraph-interval) option;
  rename -B<BCfile> to -O<BCfile>.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Replace the MEX-era argv[] rebuild + switch(argv[i][1]) dispatch with
the standard GMT pattern: iterate options directly and switch(opt->option),
reading opt->arg. Drops parse()'s argc_out/argv_out signature (now matches
other GMT modules' parse(GMT, Ctrl, options)).

-Fk's region decode no longer mutates opt->arg in place with a
strdup/strcat restore hack; uses a small local scratch buffer instead.
The netCDF History attribute is now built via GMT_Create_Cmd()/
GMT_Destroy_Cmd() instead of hand-joining a rebuilt argv[].

Verified: clean build, basic run + -Fk region/RxC parsing produce
identical output/History attribute to before. A crash in the
multi-prism (-Fk.../RxC) simulation loop was found but confirmed
pre-existing (reproduces on the prior commit too, unrelated to this
change).

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
-Fk<region>/RxC or /dx/dy forces out_maregs_nc since the whole point of
a prism grid is sampling a maregraph per prism, but nothing checked that
-T (which sizes/allocates the maregraph arrays via n_mareg/n_ptmar) was
actually given. Without it those stayed at their uninitialized stack
values, calloc'd a garbage-sized buffer, and corrupted memory once real
maregraph data got written into it (previously misread as a crash in the
kaba goto-loop itself; bisected against the prior commit and confirmed
it wasn't the parse() refactor).

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
nswing copies every grid into its own flat arrays and handles boundary
conditions itself; GMT's default grid pad is dead weight here, and the
user wants it gone outright rather than relying on gmt_M_ijp() to always
index around it. GMT_Set_Default(API, "API_PAD", "0") applied once right
after GMT_Parse_Common, before any grid is read or created.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
The parse() rewrite to switch(opt->option) dropped the -R exclusion that
used to live in the old argv-rebuild loop, so -R fell through to
default: and errored as "Unknown option". Restored as a continue filter
before the switch (not a case label) -R is already consumed by
GMT_Parse_Common into GMT->common.R and isn't one of nswing's own
options.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
When no bathy grid is given but a source (-F or -Fk) and -R are, there's
nothing to simulate -- skip the whole time-stepping machinery and just
compute the deformation directly over -R's grid geometry (region +
increment + registration, populated by GMT_Parse_Common when -R names a
grid file) and write it to -G's name. -t and -G's saving interval are no
longer required in this case.

Verified: `-Rbathy.grd -F<9 params> -Glixa_def.grd` now writes a grid
matching bathy.grd's exact geometry with sane deformation values, no
simulation run; the normal simulation path (bathy given) still requires
-t as before.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
The big one: mass() had 'for (row = row_start; row_end; row++)' -- the
loop condition was the constant row_end, so the row loop ran past the
array end and every simulation segfaulted in the first
mass_conservation() call since the Mirone conversion (2380d50).
Only the deformation-only mode worked.

Also fixed, from a cppcheck + manual audit:
- -Em always aborted: '||/&&' precedence bug skipped the workMax calloc
  and fell into the out-of-memory branch.
- work array sized with a comma expression '(a, b)' instead of MAX(a, b).
- sanitize_nestContainer typo left fluxn_d[] uninitialized (fluxm_d was
  NULLed twice).
- '!stem' (array, always false) dead check replaced by a scoped test;
  guarded the strlen(stem)-1 underflow when stem is empty.
- always-true '(nest.long_beach || nest.short_beach)' array-decay test
  replaced by the do_long_beach/do_short_beach flags.
- initialize_nestum() read nest->hdr[lev-1] before the lev > 0 guard.
- -X Manning parsing rewritten with strtok, bounded to the 10 levels the
  struct holds (was an unbounded manning[k++] overflow, plus a hidden
  t[16] overflow).
- parse() hardening: field widths on all sscanf %s, bounded strncpy for
  -S, bounds guard on the -G '.nc' extension append, short-arg guard in
  -L, Okada params zero-initialized, deformation-only mode now requires
  -G explicitly.
- error-path cleanups (matter when called repeatedly from a persistent
  GMT.jl session): destroy read grids on the error bail, only open the
  ASCII maregraph file when not writing netCDF, free mareg_names and the
  tracers container on their failure paths, close fp_oranges, call
  free_arrays if initialize_nestum fails, guard the netCDF maregs end
  block when no maregraphs were found, require a bathy grid of at least
  2 x 2 nodes.
- check_paternity(): loop-safety only (k++ and k < 10 bound). The check
  itself stays inert: reviving it rejected valid nested grids because
  check_binning()'s math predates the current half-cell nesting rules.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Replace struct srf_header (a relic of the old Surfer-binary I/O) with
the generic struct grd_header everywhere: gmtnswing_get_grid(),
deform(), kaba_source(), the hdr_b/hdr_f/hdr_mM/hdr_mN headers and the
nested-grid read block. Grids read through the API are no longer capped
at 32767 rows/columns by the short int nx/ny fields, and the (short)
casts in deformation-only mode are gone. gmtnswing_get_grid() now also
fills nm and x_inc/y_inc straight from the GMT header, and the memset
in kaba_source() computes its size in size_t.

Deformation output verified bit-identical to the previous build; plain
simulation and a real 2-level nested run pass.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…rflow

Drop the home-grown struct grd_header and use GMT's own GMT_GRID_HEADER
throughout (374 field references: nx/ny -> n_columns/n_rows, x_min etc
-> wesn[], x_inc/y_inc -> inc[], plus the ij_grd macro).
gmtnswing_get_grid() now just struct-copies *G->header. The write-only
hdr.lat_min4Coriolis copies are gone (the nestContainer keeps its own).

Also fix a long-standing heap overflow found while validating this:
open_most_nc() allocated basename with malloc(strlen(base)) - no room
for the NUL terminator, and the MOST variants strcat a further "_ha.nc"
type suffix onto it. Every 3D/MOST netCDF run wrote past the block; the
header refactor reshuffled the heap enough for the CRT to start
aborting at teardown on some -G name lengths. Now strlen(base) + 8.

Validation: single-threaded old vs new outputs are bit-identical (the
default multithreaded runs cannot be compared - moment_M/moment_N race
makes results nondeterministic run-to-run, a separate pre-existing
issue). Deformation-only output bit-identical; nested 2-level case
passes 8/8.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
-T<x/y>[+o<outmaregs>][+t<int>] is now equivalent to giving a one-point
positions file: the token before the modifiers is taken as coordinates
when it parses as exactly two /-separated numbers, otherwise it is a
file name as before. The node index is computed with the same
irint((x - wesn)/inc) math and inside-grid check that read_maregs()
uses. Useful for quick single-gauge runs without touching disk.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The #ifndef PARALLEL prologue of the cartesian moment_N() did
memset(fluxm_d, ...) - moment_M's output array - instead of its own
fluxn_d. Under DO_MULTI_THREAD moment_M and moment_N run concurrently,
so the wipe landed at a random point while moment_M was filling
fluxm_d: cartesian results were nondeterministic run to run (up to
~0.27 m differences at a gauge). In a serial build the wipe simply
zeroed the x-fluxes each step. The three sibling functions (moment_M,
moment_sp_M, moment_sp_N) and the PARALLEL-mode hoisted memsets all use
the right array, which is why geographic (-fg) runs were always
deterministic and unaffected.

Verified after the one-word fix: identical cartesian commands are now
bit-identical run to run, a serial build matches the threaded one, and
geographic output is unchanged.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Drop both the DO_MULTI_THREAD path (2 Windows threads created and
destroyed every time step, moment only) and the never-enabled PARALLEL
path (per-step slice threads whose integer-division row splits silently
skipped up to n_threads-1 rows). In their place, OpenMP parallel-for
with schedule(dynamic, 8) on the six row loops: mass, mass_sp,
moment_M, moment_N, moment_sp_M and moment_sp_N. mass is now
parallelized too, all cores participate, and threads persist in the
OpenMP pool instead of being respawned each step. About 330 lines of
scaffolding removed: MT_* wrappers, *_slice function variants and their
mid-function #ifdef body-sharing tricks, ThreadArg, the PFV function
pointers. -x<n> now maps to omp_set_num_threads() [default all cores].

dynamic scheduling matters on hybrid P/E-core CPUs, where static splits
made every per-step barrier wait for the slowest efficiency core and
ran slower than the old 2-thread code. Results are bit-identical to the
old build and across schedules and thread counts (the loops have no
reductions), so runs are finally reproducible regardless of -x.

The four moment loops also zero dd/df per cell now: the goto L121/L201
shortcuts could read them stale from the previous cell (a latent NaN
source for -S velocities, and a loop-carried dependency the parallel
loops must not have).

Timings (20-core machine, maregraphs-only runs): 765x476 geographic,
1010 steps: 5.8-7.3 s old -> 3.5 s; 1201x1201 cartesian, 300 steps:
12.0 s old -> 4.4 s. Non-Windows builds gain parallelism for the first
time (old threading was Windows-only).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Two-stage detection, applied in both the simulation path and the
deformation-only mode: first trust GMT's own grid metadata
(gmt_M_is_geographic, set by the netCDF reader when the grid carries
degree units), then fall back to the old -180/360, -90/90 range
heuristic for unit-less grids (e.g. Mirone-written ones). Either way a
warning is issued and isGeog is set.

The deformation-only mode previously had NO detection at all - it
returns before the bathymetry-reading branch where the heuristic lived,
so 'nswing -R<grid> -F... -G<name>' on a lon/lat region silently ran
the Okada deformation in Cartesian, producing wrong results with no
warning. Verified against an explicit -fg run: bit-identical now.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
… 3D cube

Previously -E's Energy/Power silently overwrote the z (sea-surface) variable
in the -G 3D netCDF cube, mislabeled with z's own name/units ("Sea surface
[meters]"), while -S always wrote velocity alongside z with no way to opt
out. Add +a to both: absent, the option's quantity is the cube's only
variable (correctly named/labeled); with +a, it's written as an extra
channel next to z, matching -S's old default behavior.

Note: -S without +a now omits the z variable (previously always included).

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
@joa-quim joa-quim requested a review from a team July 7, 2026 22:11
@Esteban82 Esteban82 added add-changelog Add PR to the changelog AI-assisted All (or most) of the code was written by Artificial Intelligence. labels Jul 8, 2026
@Esteban82

Copy link
Copy Markdown
Member

Is there any way to check this? At least to get the same figures as the ones you’ve shown?

@Esteban82 Esteban82 left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Blindly approved. Let me know if you want me to test it somehow.

@joa-quim

joa-quim commented Jul 8, 2026

Copy link
Copy Markdown
Member Author

You can get something alike with

gmt grdcut @earth_relief_02m_g -R-15/-7.5/34/39.5 -Gbat.grd
gmt nswing -Rbat.grd -F-12.49593345/35.93634937/25/58.2/90/10/215/53.75/10 -Ginit.grd
gmt nswing bat.grd init.grd -t3 -Gtsu+m,100 -N100 -V

The above created 3 grids (the -G...,100 means save only every 100'th simulation time step (which was 3 sec). Now the problem for visualization is that the grids contain also the topography and we don't see a wave of 1 or a couple meters high in presence of land topography. This is where the Mirone visualization tools come to help but we can see just the wave if we crop to an oceanic region only.

Meanwhile I hit another gmt bug. grdinterpolate fails to extract slices from the .nc cube (that is why above I used -Gtsu+m to write individual 2D grids) apparectly because it's not able to navigate well in the different variables of the cube.

gmt grdimage tsu00597.grd -R-15/-9.2/34/38 -I -png lixo
lixo

@Esteban82

Copy link
Copy Markdown
Member

I think it works fine. I tried these commands:

gmt grdcut @earth_relief_02m_g -R-15/-7.5/34/39.5 -Gbat.grd 
gmt nswing -Rbat.grd -F-12.49593345/35.93634937/25/58.2/90/10/215/53.75/10 -Ginit.grd
gmt nswing bat.grd init.grd -df -t3 -Gtsum,100 -N100 -V
gmt grdimage tsum.nc -R-15/-9.2/34/38 -I -png lixo

lixo

I got these messages in the terminal:

nswing [WARNING]: NSWING: -R region is in geographic coordinates. Assuming -fg.
nswing [INFORMATION]: netCDF grid bat.grd has a default CPT geo.
nswing [INFORMATION]: Reading grid from file bat.grd
nswing [INFORMATION]: gmt_grd_BC_set: Called with a pad < 2; skipped.
nswing [WARNING]: NSWING: Bathymetry grid metadata says geographic coordinates. Assuming -fg.
nswing [INFORMATION]: netCDF grid init.grd has no default CPT.
nswing [INFORMATION]: Reading grid from file init.grd
nswing [INFORMATION]: gmt_grd_BC_set: Called with a pad < 2; skipped.
nswing [INFORMATION]: 
NSWING: 

nswing [INFORMATION]: Layer 0  time step = 3    x_min = -15     x_max = -7.5    y_min = 34      y_max = 39.5
nswing [INFORMATION]: dtCFL = 15.7639   Courant number (sqrt(g*h)*dt / max(dx,dy)) = 0.190309
nswing [INFORMATION]: 
Using 1 OpenMP threads
nswing [INFORMATION]: 
Using DISCHARGE limit to minimize sources of instability
nswing [INFORMATION]: 
nswing [INFORMATION]:   100 %   CPU secs/ticks = 743403.000

@joa-quim

joa-quim commented Jul 8, 2026

Copy link
Copy Markdown
Member Author

You are plotting the bathymetry. To extract a layer from the .nc file you need #9095 That is why I used in my example
-Gtsu+m,100
and there is no -df option. Just -fg (or -f) to confirm grid is geographic

@Esteban82

Copy link
Copy Markdown
Member

Ok. I see. Now I got this:

lixo

I used these commands:

gmt grdcut @earth_relief_02m_g -R-15/-7.5/34/39.5 -Gbat.grd 
gmt nswing -Rbat.grd -F-12.49593345/35.93634937/25/58.2/90/10/215/53.75/10 -Ginit.grd
gmt nswing bat.grd init.grd -t3 -Gtsu+m,100 -N100 -V
gmt grdimage tsu00297.grd -R-15/-9.2/34/38 -I -png lixo

@joa-quim

joa-quim commented Jul 8, 2026

Copy link
Copy Markdown
Member Author

Yes, that is correct

@Esteban82

Copy link
Copy Markdown
Member

Great, now it can be merged.

@Esteban82

Copy link
Copy Markdown
Member

Maybe it would be good to add some references in the rst file. For example I found this one:

https://www.researchgate.net/publication/275349940_Validation_of_NSWING_a_multi-core_finite_difference_code_for_tsunami_propagation_and_run-up

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

add-changelog Add PR to the changelog AI-assisted All (or most) of the code was written by Artificial Intelligence.

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants