Skip to content

Port srf2stoch to Python - #128

Merged
lispandfound merged 13 commits into
pegasusfrom
srf2stoch_port
Sep 25, 2026
Merged

lispandfound merged 13 commits into
pegasusfrom
srf2stoch_port

Conversation

@lispandfound

@lispandfound lispandfound commented Jul 29, 2026 •

Copy link
Copy Markdown
Contributor

This PR removes the srf2stoch binary from the workflow and replaces it with a pure python implementation that does the same thing with some important improvements:

  1. Angle averaging: A number of places where rob calculates mean of angles are replaced with a circular mean.
  2. Testing: Proper tests of the stoch code asserting the important properties: cell registration, moment conservation,
  3. Overhang: srf2stoch rounds nx = (int)(len/dx + 0.5) to the nearest cell but still writes the requested dx to the header. The HF code therefore rebuilds a plane nxdx long instead of len, stretching or squeezing the slip to fit, and moment is off by nxdx/len. For example, a 3 km plane at dx=2 becomes 4 km with a third more moment. The Python port averages by fractional overlap, so moment is conserved exactly,
  4. SIGFPE: The C code would crash if the dx in the SRF and the provided stoch dx didn't place nice. The old wrapper worked around this with the min_length/2 clamp and a special case for point sources. In the Python port, these arithmetic issues are avoided.
  5. Speed and memory: srf2stoch upsamples both grids to a common multiple, so memory and time grow by nx/gcd(nstk,nx) × ny/gcd(ndip,ny). On a 1373×199 plane at 2 km (factors 69 and 10), the C binary took 1.85 s and peaked at 2.3 GB. The Python sparse kernel took 0.014 s to convert. This scales to runs we made on the alpine fault which struggled to execute on normal HPC nodes for some clusters.

Copilot AI review requested due to automatic review settings July 29, 2026 04:34
@gemini-code-assist

Copy link
Copy Markdown
Contributor

Caution

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

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

This PR ports stoch generation from the external srf2stoch binary to a pure-Python implementation, integrating SRF→Stoch conversion directly into the generate-stoch CLI and simplifying container tooling.

Changes:

  • Implement SRF→Stoch conversion in Python using NumPy/SciPy (box-averaging slip, rupture time, and slip-weighted rise).
  • Remove the srf2stoch binary dependency from both the CLI and the container build.
  • Bump source_modelling minimum/version pin to 2026.07.3 / 2026.7.3.

Reviewed changes

Copilot reviewed 3 out of 4 changed files in this pull request and generated 3 comments.

File Description
workflow/scripts/generate_stoch.py Replaces external srf2stoch invocation with in-process Python conversion and updates CLI/docs accordingly.
container/runner.def Stops building srf2stoch in the runner image.
pyproject.toml Updates source_modelling minimum version requirement.
uv.lock Updates lockfile entries (incl. source-modelling 2026.7.3 and related metadata).

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

Comment thread workflow/scripts/generate_stoch.py Outdated
Comment thread workflow/scripts/generate_stoch.py Outdated
Comment thread workflow/scripts/generate_stoch.py
@lispandfound
lispandfound changed the base branch from pegasus to nzvm/domain-fault-buffer August 31, 2026 09:09
@lispandfound
lispandfound force-pushed the srf2stoch_port branch 2 times, most recently from 52a9ee3 to b4d0555 Compare September 2, 2026 02:52
@lispandfound
lispandfound force-pushed the srf2stoch_port branch 2 times, most recently from eee3b7f to eb147f0 Compare September 22, 2026 21:51
@lispandfound
lispandfound force-pushed the srf2stoch_port branch 3 times, most recently from 65eb54c to f943533 Compare September 22, 2026 22:38

@lispandfound lispandfound left a comment

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Code review by Claude Code (run locally, since the review action is currently broken). 8 finding(s).

Comment thread workflow/scripts/generate_stoch.py Outdated
Comment thread workflow/scripts/generate_stoch.py Outdated
]
)
srf_file = srf.read_srf(srf_ffp)
dx = stoch_config.stoch_dx

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Two special cases were removed: point sources used the SRF grid (len/nstk, wid/ndip), and other planes were capped at min(stoch_dx, min_length/2). Every plane now uses the 2 km stoch_dx/stoch_dy, even when the plane is smaller than one cell.

Failure scenario: A point-source SRF, or a 0.5 km segment in a multi-fault rupture, becomes a single 2x2 km stoch cell. Moment is conserved, but HF reconstructs a plane up to 4x larger than the real source. Combined with the dip padding, the source centre also sits about dy/2 * sin(dip) deeper than it should. If this is intended, it needs a comment. If not, keep a per-plane dx for sub-cell planes.

@lispandfound lispandfound Sep 24, 2026 •

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The stoch dx behaviour is expected. The old hack to adjust the stoch dx was not desired it was an artifact of the limitations in srf2stoch. As for the point source case: yes I guess this could be an issue but I don't think it will matter, you can just adjust the stoch dx/dy if you want that.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Cell registration is a bug, which I have addressed.

Comment thread workflow/scripts/generate_stoch.py Outdated
Comment thread workflow/scripts/generate_stoch.py
Comment thread workflow/scripts/generate_stoch.py Outdated
Comment thread tests/test_generate_stoch.py Outdated
Comment thread workflow/scripts/generate_stoch.py Outdated
Comment thread workflow/scripts/generate_stoch.py Outdated
lispandfound and others added 13 commits September 25, 2026 11:43
Replaces the `srf2stoch` binary with an in-process conversion, and drops
it from the container build.

The conversion is a box average expressed as a sparse fractional-overlap
matrix: row j of the kernel holds the overlap weights between coarse bin
j and the fine cells it spans. The special case where the two grids span
exactly the same length is what srf2stoch.c implements; the sparse form
generalises it without materialising the empty overlaps, and handles the
padded case where they do not.

That padded case is why this was worth porting. The HF code demands a
single dx/dy across every SRF segment, and `nx = ceil(len / dx)` means
the coarse grid is generally *longer* than the plane it covers. The
overhang is split evenly between the two ends, so the grids stay
concentric -- the stoch format records only a centre point and an
`nx * dx` extent, so off-centre padding would sit the slip distribution
in the wrong place on the plane HF reconstructs. Edge bins then carry
weights summing to the covered fraction rather than to 1, which is what
conserves total slip*area rather than cell value.

Quantities that are not spread over a cell are handled separately.
Rupture time is a time, so it is divided by the covered fraction; every
cell is partially covered because nx and ny round up, so this never
divides by zero. Rise time is slip-averaged, and there the coverage
factor cancels, so it must *not* be divided out again.

Rake is averaged as a circular mean weighted by slip -- the arithmetic
mean of -179 and 179 degrees is 0, which is the opposite of the right
answer. A plane with no slip anywhere has no slip-weighted mean, so it
falls back to an unweighted average.

This changes every stoch file, and therefore every high-frequency
result. The previous binary's outputs were not moment-conserving for
segments whose length was not an exact multiple of the stoch dx.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>

@AndrewRidden-Harper AndrewRidden-Harper left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

very nice docstrings and comments

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants