Quickstart: QSEIS2025#

This example computes a small Green’s-function library, synthesizes displacement and writes a figure. It uses the AK135 model content included with pygrnwang, one source depth, one receiver depth and three distances. Computation is serial. The source is 10 km deep, receivers are at the surface and distances are 300, 600 and 900 km. The native library uses a 4 s sampling interval, a 0.125 Hz Nyquist limit and a 4092 s window containing 1024 samples, with the flat-Earth transformation enabled. After synthesis, the example saves 0–1020 s inclusive, giving 256 samples without changing the underlying library. The layered calculation uses the first 24 numeric model rows, down to 809.5 km. Constant Qp=600 and Qs=300 are illustrative tutorial choices, not the full AK135-F attenuation model. The source is a normalized 64 s squared half-sinusoid moment-rate pulse. The mechanism is strike 30°, dip 45°, rake 90°, scaled to M0 = 10^15 N m.

1. Prepare the environment#

Follow installation, including cloning the repository to obtain examples/. Run the following commands from the repository root in the activated environment.

2. Calculate displacement#

python examples/qseis2025.py --regional --output-dir examples/output/qseis2025-regional

On Windows with Conda, the equivalent non-interactive command is:

conda run -n pygrnwang python examples/qseis2025.py --regional --output-dir examples/output/qseis2025-regional

The script prepares a local model, writes solver input, installs the 64 s source samples, runs QSEIS2025 and converts the library. It then reads the selected source mechanism with output_type="disp" and saves displacement curves. It checks that arrays have the expected components and contain finite values. All generated files stay below the selected output directory; no files from test/ are required. The solver run takes a few minutes.

3. Inspect the result#

Look for disp.png, disp.npz, source_time_function.json and summary.json in the output directory. The saved displacement array has shape (3, 3, 256): three distances, three components and samples from 0 to 1020 s inclusive. The summary records the environment, calculation time, array dimensions and output size. The underlying library retains its full 1024-sample solver output, geometry metadata and converted binary arrays.

QSEIS2025 example displacement traces at 300, 600 and 900 km from a small AK135 Green's-function library.

Validated QSEIS2025 displacement at 300, 600 and 900 km over 0–1020 s since source origin. Its model, mechanism and numerical choices are shown in the script below.#

The vector reader uses east, north, up when rotate=True. Displacement is reported in metres for the moment specified by the script. The plotting time axis must be interpreted with the example’s reduction and sampling settings; it is not automatically a P-relative axis. See scientific conventions before changing these settings.

The source uses wavelet_type=0 with 1024 custom moment-rate samples. The helper precompensates QSEIS’s numerical damping, so the effective pulse has unit area and a 32 s centroid; see the regional tutorial. The library stores rate kernels, and the reader integrates them once for output_type="disp".

4. Add strain and stress#

Use a separate directory because the output flags change the computed library:

python examples/qseis2025.py --regional --observables all --output-dir examples/output/qseis2025-regional-tensors

The script requests displacement, strain and stress directly from the reader and saves 0–1020 s of each. The saved strain/stress arrays have shape (3, 6, 256). With geographic rotation enabled, symmetric tensors are stored as [EE, EN, EU, NN, NU, UU]; U is the same upward vertical component called Z elsewhere in the code. Strain is dimensionless and stress is in pascals for the chosen source moment. These curves are a small workflow example, not a convergence study.

5. Reuse a finished example#

For the same model, geometry and output flags:

python examples/qseis2025.py --regional --output-dir examples/output/qseis2025-regional --reuse

Use a fresh output directory after changing calculation parameters. Reuse reloads the existing model and data and writes summary-reuse.json; it does not prove that existing files match newly selected scientific settings.

The complete script#

The executable script is included directly here, so the documentation and the tested example share one source.

  1"""Build QSEIS2025 introductory or regional traces; optionally include tensors."""
  2from pathlib import Path
  3
  4from common import (MECHANISM, MOMENT_NM, REGIONAL_SAMPLING_INTERVAL_S, REGIONAL_STF, finish, parser_for, prepare,
  5                    require_library_settings, save_waveforms)
  6from pygrnwang.create_qseis2025_bulk import (
  7    pre_process_qseis2025, create_grnlib_qseis2025_sequential)
  8from pygrnwang.read_qseis2025 import get_outfile_name_list, seek_qseis2025
  9from source_time_function import prepare_qseis_stf, validate_qseis_stf
 10from spectral_settings import qseis_spectral_settings
 11
 12
 13def main():
 14    parser = parser_for("qseis2025")
 15    parser.add_argument("--regional", action="store_true",
 16                        help="Use 300/600/900 km, a 64 s wavelet and Earth flattening")
 17    parser.add_argument("--point-source", action="store_true",
 18                        help="Regional control run: disable the default Gaussian spatial source smoothing")
 19    parser.set_defaults(output_dir=None)
 20    parser.add_argument("--observables", choices=("disp", "all"), default="disp",
 21                        help="all also computes strain and stress")
 22    args = parser.parse_args()
 23    if args.point_source and not args.regional:
 24        parser.error("--point-source requires --regional")
 25    source_radius_ratio = 0.0 if args.point_source else 0.05
 26    if args.output_dir is None:
 27        directory = "qseis2025-regional" if args.regional else "qseis2025"
 28        if args.point_source:
 29            directory += "-point-source"
 30        args.output_dir = Path(__file__).resolve().parent / "output" / directory
 31    output, library, model, report, started = prepare(args, "QSEIS2025")
 32    dt, window = (REGIONAL_SAMPLING_INTERVAL_S, 4092.0) if args.regional else (0.5, 127.5)
 33    output_end = 1020.0 if args.regional else 100.0
 34    output_samples = int(round(output_end / dt)) + 1  # Include the final sample.
 35    native_samples = int(round(window / dt)) + 1
 36    distances = [300.0, 600.0, 900.0] if args.regional else [30.0, 60.0, 90.0]
 37    wavelet_duration = 16 if args.regional else 4
 38    wavelet_type = 0 if args.regional else 2
 39    source_time_function = None
 40    if not args.reuse:
 41        pre_process_qseis2025(
 42            processes_num=1, path_green=library, event_depth_list=[10.0],
 43            receiver_depth_list=[0.0], dist_range=[distances[0], distances[-1]],
 44            delta_dist=distances[0],
 45            N_each_group=3, time_window=window, sampling_interval=dt, source_radius_ratio=source_radius_ratio,
 46            output_observables=([1, 0, 1, 1, 0] if args.observables == "all"
 47                                else [1, 0, 0, 0, 0]),
 48            wavelet_type=wavelet_type, wavelet_duration=wavelet_duration, time_reduction_velo=0,
 49            flat_earth_transform=args.regional, path_nd=model, earth_model_layer_num=24,
 50        )
 51        if args.regional:
 52            source_time_function = prepare_qseis_stf(
 53                library, duration_s=64.0, samples=1024)
 54        create_grnlib_qseis2025_sequential(library, remove_pd=False)
 55    require_library_settings(
 56        library, event_depth_list=[10.0], receiver_depth_list=[0.0],
 57        grn_dist_range=[distances[0], distances[-1]], grn_delta_dist=distances[0],
 58        sampling_interval=dt, time_window=window, sampling_num=native_samples,
 59        wavelet_type=wavelet_type, wavelet_duration=wavelet_duration, time_reduction_velo=0,
 60        flat_earth_transform=args.regional, earth_model_layer_num=24,
 61        slowness_window=None, wavenumber_sampling_rate=12, anti_alias=0.01,
 62        free_surface=0,
 63    )
 64    native_input = Path(library) / "10.00" / "0.00" / "0_0" / "grn.inp"
 65    lines = native_input.read_text(encoding="utf-8-sig").splitlines()
 66    headings = [i for i, line in enumerate(lines) if "WAVENUMBER INTEGRATION PARAMETERS" in line]
 67    if len(headings) != 1:
 68        raise ValueError("Expected one native wavenumber section")
 69    data = [line.split("#", 1)[0].strip() for line in lines[headings[0] + 1:]]
 70    data = [line for line in data if line]
 71    controls = [float(value) for value in data[1].split()]
 72    if controls != [1e-6, source_radius_ratio]:
 73        raise ValueError("Existing native spatial-source settings differ; use a fresh --output-dir")
 74    if args.regional and args.reuse:
 75        source_time_function = validate_qseis_stf(library)
 76    observables = ("disp", "strain", "stress") if args.observables == "all" else ("disp",)
 77    if args.reuse:
 78        # Backend metadata does not store observable flags. Check the requested
 79        # binary outputs before reading a displacement-only library as tensors.
 80        native_dir = Path(library) / "10.00" / "0.00" / "0_0"
 81        required = []
 82        for observable in observables:
 83            psv, sh = get_outfile_name_list(observable)
 84            required.extend(native_dir / ("grn_%s.bin" % name) for name in psv + sh)
 85        if any(not path.is_file() for path in required):
 86            raise ValueError("Requested outputs are absent. Recalculate in a fresh "
 87                             "--output-dir without --reuse using --observables all.")
 88    for observable in observables:
 89        # Regional type-0 kernels are rates; the reader integrates them before cropping.
 90        arrays = [MOMENT_NM * seek_qseis2025(
 91            path_green=library, event_depth_km=10.0, receiver_depth_km=0.0,
 92            az_deg=30.0, dist_km=distance, focal_mechanism=MECHANISM,
 93            srate=1 / dt, output_type=observable, rotate=True,
 94            before_p=None, shift=False, pad_zeros=False,
 95        )[:, :output_samples] for distance in distances]
 96        labels = ["E", "N", "U"] if observable == "disp" else ["EE", "EN", "EU", "NN", "NU", "UU"]
 97        unit = {"disp": "m", "strain": "1", "stress": "Pa"}[observable]
 98        save_waveforms(output, report, observable, arrays, distances, dt, labels, unit,
 99                       expected_samples=output_samples, time_limits=(0.0, output_end))
100    report.update(sampling_interval_s=dt, max_frequency_hz=0.5 / dt, time_window_s=window,
101                  output_time_range_s=[0.0, output_end], distances_km=distances,
102                  native_samples=native_samples, earth_model_numeric_rows=24,
103                  wavelet_type=wavelet_type, wavelet_duration_samples=wavelet_duration,
104                  wavelet_duration_s=wavelet_duration * dt,
105                  flat_earth_transform=args.regional, regional=args.regional, source_radius_ratio=source_radius_ratio, point_source=args.point_source)
106    if args.regional:
107        report["source_time_function"] = source_time_function
108        report["physical_source_time_function"] = dict(REGIONAL_STF)
109        report["spectral_settings"] = qseis_spectral_settings(library)
110    finish(output, report, started)
111
112
113if __name__ == "__main__":
114    main()

Shorter introductory calculation#

Without --regional, the same script runs a lighter calculation at 30, 60 and 90 km with 0.5 s sampling. It exports 0–100 s, giving 201 samples, from a 256-sample native library and completes in seconds:

python examples/qseis2025.py --output-dir examples/output/qseis2025

See the QSEIS2025 tutorial for its figures and the backend comparison for cross-backend results and limitations.

Continue#