Skip to content

Revert the boundary block constrain that blew up the Catalina benchmark 1 #9180

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
joa-quim merged 4 commits into from
Sep 5, 2026
+95 −14
Merged

Revert the boundary block constrain that blew up the Catalina benchmark 1 #9180

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
7fc4aba
Start fixing a nswing regression that caused a numeric instability
joa-quim Sep 4, 2026
74c21eb
Move location of the check paternity test
joa-quim Sep 4, 2026
fef6558
Last fix. Turn the "Margin block" into a warning instead of erroring.
joa-quim Sep 4, 2026
63c7c01
Add -W option to allow communicate to external precesses (iGMT in thi…
joa-quim Sep 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/nswing/benchmark1_catalina.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -241,7 +241,7 @@ function run_nswing(bat::GMTgrid, src::GMTgrid, nests::Vector{<:GMTgrid}, pato::
push!(args, "-t$(dt)")
#push!(args, "-T$(maregs_pos[1])/$(maregs_pos[2])+o$(outfile)+t$(mareg_interval)")

cmd = "nswing -V " * join(args, " ")
cmd = "nswing -v " * join(args, " ")
println(cmd)
gmt(cmd, bat, src, nests...)
return nothing
Expand Down
107 changes: 94 additions & 13 deletions src/nswing/nswing.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,43 @@ union {uint64_t i; double d;} loc_nan = {0x7ff8000000000000};

static double EPS4 = EPS4_; /* Kinda trick to be able to change EPS4 via a command line option */

/* ------------------------------------------------------------------------------------------ *
* Progress reporting to an external caller (iGMT, GMT.jl, ...) that drives us in-process.
* Neither stdout nor a file is involved: the caller either registers a callback that we invoke
* on every waitbar tick (push), or polls nswing_get_progress() from another thread (pull).
* Both are updated from the serial part of the main loop only, never from inside an OpenMP
* region, so the callback runs on the same thread that entered GMT_nswing.
* ------------------------------------------------------------------------------------------ */
typedef void (*nswing_progress_fn)(int percent, int cycle, int n_cycles, double model_time, void *data);

/* A caller that drives us OUT-OF-PROCESS (iGMT runs a detached `gmt nswing`, so its window never
* blocks) cannot see the statics or the callback above. For those, -W<file> is rewritten in place
* once per tick. Losing progress must never fail a run, so a file that cannot be opened is only a
* warning. */


static nswing_progress_fn nswing_progress_cb = NULL;
static void *nswing_progress_data = NULL;
static int nswing_prog_percent = 0, nswing_prog_cycle = 0, nswing_prog_n_cycles = 0;
static double nswing_prog_time = 0.0;

EXTERN_MSC void nswing_set_progress_callback(nswing_progress_fn fn, void *data) {
nswing_progress_cb = fn; nswing_progress_data = data;
}

EXTERN_MSC void nswing_get_progress(int *percent, int *cycle, int *n_cycles, double *model_time) {
if (percent) *percent = nswing_prog_percent;
if (cycle) *cycle = nswing_prog_cycle;
if (n_cycles) *n_cycles = nswing_prog_n_cycles;
if (model_time) *model_time = nswing_prog_time;
}

GMT_LOCAL void nswing_report_progress(int percent, int cycle, int n_cycles, double model_time) {
nswing_prog_percent = percent; nswing_prog_cycle = cycle;
nswing_prog_n_cycles = n_cycles; nswing_prog_time = model_time;
if (nswing_progress_cb) nswing_progress_cb(percent, cycle, n_cycles, model_time, nswing_progress_data);
}

#define MAXRUNUP -50 /* Do not waste time computing flood above this altitude */
#define V_LIMIT 20 /* Upper limit of maximum velocity */

Expand Down Expand Up @@ -405,6 +442,9 @@ GMT_LOCAL int usage(struct GMTAPI_CTRL *API, int level) {
GMT_Usage(API, -2, "Manning friction coefficients. If only one is provided, use it for all nesting levels; "
"otherwise specify one per level, comma separated. Append +<depth> to apply Manning only at depths "
"shallower than <depth> (positive up).");
GMT_Usage(API, 1, "\n-W<file> Report the run progress to another process.");
GMT_Usage(API, -2, "The file is rewritten in place once per tick; its payload is: "
"percent, cycle, n_cycles, model time.");
GMT_Usage(API, 1, "\n-x<n>");
GMT_Usage(API, -2, "Number of cores to use in a parallel run [Default is the max in the machine].");
GMT_Option(API, "V,f,.");
Expand All @@ -424,7 +464,7 @@ struct NSWING_CTRL {
bool verbose, out_velocity, out_velocity_x, out_velocity_y, out_velocity_r, out_maregs_velocity;
bool cumpt, do_2Dgrids, do_maxs, mareg_xy;
bool append_z;
char *bathy, *fonte, *fname_sww, *basename_most, *bnc_file;
char *bathy, *fonte, *fname_sww, *basename_most, *bnc_file, *prog_file;
char *nesteds[10];
char hcum[256];
char maregs[256];
Expand Down Expand Up @@ -484,6 +524,7 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC
char *fname3D = NULL;
char *fonte = NULL;
char *bnc_file = NULL;
char *prog_file = NULL;
char fname_mask_lbeach[256] = "";
char fname_mask_sbeach[256] = "";
char tracers_infile[256] = "", tracers_outfile[256] = "";
Expand Down Expand Up @@ -842,6 +883,14 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC
case 'U':
nest.do_upscale = true;
break;
case 'W': /* Report our advance to an external caller: -W<file> */
if (opt->arg[0] == '\0') {
GMT_Report(GMT->parent, GMT_MSG_ERROR, "NSWING: Error, -W option. Must provide a file name\n");
error++;
break;
}
prog_file = opt->arg;
break;
case 'v': /* Undocumented: prints only the NSWING setup/diagnostics block */
verbose = true;
break;
Expand Down Expand Up @@ -875,8 +924,6 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC
}
}

if (gmt_M_is_verbose(GMT, GMT_MSG_INFORMATION)) verbose = true; /* Global -V (parsed by GMT_Parse_Common) also triggers the -v printout */

if (GMT->common.R.active[RSET]) { /* -R was parsed by GMT_Parse_Common */
got_R = true;
dfXmin = GMT->common.R.wesn[XLO]; dfXmax = GMT->common.R.wesn[XHI];
Expand Down Expand Up @@ -1037,6 +1084,7 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC
Ctrl->fname_sww = fname_sww;
Ctrl->basename_most = basename_most;
Ctrl->bnc_file = bnc_file;
Ctrl->prog_file = prog_file;
strcpy(Ctrl->hcum, hcum);
strcpy(Ctrl->maregs, maregs);
strcpy(Ctrl->stem, stem);
Expand Down Expand Up @@ -1107,6 +1155,8 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) {
char *fname3D = NULL; /* Name pointer for the 3D netCDF file */
char *fonte = NULL; /* Name pointer for tsunami source file */
char *bnc_file = NULL; /* Name pointer for a boundary condition file */
char *prog_file = NULL; /* Name pointer for the -W progress file */
FILE *fp_prog = NULL; /* ... and its handle, kept open for the whole run */
char fname_mask_lbeach[256] = ""; /* Name pointer for the "long_beach" mask grid */
char fname_mask_sbeach[256] = ""; /* Name pointer for the "short_beach" mask grid */
char tracers_infile[256] = "", tracers_outfile[256] = ""; /* Names for in and out tracers files */
Expand Down Expand Up @@ -1243,6 +1293,7 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) {
fname_sww = Ctrl->fname_sww;
basename_most = Ctrl->basename_most;
bnc_file = Ctrl->bnc_file;
prog_file = Ctrl->prog_file;
strcpy(hcum, Ctrl->hcum);
strcpy(maregs, Ctrl->maregs);
strcpy(stem, Ctrl->stem);
Expand Down Expand Up @@ -1671,7 +1722,8 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) {
if (initialize_nestum(API, &nest, isGeog, k))
{free_arrays(&nest, isGeog, num_of_nestGrids); Return(-1);}
}
/* Check if nesting grids fit nicely within each others. Maybe too late? */
/* Check if nesting grids fit nicely within each others. Only here because check_paternity's
* loop is gated on nest->level[], which is set by initialize_nestum just above. */
if (check_paternity(API, &nest)) Return(-1);
nest.time_h = time_h;
/* Resample eta(s) in descendent grids to avoid initial jumps at borders */
Expand Down Expand Up @@ -1771,7 +1823,10 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) {

/* --------------------------------------------------------------------------------------- */
if (verbose) {
API->GMT->current.setting.verbose = GMT_MSG_INFORMATION; /* So that next messages are printed */
unsigned int verbose_bak = API->GMT->current.setting.verbose;
/* Raise the message level just for this block so -v prints the same setup info that -V used to,
* then put it back: leaving it raised would also turn on GMT's own chatter (grid reads/writes). */
API->GMT->current.setting.verbose = GMT_MSG_INFORMATION;
GMT_Report(API, GMT_MSG_INFORMATION, "\nNSWING: \n\n");
GMT_Report(API, GMT_MSG_INFORMATION, "Layer 0 time step = %g\tx_min = %g\tx_max = %g\ty_min = %g\ty_max = %g\n",
dt, hdr_b.wesn[XLO], hdr_b.wesn[XHI], hdr_b.wesn[YLO], hdr_b.wesn[YHI]);
Expand Down Expand Up @@ -1815,6 +1870,7 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) {
GMT_Report(API, GMT_MSG_INFORMATION, "\nUsing DISCHARGE limit to minimize sources of instability\n");
#endif
GMT_Report(API, GMT_MSG_INFORMATION, "\n");
API->GMT->current.setting.verbose = verbose_bak;
}
/* --------------------------------------------------------------------------------------- */

Expand All @@ -1839,6 +1895,11 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) {

one_100 = (double)(n_of_cycles) / 100.0;

nswing_report_progress(0, 0, n_of_cycles, 0.0); /* Reset, so a second run in the same session does not start at 100 */

if (prog_file && (fp_prog = fopen(prog_file, "w")) == NULL) /* Not fatal: losing progress must not kill a run */
GMT_Report(API, GMT_MSG_WARNING, "NSWING: Could not create progress file %s. Proceeding without it.\n", prog_file);

LoopKabas: /* When computing a grid of Kabas we use a GOTO to simulate a loop. Sorry but have to. */
/* --------------------------------------------------------------------------------------- */
/* Begin main iteration */
Expand All @@ -1849,7 +1910,19 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) {
prc = (double)iprc / 100.;
iprc++;
prc = (double)(k+1) / (double)n_of_cycles;
GMT_Report(API, GMT_MSG_INFORMATION, "\t%d %%\r", iprc);
/* GMT_Message ignores the message level, so -v shows the progress without -V's chatter.
* Flush it: when our stream is a pipe (an external caller capturing us) the CRT switches
* to full buffering and this line, having no '\n', would only surface when we exit. */
if (verbose || gmt_M_is_verbose(API->GMT, GMT_MSG_INFORMATION)) {
GMT_Message(API, GMT_TIME_NONE, "\t%d %%\r", iprc);
fflush(API->GMT->session.std[GMT_ERR]);
}
nswing_report_progress(iprc, k+1, n_of_cycles, nest.time_h);
if (fp_prog) { /* Rewrite the progress file in place so pollers always read one short line */
rewind(fp_prog);
fprintf(fp_prog, "%d %d %d %.3f\n", iprc, k+1, n_of_cycles, nest.time_h);
fflush(fp_prog);
}
}

/* ------------------------------------------------------------------------------------ */
Expand Down Expand Up @@ -2303,7 +2376,17 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) {
free(oranges);
}

GMT_Report(API, GMT_MSG_INFORMATION, "\t100 %%\tCPU secs/ticks = %.3f\n", (double)(clock() - tic));
if (verbose || gmt_M_is_verbose(API->GMT, GMT_MSG_INFORMATION)) {
GMT_Message(API, GMT_TIME_NONE, "\t100 %%\tCPU secs/ticks = %.3f\n", (double)(clock() - tic));
fflush(API->GMT->session.std[GMT_ERR]);
}
nswing_report_progress(100, n_of_cycles, n_of_cycles, nest.time_h);
if (fp_prog) { /* Final 100% so a poller sees completion without having to watch the process */
rewind(fp_prog);
fprintf(fp_prog, "%d %d %d %.3f\n", 100, n_of_cycles, n_of_cycles, nest.time_h);
fclose(fp_prog);
fp_prog = NULL;
}

if (cumpt) {
if (fp) fclose(fp); /* Not opened when maregraphs went to netCDF */
Expand Down Expand Up @@ -2412,16 +2495,14 @@ int check_paternity(void *API, struct nestContainer *nest) {
error++;
}

if (k == 1) { /* First nested grid must sit at least 2 bathymetry cells inside the bathymetry grid, on every side */
if (k == 1) { /* Only warn: a tight fit is legal (benchmark 1 uses it) but is a classic trouble source */
double margin_w = (nest->hdr[0].wesn[XLO] - nest->hdr[1].wesn[XLO]) / nest->hdr[0].inc[GMT_X];
double margin_e = (nest->hdr[1].wesn[XHI] - nest->hdr[0].wesn[XHI]) / nest->hdr[0].inc[GMT_X];
double margin_s = (nest->hdr[0].wesn[YLO] - nest->hdr[1].wesn[YLO]) / nest->hdr[0].inc[GMT_Y];
double margin_n = (nest->hdr[1].wesn[YHI] - nest->hdr[0].wesn[YHI]) / nest->hdr[0].inc[GMT_Y];
if (margin_w > -2 || margin_e > -2 || margin_s > -2 || margin_n > -2) {
GMT_Report(API, GMT_MSG_ERROR, "NSWING: Error, bad input. The first nested grid must be at least 2 rows/columns "
"shorter than the bathymetry grid on every side (West, East, South, North).\n");
error++;
}
if (margin_w > -2 || margin_e > -2 || margin_s > -2 || margin_n > -2)
GMT_Report(API, GMT_MSG_WARNING, "NSWING: First nested grid is less than 2 bathymetry cells inside the "
"bathymetry grid. Legal (benchmark 1 is such a case) but may cause trouble.\n");
}

if (error) /* Abort since any further info would be false/useless */
Expand Down
Loading