Actual source code: ex4fwd.c
1: static char help[] = "2D Shallow water equations forward model with MMS verification.\n"
2: "Implements 2D shallow water equations with 3 DOF per grid point (h, hu, hv).\n\n"
3: "Usage:\n"
4: " ./ex4fwd -steps 500 -nx 80 -ny 80\n"
5: " ./ex4fwd -verify_mms -steps 100 -nx 40 -ny 40 -verification_freq 10\n"
6: " ./ex4fwd -test_mms_spatial_order -steps 5 -dt 1e-4\n"
7: " ./ex4fwd -test_mms_spatial_order -conv_nx_coarse 20 -conv_ny_coarse 20 -conv_refine 2 -steps 5 -dt 1e-4\n\n";
9: #include <petscdmda.h>
10: #include <petscts.h>
11: #include "ex4.h"
13: /* Default parameter values */
14: #define DEFAULT_NX 40
15: #define DEFAULT_NY 40
16: #define DEFAULT_STEPS 100
17: #define DEFAULT_G 9.81
18: #define DEFAULT_DT 0.02
19: #define DEFAULT_LX 80.0
20: #define DEFAULT_LY 80.0
21: #define DEFAULT_H0 1.5
22: #define DEFAULT_AX 0.2
23: #define DEFAULT_AY 0.2
24: #define DEFAULT_PROGRESS_FREQ 10
25: #define DEFAULT_VERIFICATION_FREQ 10
26: #define DEFAULT_CONV_NX_COARSE 12
27: #define DEFAULT_CONV_NY_COARSE 12
28: #define DEFAULT_CONV_REFINE 2
30: /* Minimum acceptable observed convergence order for the MMS spatial-order check.
31: First-order Rusanov has asymptotic order 1; in the pre-asymptotic regime exercised by
32: the test grids the observed values are higher (~1.8). A floor well below 1 catches a
33: clear degradation without flagging routine variation. */
34: #define MIN_OBSERVED_ORDER 0.8
36: static PetscErrorCode ComputeManufacturedError(Vec numerical, DM da, PetscReal time, PetscReal Lx, PetscReal Ly, PetscReal h0, PetscReal A, PetscReal *L1_error, PetscReal *L2_error, PetscReal *Linf_error)
37: {
38: const PetscScalar ***x_num;
39: PetscInt xs, ys, xm, ym;
40: PetscInt nx, ny;
41: PetscReal dx, dy, dA;
42: PetscReal L1_local = 0.0, L2_local = 0.0, Linf_local = 0.0;
43: PetscReal sums[2];
45: PetscFunctionBeginUser;
46: PetscCall(DMDAGetInfo(da, NULL, &nx, &ny, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
47: dx = Lx / nx;
48: dy = Ly / ny;
49: dA = dx * dy;
50: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
51: PetscCall(DMDAVecGetArrayDOFRead(da, numerical, (void *)&x_num));
53: for (PetscInt j = ys; j < ys + ym; j++) {
54: for (PetscInt i = xs; i < xs + xm; i++) {
55: PetscReal x = ((PetscReal)i + 0.5) * dx;
56: PetscReal y = ((PetscReal)j + 0.5) * dy;
57: PetscReal h_exact, hu_exact, hv_exact;
58: PetscReal error;
60: ManufacturedSolution2D(Lx, Ly, x, y, time, h0, A, &h_exact, &hu_exact, &hv_exact);
61: error = PetscAbsReal(PetscRealPart(x_num[j][i][0]) - h_exact);
62: L1_local += error * dA;
63: L2_local += error * error * dA;
64: Linf_local = PetscMax(Linf_local, error);
65: }
66: }
68: PetscCall(DMDAVecRestoreArrayDOFRead(da, numerical, (void *)&x_num));
69: sums[0] = L1_local;
70: sums[1] = L2_local;
71: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, sums, 2, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)da)));
72: *L1_error = sums[0];
73: *L2_error = PetscSqrtReal(sums[1]);
74: PetscCallMPI(MPIU_Allreduce(&Linf_local, Linf_error, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)da)));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode RunManufacturedCase(const ShallowWater2DConfig *cfg, PetscInt steps, PetscReal *L1_err, PetscReal *L2_err, PetscReal *Linf_err)
79: {
80: DM da_state;
81: ShallowWater2DCtx *sw_ctx;
82: Vec x_numerical;
84: PetscFunctionBeginUser;
85: PetscCall(SetupForwardProblem(cfg, &da_state, &sw_ctx, &x_numerical));
86: PetscCall(SetInitialCondition(da_state, x_numerical, sw_ctx, PETSC_TRUE));
87: for (PetscInt step = 0; step < steps; step++) PetscCall(ShallowWaterStep2DVec(sw_ctx, step * cfg->dt, x_numerical));
88: PetscCall(ComputeManufacturedError(x_numerical, da_state, steps * cfg->dt, cfg->Lx, cfg->Ly, cfg->h0, cfg->Ax, L1_err, L2_err, Linf_err));
89: PetscCall(VecDestroy(&x_numerical));
90: PetscCall(ShallowWater2DContextDestroy(&sw_ctx));
91: PetscCall(DMDestroy(&da_state));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: int main(int argc, char **argv)
96: {
97: PetscInt nx = DEFAULT_NX, ny = DEFAULT_NY, steps = DEFAULT_STEPS, progress_freq = DEFAULT_PROGRESS_FREQ, verification_freq = DEFAULT_VERIFICATION_FREQ;
98: PetscInt conv_nx_coarse = DEFAULT_CONV_NX_COARSE, conv_ny_coarse = DEFAULT_CONV_NY_COARSE, conv_refine = DEFAULT_CONV_REFINE;
99: PetscReal g = DEFAULT_G, dt = DEFAULT_DT, Lx = DEFAULT_LX, Ly = DEFAULT_LY, h0 = DEFAULT_H0, Ax = DEFAULT_AX, Ay = DEFAULT_AY;
100: PetscBool verify_mms = PETSC_FALSE, test_mms_spatial_order = PETSC_FALSE;
101: ShallowWater2DConfig cfg;
103: PetscFunctionBeginUser;
104: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
105: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "2D Shallow Water Forward Model", NULL);
106: PetscCall(PetscOptionsInt("-nx", "Number of grid points in x", "", nx, &nx, NULL));
107: PetscCall(PetscOptionsInt("-ny", "Number of grid points in y", "", ny, &ny, NULL));
108: PetscCall(PetscOptionsInt("-steps", "Number of time steps", "", steps, &steps, NULL));
109: PetscCall(PetscOptionsReal("-g", "Gravitational constant", "", g, &g, NULL));
110: PetscCall(PetscOptionsReal("-dt", "Time step size", "", dt, &dt, NULL));
111: PetscCall(PetscOptionsReal("-Lx", "Domain length in x", "", Lx, &Lx, NULL));
112: PetscCall(PetscOptionsReal("-Ly", "Domain length in y", "", Ly, &Ly, NULL));
113: PetscCall(PetscOptionsReal("-h0", "Mean water height", "", h0, &h0, NULL));
114: PetscCall(PetscOptionsReal("-Ax", "Wave amplitude in x", "", Ax, &Ax, NULL));
115: PetscCall(PetscOptionsReal("-Ay", "Wave amplitude in y", "", Ay, &Ay, NULL));
116: PetscCall(PetscOptionsInt("-progress_freq", "Print progress every N steps (0 = only last)", "", progress_freq, &progress_freq, NULL));
117: PetscCall(PetscOptionsInt("-verification_freq", "Frequency for MMS error output", "", verification_freq, &verification_freq, NULL));
118: PetscCall(PetscOptionsBool("-verify_mms", "Enable manufactured-solution verification forcing", "", verify_mms, &verify_mms, NULL));
119: PetscCall(PetscOptionsBool("-test_mms_spatial_order", "Run a three-grid manufactured-solution spatial-order check and exit", "", test_mms_spatial_order, &test_mms_spatial_order, NULL));
120: PetscCall(PetscOptionsInt("-conv_nx_coarse", "Coarse-grid nx for manufactured-solution spatial-order check", "", conv_nx_coarse, &conv_nx_coarse, NULL));
121: PetscCall(PetscOptionsInt("-conv_ny_coarse", "Coarse-grid ny for manufactured-solution spatial-order check", "", conv_ny_coarse, &conv_ny_coarse, NULL));
122: PetscCall(PetscOptionsInt("-conv_refine", "Grid refinement factor for manufactured-solution spatial-order check", "", conv_refine, &conv_refine, NULL));
123: PetscOptionsEnd();
125: PetscCheck(nx > 0 && ny > 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Grid dimensions must be positive");
126: PetscCheck(steps >= 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of steps must be non-negative");
127: PetscCheck(dt > 0.0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Time step must be positive");
128: PetscCheck(!test_mms_spatial_order || steps > 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "-test_mms_spatial_order requires -steps > 0 (got %" PetscInt_FMT ")", steps);
129: PetscCheck(!test_mms_spatial_order || conv_refine >= 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "-conv_refine must be >= 2 for spatial-order check (got %" PetscInt_FMT ")", conv_refine);
130: PetscCheck(verification_freq > 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "-verification_freq must be positive (got %" PetscInt_FMT ")", verification_freq);
132: cfg.Lx = Lx;
133: cfg.Ly = Ly;
134: cfg.g = g;
135: cfg.dt = dt;
136: cfg.h0 = h0;
137: cfg.Ax = Ax;
138: cfg.Ay = Ay;
139: cfg.verify_mms = verify_mms;
141: if (test_mms_spatial_order) {
142: PetscInt medium_nx = conv_refine * conv_nx_coarse, medium_ny = conv_refine * conv_ny_coarse, fine_nx = conv_refine * medium_nx, fine_ny = conv_refine * medium_ny;
143: PetscReal coarse_L1, coarse_L2, coarse_Linf, medium_L1, medium_L2, medium_Linf, fine_L1, fine_L2, fine_Linf;
144: PetscReal order_cm_L1, order_cm_L2, order_cm_Linf;
145: PetscReal order_mf_L1, order_mf_L2, order_mf_Linf;
147: /* MMS spatial-order check: each refinement reuses cfg with Ay=Ax=A and verify_mms=TRUE
148: so the manufactured source is symmetric in x/y. */
149: cfg.Ay = Ax;
150: cfg.verify_mms = PETSC_TRUE;
151: cfg.nx = conv_nx_coarse;
152: cfg.ny = conv_ny_coarse;
153: PetscCall(RunManufacturedCase(&cfg, steps, &coarse_L1, &coarse_L2, &coarse_Linf));
154: cfg.nx = medium_nx;
155: cfg.ny = medium_ny;
156: PetscCall(RunManufacturedCase(&cfg, steps, &medium_L1, &medium_L2, &medium_Linf));
157: cfg.nx = fine_nx;
158: cfg.ny = fine_ny;
159: PetscCall(RunManufacturedCase(&cfg, steps, &fine_L1, &fine_L2, &fine_Linf));
160: PetscCheck(coarse_L1 > 0.0 && medium_L1 > 0.0 && fine_L1 > 0.0 && coarse_L2 > 0.0 && medium_L2 > 0.0 && fine_L2 > 0.0 && coarse_Linf > 0.0 && medium_Linf > 0.0 && fine_Linf > 0.0, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "MMS error norm collapsed to zero on at least one grid (truncation below print precision); increase -steps or -dt to obtain a measurable error before the order check");
161: order_cm_L1 = PetscLogReal(coarse_L1 / medium_L1) / PetscLogReal((PetscReal)conv_refine);
162: order_cm_L2 = PetscLogReal(coarse_L2 / medium_L2) / PetscLogReal((PetscReal)conv_refine);
163: order_cm_Linf = PetscLogReal(coarse_Linf / medium_Linf) / PetscLogReal((PetscReal)conv_refine);
164: order_mf_L1 = PetscLogReal(medium_L1 / fine_L1) / PetscLogReal((PetscReal)conv_refine);
165: order_mf_L2 = PetscLogReal(medium_L2 / fine_L2) / PetscLogReal((PetscReal)conv_refine);
166: order_mf_Linf = PetscLogReal(medium_Linf / fine_Linf) / PetscLogReal((PetscReal)conv_refine);
167: PetscCall(PetscPrintf(PETSC_COMM_WORLD,
168: "MMS spatial-order check\n"
169: " coarse grid : %" PetscInt_FMT " x %" PetscInt_FMT ", dt=%.5e, L1=%.5e, L2=%.5e, Linf=%.5e\n"
170: " medium grid : %" PetscInt_FMT " x %" PetscInt_FMT ", dt=%.5e, L1=%.5e, L2=%.5e, Linf=%.5e\n"
171: " fine grid : %" PetscInt_FMT " x %" PetscInt_FMT ", dt=%.5e, L1=%.5e, L2=%.5e, Linf=%.5e\n"
172: " observed p (coarse->medium) : L1=%.2f, L2=%.2f, Linf=%.2f\n"
173: " observed p (medium->fine) : L1=%.2f, L2=%.2f, Linf=%.2f\n",
174: conv_nx_coarse, conv_ny_coarse, (double)dt, (double)coarse_L1, (double)coarse_L2, (double)coarse_Linf, medium_nx, medium_ny, (double)dt, (double)medium_L1, (double)medium_L2, (double)medium_Linf, fine_nx, fine_ny, (double)dt, (double)fine_L1, (double)fine_L2, (double)fine_Linf, (double)order_cm_L1, (double)order_cm_L2, (double)order_cm_Linf, (double)order_mf_L1, (double)order_mf_L2, (double)order_mf_Linf));
175: /* Guard against silent regression: the printed rates are rounded for portability and would
176: hide a degraded order. Assert each observed rate exceeds MIN_OBSERVED_ORDER explicitly. */
177: PetscCheck(order_cm_L1 >= MIN_OBSERVED_ORDER && order_cm_L2 >= MIN_OBSERVED_ORDER && order_cm_Linf >= MIN_OBSERVED_ORDER && order_mf_L1 >= MIN_OBSERVED_ORDER && order_mf_L2 >= MIN_OBSERVED_ORDER && order_mf_Linf >= MIN_OBSERVED_ORDER, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MMS spatial-order regression: observed orders fell below %g (coarse->medium L1=%g L2=%g Linf=%g; medium->fine L1=%g L2=%g Linf=%g)", (double)MIN_OBSERVED_ORDER, (double)order_cm_L1, (double)order_cm_L2, (double)order_cm_Linf, (double)order_mf_L1, (double)order_mf_L2, (double)order_mf_Linf);
178: } else {
179: DM da_state;
180: ShallowWater2DCtx *sw_ctx;
181: Vec x_numerical;
183: cfg.nx = nx;
184: cfg.ny = ny;
185: PetscCall(SetupForwardProblem(&cfg, &da_state, &sw_ctx, &x_numerical));
186: PetscCall(SetInitialCondition(da_state, x_numerical, sw_ctx, verify_mms));
188: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "2D Shallow Water Forward Example\n==============================\n"));
190: for (PetscInt step = 1; step <= steps; step++) {
191: PetscReal time = step * dt;
192: PetscCall(ShallowWaterStep2DVec(sw_ctx, (step - 1) * dt, x_numerical));
193: if (verify_mms && (step % verification_freq == 0 || step == steps)) {
194: PetscReal L1_err, L2_err, Linf_err;
195: PetscCall(ComputeManufacturedError(x_numerical, da_state, time, Lx, Ly, h0, Ax, &L1_err, &L2_err, &Linf_err));
196: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Step %4" PetscInt_FMT ", time %6.3f L1=%.5e L2=%.5e Linf=%.5e\n", step, (double)time, (double)L1_err, (double)L2_err, (double)Linf_err));
197: } else if (step == steps || (progress_freq > 0 && step % progress_freq == 0)) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Step %4" PetscInt_FMT ", time %6.3f Forward step complete\n", step, (double)time));
198: }
200: PetscCall(PetscPrintf(PETSC_COMM_WORLD, verify_mms ? "\nMMS forward run complete.\n" : "\nForward simulation complete.\n"));
201: PetscCall(VecDestroy(&x_numerical));
202: PetscCall(ShallowWater2DContextDestroy(&sw_ctx));
203: PetscCall(DMDestroy(&da_state));
204: }
205: PetscCall(PetscFinalize());
206: return 0;
207: }
209: /*TEST
211: test:
212: suffix: verification
213: requires: !complex
214: args: -verify_mms -steps 50 -nx 40 -ny 40 -verification_freq 10
215: output_file: output/ex4fwd_verification.out
217: test:
218: suffix: verification_parallel
219: requires: !complex
220: nsize: 2
221: args: -verify_mms -steps 20 -nx 30 -ny 30 -verification_freq 5
222: output_file: output/ex4fwd_verification_parallel.out
224: test:
225: suffix: mms_spatial
226: requires: !complex !single
227: args: -test_mms_spatial_order -steps 5 -dt 1e-4
228: filter: grep -E "MMS spatial-order check|observed p" | sed -E "s/=([0-9]+)\.([0-9])[0-9]/=\\1.\\2/g"
230: TEST*/