Fix velocity model node sampling

[lispandfound]

The velocity model has been incorrectly sampling at the midpoint between each velocity layer (50m, 150m, 250m, etc in a 100m simulation). The correct sampling should be based on the boundaries for most points (e.g. 0m, 100m, 200m, 300m) with each velocity a tributary for the model region [0-h/2m], [h/2-3h/2], [3h/2-5h/2], and so on.

We make the modelling decision to sample the top bounding point at h/4m (e.g. 25m in a 100m simulation) because we believe the velocities at this depth represent the region of 0-h/2m better than the velocities directly at the surface. As a consequence, there is no longer a fixed dz in the velocity model anymore. I have updated the most important places for this including the velocity model compression tool which now reports the correct z values for each gridpoint.

[gemini-code-assist]

Caution

The consumer version of Gemini Code Assist on GitHub has been sunset. All code review activity has officially ceased.

[lispandfound]

Ruff checks aren't passing. I'll make a separate auto-generated PR to resolve those separately.

[Copilot]

Pull request overview

Adjusts vertical (z) sampling in the velocity model to avoid midpoint-based sampling and better align sampling points with layer boundaries, including a special-case near-surface sampling offset. This also updates downstream utilities that depend on z spacing/coordinates.

Changes:

• Updates global mesh z-coordinate generation to sample at boundaries (with a top-node offset of h_depth/4).
• Updates VS averaging to use trapezoidal integration over (potentially) non-uniform z spacing.
• Updates 1D profile generation and compressed VM output to reflect the new z sampling scheme.
• Removes the .github/workflows/pr-review.yml workflow.

Reviewed changes

Copilot reviewed 5 out of 5 changed files in this pull request and generated 2 comments.

Show a summary per file

File	Description
velocity_modelling/geometry.py	Changes how global_mesh.z is constructed to align sampling with layer boundaries and apply a top-node offset.
velocity_modelling/threshold.py	Replaces constant-dz harmonic-mean accumulation with trapezoidal integration for VS metrics.
velocity_modelling/scripts/generate_1d_profiles.py	Updates depth extents handling and adjusts 1D profile layer thickness logic for the new sampling approach.
velocity_modelling/tools/compress_vm.py	Updates compressed VM depth coordinate construction to reflect the new top-node offset and formatting cleanups.
.github/workflows/pr-review.yml	Removes the AI Code Review workflow.

Comments suppressed due to low confidence (1)

.github/workflows/pr-review.yml:1

• This PR removes .github/workflows/pr-review.yml (AI Code Review) but the PR description is focused on velocity model z-sampling changes. Either restore this workflow or update the PR description/rationale so the removal is intentional and doesn’t surprise maintainers relying on it.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

Pull request overview

Copilot reviewed 5 out of 5 changed files in this pull request and generated 1 comment.

Comments suppressed due to low confidence (1)

velocity_modelling/geometry.py:1130

• This new z sampling changes the meaning of extent_zmin/extent_zmax without updating how nz is computed/validated. With the current nz = int((zmax - zmin)/h_depth + 0.5), h_depth * np.arange(nz) + zmin ends at zmin + (nz-1)h_depth (typically zmax - h_depth when the range is an exact multiple), so the deepest z-plane no longer reaches the configured extent_zmax (e.g., VS30 grid ends at 29.5 m instead of ~30 m). This likely impacts any code assuming the last z sample corresponds to the requested max depth.

    global_mesh.z = -1000.0 * (h_depth * np.arange(nz) + zmin)
    global_mesh.z[0] -= h_depth / 4 * 1000.0

[Copilot]

Pull request overview

Copilot reviewed 5 out of 5 changed files in this pull request and generated 1 comment.

Comments suppressed due to low confidence (1)

velocity_modelling/scripts/generate_1d_profiles.py:220

• dep_top is never updated inside the loop, so thickness = abs(dep_bot - dep_top) is computed relative to the surface for every layer after the first. This produces cumulative thicknesses (increasing each row) instead of per-layer thicknesses.

                thickness = abs(dep_bot - dep_top)
                f.write(
                    f"{thickness / 1000:.3f} \t {qualities_vector.vp[i]:.3f} \t "
                    f"{vs:.3f} \t {qualities_vector.rho[i]:.3f} \t "
                    f"{qp:.3f} \t {qs:.3f}\n"
                )

[lispandfound]

Tracing the logic for this +1 padding:

1. In EMOD3D we specify model_style=1. This causes the code (genmodel.c:2908) to use the path reedmedslice which reads the media from the files we supply.
2. Inside reedmedslice there is a loop (genmodel.c:1008) that iterates over the velocity model to supply material values to the finite difference code. When free surface (fs) is enabled (genmodel.c:1005) it shifts it's entire search up one z-slice meaning the last layer of the velocity model read is nz - 2 not nz - 1.

shft = 0;
if(fs)      /* shift model down one grid point for free surface */
   shft = 1;
for(iz=nz-1;iz>=1;iz--) // note index pointer
   {
   for(ix=0;ix<nx;ix++)
      {
      i = iz*nx + ix;
      ip = (iz-shft)*nx + ix; // shft = 1 means that that we are never reading the nz - 1 layer. When iz = nz - 1, iz - shft = *nz - 2*.
      lam2mu[i] = a[ip]*a[ip]*rho[ip];   // a/b/rho are the raw pmodfile/smodfile/dmodfile buffers
      ...

3. This means that, if nz = 3 and we supply velocities for 25m, 100m, 200m, the last row (200m) is not read, instead we get velocities read from the 25m layer and 100m layer. What happens to the other layer? We lost it when EMOD3D adds in the ghost layer for the free surface boundary condition. The nz we supply to EMOD3D includes the ghost layer
4. So the fix is to here add one to nz, and then to add 1 to nz downstream in create_e3d_par.py in the workflow. This is the aim of 29f2347. I've included a description in the commit (not that anybody ever reads those).

[lispandfound]

Regrettably I have had to remove several tests that won't work because they are golden tests against NZVM output. This strategy made sense back when "behave like the old C binary exactly" was the goal. But we have moved on from that understanding now and these tests are no longer useful. This should be a lesson for the SW team more broadly: structuring code in such a way that each piece is independently testable from properties makes it easier to validate correctness and easier to refactor tests. As it is I see no way to recover the tests that previously existed other than to create another set of golden tests that would then pin our understanding to this PR (which is inevitably also wrong as all things are in the limit).

[AndrewRidden-Harper]

Note

This review was written by Claude Code (Opus 5), run locally by @AndrewRidden-Harper against a full nzcvm_data checkout. The end-to-end runs, the comparisons against main, and the verification of each suggested fix below were all actually executed rather than inferred from reading the diff.

Reviewed the sampling change and ran it end-to-end against nzcvm_data on a local checkout, comparing against main.

The design is right. I checked the padding-row reasoning against EMOD3D/Emod3d/V3.0.13/genmodel.c: with fs=1, reedmedslice (L1008-L1013) reads file row iz-shft for grid point iz, and adj_medsten_fs (L104) applies the same +fs index shift. So grid point iz is physical depth (iz-1)*h, file row j lands at depth j*h, and row nz-1 genuinely is never read. nz+1 plus the h/4 top node is self-consistent, and it lines up with the nz+1 and file-size check in the companion workflow PR - my 2x2x6 run produced exactly the 96-byte files that check expects. compress_vm is correct too: depths come out [0.25, 1, 2, 3, 4, 5] km and round-trip through netCDF intact.

Two blockers in the implementation, though:

1. generate-3d-model raises on every config - the +1 on nz isn't matched by the validator in gen_full_model_grid_great_circle.
2. Vs30/Vs500 silently return empty for every station - get_depth_parameters still carries the half-spacing zmin offsets that existed only to cancel the removed 0.5*h_depth term, so the top node lands above ground and nan propagates through np.trapezoid.

Neither is visible to CI as configured - see the Jenkinsfile note. I verified fixes for both, plus the three smaller issues; details inline.