Green Function DB computation using HDF5#878
Closed
lsawade wants to merge 21 commits into
Closed
Conversation
Implement configuration, force detection, and Butterworth-filtered STF for the Green function (GF) database (Sawade et al. 2025). This enables reciprocal simulations where stations are placed as force sources and displacement fields are recorded at earthquake element locations. Stage 1 — Configuration: - Add GF parameters to Par_file (GF_DATABASE_ENABLED, GF_DATABASE_PATH, GF_SUBSAMPLE_STEP, GF_BUFFER_SIZE, GF_NEIGHBOR_SHELLS) - Read, validate, and broadcast parameters across MPI ranks - Add GF constants to constants.h.in Stage 2 — Force Detection: - Add green_function_par module for GF runtime variables - Auto-detect force component (N/E/Z) from FORCESOLUTION direction vectors - Parse station identity (NET.STA) from FORCESOLUTION header - Validate NSOURCES=1 and USE_FORCE_POINT_SOURCE when GF enabled - Fix backward-compatible FORCESOLUTION header parsing in get_force.f90 and get_event_info.f90 (skip unused integer ID without breaking e_n) Stage 3 — STF Generation: - Implement Butterworth lowpass filter via bilinear transform (order 4, cascaded second-order sections, zero-phase forward+backward filtering) - Precompute Gaussian STF filtered at Nyquist of subsampled output - Override get_stf_viscoelastic with precomputed gf_stf(it) array - Write STF to OUTPUT_FILES/gf_stf.txt for cross-validation - Set t0 buffer to 15*hdur to accommodate filter ringing - Add unit test (test_gf_stf.f90) and cross-validation against scipy (relative error 1e-10)
…elements - New green_function_locate.F90: gf_read_locations() parses lat/lon/depth_km from DATA/GF_LOCATIONS; gf_locate_elements() converts geographic coords to Cartesian (with ellipticity/topography), uses KD-tree + locate_MPI_slice to find containing elements across MPI ranks, then deduplicates with O(n) boolean mask per rank. - Errors if GF_LOCATIONS is missing or empty when GF_DATABASE_ENABLED=.true. - Hooked into setup_sources_receivers.f90 before xadj/adjncy deallocation (required for Stage 5 neighbor expansion). - Updated green_function_par.F90 with location arrays and rules.mk with build target.
green_function_expand.F90 expands tagged elements by GF_NEIGHBOR_SHELLS layers using xadj/adjncy CSR adjacency (local) and ibool_interfaces exchange (cross-MPI). Shared gf_report_tagged_elements() gathers and prints per-rank element counts and IDs on rank 0.
Compute 63-bit Morton codes from element center GLL node coordinates
(21-bit quantization, bit-interleaved via spread_bits). Create per-element
directories under GF_DATABASE_PATH/elements/{morton_hex}/. Write
centroids.bin (binary, 32 bytes/record: int64 morton + 3x float64
centroid) and manifest.csv gathered from all MPI ranks via send_cr/recv_cr.
New files:
- green_function_morton.F90: gf_spread_bits, gf_morton_encode,
gf_compute_morton_codes, gf_create_directories, gf_write_manifest
- prepare_green_function_storage.f90: central wrapper consolidating all
GF preparation (STF + Morton + dirs + manifest) into a single call
from prepare_timerun
- utils/green_function/verify_morton.py: cross-validates Fortran Morton codes
against Python reimplementation
- utils/green_function/decode_morton.py: decodes Morton hex back to approximate
coordinates (lat/lon/r)
Modified files:
- green_function_par.F90: add gf_morton_codes, gf_morton_hex, gf_center_xyz
- prepare_timerun.F90: replace two GF calls with single
call prepare_green_function_storage()
- green_function_stf.F90: remove dead prepare_greenfunction_stf wrapper
- rules.mk: add green_function_morton and prepare_green_function_storage
Tested with 98 elements on 4 MPI procs. All 98 Morton codes verified
bit-for-bit against Python implementation. Max quantization error 1.4e-6.
Add green_function_io.F90 with serial HDF5 file creation, chunked 6D datasets (force_comp, disp_comp, NGLLX, NGLLY, NGLLZ, nt_sub), and buffered hyperslab writing. Each MPI rank independently manages its local tagged elements. Files are create-or-open so three sequential force component runs (N/E/Z) fill the same dataset without overwriting. Buffer layout (3,5,5,5,bufsize,nelem) uses element-last ordering for contiguous memory access per element. Uses c_loc for rank-agnostic h5dwrite_f calls. Fill value 0 ensures unwritten components are zero. Hooked into iterate_time.F90 (subsampled extraction) and finalize_simulation.F90 (final flush + cleanup). There is still the issue of the metadata, that is topography etc.
Add per-element coordinates.h5 (GLL node coords), shared mesh_info.h5 (simulation params, topography, ellipticity), and per-station metadata (lat/lon/depth, STF). Track completion via computed_N/E/Z/ALL attributes on element station files to enable skip-on-rerun (GF_OVERWRITE parameter). Reverse rotation for reciprocity, multiply displacement by scale_displ for float32 precision, and fix station name parsing from FORCESOLUTION. Add Python post-processing: gf_build_manifest.py (centroids + completion summary) and read_centroids.py (binary reader example).
…ude the validation script plot.
Contributor
Author
|
@danielpeter I didn't try to compile without HDF5, so the CI errors are likely about that. Compiled fine on cluster and mac but with hdf5 enabled. |
lsawade
commented
Jun 26, 2026
| @@ -0,0 +1,11 @@ | |||
| FORCE IU.SJG | |||
| @@ -13,10 +13,10 @@ SAVE_FORWARD = .false. # save last frame of forward simulat | |||
| NCHUNKS = 1 | |||
| @@ -442,7 +442,7 @@ endif | |||
| @COND_ASDF_FALSE@ASDF = no | |||
Member
|
hi Lucas, that will be a great contribution! can you try fix the compilation errors before merging? that is, add a stub function for i'll wait with merging to give you time to add more commits to this PR. |
Member
|
i'll close and wait for the updated one :) |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This PR implements the writing of an database that is element and station-wise!
The most involved edits were adding signal processing for pre simulation low pass filtering of the source time function, to enable subsampling and the two
snakemakeworkflows to run the forward and reciprocal workflows and run an interpolation to test the output elements.The workflow works as follows:
EXAMPLES/green_function_database/regional/db_base/DATA/FORCESOLUTIONNote that the component is set to a force in one direction. if more than one component are non-zero the code will error outEXAMPLES/green_function_database/regional/db_base/DATA/GF_LOCATIONSOf course there are more details to it. But that could be added through the documentation in a follow up PR.
@danielpeter Note how I'm currently setting the
GRAVITY = .false.because of the mismatch. You will see that the traces with gravity off only mismatch at the end, which is an interpolation artifact.To run the snakemake workflow see
EXAMPLES/green_function_database/README.mdIt's pretty straight forward.