Prepare an Earth model#
Use the same physical model for Green’s functions, material normalization and travel times. A model filename in one step does not override all other steps.
Propagation input#
Preprocessors accept path_nd containing six numeric columns:
depth_km vp_km_s vs_km_s density_g_cm3 Qp Qs
Repeated depths represent discontinuities: rows above/below give the material
on each side. Preserve depth order and named boundaries such as mantle,
outer-core and inner-core where the TauP format uses them. Do not add
arbitrary prose, extra columns or headings: the conversion helpers distinguish
numeric rows from one-word boundary labels.
The examples use AK135 elastic velocities/density included in the package and
append illustrative constant Qp=600 and Qs=300. This is not the
original AK135-F attenuation model. The QSEIS and EDGRN examples use the
first 24 numeric rows (ending at 809.5 km); spherical examples retain the
full Earth. No model download or developer-specific absolute path is needed.
Their shared helper is:
"""Shared paths, model preparation, and output checks for executable tutorials."""
import argparse
import importlib.metadata
import json
from pathlib import Path
import platform
import sys
import time
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from pygrnwang.ak135fc import s as AK135_ELASTIC_MODEL
MOMENT_NM = 1e15
MECHANISM = [30.0, 45.0, 90.0]
REGIONAL_SAMPLING_INTERVAL_S = 4.0
REGIONAL_MAX_FREQUENCY_HZ = 0.5 / REGIONAL_SAMPLING_INTERVAL_S
REGIONAL_SOURCE_DURATION_S = 64.0
REGIONAL_STF = {"shape": "normalized_sin_squared_moment_rate",
"duration_s": REGIONAL_SOURCE_DURATION_S,
"physical_integral": 1.0, "centroid_s": 32.0}
def parser_for(name):
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--output-dir", type=Path,
default=Path(__file__).resolve().parent / "output" / name,
help="Directory for the model, library, figures and summary.")
parser.add_argument("--reuse", action="store_true",
help="Read the existing library without rerunning a solver.")
return parser
def prepare(args, backend, extra=None):
"""Resolve paths before solvers change cwd, and write a six-column ND model."""
output = args.output_dir.expanduser().resolve()
output.mkdir(parents=True, exist_ok=True)
if not args.reuse and (output / "library" / "green_lib_info.json").exists():
raise FileExistsError("This tutorial directory already contains a library. "
"Use --reuse for matching settings, or a fresh --output-dir.")
model = output / "ak135_tutorial.nd"
# AK135 elastic velocities/density are bundled; constant Q is an explicit
# tutorial choice, not the frequency-dependent AK135-F attenuation model.
model_text = "\n".join(line + " 600.0 300.0" if len(line.split()) > 1 else line
for line in AK135_ELASTIC_MODEL.splitlines()) + "\n"
if args.reuse:
if not model.is_file() or model.read_text(encoding="utf-8") != model_text:
raise ValueError("--reuse requires the model from a completed tutorial run")
else:
model.write_text(model_text, encoding="utf-8")
library = output / "library"
library.mkdir(exist_ok=True)
report = {"backend": backend, "started_utc": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()),
"platform": platform.platform(), "python": sys.version.split()[0],
"dependencies": {name: importlib.metadata.version(name)
for name in ("numpy", "scipy", "pandas", "matplotlib", "obspy")},
"source_depth_km": 10.0, "receiver_depth_km": 0.0,
"moment_nm": MOMENT_NM, "strike_dip_rake_deg": MECHANISM,
"azimuth_deg": 30.0, "reuse": args.reuse,
"model": "Bundled AK135 elastic properties; illustrative constant Qp=600, Qs=300",
"outputs": {}}
report.update(extra or {})
return output, str(library), str(model), report, time.perf_counter()
def require_library_settings(library, **expected):
"""Reject an older library when a tutorial's numerical settings change."""
info = json.loads((Path(library) / "green_lib_info.json").read_text(encoding="utf-8"))
mismatches = [name for name, value in expected.items() if info.get(name) != value]
if mismatches:
raise ValueError("Existing library settings differ: %s. Recalculate in a fresh "
"--output-dir without --reuse." % ", ".join(mismatches))
return info
def save_waveforms(output, report, name, arrays, distances, dt, labels,
unit, time_label="Time since origin (s)", start_times=None,
expected_samples=256, time_limits=None):
"""Check shape/finiteness/nonzero output and save physical-unit arrays/plots."""
values = np.asarray(arrays)
if values.ndim != 3 or values.shape[:2] != (len(distances), len(labels)):
raise AssertionError("Unexpected waveform shape: %s" % (values.shape,))
if values.shape[2] != expected_samples or not np.isfinite(values).all() or not np.all(np.any(values != 0, axis=(1, 2))):
raise AssertionError("Each distance must have %d finite samples and a nonzero waveform"
% expected_samples)
starts = np.zeros(len(distances)) if start_times is None else np.asarray(start_times)
times = starts[:, None] + np.arange(values.shape[2])[None, :] * dt
np.savez_compressed(output / (name + ".npz"), values=values, time_s=times,
distance_km=distances, components=labels, unit=unit)
fig, axes = plt.subplots(len(labels), 1, figsize=(8, max(4, len(labels) * 1.45)),
sharex=True, constrained_layout=True, squeeze=False)
for component, label in enumerate(labels):
axis = axes[component, 0]
for index, distance in enumerate(distances):
axis.plot(times[index], values[index, component], lw=1.1,
label="%g km" % distance)
axis.set_ylabel("%s (%s)" % (label, unit))
axis.ticklabel_format(axis="y", style="sci", scilimits=(-2, 2))
axis.grid(alpha=0.2)
axes[0, 0].legend(ncol=len(distances), fontsize=8)
axes[0, 0].set_title("%s: %s, M0 = 10^15 N m"
% (report.get("display_name", report["backend"]), name))
axes[-1, 0].set_xlabel(time_label)
if time_limits is not None:
axes[-1, 0].set_xlim(*time_limits)
fig.savefig(output / (name + ".png"), dpi=140)
plt.close(fig)
report["outputs"][name] = {"shape": list(values.shape), "components": labels,
"unit": unit, "finite": True,
"peak_absolute": float(np.max(np.abs(values)))}
def finish(output, report, started):
report["elapsed_seconds"] = round(time.perf_counter() - started, 3)
files = [p for p in output.rglob("*") if p.is_file() and not p.name.startswith("summary")]
report["file_count"] = len(files)
report["output_bytes"] = sum(p.stat().st_size for p in files)
summary_name = "summary-reuse.json" if report["reuse"] else "summary.json"
(output / summary_name).write_text(json.dumps(report, indent=2) + "\n", encoding="utf-8")
print(json.dumps(report, indent=2))
print("Results: %s" % output)
Resolve path_nd and path_green to absolute paths before launching a
backend. Some low-level calls change the current working directory, so
relative paths can fail after the first job starts.
Extent and row selection#
QSEIS uses a layered half-space with optional flat-Earth transformation. Spherical SPGRN/QSSP require an appropriate radially layered Earth model. A shallow crustal table is not automatically a full spherical model. Keep the complete model in the spherical tutorials.
earth_model_layer_num controls the number of numeric model rows, not the
number of distinct geological layers. Leave it None for the supplied
model. If changed, inspect the generated input and confirm its declared
row count and final model boundary.
Zero shear velocity represents fluid. EDCMP moment normalization divides by shear modulus and therefore requires a suitable solid source layer. Strain-to-stress conversion also needs the intended receiver-layer material.
Material values#
read_material_nd(model_name, depth) returns
[depth, vp, vs, density] at or just below the requested depth. It recognizes
ak135fc or a four-column model filename, a different model-name
convention from TauP. Pass the generated noQ.nd, not the six-column
propagation input: read_nd defaults to four-column reshaping and cannot
infer the intended columns.
For solid material in these units,
give Pa. Use source-layer material for EDCMP normalization and receiver-layer material for local strain-to-stress conversion. At an interface, choose the physically relevant side explicitly.
Travel-time input#
Preprocessing writes noQ.nd by removing attenuation columns. Java reads
this model directly; ObsPy builds an adjacent .npz, requiring a writable
model directory. taup_create_npz_file returns a path usable by the active
backend: .nd for Java, .npz for ObsPy.
Prepare a custom model before starting workers. Changing the propagation model requires rebuilding travel-time tables as well; reuse flags do not compare model contents. See TauP and resuming.
Before scaling up#
Run the original backend tutorial and retain its summary.
Substitute your model in a new output directory and inspect
grn.inporspec.inp.Calculate a few P/S times at the intended depths and distances; check finite arrivals for the phases the reader will use.
Build one source/receiver pair and check units, shape and a known arrival or limiting case before expanding the grid.