diff --git a/.gitignore b/.gitignore index 04984d8..06b46f0 100644 --- a/.gitignore +++ b/.gitignore @@ -54,6 +54,7 @@ __pycache__/ # Output files produced by the program OUTPUT_FILES/ +figures/ tests/*/*.png tests/*/*.h5 diff --git a/examples/00_checkerboard_iso/run_this_example.sh b/examples/00_checkerboard_iso/run_this_example.sh index 055273e..b55da2a 100644 --- a/examples/00_checkerboard_iso/run_this_example.sh +++ b/examples/00_checkerboard_iso/run_this_example.sh @@ -7,7 +7,7 @@ cp src_rec_file_ph.csv OUTPUT_FILES/src_rec_file_forward_RL_PH.csv # cp src_rec_file_ph.csv OUTPUT_FILES/src_rec_file_forward_LV_PH.csv # create 2x3x2 checkers and forward simulate surface traveltimes -mpirun -np $NPROC ../../bin/surfatt_cb_fwd -i input_params.yml -n 2/3/2 -m 0.2 -p 0.08 +mpirun -np $NPROC ../../bin/SURFATT_cb_fwd -i input_params.yml -n 2/3/2 -m 0.2 -p 0.08 # inversion -mpirun -np $NPROC ../../bin/surfatt_tomo -i input_params.yml +mpirun -np $NPROC ../../bin/SURFATT_tomo -i input_params.yml diff --git a/include/h5io.h b/include/h5io.h index e408da3..e23d042 100644 --- a/include/h5io.h +++ b/include/h5io.h @@ -94,6 +94,9 @@ class H5IO { std::vector read_vector(const std::string &name) const { H5::DataSet ds = file_.openDataSet(name); H5::DataSpace sp = ds.getSpace(); + if (sp.getSimpleExtentNdims() != 1) + throw std::runtime_error( + "H5IO::read_vector: dataset '" + name + "' is not 1-D"); hsize_t n = 0; sp.getSimpleExtentDims(&n, nullptr); std::vector v(n); diff --git a/include/inversion1d.h b/include/inversion1d.h index 1f101b1..0daa384 100644 --- a/include/inversion1d.h +++ b/include/inversion1d.h @@ -7,6 +7,8 @@ #include #include +namespace surfker { struct DepthKernel1D; } + class Inversion1D { public: explicit Inversion1D(WaveType wavetype); @@ -26,4 +28,23 @@ class Inversion1D { int niter; std::vector misfits; + + // Validate depth grid, initial Vs and the active src_rec tables before + // the first forward call; logs every problem found and aborts if any. + void check_inputs(const Eigen::VectorX& zarr) const; + + // Abort with a diagnostic if disper() failed (returned 0 / non-finite). + void check_pred_vel(const Eigen::VectorX& pred_vel, + const Eigen::VectorX& periods, + const Eigen::VectorX& zarr, + int iter, WaveType wt, SurfType tp) const; + + // Abort with a diagnostic if the depth kernels contain inf/NaN. + void check_kernels(const surfker::DepthKernel1D& kernels, + const Eigen::VectorX& periods, + const Eigen::VectorX& zarr, + int iter, WaveType wt, SurfType tp) const; + + // Log the current Vs profile and its first non-physical node, if any. + void log_vs_profile(const Eigen::VectorX& zarr, int n_update) const; }; diff --git a/include/src_rec.h b/include/src_rec.h index 63802c6..0c43c57 100644 --- a/include/src_rec.h +++ b/include/src_rec.h @@ -72,7 +72,13 @@ class SrcRec { // Rank 0 reads the file; each field gets its own per-node MPI shared-memory // window via Parallel::alloc_shared (1-D for numeric, 2-D for strings). // Must call release_shm() before MPI_Finalize(). - void load(const std::string& filepath); + // The period column is always validated (finite, > 0); with check_obs the + // observed tt and vel columns are too. Aborts all ranks on invalid rows. + void load(const std::string& filepath, bool check_obs = false); + + // Log an error listing rows whose value in `col` is non-finite or <= 0. + // Returns the number of such rows (0 = column is valid). + int check_positive(const real_t* col, const std::string& col_name) const; // gather synthetic travel times void gather_syn_tt(); @@ -134,6 +140,7 @@ class SrcRec { void get_events(); void get_periods(); + std::string filepath_; // CSV path, for diagnostics int nsrc_total_ = 0; // Number of rows after load() int n_obs_ = 0; diff --git a/scripts/plot_1d_model_and_data.py b/scripts/plot_1d_model_and_data.py new file mode 100644 index 0000000..9fa6403 --- /dev/null +++ b/scripts/plot_1d_model_and_data.py @@ -0,0 +1,339 @@ +#!/usr/bin/env python +""" +Plot 1-D summaries of a SurfATT inversion. + +1. Horizontally averaged 1-D Vs profile (depth vs Vs) of the initial and final + models, with the lateral standard deviation at each depth. +2. For each src_rec file: number of data at each period, and a period-depth + map of the 1-D sensitivity kernels computed on the averaged initial model. + +Kernels are computed with bin/SURFATT_kernel1d, which uses the same surfker +engine as SURFATT_tomo (same layering, Earth flattening, fundamental mode and +empirical Vp/rho relations). Build it with: + cd build && cmake .. && make -j SURFATT_kernel1d + +Edit the parameters below and run it from the project folder (the parent of +OUTPUT_FILES), e.g.: + cd examples/00_checkerboard_iso + python ../../scripts/plot_1d_model_and_data.py + +Requirements: numpy, pandas, h5py, matplotlib +""" + +import os +import shutil +import subprocess +import tempfile + +import h5py +import numpy as np +import pandas as pd +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from matplotlib.colors import LinearSegmentedColormap + +SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) + +# ============================================================================= +# User parameters (relative paths are relative to the current working directory) +# ============================================================================= +OUTPUT_DIR = "OUTPUT_FILES" +INITIAL_MODEL = os.path.join(OUTPUT_DIR, "initial_model.h5") +FINAL_MODEL = os.path.join(OUTPUT_DIR, "final_model.h5") + +# Dataset to average (e.g. "vs", "vsv", "vsh"); kernels use the initial profile +MODEL_KEY = "vs" + +# Data files used in the inversion. Keys are _; set a value to +# None (or remove the entry) to skip that data type. +SRC_REC_FILES = { + "RL_PH": os.path.join(OUTPUT_DIR, "src_rec_file_forward_RL_PH.csv"), + "RL_GR": None, + "LV_PH": None, + "LV_GR": None, +} + +# Horizontal region for averaging: [lon_min, lon_max, lat_min, lat_max]. +# None uses all grid points (including the margin grid). +REGION = None + +# Path to the kernel tool (default: bin/ of this repository); if not found, +# SURFATT_kernel1d is searched in PATH +KERNEL_BIN = os.path.join(SCRIPT_DIR, "..", "bin", "SURFATT_kernel1d") +# Launcher prefix, e.g. ["mpirun", "-np", "1"] if the bare MPI binary fails +MPI_LAUNCHER = [] + +# "total": Vs kernel as used in the vs-only inversion (Rayleigh: Vp and rho +# scaled from Vs; Love: same as "vs", as in SURFATT_tomo) +# "vs": partial d(vel)/d(Vs) with Vp and rho fixed +KERNEL_TYPE = "total" +# Normalize each period's kernel by its maximum absolute value +NORMALIZE_KERNEL = True + +# Figure output +FIG_DIR = "figures" +FIG_FORMAT = "png" +DPI = 300 + +# ============================================================================= +# Plot style +# ============================================================================= +INK = "#0b0b0b" # primary text +INK_2 = "#52514e" # secondary text, axes +GRID = "#e2e1dc" # recessive grid +MODEL_STYLE = { # line color and style for each model + "Initial": {"color": "#eb6834", "ls": "--", "zorder": 3}, + "Final": {"color": "#2a78d6", "ls": "-", "zorder": 2}, +} +BAR_COLOR = "#2a78d6" +# Sequential one-hue ramp for non-negative kernels +KERNEL_CMAP = LinearSegmentedColormap.from_list( + "kernel_blue", + ["#fcfcfb", "#cde2fb", "#9ec5f4", "#6da7ec", "#3987e5", "#256abf", "#184f95", "#0d366b"], +) +# Diverging ramp (red - gray - blue) used when kernels have negative lobes +KERNEL_CMAP_DIV = LinearSegmentedColormap.from_list( + "kernel_div", + ["#9e2322", "#e34948", "#f4a8a0", "#f0efec", "#9ec5f4", "#3987e5", "#104281"], +) + +plt.rcParams.update({ + "font.size": 10, + "axes.edgecolor": INK_2, + "axes.labelcolor": INK, + "axes.titlecolor": INK, + "xtick.color": INK_2, + "ytick.color": INK_2, + "axes.spines.top": False, + "axes.spines.right": False, + "axes.grid": True, + "grid.color": GRID, + "grid.linewidth": 0.6, + "legend.frameon": False, +}) + +BASE_DIR = os.getcwd() + + +def resolve(path): + """Return an absolute path, interpreting relative paths from the working directory.""" + return path if os.path.isabs(path) else os.path.join(BASE_DIR, path) + + +# ============================================================================= +# Data loading +# ============================================================================= +def load_model(path, key): + """Read lon (x), lat (y), depth (z) and a model field of shape (nx, ny, nz).""" + with h5py.File(path, "r") as f: + if key not in f: + raise KeyError(f"Dataset '{key}' not found in {path}; available: {list(f.keys())}") + x, y, z = f["x"][:], f["y"][:], f["z"][:] + field = f[key][:] + if field.shape != (len(x), len(y), len(z)): + raise ValueError(f"{path}: {key} shape {field.shape} != (nx, ny, nz) = " + f"{(len(x), len(y), len(z))}") + return x, y, z, field + + +def horizontal_stats(x, y, field, region=None): + """Mean and standard deviation over horizontal grid points at each depth.""" + if region is not None: + lon_min, lon_max, lat_min, lat_max = region + ix = (x >= lon_min) & (x <= lon_max) + iy = (y >= lat_min) & (y <= lat_max) + if not ix.any() or not iy.any(): + raise ValueError(f"REGION {region} contains no grid points") + field = field[np.ix_(ix, iy)] + return field.mean(axis=(0, 1)), field.std(axis=(0, 1)) + + +def load_src_rec(path): + """Number of data and mean observed velocity at each period.""" + df = pd.read_csv(path) + grp = df.groupby("period") + return pd.DataFrame({"count": grp.size(), "mean_vel": grp["vel"].mean()}) + + +# ============================================================================= +# 1-D kernels via SURFATT_kernel1d +# ============================================================================= +def find_kernel_bin(): + """Locate the SURFATT_kernel1d executable.""" + path = resolve(KERNEL_BIN) + if os.path.isfile(path): + return path + found = shutil.which(os.path.basename(KERNEL_BIN)) + if found: + return found + raise FileNotFoundError( + f"SURFATT_kernel1d not found at {path} or in PATH. Build it with:\n" + " cd build && cmake .. && make -j SURFATT_kernel1d") + + +def compute_kernel_1d(z, vs, periods, data_type): + """Dispersion and depth kernels of a 1-D model for one data type (e.g. RL_PH).""" + wave, vtype = data_type.split("_") + with tempfile.TemporaryDirectory() as tmp: + fin = os.path.join(tmp, "model1d.h5") + fout = os.path.join(tmp, "kernel1d.h5") + with h5py.File(fin, "w") as f: + f["z"] = np.asarray(z, dtype=float) + f["vs"] = np.asarray(vs, dtype=float) + f["periods"] = np.asarray(periods, dtype=float) + cmd = MPI_LAUNCHER + [find_kernel_bin(), "-i", fin, "-o", fout, "-w", wave, "-t", vtype] + res = subprocess.run(cmd, capture_output=True, text=True) + if res.returncode != 0 or not os.path.isfile(fout): + raise RuntimeError(f"Command failed: {' '.join(cmd)}\n{res.stdout}\n{res.stderr}") + with h5py.File(fout, "r") as f: + ker = {k: f[k][:] for k in f.keys()} + bad = [k for k, v in ker.items() if not np.isfinite(v).all()] + if bad: + raise ValueError(f"{data_type}: NaN/Inf in {bad} from SURFATT_kernel1d") + return ker + + +# ============================================================================= +# Plotting +# ============================================================================= +def plot_average_model(z, stats, fname): + """Mean +/- std profiles and std profiles of the averaged models.""" + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 6), sharey=True, + gridspec_kw={"width_ratios": [2, 1]}) + for label, (mean, std) in stats.items(): + st = MODEL_STYLE[label] + ax1.fill_betweenx(z, mean - std, mean + std, color=st["color"], alpha=0.18, lw=0) + ax1.plot(mean, z, color=st["color"], ls=st["ls"], lw=1.8, zorder=st["zorder"], + label=f"{label} mean ± 1σ") + ax2.plot(std, z, color=st["color"], ls=st["ls"], lw=1.8, zorder=st["zorder"], label=label) + + ax1.set_xlabel(f"{MODEL_KEY} (km/s)") + ax1.set_ylabel("Depth (km)") + ax1.set_title("(a) Average 1-D model", loc="left") + ax1.legend(loc="lower left") + ax2.set_xlabel(f"Std of {MODEL_KEY} (km/s)") + ax2.set_title("(b) Lateral std", loc="left") + ax1.set_ylim(z.max(), z.min()) + + region = "all grid points" if REGION is None else f"region {REGION}" + fig.suptitle(f"Horizontally averaged model ({region})", color=INK) + fig.tight_layout() + fig.savefig(fname, dpi=DPI) + plt.close(fig) + print(f"Saved {fname}") + + +def period_edges(periods): + """Cell edges halfway between neighbouring periods.""" + if len(periods) == 1: + return np.array([periods[0] - 0.5, periods[0] + 0.5]) + mid = 0.5 * (periods[1:] + periods[:-1]) + return np.concatenate(([2 * periods[0] - mid[0]], mid, [2 * periods[-1] - mid[-1]])) + + +def plot_data_and_kernels(data_type, table, z, ker, fname): + """Data count per period (top) and period-depth kernel map (bottom) sharing the period axis.""" + periods = table.index.values + edges = period_edges(periods) + key = "sen_vs_total" if KERNEL_TYPE == "total" else "sen_vs" + vel_sym = "c" if data_type.endswith("PH") else "U" + + # kernel matrix of shape (nz, nper) + kmat = ker[key].T + if NORMALIZE_KERNEL: + kmat = kmat / np.maximum(np.abs(kmat).max(axis=0), 1e-30) + + fig = plt.figure(figsize=(7.5, 7.5)) + gs = fig.add_gridspec(2, 2, height_ratios=[1, 2.6], width_ratios=[1, 0.03], + hspace=0.2, wspace=0.03) + ax_cnt = fig.add_subplot(gs[0, 0]) + ax_ker = fig.add_subplot(gs[1, 0], sharex=ax_cnt) + cax = fig.add_subplot(gs[1, 1]) + + # (a) number of data per period + ax_cnt.bar(periods, table["count"].values, width=0.8 * np.diff(edges), color=BAR_COLOR) + ax_cnt.set_ylabel("Number of data") + ax_cnt.set_title(f"(a) Data count (total {int(table['count'].sum())})", loc="left") + ax_cnt.grid(axis="x", visible=False) + ax_cnt.tick_params(labelbottom=False) + + # (b) kernel map; diverging colors centred at zero if kernels have negative lobes + kmin, kmax = float(kmat.min()), float(kmat.max()) + if kmin < -0.01 * max(abs(kmax), 1e-30): + vabs = max(abs(kmin), abs(kmax)) + cmap, vmin, vmax = KERNEL_CMAP_DIV, -vabs, vabs + else: + cmap, vmin, vmax = KERNEL_CMAP, 0.0, kmax + pcm = ax_ker.pcolormesh(periods, z, kmat, cmap=cmap, vmin=vmin, vmax=vmax, + shading="nearest", rasterized=True) + ax_ker.set_xlim(edges[0], edges[-1]) + ax_ker.set_ylim(z.max(), z.min()) + ax_ker.grid(False) + ax_ker.set_xlabel("Period (s)") + ax_ker.set_ylabel("Depth (km)") + ax_ker.set_title("(b) 1-D sensitivity, initial model", loc="left") + + prefix = "Normalized " if NORMALIZE_KERNEL else "" + cbar = fig.colorbar(pcm, cax=cax) + scaled = KERNEL_TYPE == "total" and data_type.startswith("RL") + cbar.set_label(f"{prefix}∂{vel_sym}/∂{MODEL_KEY}" + (" (Vp, ρ scaled)" if scaled else "")) + cbar.outline.set_visible(False) + cax.grid(False) + + fig.suptitle(f"{data_type}: data count and 1-D sensitivity", color=INK, y=0.95) + fig.savefig(fname, dpi=DPI, bbox_inches="tight") + plt.close(fig) + print(f"Saved {fname}") + + +# ============================================================================= +# Main +# ============================================================================= +def main(): + fig_dir = resolve(FIG_DIR) + os.makedirs(fig_dir, exist_ok=True) + + # ---- averaged 1-D models ---- + profiles, stats, z = {}, {}, None + for label, path in (("Initial", INITIAL_MODEL), ("Final", FINAL_MODEL)): + path = resolve(path) + if not os.path.isfile(path): + print(f"Warning: {path} not found; {label.lower()} model skipped") + continue + x, y, zz, field = load_model(path, MODEL_KEY) + if z is not None and not np.allclose(z, zz): + raise ValueError("Initial and final models have different depth grids") + z = zz + stats[label] = horizontal_stats(x, y, field, REGION) + profiles[label] = stats[label][0] + if not stats: + raise FileNotFoundError("Neither the initial nor the final model was found") + + plot_average_model(z, stats, os.path.join(fig_dir, f"model_1d_average.{FIG_FORMAT}")) + + # ---- data count and kernels (on the averaged initial model) ---- + if "Initial" not in profiles: + print("Warning: initial model not found; data count and sensitivity figures skipped") + return + for data_type, path in SRC_REC_FILES.items(): + if path is None: + continue + path = resolve(path) + if not os.path.isfile(path): + print(f"Warning: {path} not found; {data_type} skipped") + continue + table = load_src_rec(path) + ker = compute_kernel_1d(z, profiles["Initial"], table.index.values, data_type) + table["vel_initial"] = ker["vel"] + + print(f"\n{data_type} ({os.path.relpath(path, BASE_DIR)}):") + print(" mean_vel: mean observed velocity; vel_initial: predicted on the averaged initial model") + print(table.to_string(float_format=lambda v: f"{v:.4f}")) + + plot_data_and_kernels(data_type, table, z, ker, + os.path.join(fig_dir, f"data_sensitivity_{data_type}.{FIG_FORMAT}")) + +if __name__ == "__main__": + main() diff --git a/src/SURFATT_kernel1d.cxx b/src/SURFATT_kernel1d.cxx new file mode 100644 index 0000000..4b31e98 --- /dev/null +++ b/src/SURFATT_kernel1d.cxx @@ -0,0 +1,240 @@ +// SURFATT_kernel1d +// +// Compute surface-wave dispersion and 1-D depth sensitivity kernels for a +// single 1-D Vs profile with the same surfker engine used by SURFATT_tomo +// (same layering, Earth flattening, fundamental mode and Brocher vp/rho). +// +// Input HDF5 datasets (1-D): +// z depth nodes (km), strictly increasing, shape (nz) +// vs Vs at the depth nodes (km/s), > 0, shape (nz) +// periods periods (s), > 0 and strictly increasing, shape (nper) +// +// Output HDF5 datasets: +// z, periods +// vel phase or group velocity (km/s), shape (nper) +// vp, rho empirical Vp (km/s) and density (g/cm^3), shape (nz) +// sen_vs d(vel)/d(vs), shape (nper, nz) +// sen_vp d(vel)/d(vp), shape (nper, nz); zeros for Love waves +// sen_rho d(vel)/d(rho), shape (nper, nz) +// sen_vs_total vs kernel as combined by the vs-only inversion of SURFATT_tomo, +// shape (nper, nz): for Rayleigh waves vp and rho are scaled +// from vs by the empirical relations; for Love waves it equals +// sen_vs +// +// A period without a root, or any NaN/Inf in the velocities or kernels, is +// reported as an error. + +#include "surfker/surfker.hpp" +#include "h5io.h" +#include "argparser.h" +#include "logger.h" +#include "parallel.h" +#include "utils.h" + +#include +#include +#include +#include + +namespace { + +struct Kernel1DArgs { + std::string fname; + std::string outfname; + WaveType wave_type = WaveType::RL; + SurfType surf_type = SurfType::PH; +}; + +Kernel1DArgs argparse_kernel1d(int argc, char* argv[]) { + ArgList al(argc, argv); + if (al.empty() || al.has("-h")) { + std::cout << + "Usage: SURFATT_kernel1d -i model_file -o out_file [-w RL|LV] [-t PH|GR] [-h]\n\n" + "Compute dispersion and 1-D depth sensitivity kernels of a 1-D Vs model\n" + "with the same surfker engine used by SURFATT_tomo.\n\n" + "required arguments:\n" + " -i model_file HDF5 file with datasets z (km), vs (km/s) and periods (s)\n" + " -o out_file Output HDF5 file with vel, sen_vs, sen_vp, sen_rho\n" + " and sen_vs_total\n\n" + "optional arguments:\n" + " -w RL|LV Wave type: Rayleigh or Love (default: RL)\n" + " -t PH|GR Velocity type: phase or group (default: PH)\n" + " -h Print help message\n"; + std::exit(0); + } + Kernel1DArgs out; + out.fname = al.require("-i"); + out.outfname = al.require("-o"); + if (auto v = al.get("-w")) { + if (*v == "RL") out.wave_type = WaveType::RL; + else if (*v == "LV") out.wave_type = WaveType::LV; + else throw std::runtime_error("-w must be RL or LV, got \"" + *v + "\""); + } + if (auto v = al.get("-t")) { + if (*v == "PH") out.surf_type = SurfType::PH; + else if (*v == "GR") out.surf_type = SurfType::GR; + else throw std::runtime_error("-t must be PH or GR, got \"" + *v + "\""); + } + return out; +} + +Eigen::VectorX read_eigen_vector(const H5IO &file, const std::string &name) { + if (!file.exists(name)) { + throw std::runtime_error("Required dataset '" + name + "' was not found"); + } + auto v = file.read_vector(name); + return Eigen::Map>(v.data(), static_cast(v.size())); +} + +// SURFATT_tomo always uses a strictly increasing depth grid and a sorted, +// unique period list; surfker relies on both (layer thicknesses from z, +// root bracketing from the previous period). +void require_strictly_increasing(const Eigen::VectorX &v, const std::string &name) { + for (Eigen::Index i = 0; i < v.size(); ++i) { + if (!std::isfinite(v(i))) { + throw std::runtime_error(fmt::format("{} has a non-finite value at index {}", name, i)); + } + if (i > 0 && v(i) <= v(i - 1)) { + throw std::runtime_error(fmt::format( + "{} must be strictly increasing (index {}: {} <= {})", name, i, v(i), v(i - 1))); + } + } +} + +// Throw if any entry of a kernel matrix is NaN/Inf. +void require_finite(const Eigen::MatrixX &M, const std::string &name, + const Eigen::VectorX &periods) { + for (Eigen::Index i = 0; i < M.rows(); ++i) { + if (!M.row(i).allFinite()) { + throw std::runtime_error(fmt::format( + "NaN/Inf in {} at period {:.3f} s", name, periods(i))); + } + } +} + +int kernel1d(const Kernel1DArgs &args, ATTLogger &logger) { + const std::string tag = waveTypeStr[static_cast(args.wave_type)] + "_" + + surfTypeStr[static_cast(args.surf_type)]; + + // read 1-D model + H5IO fin(args.fname, H5IO::RDONLY); + const Eigen::VectorX z = read_eigen_vector(fin, "z"); + const Eigen::VectorX vs = read_eigen_vector(fin, "vs"); + const Eigen::VectorX periods = read_eigen_vector(fin, "periods"); + const int nz = static_cast(z.size()); + const int nper = static_cast(periods.size()); + if (vs.size() != nz || nz < 2) { + throw std::runtime_error("z and vs must have the same length (>= 2)"); + } + if (nper < 1) { + throw std::runtime_error("periods must not be empty"); + } + require_strictly_increasing(z, "z"); + require_strictly_increasing(periods, "periods"); + if (periods(0) <= _0_CR) { + throw std::runtime_error("periods must be positive"); + } + if (!vs.allFinite() || (vs.array() <= _0_CR).any()) { + throw std::runtime_error("vs must be finite and positive"); + } + logger.Info(fmt::format("Computing {} dispersion and kernels: {} depth nodes, {} periods", + tag, nz, nper), MODULE_MAIN); + + // dispersion and kernels, vp and rho from empirical relations + auto req = surfker::build_disp_req(z, vs, periods, IFLSPH, iwave_of(args.wave_type), + IMODE, static_cast(args.surf_type)); + const Eigen::VectorX vel = surfker::surfdisp(req); + for (int i = 0; i < nper; ++i) { + if (!std::isfinite(vel(i))) { + throw std::runtime_error(fmt::format( + "NaN/Inf {} velocity at period {:.3f} s", tag, periods(i))); + } + // surfker zero-fills vel from the first failed period on, and its + // kernels there are NaN + if (vel(i) <= _0_CR) { + throw std::runtime_error(fmt::format( + "No {} root found in the fundamental mode at period {:.3f} s", tag, periods(i))); + } + } + surfker::DepthKernel1D K = surfker::depthkernel1d(req); + + // Love-wave kernels have no vp sensitivity + const auto full_or_zero = [&](const Eigen::MatrixX &M) { + if (M.rows() == nper && M.cols() == nz) return Eigen::MatrixX(M); + return Eigen::MatrixX(Eigen::MatrixX::Zero(nper, nz)); + }; + const Eigen::MatrixX sen_vs = full_or_zero(K.sen_vs); + const Eigen::MatrixX sen_vp = full_or_zero(K.sen_vp); + const Eigen::MatrixX sen_rho = full_or_zero(K.sen_rho); + + // vs kernel as combined in preproc::combine_kernels (vs-only parametrisation): + // Rayleigh: K_vs + K_vp * d(vp)/d(vs) + K_rho * d(rho)/d(vp) * d(vp)/d(vs) + // Love: K_vs (SURFATT_tomo keeps no vp or rho kernel for Love waves) + Eigen::MatrixX sen_vs_total = sen_vs; + if (args.wave_type == WaveType::RL) { + for (int k = 0; k < nz; ++k) { + const real_t dab = dalpha_dbeta(req.vs_km_s(k)); + const real_t dra = drho_dalpha(req.vp_km_s(k)); + sen_vs_total.col(k) += sen_vp.col(k) * dab + sen_rho.col(k) * dra * dab; + } + } + require_finite(sen_vs, "sen_vs", periods); + require_finite(sen_vp, "sen_vp", periods); + require_finite(sen_rho, "sen_rho", periods); + + // write results + H5IO fout(args.outfname, H5IO::TRUNC); + fout.write_vector("z", z); + fout.write_vector("periods", periods); + fout.write_vector("vel", vel); + fout.write_vector("vp", req.vp_km_s.head(nz)); + fout.write_vector("rho", req.rho_g_cm3.head(nz)); + fout.write_matrix("sen_vs", sen_vs); + fout.write_matrix("sen_vp", sen_vp); + fout.write_matrix("sen_rho", sen_rho); + fout.write_matrix("sen_vs_total", sen_vs_total); + logger.Info(fmt::format("Kernels written to {}", args.outfname), MODULE_MAIN); + return EXIT_SUCCESS; +} + +} // namespace + + +int main(int argc, char* argv[]) { + auto args = argparse_kernel1d(argc, argv); + + // initialise MPI + Parallel::init(); + auto &mpi = Parallel::mpi(); + + // logger + ATTLogger::init("", /*log_level=*/2, /*console_only=*/true); + auto &logger = ATTLogger::logger(); + + if (mpi.size() > 1) { + logger.Error("SURFATT_kernel1d is not designed for parallel execution. Please run with a single process.", MODULE_MAIN); + mpi.finalize(); + return EXIT_FAILURE; + } + + // Prevent the HDF5 library from printing its own diagnostic stack. Errors + // are caught below and reported once through the application logger. + H5::Exception::dontPrint(); + int status = EXIT_FAILURE; + try { + status = kernel1d(args, logger); + } catch (const H5::Exception &e) { + logger.Error(fmt::format( + "HDF5 error while processing '{}' or '{}': {}", + args.fname, args.outfname, e.getDetailMsg()), MODULE_MAIN); + } catch (const std::exception &e) { + logger.Error(fmt::format( + "Failed to compute 1-D kernels for '{}': {}", args.fname, e.what()), MODULE_MAIN); + } catch (...) { + logger.Error(fmt::format( + "Unknown error while processing '{}'", args.fname), MODULE_MAIN); + } + + mpi.finalize(); + return status; +} diff --git a/src/SURFATT_tomo.cxx b/src/SURFATT_tomo.cxx index f93cbdf..b013053 100644 --- a/src/SURFATT_tomo.cxx +++ b/src/SURFATT_tomo.cxx @@ -34,9 +34,10 @@ int main(int argc, char* argv[]) false ); - // load source-receiver tables into shared memory + // load source-receiver tables into shared memory; tt/vel are observations + // (and so validated) only when inverting for (auto [wt, vt] : IP.data().active_data) - SrcRec::SR(wt, vt).load(IP.data().file_of(wt, vt)); + SrcRec::SR(wt, vt).load(IP.data().file_of(wt, vt), run_mode == INVERSION_MODE); SrcRec::build_stas(); // build model grid diff --git a/src/inversion1d.cpp b/src/inversion1d.cpp index 608bd53..393095c 100644 --- a/src/inversion1d.cpp +++ b/src/inversion1d.cpp @@ -4,7 +4,39 @@ #include "utils.h" #include "input_params.h" #include "src_rec.h" +#include "parallel.h" +#include +#include +#include #include +#include + +namespace { + +// inv1d runs redundantly on every rank with identical inputs, so all ranks +// reach the same failure. Only rank 0 logs; the barrier keeps the other ranks' +// MPI_Abort from killing rank 0 before its diagnostics are flushed. +[[noreturn]] void abort_all_ranks() { + auto& mpi = Parallel::mpi(); + mpi.barrier(); + mpi.abort(EXIT_FAILURE); +} + +// Index of the first row containing inf/NaN, or -1 if all entries are finite. +// (Unlike lpNorm/maxCoeff, allFinite() is guaranteed to see every NaN.) +template +int first_nonfinite_row(const Eigen::DenseBase& m) { + for (Eigen::Index i = 0; i < m.rows(); ++i) { + if (!m.row(i).allFinite()) return static_cast(i); + } + return -1; +} + +std::string data_tag(WaveType wt, SurfType tp) { + return waveTypeStr[static_cast(wt)] + (tp == SurfType::PH ? "_PH" : "_GR"); +} + +} // namespace Inversion1D::Inversion1D(WaveType wavetype) : wavetype_(wavetype) { @@ -13,6 +45,158 @@ Inversion1D::Inversion1D(WaveType wavetype) vs1d.resize(0); } +void Inversion1D::check_inputs(const Eigen::VectorX& zarr) const { + auto& IP = InputParams::IP(); + auto& logger = ATTLogger::logger(); + int n_err = 0; + auto error = [&](const std::string& msg) { + logger.Error(msg, MODULE_INV1D); + ++n_err; + }; + + // 1. Depth grid: finite and strictly increasing (layer thickness > 0) + const int nz = static_cast(zarr.size()); + if (nz < 2) { + error(fmt::format("Depth grid has {} node(s); need >= 2. Check domain.depth_min_max / interval.", nz)); + } + for (int k = 0; k < nz; ++k) { + if (!std::isfinite(zarr(k))) { + error(fmt::format("Depth grid node {} is non-finite ({})", k, zarr(k))); + } else if (k > 0 && !(zarr(k) > zarr(k - 1))) { + error(fmt::format("Depth grid not strictly increasing at node {}: z={} <= z[{}]={}", + k, zarr(k), k - 1, zarr(k - 1))); + } + } + + // 2. Initial Vs: finite and positive at every node + for (int k = 0; k < vs1d.size(); ++k) { + if (!std::isfinite(vs1d(k)) || vs1d(k) <= _0_CR) { + error(fmt::format("Initial Vs invalid at node {} (z={} km): Vs={}. Check model.vel_range.", + k, k < nz ? zarr(k) : NAN, vs1d(k))); + } + } + + // 3. Active src_rec tables for this wave type + for (auto [wt, tp] : IP.data().active_data) { + if (wt != wavetype_) continue; + auto& sr = SrcRec::SR(wt, tp); + const std::string tag = data_tag(wt, tp); + const std::string& file = IP.data().file_of(wt, tp); + + // Periods are validated by SrcRec::load in every run mode; vel only + // when inverting, but inv1d fits vel in forward-only runs too. + if (sr.check_positive(sr.vel, "vel (or dist/tt)") > 0) ++n_err; + + const auto& pinfo = sr.periods_info; + if (pinfo.nperiod <= 0) { + error(fmt::format("[{}] {}: no periods found (empty table?)", tag, file)); + continue; + } + for (int ip = 0; ip < pinfo.nperiod; ++ip) { + if (!std::isfinite(pinfo.meanvel(ip)) || pinfo.meanvel(ip) <= _0_CR) { + error(fmt::format("[{}] {}: mean velocity at period {} s is invalid ({})", + tag, file, pinfo.periods(ip), pinfo.meanvel(ip))); + } + // Distinct values closer than real_t_equal's tolerance are split + // into separate periods by get_periods(); usually a formatting issue. + if (ip > 0 && std::isfinite(pinfo.periods(ip)) && + real_t_equal(pinfo.periods(ip), pinfo.periods(ip - 1))) { + logger.Warn(fmt::format("[{}] {}: near-duplicate periods {} and {} are treated as different periods", + tag, file, pinfo.periods(ip - 1), pinfo.periods(ip)), MODULE_INV1D); + } + } + std::string plist; + for (int ip = 0; ip < pinfo.nperiod; ++ip) plist += fmt::format(" {}", pinfo.periods(ip)); + logger.Debug(fmt::format(" [{}] {} rows, {} periods:{}", tag, sr.n_obs(), pinfo.nperiod, plist), + MODULE_INV1D); + } + + if (n_err > 0) { + logger.Error(fmt::format("1D inversion input check failed with {} error(s); see messages above.", n_err), + MODULE_INV1D); + abort_all_ranks(); + } +} + +void Inversion1D::log_vs_profile(const Eigen::VectorX& zarr, int n_update) const { + auto& logger = ATTLogger::logger(); + // minCoeff/maxCoeff are unreliable with NaN, so take the range over finite + // nodes and count the rest separately. + std::string prof; + real_t vmin = std::numeric_limits::infinity(); + real_t vmax = -vmin; + int n_nonfinite = 0; + for (int k = 0; k < vs1d.size(); ++k) { + prof += fmt::format(" {:.2f}:{:.4f}", zarr(k), vs1d(k)); + if (!std::isfinite(vs1d(k))) { + ++n_nonfinite; + } else { + vmin = std::min(vmin, vs1d(k)); + vmax = std::max(vmax, vs1d(k)); + } + } + logger.Error(fmt::format(" Vs min={:.4f}, max={:.4f} km/s over finite nodes, {} non-finite node(s); " + "profile (z:Vs):{}", vmin, vmax, n_nonfinite, prof), MODULE_INV1D); + + for (int k = 0; k < vs1d.size(); ++k) { + if (!std::isfinite(vs1d(k)) || vs1d(k) <= _0_CR) { + logger.Error(fmt::format(" First non-physical Vs: node {} (z={} km) is {} after {} update(s).", + k, zarr(k), vs1d(k), n_update), MODULE_INV1D); + break; + } + } +} + +void Inversion1D::check_pred_vel(const Eigen::VectorX& pred_vel, + const Eigen::VectorX& periods, + const Eigen::VectorX& zarr, + int iter, WaveType wt, SurfType tp) const { + int first_bad = -1; + for (int ip = 0; ip < pred_vel.size(); ++ip) { + if (!std::isfinite(pred_vel(ip)) || pred_vel(ip) <= _0_CR) { first_bad = ip; break; } + } + if (first_bad < 0) return; + + auto& logger = ATTLogger::logger(); + const real_t bad = pred_vel(first_bad); + if (bad == _0_CR) { + // disper() writes exactly 0 from the first period whose root search fails + logger.Error(fmt::format( + "[{}] iter {}: dispersion solver found no fundamental-mode root from period {} s " + "(index {} of {}); predicted velocities from here on are 0.", + data_tag(wt, tp), iter, periods(first_bad), first_bad, periods.size()), MODULE_INV1D); + } else { + logger.Error(fmt::format( + "[{}] iter {}: predicted velocity at period {} s (index {} of {}) is non-physical ({}).", + data_tag(wt, tp), iter, periods(first_bad), first_bad, periods.size(), bad), MODULE_INV1D); + } + log_vs_profile(zarr, iter); + abort_all_ranks(); +} + +void Inversion1D::check_kernels(const surfker::DepthKernel1D& kernels, + const Eigen::VectorX& periods, + const Eigen::VectorX& zarr, + int iter, WaveType wt, SurfType tp) const { + // Empty matrices (sen_vp for Love waves) have no rows and pass trivially. + int first_bad = -1; + for (const auto* m : {&kernels.sen_vs, &kernels.sen_vp, &kernels.sen_rho}) { + const int ib = first_nonfinite_row(*m); + if (ib >= 0 && (first_bad < 0 || ib < first_bad)) first_bad = ib; + } + if (first_bad < 0) return; + + ATTLogger::logger().Error(fmt::format( + "[{}] iter {}: depth kernel is non-finite at period {} s (index {} of {}){}", + data_tag(wt, tp), iter, periods(first_bad), first_bad, periods.size(), + tp == SurfType::GR + ? "; group kernels re-run disper() at periods perturbed by +-0.5%, " + "which can fail even when the nominal period succeeds." + : "."), MODULE_INV1D); + log_vs_profile(zarr, iter); + abort_all_ranks(); +} + Eigen::VectorX Inversion1D::inv1d( Eigen::VectorX zarr, Eigen::VectorX init_vs @@ -27,12 +211,6 @@ Eigen::VectorX Inversion1D::inv1d( int nz = static_cast(zarr.size()); real_t step_length = IP.inversion().step_length; - real_t sigma = _0_CR; - if (IP.postproc().smooth_method == 0) { - sigma = IP.postproc().sigma[1]; - } else if (IP.postproc().smooth_method == 1) { - sigma = 0.68 * (zarr(zarr.size() - 1) - zarr(0)) / IP.postproc().n_inv_grid[2]; - } int n_active_for_wave = 0; for (const auto& [wt, tp] : IP.data().active_data) { (void)tp; @@ -50,6 +228,16 @@ Eigen::VectorX Inversion1D::inv1d( waveTypeStr[static_cast(wavetype_)]), MODULE_INV1D ); + + // Must precede anything that indexes zarr / vs1d (sigma, min/max below). + check_inputs(zarr); + + real_t sigma = _0_CR; + if (IP.postproc().smooth_method == 0) { + sigma = IP.postproc().sigma[1]; + } else if (IP.postproc().smooth_method == 1) { + sigma = 0.68 * (zarr(zarr.size() - 1) - zarr(0)) / IP.postproc().n_inv_grid[2]; + } logger.Debug( fmt::format(" nz={}, sigma={:.3f}, initial step_length={:.3e}, max_iter={}", nz, sigma, step_length, MAX_ITER_1D), @@ -82,6 +270,7 @@ Eigen::VectorX Inversion1D::inv1d( ); Eigen::VectorX pred_vel = surfker::surfdisp(req); + check_pred_vel(pred_vel, sr.periods_info.periods, zarr, iter, wt, tp); real_t misfit = 0.5 * (pred_vel - sr.periods_info.meanvel).array().square().sum(); logger.Debug( fmt::format(" iter {:3d} | {}_{} misfit={:.6e} (weight={:.3f})", @@ -92,6 +281,7 @@ Eigen::VectorX Inversion1D::inv1d( misfit_total += misfit * IP.data().weights[itype]; surfker::DepthKernel1D kernels = surfker::depthkernel1d(req); + check_kernels(kernels, sr.periods_info.periods, zarr, iter, wt, tp); update.setZero(); auto vp = vs2vp(vs1d); @@ -120,7 +310,19 @@ Eigen::VectorX Inversion1D::inv1d( if (iter > 0 && misfits[iter] > misfits[iter - 1]) { step_length *= IP.inversion().maxshrink; } - update_total = step_length * update_total / update_total.lpNorm(); + const int bad_node = first_nonfinite_row(update_total); + if (bad_node >= 0) { + // data, predictions and kernels are checked above, so this points + // at the smoothing (e.g. sigma <= 0 gives exp(0/0) in gaussian_smooth_1d) + logger.Error(fmt::format("iter {}: model update is non-finite at node {} (z={} km, value {}); " + "check the smoothing sigma.", iter, bad_node, zarr(bad_node), update_total(bad_node)), + MODULE_INV1D); + abort_all_ranks(); + } + const real_t update_norm = update_total.lpNorm(); + if (update_norm > _0_CR) { + update_total = step_length * update_total / update_norm; + } // zero update (perfect fit) would otherwise give 0/0 = NaN logger.Debug( fmt::format("Iteration {}: misfit = {:.6e}, step_length = {:.3e}", iter, misfits.back(), step_length), @@ -143,6 +345,14 @@ Eigen::VectorX Inversion1D::inv1d( } niter = iter; + // The last update is applied after the final forward check, and the result + // may seed the next wave type's inversion; make sure it is still physical. + if (!vs1d.allFinite() || (vs1d.array() <= _0_CR).any()) { + logger.Error("1D inversion produced a non-physical Vs model in its final update.", MODULE_INV1D); + log_vs_profile(zarr, std::min(niter + 1, MAX_ITER_1D)); + abort_all_ranks(); + } + if (niter < MAX_ITER_1D) { logger.Info( fmt::format("1D inversion converged after {} iterations, final misfit={:.6e}", diff --git a/src/src_rec.cpp b/src/src_rec.cpp index 3c2a221..0b16223 100644 --- a/src/src_rec.cpp +++ b/src/src_rec.cpp @@ -3,6 +3,7 @@ #include "utils.h" #include "config.h" +#include #include #include #include @@ -17,11 +18,12 @@ // - Only rank 0 performs file I/O (single reader, deterministic parsing). // - Other ranks obtain the same data via broadcast + shared-memory sync, // avoiding duplicated memory footprint per rank. -void SrcRec::load(const std::string& filepath) +void SrcRec::load(const std::string& filepath, bool check_obs) { auto& mpi = Parallel::mpi(); // auto& IP = InputParams::IP(); auto& logger = ATTLogger::logger(); + filepath_ = filepath; std::vector v_stla, v_stlo, v_evla, v_evlo, v_dist, v_period, v_tt, v_vel, v_weight; @@ -121,6 +123,22 @@ void SrcRec::load(const std::string& filepath) mpi.sync_from_main_rank(vel, n_obs_); mpi.sync_from_main_rank(weight, n_obs_); + // Validate before get_periods(): its std::sort needs NaN-free periods. + // tt/vel are only checked when they are observations; forward runs may + // carry placeholders there. Every rank sees the same shared data, so all + // reach the same verdict; the barrier lets rank 0 flush its log first. + int n_bad_cols = (check_positive(period_all, "period") > 0); + if (check_obs) { + n_bad_cols += (check_positive(tt, "tt") > 0); + n_bad_cols += (check_positive(vel, "vel") > 0); + } + if (n_bad_cols > 0) { + logger.Error(fmt::format("SrcRec::load: {} column(s) of {} contain invalid values; see messages above.", + n_bad_cols, filepath), MODULE_SRCREC); + mpi.barrier(); + mpi.abort(EXIT_FAILURE); + } + // Build derived metadata used by inversion/forward steps: // - event-wise receiver index lists // - per-period statistics (counts and mean velocities) @@ -128,6 +146,25 @@ void SrcRec::load(const std::string& filepath) get_events(); } +// Rows are reported as CSV line numbers: data row index + 2 (header is line 1). +int SrcRec::check_positive(const real_t* col, const std::string& col_name) const { + constexpr int MAX_BAD_ROWS_SHOWN = 10; + std::vector rows; + for (int i = 0; i < n_obs_; ++i) { + if (!std::isfinite(col[i]) || col[i] <= _0_CR) rows.push_back(i); + } + if (rows.empty()) return 0; + + std::string msg = fmt::format("{}: '{}' must be finite and > 0; {} row(s):", + filepath_, col_name, rows.size()); + for (size_t i = 0; i < rows.size() && i < MAX_BAD_ROWS_SHOWN; ++i) { + msg += fmt::format(" line {} ({})", rows[i] + 2, col[rows[i]]); + } + if (rows.size() > MAX_BAD_ROWS_SHOWN) msg += " ..."; + ATTLogger::logger().Error(msg, MODULE_SRCREC); + return static_cast(rows.size()); +} + // Build an event -> observation-index mapping from evtname and distribute // event payloads to destination ranks. // diff --git a/src/surf_grid.cpp b/src/surf_grid.cpp index 1267bd0..6854504 100644 --- a/src/surf_grid.cpp +++ b/src/surf_grid.cpp @@ -1,5 +1,7 @@ #include "surf_grid.h" +#include + namespace { inline int kernel_idx4(const int ix, const int iy, const int iz, const int iper, @@ -7,6 +9,58 @@ inline int kernel_idx4(const int ix, const int iy, const int iz, const int iper, return (((ix * ngrid_j) + iy) * ngrid_k + iz) * nperiod_ + iper; } +// Grid points on this rank where the dispersion solver failed; the 1-D model +// of the first one is kept so that the failure can be reproduced offline. +struct DisperFailure { + int count = 0; + int ix = -1, iy = -1; + Eigen::VectorX vs, vp, rho; + + void add(int ix_glob, int iy_glob, const Eigen::VectorX &vs1d, + const Eigen::VectorX &vp1d, const Eigen::VectorX &rho1d) { + if (count++ > 0) return; + ix = ix_glob; iy = iy_glob; + vs = vs1d; vp = vp1d; rho = rho1d; + } +}; + +// Collective. disper() returns 0 from the first period without a +// fundamental-mode root; such velocities must not reach the eikonal solver, +// and non-finite kernels would turn the model into NaN. If any rank failed, write the +// first failed 1-D model of each rank to +// /disper_fail__rank.txt and abort all ranks. +void abort_on_disper_failure(const DisperFailure &fail, const std::string &what, + const std::string &type_name, + const Eigen::VectorX &periods) { + auto &mpi = Parallel::mpi(); + int n_fail = 0; + mpi.sum_all_all(fail.count, n_fail); + if (n_fail == 0) return; + + auto &mg = ModelGrid::MG(); + const std::string &out_dir = InputParams::IP().output().output_path; + if (fail.count > 0) { + std::ofstream f(fmt::format("{}/disper_fail_{}_rank{:04d}.txt", out_dir, type_name, mpi.rank())); + f << fmt::format("# {} {}: {} grid point(s) failed on this rank, first one at " + "ix={} iy={} (x={:.4f} y={:.4f})\n", + type_name, what, fail.count, fail.ix, fail.iy, + mg.xgrids(fail.ix), mg.ygrids(fail.iy)); + f << "# periods_s:"; + for (int i = 0; i < periods.size(); ++i) f << fmt::format(" {:.8g}", periods(i)); + f << "\n# depth_km vs_km_s vp_km_s rho_g_cm3\n"; + for (int k = 0; k < fail.vs.size(); ++k) + f << fmt::format("{:.8g} {:.8g} {:.8g} {:.8g}\n", + mg.zgrids(k), fail.vs(k), fail.vp(k), fail.rho(k)); + } + ATTLogger::logger().Error(fmt::format( + "{}: {} at {} of {} surface grid points. 1-D models of the failed points " + "are written to {}/disper_fail_{}_rank*.txt", + type_name, what, n_fail, ngrid_i * ngrid_j, out_dir, type_name), MODULE_GRID); + // Let rank 0 flush its log before any rank's MPI_Abort tears the job down. + mpi.barrier(); + mpi.abort(EXIT_FAILURE); +} + } SurfGrid::SurfGrid(WaveType wt, SurfType vt){ @@ -224,6 +278,7 @@ void SurfGrid::fwdsurf(){ const int n_elem = ngrid_i * ngrid_j * nperiod_; std::vector tmp_svel(n_elem, _0_CR); + DisperFailure fail; for (int ix = 0; ix < dcp.loc_nx(); ++ix) { for (int iy = 0; iy < dcp.loc_ny(); ++iy) { const int ix_glob = dcp.loc_I_start() + ix; @@ -248,12 +303,17 @@ void SurfGrid::fwdsurf(){ auto req = surfker::build_disp_req(mg.zgrids, vs1d, vp1d, rho1d, periods, IFLSPH, iwave_of(wt_), IMODE, itype_); Eigen::VectorX svel_point = surfker::surfdisp(req); + if (!svel_point.allFinite() || (svel_point.array() <= _0_CR).any()) { + fail.add(ix_glob, iy_glob, vs1d, vp1d, rho1d); + } for (int iper = 0; iper < nperiod_; ++iper) { const int idx = surf_idx(ix_glob, iy_glob, iper); tmp_svel[idx] = svel_point(iper); } } } + abort_on_disper_failure(fail, "dispersion solver failed (zero or non-finite velocity)", + type_name(), periods); // Reduce into a private buffer first to avoid the shared-memory double-write // problem: svel is MPI shared memory; all node-local ranks point to the same // physical address, so MPI_Allreduce writing directly to svel would accumulate @@ -282,6 +342,7 @@ void SurfGrid::compute_dispersion_kernel() { using MatRM = Eigen::Matrix; + DisperFailure fail; for (int ix = 0; ix < dcp.loc_nx(); ++ix) { for (int iy = 0; iy < dcp.loc_ny(); ++iy) { Eigen::VectorX vs1d(ngrid_k); @@ -310,7 +371,13 @@ void SurfGrid::compute_dispersion_kernel() { } else { kernels = surfker::depthkernelHTI1d(req); } - + // Empty matrices (e.g. sen_vp for Love waves) pass trivially. + if (!kernels.sen_vs.allFinite() || !kernels.sen_vp.allFinite() || + !kernels.sen_rho.allFinite() || !kernels.sen_gc.allFinite() || + !kernels.sen_gs.allFinite()) { + fail.add(dcp.loc_I_start() + ix, dcp.loc_J_start() + iy, vs1d, vp1d, rho1d); + } + // Copy the kernels for this grid point into the corresponding location in the global sensitivity arrays. const int id0 = kernel_idx4(ix, iy, 0, 0, dcp.loc_ny(), ngrid_k, nperiod_); Eigen::Map vs_block(sen_vs_loc.data() + id0, ngrid_k, nperiod_); @@ -329,6 +396,10 @@ void SurfGrid::compute_dispersion_kernel() { } } } + // Group kernels re-run disper() at periods perturbed by +-0.5%, which can + // fail even where fwdsurf() succeeded. + abort_on_disper_failure(fail, "depth kernels are non-finite (dispersion solver failed)", + type_name(), periods); mpi.barrier(); } diff --git a/src/surfker/disper.cpp b/src/surfker/disper.cpp index 90983c0..79c84a5 100644 --- a/src/surfker/disper.cpp +++ b/src/surfker/disper.cpp @@ -564,7 +564,19 @@ static void getsol(double t1, double &c1, double clow, double dc, double cm, else idir = -1; - while (true) { + /* A non-finite model makes every exit test below compare false, so the + * search would never end. Cap the number of steps well above what a + * healthy search needs and report a failure instead. */ + long max_steps = 1000000; + if (std::isfinite(betmx) && dc > 0.0) + max_steps = std::min(max_steps, 10L * (long)((betmx + dc) / dc) + 1000L); + + for (long nstep = 0; ; ++nstep) { + if (nstep >= max_steps) { + iret = -1; + return; + } + double c2; if (idir > 0) c2 = c1 + dc; @@ -690,7 +702,10 @@ std::vector disper(const float *thkm, const float *vpm, const float *vsm for (int i = 0; i < kmax; i++) { c[i] = 0.0; cb[i] = 0.0; } - int ift = 999; + /* First period index at which a lower mode failed; kmax = none yet. + * (Fortran used 999, which aborts the fundamental mode at k=999 + * whenever kmax >= 1000.) */ + int ift = kmax; double del1st = 0.0; /* saved state for getsol direction logic */ for (int iq = 1; iq <= mode; iq++) {