From 40621205e6114025f1e23217dfaad144938be3cc Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 31 Aug 2026 13:33:21 +1200 Subject: [PATCH] Check the velocity model matches the domain before writing e3d.par The downstream half of #131. That change made `create-e3d-par` write `nz + 1`, because `velocity_modelling` generates one more layer than the domain implies: EMOD3D shifts the model down a gridpoint for the free surface and so never reads the last layer (`genmodel.c`). This adds the check that the two halves actually agree. A velocity model whose grid does not match the domain in `e3d.par` does not fail loudly. EMOD3D reads whatever is there and runs on a silently misaligned model. The sizes are fully determined by nx, ny and nz, so compare the three binaries against what the domain implies before templating the parameter file. `padded_nz` is now the single definition of that padded count, used both for the `e3d.par` value and for the check. They must not compute it separately: a check that drifts from the value it is checking is worse than no check at all. A size mismatch is an error. A missing or unreadable file is not: `create-e3d-par` is routinely run in a container where the velocity model paths are only being templated and nothing is on disk yet, so raising there would break workflows that are behaving correctly. That is handled per file rather than around the loop, so an absent first file does not stop the remaining two from being checked -- otherwise a genuinely mis-sized model can hide behind a missing one. Bumps the `velocity-modelling` floor to 2026.8.1, the release that generates the padding row. #131 made the workflow depend on that behaviour without raising the floor; the lockfile happened to resolve high enough, but `pip install workflow` -- which is what the container does -- reads the floor, not the lock. Co-Authored-By: Claude Opus 5 --- pyproject.toml | 2 +- uv.lock | 2 +- workflow/scripts/create_e3d_par.py | 99 +++++++++++++++++++++++++++++- 3 files changed, 99 insertions(+), 4 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index bd8f48ef..15d482eb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,7 +12,7 @@ dynamic = ["version"] dependencies = [ # UCGMSim Dependencies "im-calculation>=2025.12.5", - "velocity-modelling>=2026.2.1", + "velocity-modelling>=2026.8.1", "nshmdb>=2026.09.1", "oq_wrapper>=2025.12.3", "qcore-utils>=2025.12.2", diff --git a/uv.lock b/uv.lock index 104bbcab..f165710f 100644 --- a/uv.lock +++ b/uv.lock @@ -3548,7 +3548,7 @@ requires-dist = [ { name = "typer" }, { name = "types-geopandas", marker = "extra == 'types'" }, { name = "types-requests", marker = "extra == 'types'" }, - { name = "velocity-modelling", specifier = ">=2026.2.1" }, + { name = "velocity-modelling", specifier = ">=2026.8.1" }, { name = "xarray", extras = ["io"] }, ] provides-extras = ["test", "types", "dev"] diff --git a/workflow/scripts/create_e3d_par.py b/workflow/scripts/create_e3d_par.py index 1d5a5f34..997352a2 100644 --- a/workflow/scripts/create_e3d_par.py +++ b/workflow/scripts/create_e3d_par.py @@ -47,6 +47,33 @@ app = typer.Typer() +def padded_nz(domain_parameters: DomainParameters, resolution: Resolution) -> int: + """The number of z gridpoints EMOD3D reads, including the padding row. + + `velocity_modelling` generates one more layer than the domain implies, + because EMOD3D shifts the model down a gridpoint for the free surface and + so never reads the last layer (`genmodel.c`). EMOD3D has to be told the + padded count, and the velocity model files are sized by it. + + Both the `e3d.par` `nz` and `check_domain_against_velocity_model` come + through here. They must not compute it separately: a check that drifts + from the value it is checking is worse than no check at all. + + Parameters + ---------- + domain_parameters : DomainParameters + The realisation domain parameters. + resolution : Resolution + The simulation resolution. + + Returns + ------- + int + The z gridpoint count including the free-surface padding row. + """ + return domain_parameters.nz(resolution.resolution) + 1 + + def emod3d_domain_parameters( resolution: Resolution, domain_parameters: DomainParameters, @@ -68,8 +95,7 @@ def emod3d_domain_parameters( nx = domain_parameters.nx(resolution.resolution) ny = domain_parameters.ny(resolution.resolution) - # nz + 1 for consistency with the velocity model - nz = domain_parameters.nz(resolution.resolution) + 1 + nz = padded_nz(domain_parameters, resolution) return { "nx": nx, "ny": ny, @@ -232,6 +258,70 @@ def format_as_emod3d_value(value: float | str | Path) -> str: return str(value) +def check_domain_against_velocity_model( + domain_parameters: DomainParameters, + resolution: Resolution, + parameters: EMOD3DParameters, + velocity_model_ffp: Path, +) -> None: + """Validate that generated velocity model files match the expected domain size. + + Computes the expected file size for the pmodfile, smodfile, and dmodfile + from the domain's nx, ny, nz (with the EMOD3D free-surface padding row added + to nz), and compares it against the actual size of each file on disk. A + mismatch is an error, but a missing file or unreadable file is not + considered an error because this code is often run in a container with the + paths only used for templating and hence it would fail on many workflows. + + Parameters + ---------- + domain_parameters : DomainParameters + Domain parameters used to compute the expected nx, ny, nz grid + dimensions at the given resolution. + resolution : Resolution + Resolution at which to evaluate the domain's grid dimensions. + parameters : EMOD3DParameters + EMOD3D parameters providing the pmodfile, smodfile, and dmodfile + filenames to check. + velocity_model_ffp : Path + Directory containing the velocity model files to validate. + + Returns + ------- + None + + Raises + ------ + RuntimeError + If a velocity model file's size does not match the expected size + computed from the domain parameters. + + """ + nx = domain_parameters.nx(resolution.resolution) + ny = domain_parameters.ny(resolution.resolution) + nz = padded_nz(domain_parameters, resolution) + + expected_file_size = nx * ny * nz * np.float32().nbytes + for filename in [parameters.pmodfile, parameters.smodfile, parameters.dmodfile]: + velocity_model_file = velocity_model_ffp / filename + try: + file_size = velocity_model_file.stat().st_size + except OSError as e: + # Unreadable is not a mismatch. This stage is routinely run in a + # container where the velocity model paths are only being + # templated and nothing is on disk yet. Handled per file rather + # than around the loop, so one absent file does not stop the + # others from being checked. + print( + f"WARNING: could not validate domain parameters against velocity model supplied:\n{e}" + ) + continue + if file_size != expected_file_size: + raise RuntimeError( + f"Velocity model file {velocity_model_file} does not have the expected size (expected: {expected_file_size}, found: {file_size})" + ) + + @cli.from_docstring(app) def create_e3d_par( realisation_ffp: Path, @@ -270,6 +360,10 @@ def create_e3d_par( emod3d_parameters = EMOD3DParameters.read_from_realisation_or_defaults( realisation_ffp, metadata.defaults_version ) + check_domain_against_velocity_model( + domain_parameters, resolution, emod3d_parameters, velocity_model_ffp + ) + e3d_par_values = ( emod3d_parameters.to_dict() | emod3d_domain_parameters(resolution, domain_parameters) @@ -291,4 +385,5 @@ def create_e3d_par( for key, value in e3d_par_values.items() ) ) + realisations.append_log_entry(realisation_ffp)