diff --git a/CMakeLists.txt b/CMakeLists.txt index 07beb8f5..e2d3c1c7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -48,7 +48,7 @@ add_library(wtm src/update_effective_storativity.cpp ) target_compile_options(wtm PRIVATE -Wfloat-conversion -Wall -Wextra -pedantic -Wshadow) -target_link_libraries(wtm PUBLIC richdem OpenMP::OpenMP_CXX fmt::fmt PkgConfig::PETSC ) +target_link_libraries(wtm PUBLIC richdem OpenMP::OpenMP_CXX fmt::fmt PkgConfig::PETSC) target_compile_features(wtm PUBLIC cxx_std_20) add_executable(wtm.x diff --git a/src/CreateSNES.cpp b/src/CreateSNES.cpp index 3b5eb63f..396d7630 100644 --- a/src/CreateSNES.cpp +++ b/src/CreateSNES.cpp @@ -27,9 +27,14 @@ void InitialiseSNES(AppCtx& user_context, Parameters& params) { DMSetUp(user_context.da); user_context.make_global_vectors(); + user_context.make_local_vectors(); DMSetApplicationContext(user_context.da, &user_context); SNESSetDM(user_context.snes, user_context.da); + // Anderson mixing converges reliably without a Jacobian for this nonlinear problem. + // m=1 (1 history vector) is sufficient and avoids the instability seen with m>1. + // Override with -snes_type or -snes_anderson_m at runtime if needed. + SNESSetType(user_context.snes, SNESANDERSON); SNESSetFromOptions(user_context.snes); } diff --git a/src/CreateSNES.hpp b/src/CreateSNES.hpp index a7472f25..2c919761 100644 --- a/src/CreateSNES.hpp +++ b/src/CreateSNES.hpp @@ -18,10 +18,15 @@ struct AppCtx { Vec mask = nullptr; Vec topo_vec = nullptr; Vec rech_vec = nullptr; - Vec T_vec = nullptr; Vec porosity_vec = nullptr; Vec starting_wtd = nullptr; + // Local ghost vectors for fields accessed at neighbor indices in FormFunctionLocal + Vec topo_local = nullptr; + Vec fdepth_local = nullptr; + Vec ksat_local = nullptr; + Vec T_local = nullptr; // scratch: 1/T, computed over ghost range each F eval + // Extract global vectors from DM; then duplicate for remaining // vectors that are the same types void make_global_vectors() { @@ -33,10 +38,16 @@ struct AppCtx { VecDuplicate(x, &mask); VecDuplicate(x, &topo_vec); VecDuplicate(x, &rech_vec); - VecDuplicate(x, &T_vec); VecDuplicate(x, &porosity_vec); VecDuplicate(x, &starting_wtd); } + + void make_local_vectors() { + DMCreateLocalVector(da, &topo_local); + DMCreateLocalVector(da, &fdepth_local); + DMCreateLocalVector(da, &ksat_local); + DMCreateLocalVector(da, &T_local); + } }; void InitialiseSNES(AppCtx& user_context, Parameters& params); diff --git a/src/DMDA_array_pack.cpp b/src/DMDA_array_pack.cpp index 07cf852a..da7d8807 100644 --- a/src/DMDA_array_pack.cpp +++ b/src/DMDA_array_pack.cpp @@ -12,10 +12,39 @@ void populate_DMDA_array_pack(AppCtx& user_context, ArrayPack& arp, DMDA_Array_P for (auto i = xs; i < xs + xm; i++) { dmdapack.cellsize_EW_squared[j][i] = arp.cellsize_e_w_metres[j] * arp.cellsize_e_w_metres[j]; dmdapack.mask[j][i] = arp.land_mask(i, j); - dmdapack.fdepth_vec[j][i] = arp.fdepth(i, j); - dmdapack.ksat_vec[j][i] = arp.ksat(i, j); - dmdapack.topo_vec[j][i] = arp.topo(i, j); dmdapack.porosity_vec[j][i] = arp.porosity(i, j); } } } + +// Populate the global topo/fdepth/ksat vecs from arp and scatter to local ghost vectors. +// Must be called while these global vecs are NOT under DMDAVecGetArray (i.e., before or after +// DMDA_Array_Pack holds them — which it does NOT, by design). +void scatter_static_fields(AppCtx& user_context, ArrayPack& arp) { + const auto [xs, ys, xm, ym] = get_corners(user_context.da); + PetscScalar **topo_arr, **fdepth_arr, **ksat_arr; + + DMDAVecGetArray(user_context.da, user_context.topo_vec, &topo_arr); + DMDAVecGetArray(user_context.da, user_context.fdepth_vec, &fdepth_arr); + DMDAVecGetArray(user_context.da, user_context.ksat_vec, &ksat_arr); + for (auto j = ys; j < ys + ym; j++) { + for (auto i = xs; i < xs + xm; i++) { + topo_arr[j][i] = arp.topo(i, j); + fdepth_arr[j][i] = arp.fdepth(i, j); + ksat_arr[j][i] = arp.ksat(i, j); + } + } + DMDAVecRestoreArray(user_context.da, user_context.topo_vec, &topo_arr); + DMDAVecRestoreArray(user_context.da, user_context.fdepth_vec, &fdepth_arr); + DMDAVecRestoreArray(user_context.da, user_context.ksat_vec, &ksat_arr); + + // The DMDA's internal PetscSF is shared across all GlobalToLocal operations on the same DM. + // Overlapping Begin calls (Begin A, Begin B, End A, End B) confuse the SF state machine; + // each pair must be completed sequentially. + DMGlobalToLocalBegin(user_context.da, user_context.topo_vec, INSERT_VALUES, user_context.topo_local); + DMGlobalToLocalEnd(user_context.da, user_context.topo_vec, INSERT_VALUES, user_context.topo_local); + DMGlobalToLocalBegin(user_context.da, user_context.fdepth_vec, INSERT_VALUES, user_context.fdepth_local); + DMGlobalToLocalEnd(user_context.da, user_context.fdepth_vec, INSERT_VALUES, user_context.fdepth_local); + DMGlobalToLocalBegin(user_context.da, user_context.ksat_vec, INSERT_VALUES, user_context.ksat_local); + DMGlobalToLocalEnd(user_context.da, user_context.ksat_vec, INSERT_VALUES, user_context.ksat_local); +} diff --git a/src/DMDA_array_pack.hpp b/src/DMDA_array_pack.hpp index c054121a..7f028b22 100644 --- a/src/DMDA_array_pack.hpp +++ b/src/DMDA_array_pack.hpp @@ -3,27 +3,23 @@ struct DMDA_Array_Pack { PetscScalar** x = nullptr; PetscScalar** cellsize_EW_squared = nullptr; - PetscScalar** fdepth_vec = nullptr; - PetscScalar** ksat_vec = nullptr; PetscScalar** mask = nullptr; - PetscScalar** topo_vec = nullptr; PetscScalar** rech_vec = nullptr; - PetscScalar** T_vec = nullptr; PetscScalar** porosity_vec = nullptr; PetscScalar** starting_wtd = nullptr; const AppCtx* context = nullptr; + // topo_vec, fdepth_vec, ksat_vec are intentionally NOT held here. + // They are scattered to AppCtx local ghost vectors before each solve so that + // FormFunctionLocal can safely access neighbor indices across MPI boundaries. + DMDA_Array_Pack(const AppCtx& user) { assert(!context); // Make sure we're not already initialized context = &user; DMDAVecGetArray(user.da, user.x, &x); DMDAVecGetArray(user.da, user.cellsize_EW_squared, &cellsize_EW_squared); - DMDAVecGetArray(user.da, user.fdepth_vec, &fdepth_vec); - DMDAVecGetArray(user.da, user.ksat_vec, &ksat_vec); DMDAVecGetArray(user.da, user.mask, &mask); - DMDAVecGetArray(user.da, user.topo_vec, &topo_vec); DMDAVecGetArray(user.da, user.rech_vec, &rech_vec); - DMDAVecGetArray(user.da, user.T_vec, &T_vec); DMDAVecGetArray(user.da, user.porosity_vec, &porosity_vec); DMDAVecGetArray(user.da, user.starting_wtd, &starting_wtd); } @@ -32,12 +28,8 @@ struct DMDA_Array_Pack { assert(context); // Make sure we are already initialized DMDAVecRestoreArray(context->da, context->x, &x); DMDAVecRestoreArray(context->da, context->cellsize_EW_squared, &cellsize_EW_squared); - DMDAVecRestoreArray(context->da, context->fdepth_vec, &fdepth_vec); - DMDAVecRestoreArray(context->da, context->ksat_vec, &ksat_vec); DMDAVecRestoreArray(context->da, context->mask, &mask); - DMDAVecRestoreArray(context->da, context->topo_vec, &topo_vec); DMDAVecRestoreArray(context->da, context->rech_vec, &rech_vec); - DMDAVecRestoreArray(context->da, context->T_vec, &T_vec); DMDAVecRestoreArray(context->da, context->porosity_vec, &porosity_vec); DMDAVecRestoreArray(context->da, context->starting_wtd, &starting_wtd); context = nullptr; @@ -45,3 +37,4 @@ struct DMDA_Array_Pack { }; void populate_DMDA_array_pack(AppCtx& user_context, ArrayPack& arp, DMDA_Array_Pack& dmdapack); +void scatter_static_fields(AppCtx& user_context, ArrayPack& arp); diff --git a/src/WTM.cpp b/src/WTM.cpp index d9b9b4f6..97c59a90 100644 --- a/src/WTM.cpp +++ b/src/WTM.cpp @@ -104,8 +104,13 @@ void update( if ((params.cycles_done % params.cycles_to_save) == 0) { // Save the output every "cycles_to_save" iterations, under a new filename // so we can compare how the water table has changed through time. - arp.wtd.setNoData(-9999); - arp.wtd.saveGDAL(fmt::format("{}{:09}.tif", params.outfile_prefix, params.cycles_done)); + // wtd is fully assembled on all ranks by FanDarcyGroundwater::update; rank 0 writes. + PetscMPIInt rank; + MPI_Comm_rank(PETSC_COMM_WORLD, &rank); + if (rank == 0) { + arp.wtd.setNoData(-9999); + arp.wtd.saveGDAL(fmt::format("{}{:09}.tif", params.outfile_prefix, params.cycles_done)); + } } arp.wtd_old = arp.wtd; // These are used to see how much change occurs @@ -228,9 +233,13 @@ void finalise(Parameters& params, ArrayPack& arp, AppCtx& user_context) { std::ofstream textfile(params.textfilename, std::ios_base::app); textfile << "p done with processing" << std::endl; - // save the final answer for water table depth. - arp.wtd.setNoData(-9999); - arp.wtd.saveGDAL(fmt::format("{}{:09}.tif", params.outfile_prefix, params.cycles_done)); + // Save the final answer. wtd is assembled on all ranks; only rank 0 writes to avoid conflicts. + PetscMPIInt rank; + MPI_Comm_rank(PETSC_COMM_WORLD, &rank); + if (rank == 0) { + arp.wtd.setNoData(-9999); + arp.wtd.saveGDAL(fmt::format("{}{:09}.tif", params.outfile_prefix, params.cycles_done)); + } textfile.close(); @@ -244,9 +253,12 @@ void finalise(Parameters& params, ArrayPack& arp, AppCtx& user_context) { VecDestroy(&user_context.mask); VecDestroy(&user_context.topo_vec); VecDestroy(&user_context.rech_vec); - VecDestroy(&user_context.T_vec); VecDestroy(&user_context.porosity_vec); VecDestroy(&user_context.starting_wtd); + VecDestroy(&user_context.topo_local); + VecDestroy(&user_context.fdepth_local); + VecDestroy(&user_context.ksat_local); + VecDestroy(&user_context.T_local); } int main(int argc, char** argv) { @@ -269,6 +281,9 @@ int main(int argc, char** argv) { DMDA_Array_Pack dmdapack(user_context); // this needs to come after initialise populate_DMDA_array_pack(user_context, arp, dmdapack); + // Scatter topo/fdepth/ksat to local ghost vectors. These global vecs are not held by + // dmdapack so there is no GetArray lock conflict. + scatter_static_fields(user_context, arp); run(params, arp, user_context, dmdapack); diff --git a/src/transient_groundwater.cpp b/src/transient_groundwater.cpp index 467f1dee..bd6620c7 100644 --- a/src/transient_groundwater.cpp +++ b/src/transient_groundwater.cpp @@ -150,6 +150,22 @@ int update(Parameters& params, ArrayPack& arp, AppCtx& user_context, DMDA_Array_ } } + // Assemble the full wtd field across all MPI ranks. Each rank has updated + // only its owned cells; non-owned entries are left at their previous values. + // Build a buffer with owned cells non-zero and everything else zeroed, then + // sum across ranks so that every rank ends up with the complete correct field. + { + const int total = params.ncells_x * params.ncells_y; + std::vector owned_only(total, 0.0); + for (int j = ys; j < ys + ym; j++) + for (int i = xs; i < xs + xm; i++) + owned_only[j * params.ncells_x + i] = arp.wtd(i, j); + MPI_Allreduce(MPI_IN_PLACE, owned_only.data(), total, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); + for (int j = 0; j < params.ncells_y; j++) + for (int i = 0; i < params.ncells_x; i++) + arp.wtd(i, j) = owned_only[j * params.ncells_x + i]; + } + return 0; } @@ -230,21 +246,23 @@ static PetscErrorCode FormFunctionLocal(DMDALocalInfo* info, PetscScalar** x, Pe **my_porosity; /* - Compute function over the locally owned part of the grid - */ + Compute function over the locally owned part of the grid. + topo/fdepth/ksat/T use local ghost vectors so neighbor accesses [j][i±1] are valid under MPI. + */ PetscCall(DMDAVecGetArray(da, user_context->mask, &my_mask)); PetscCall(DMDAVecGetArray(da, user_context->cellsize_EW_squared, &cellsize_ew_sq)); - PetscCall(DMDAVecGetArray(da, user_context->fdepth_vec, &my_fdepth)); - PetscCall(DMDAVecGetArray(da, user_context->ksat_vec, &my_ksat)); - PetscCall(DMDAVecGetArray(da, user_context->topo_vec, &my_topo)); + PetscCall(DMDAVecGetArray(da, user_context->fdepth_local, &my_fdepth)); + PetscCall(DMDAVecGetArray(da, user_context->ksat_local, &my_ksat)); + PetscCall(DMDAVecGetArray(da, user_context->topo_local, &my_topo)); PetscCall(DMDAVecGetArray(da, user_context->rech_vec, &my_rech)); - PetscCall(DMDAVecGetArray(da, user_context->T_vec, &my_T)); + PetscCall(DMDAVecGetArray(da, user_context->T_local, &my_T)); PetscCall(DMDAVecGetArray(da, user_context->porosity_vec, &my_porosity)); PetscCall(DMDAVecGetArray(da, user_context->starting_wtd, &my_starting_wtd)); + // Compute 1/T over the full ghost range so neighbor lookups in the owned-range loop below are valid. #pragma omp parallel for default(none) shared(info, my_T, x, my_topo, my_fdepth, my_ksat) collapse(2) - for (auto j = info->ys; j < info->ys + info->ym; j++) { - for (auto i = info->xs; i < info->xs + info->xm; i++) { + for (auto j = info->gys; j < info->gys + info->gym; j++) { + for (auto i = info->gxs; i < info->gxs + info->gxm; i++) { my_T[j][i] = 1. / depthIntegratedTransmissivity(x[j][i] - my_topo[j][i], my_fdepth[j][i], my_ksat[j][i]); } } @@ -283,11 +301,11 @@ static PetscErrorCode FormFunctionLocal(DMDALocalInfo* info, PetscScalar** x, Pe PetscCall(DMDAVecRestoreArray(da, user_context->mask, &my_mask)); PetscCall(DMDAVecRestoreArray(da, user_context->cellsize_EW_squared, &cellsize_ew_sq)); - PetscCall(DMDAVecRestoreArray(da, user_context->fdepth_vec, &my_fdepth)); - PetscCall(DMDAVecRestoreArray(da, user_context->ksat_vec, &my_ksat)); - PetscCall(DMDAVecRestoreArray(da, user_context->topo_vec, &my_topo)); + PetscCall(DMDAVecRestoreArray(da, user_context->fdepth_local, &my_fdepth)); + PetscCall(DMDAVecRestoreArray(da, user_context->ksat_local, &my_ksat)); + PetscCall(DMDAVecRestoreArray(da, user_context->topo_local, &my_topo)); PetscCall(DMDAVecRestoreArray(da, user_context->rech_vec, &my_rech)); - PetscCall(DMDAVecRestoreArray(da, user_context->T_vec, &my_T)); + PetscCall(DMDAVecRestoreArray(da, user_context->T_local, &my_T)); PetscCall(DMDAVecRestoreArray(da, user_context->porosity_vec, &my_porosity)); PetscCall(DMDAVecRestoreArray(da, user_context->starting_wtd, &my_starting_wtd)); diff --git a/tests/ghost_cell/check_results.py b/tests/ghost_cell/check_results.py new file mode 100755 index 00000000..63c3b79d --- /dev/null +++ b/tests/ghost_cell/check_results.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python3 +""" +Compare 1-process and 2-process WTM output TIFs. + +With the ghost-cell fix both runs must agree at all interior cells. +Without it the two halves of the domain are hydrologically decoupled at the +MPI processor boundary (x = NX//2), producing a clearly visible (~O(1) m) error. +""" + +import sys +import os +import glob +import numpy as np +import rasterio + +NX = 12 + + +def last_tif(prefix, outdir): + pattern = os.path.join(outdir, f"{prefix}*.tif") + tifs = sorted(glob.glob(pattern)) + if not tifs: + raise FileNotFoundError(f"No TIF files matching {pattern}") + return tifs[-1] + + +def load(path): + with rasterio.open(path) as src: + data = src.read(1).astype(np.float64) + return data + + +def main(): + dir1p = "out_1p" + dir2p = "out_2p" + + f1 = last_tif("out_", dir1p) + f2 = last_tif("out_", dir2p) + print(f"1-process output : {f1}") + print(f"2-process output : {f2}") + + w1 = load(f1) + w2 = load(f2) + + if w1.shape != w2.shape: + print(f"FAIL: shape mismatch {w1.shape} vs {w2.shape}", file=sys.stderr) + sys.exit(1) + + diff = np.abs(w1 - w2) + + # Interior land cells only (exclude ocean edges at row 0, row NY-1, col 0, col NX-1) + interior = diff[1:-1, 1:-1] + + max_diff = interior.max() + mean_diff = interior.mean() + boundary_diff = diff[1:-1, NX // 2 - 1 : NX // 2 + 1].max() + + print(f"Max |Δwtd| interior : {max_diff:.6f} m") + print(f"Mean |Δwtd| interior : {mean_diff:.6f} m") + print(f"Max |Δwtd| at MPI bound : {boundary_diff:.6f} m") + + # Threshold: numerical summation-order differences are O(1e-10) m. + # The ghost-cell bug produces O(1) m errors near the boundary. + TOLERANCE = 1e-4 # generous; anything above ~0.01 m signals the bug + + if max_diff > TOLERANCE: + print( + f"\nFAIL: max difference {max_diff:.4f} m exceeds tolerance {TOLERANCE} m.\n" + "Ghost-cell error is present — the MPI boundary suppresses inter-rank flux.", + file=sys.stderr, + ) + sys.exit(1) + else: + print(f"\nPASS: 1-process and 2-process outputs agree to within {TOLERANCE} m.") + sys.exit(0) + + +if __name__ == "__main__": + main() diff --git a/tests/ghost_cell/make_inputs.py b/tests/ghost_cell/make_inputs.py new file mode 100755 index 00000000..3cea1663 --- /dev/null +++ b/tests/ghost_cell/make_inputs.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python3 +""" +Generate synthetic GeoTIFF inputs for the ghost-cell MPI validation test. + +Domain: 102 x 3 cells (100 interior x-cells, 1 interior y-row). +Left half (x=1..50): ksat = 1e-4 m/s +Right half (x=51..100): ksat = 1e-3 m/s (10x higher) + +The ksat discontinuity at x=50/51 coincides with the MPI processor boundary +when the job is split with -da_processors_x 2 -da_processors_y 1. With the +ghost-cell bug the two halves are hydrologically decoupled (no flux across the +boundary); with the fix, the 2-process result matches the 1-process result. +""" + +import numpy as np +import os +import rasterio +from rasterio.transform import from_bounds + +NX = 12 # total columns (10 interior + 2 ocean-edge columns) +NY = 3 # total rows (1 interior row + 2 ocean-edge rows) + +REGION = "ghost_cell_test" +TIME = "t0" +OUTDIR = os.path.join(os.path.dirname(__file__), "inputs") +os.makedirs(OUTDIR, exist_ok=True) + +# Arbitrary geographic extent; richdem only needs dimensions + data values. +transform = from_bounds(0, 0, NX, NY, NX, NY) +CRS = "EPSG:4326" + + +def write_tif(path, data, dtype="float32"): + with rasterio.open( + path, "w", + driver="GTiff", + height=NY, width=NX, + count=1, dtype=dtype, + crs=CRS, + transform=transform, + ) as dst: + dst.write(data.astype(dtype), 1) + + +# Flat topography at 100 m. Edge cells will be forced to 0 (ocean) by the +# model's land_mask.setEdges(0) + wtd=topo=0 for mask==0 logic. +topo = np.full((NY, NX), 100.0, dtype=np.float32) + +# Flat slope → fdepth = fdepth_a (set to 200 in the config) +slope = np.zeros((NY, NX), dtype=np.float32) + +# Land mask: 1 inside, 0 at edges. The model also calls setEdges(0) itself, +# so this is just for the file to have valid dimensions. +mask = np.ones((NY, NX), dtype=np.float32) +mask[0, :] = 0 +mask[-1, :] = 0 +mask[:, 0] = 0 +mask[:, -1] = 0 + +precip = np.full((NY, NX), 0.3, dtype=np.float32) # m/yr +evap = np.zeros((NY, NX), dtype=np.float32) # m/yr +open_water_evap = np.full((NY, NX), 0.4, dtype=np.float32) # m/yr +winter_temp = np.zeros((NY, NX), dtype=np.float32) # deg C (>-5 → no permafrost) + +# Heterogeneous ksat: factor-of-10 jump at the midpoint. +# This asymmetry means there is non-zero flux at the MPI boundary in the +# correct solution; the ghost-cell bug suppresses it, causing a clear error. +ksat = np.full((NY, NX), 1e-4, dtype=np.float32) +ksat[:, NX // 2:] = 1e-3 # right half: 10× higher + +porosity = np.full((NY, NX), 0.25, dtype=np.float32) + +# Write all files +files = { + f"{REGION}_{TIME}_topography.tif": topo, + f"{REGION}_{TIME}_slope.tif": slope, + f"{REGION}_{TIME}_mask.tif": mask, + f"{REGION}_{TIME}_precipitation.tif": precip, + f"{REGION}_{TIME}_evaporation.tif": evap, + f"{REGION}_{TIME}_open_water_evaporation.tif": open_water_evap, + f"{REGION}_{TIME}_winter_temperature.tif": winter_temp, + f"{REGION}_horizontal_ksat.tif": ksat, + f"{REGION}_porosity.tif": porosity, +} + +for fname, arr in files.items(): + path = os.path.join(OUTDIR, fname) + write_tif(path, arr) + print(f" wrote {path}") + +print("Done.") diff --git a/tests/ghost_cell/run_test.sh b/tests/ghost_cell/run_test.sh new file mode 100755 index 00000000..b7896aac --- /dev/null +++ b/tests/ghost_cell/run_test.sh @@ -0,0 +1,70 @@ +#!/usr/bin/env bash +# Ghost-cell MPI validation test. +# +# Runs the WTM on a 102x3 heterogeneous-ksat domain with 1 and then 2 MPI +# processes, then compares the output TIFs. With the ghost-cell fix the two +# runs agree; without it they diverge at the MPI processor boundary. +# +# Usage: +# cd tests/ghost_cell +# ./run_test.sh [path/to/wtm.x] +# +# Default binary: ../../build/wtm.x + +set -euo pipefail +cd "$(dirname "$0")" + +WTM=${1:-../../build/wtm.x} + +if [[ ! -x "$WTM" ]]; then + echo "ERROR: WTM binary not found at $WTM" >&2 + echo "Build the project first (cmake --build build) or supply the path as \$1." >&2 + exit 1 +fi + +echo "=== Ghost-cell MPI validation test ===" +echo "WTM binary: $WTM" +echo + +# 1. Generate synthetic inputs +echo "--- Generating inputs ---" +python3 make_inputs.py + +# 2. One-process reference run +echo +echo "--- 1-process reference run ---" +rm -rf out_1p run_1p.txt +mkdir -p out_1p +# outfile_prefix and textfilename are relative to CWD when wtm.x is called. +# Override them with sed-generated temp configs so both runs share the same +# base config without modifying it. +CFG_1P=$(mktemp /tmp/ghost_cell_1p_XXXXXX.cfg) +sed 's|^outfile_prefix.*|outfile_prefix out_1p/out_|; + s|^textfilename.*|textfilename run_1p.txt|' ghost_cell.cfg > "$CFG_1P" +mpirun -n 1 "$WTM" "$CFG_1P" \ + -snes_mf \ + -snes_stol 1e-6 \ + 2>&1 | grep -E 'SNES|converged|norm|Error|error' || true +echo "1-process run complete." +rm -f "$CFG_1P" + +# 3. Two-process run (processor boundary at the ksat discontinuity) +echo +echo "--- 2-process run (split at ksat boundary) ---" +rm -rf out_2p run_2p.txt +mkdir -p out_2p +CFG_2P=$(mktemp /tmp/ghost_cell_2p_XXXXXX.cfg) +sed 's|^outfile_prefix.*|outfile_prefix out_2p/out_|; + s|^textfilename.*|textfilename run_2p.txt|' ghost_cell.cfg > "$CFG_2P" +mpirun -n 2 "$WTM" "$CFG_2P" \ + -snes_mf \ + -snes_stol 1e-6 \ + -da_processors_x 2 -da_processors_y 1 \ + 2>&1 | grep -E 'SNES|converged|norm|Error|error' || true +echo "2-process run complete." +rm -f "$CFG_2P" + +# 4. Compare outputs +echo +echo "--- Comparing outputs ---" +python3 check_results.py