Quickstart: static Coulomb stress

This calculation uses one 1 km × 1 km source patch with 1 m slip and receivers at 30, 60 and 90 km epicentral distance. The patch is 10 km deep; receivers are 5 km deep. Both mechanisms are strike 30°, dip 45°, rake 90°. Friction is 0.4 and B_pore=0.

The script copies examples/wenchuan/input/model.nd and uses its first 24 numeric rows for the static backend. This illustrative geometry is not a Wenchuan reproduction. The coarse library checks the workflow; it is not a convergence study.

1. Run the calculation

Follow installation, then run from the DynCFS repository root:

python docs/examples/quickstart.py

On Windows with Conda:

conda run -n cfs python docs/examples/quickstart.py

Files are written below docs/_build/quickstart/. The script requires a new or empty output directory. To repeat:

python docs/examples/quickstart.py --output-dir docs/_build/quickstart-repeat

2. Inspect the inputs and library

The script writes source/receiver CSV files and an INI with absolute paths, then calls create_static_lib and compute_static_cfs.

The library has source depths 10 and 11 km, receiver depth 5 km, and distances 1–121 km at 10 km spacing. EDGRN needs at least two source depths. Queries stay inside the library range.

To prepare inputs and run the equivalent CLI steps:

python docs/examples/quickstart.py --prepare-only --output-dir docs/_build/prepared
python -m dyncfs.main --config docs/_build/prepared/quickstart.ini --create-static-lib --compute-static-cfs

The CLI writes numerical results; plotting and the summary are steps in the complete Python script.

3. Read the results

File below the output directory	Meaning
quickstart.ini	Configuration with resolved paths
input/	Source, receivers and copied Earth model
grn_s/	Static library and metadata
results/static/stress_tensor_plane1.npy	Shape (3, 6), Pa, NED order
results/static/cfs_static_plane1.csv	Three rows, one column, Pa
static_cfs.png	Normal, shear and Coulomb stress in kPa
summary.json	Environment, timing, shapes and checks

Local quickstart output. CSV values remain in pascals; this figure displays kilopascals.

The script checks dimensions, finite nonzero values and CFS = shear_stress + 0.4 * normal_stress. See the validation record for what was run.

from pathlib import Path
import numpy as np

result = Path("docs/_build/quickstart/results/static")
stress = np.load(result / "stress_tensor_plane1.npy")
cfs = np.loadtxt(result / "cfs_static_plane1.csv", delimiter=",", ndmin=1)
print(stress.shape)  # (3, 6)
print(cfs / 1e6)     # MPa

4. Continue to dynamic stress

Follow the dynamic workflow and provide a physically appropriate STF. Dynamic synthesis scales the STF to moment; static synthesis uses area and slip. Verify their consistency before comparing them.

The introductory INI includes time fields required by the parser, but its dynamic settings are not a validated dynamic tutorial.

Complete script and configuration

"""Small static DynCFS calculation; run from a source checkout."""
import argparse
import configparser
import json
import os
from pathlib import Path
import platform
import shutil
import sys
import time

ROOT = Path(__file__).resolve().parents[2]
sys.path.insert(0, str(ROOT))
os.environ.setdefault("MPLBACKEND", "Agg")

import matplotlib.pyplot as plt
import numpy as np
import pygrnwang

from dyncfs import __version__
from dyncfs.configuration import CfsConfig
from dyncfs.cfs_static import create_static_lib, compute_static_cfs
from pygrnwang.geo import d2km


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--output-dir", type=Path, default=ROOT / "docs/_build/quickstart")
    parser.add_argument("--prepare-only", action="store_true", help="Write inputs without running solvers")
    args = parser.parse_args()
    output = args.output_dir.resolve()
    if output.exists() and any(output.iterdir()):
        parser.error("Use a new or empty output directory to keep calculations separate.")
    started = time.perf_counter()
    input_dir = output / "input"
    input_dir.mkdir(parents=True, exist_ok=True)
    shutil.copyfile(ROOT / "examples/wenchuan/input/model.nd", input_dir / "model.nd")
    distances = np.array([30.0, 60.0, 90.0])
    # Patch center; static synthesis uses area and slip. STF is included for the file schema.
    source = [[0, 0, 10, 30, 45, 90, 1, 1, 1, 3.112616e16, 0, 0.5, 1, 0.5, 0]]
    receivers = np.array([[0, d / d2km, 5, 30, 45, 90] for d in distances])
    np.savetxt(input_dir / "source_plane1.csv", source, delimiter=",")
    np.savetxt(input_dir / "obs_plane1.csv", receivers, delimiter=",")
    ini = configparser.ConfigParser()
    ini.read(ROOT / "docs/examples/quickstart.ini", encoding="utf-8")
    ini["path"]["path_input"] = input_dir.as_posix()
    ini["path"]["path_output"] = output.as_posix()
    config_path = output / "quickstart.ini"
    with config_path.open("w", encoding="utf-8") as stream:
        ini.write(stream)
    print(f"Configuration: {config_path}")
    if args.prepare_only:
        return
    config = CfsConfig()
    config.read_config(str(config_path))
    create_static_lib(config)
    compute_static_cfs(config)
    result_dir = Path(config.path_output_results_static)
    stress = np.load(result_dir / "stress_tensor_plane1.npy")
    names = ["normal_stress_static", "shear_stress_static", "cfs_static"]
    values = np.column_stack([
        np.loadtxt(result_dir / f"{name}_plane1.csv", delimiter=",", ndmin=1)
        for name in names
    ])
    if stress.shape != (3, 6) or values.shape != (3, 3):
        raise AssertionError(f"Unexpected output shapes: {stress.shape}, {values.shape}")
    if not np.isfinite(stress).all() or not np.isfinite(values).all() or not np.any(values):
        raise AssertionError("Results must be finite and nonzero")
    np.testing.assert_allclose(values[:, 2], values[:, 1] + config.mu_f * values[:, 0],
                               rtol=1e-12, atol=1e-10)
    fig, ax = plt.subplots(figsize=(8, 4.5), constrained_layout=True)
    for index, label in enumerate(["Normal stress", "Shear stress", "Coulomb stress"]):
        ax.plot(distances, values[:, index] / 1000, "o-", label=label)
    ax.axhline(0, color="0.5", lw=0.8)
    ax.set(xlabel="Epicentral distance (km)", ylabel="Stress change (kPa)",
           title="DynCFS · a 1 km² patch with 1 m slip")
    ax.grid(alpha=0.2)
    ax.legend()
    fig.savefig(output / "static_cfs.png", dpi=160)
    plt.close(fig)
    report = {
        "dyncfs": __version__, "pygrnwang_source_version": pygrnwang.__version__,
        "pygrnwang_source": str(Path(pygrnwang.__file__).resolve()),
        "python": platform.python_version(), "platform": platform.platform(),
        "backend": "EDGRN2 + EDCMP2", "elapsed_seconds": time.perf_counter() - started,
        "stress_shape": list(stress.shape), "resolved_stress_shape": list(values.shape),
        "distance_km": distances.tolist(), "stress_unit": "Pa",
        "normal_shear_cfs_pa": values.tolist(),
        "checks": ["finite nonzero outputs", "expected array shapes", "CFS = shear + 0.4 * normal"],
    }
    (output / "summary.json").write_text(json.dumps(report, indent=2) + "\n", encoding="utf-8")
    print(json.dumps(report, indent=2))


if __name__ == "__main__":
    main()

[path]
# The preparation script replaces these with absolute paths.
path_input = docs/_build/quickstart/input
path_output = docs/_build/quickstart

[input_addition]
optimal_type = 0
tectonic_stress_type = 1
tectonic_stress = [-4, 0, 0, -6, 0, -2]
mu_f = 0.4
B_pore = 0
source_inds = [1]
source_shapes = [[1, 1]]
source_ref = [0, 0]
obs_inds = [1]
obs_shapes = [[3, 1]]
obs_ref = [0, 0]
earth_model_layer_num = 24
use_spherical = False
slip_thresh = 0
cut_stf = 0
correct_zero_freq = False

[fixed_obs_depth]
fixed_obs_depth = 0
obs_lat_range = [0, 0]
obs_lon_range = [0.3, 0.8]
obs_delta_lat = 0.1
obs_delta_lon = 0.1
receiver_mechanism = None

[grn_region]
# EDGRN needs at least two source depths.
grn_source_depth_range = [10, 11]
grn_delta_source_depth = 1
grn_obs_depth_range = [5, 5]
grn_delta_obs_depth = 1
grn_dist_unit = km
grn_dist_range = [1, 121]
grn_delta_dist = 10

[time_window]
sampling_interval_stf = 0.5
sampling_interval_cfs = 0.5
sampling_num = 256
max_frequency = 1

[parallel]
processes_num = 1
check_finished = False

[default_config]
default_config = True