Rewrite hf-sim on the hf-simulation package

pyproject.toml

@@ -16,6 +16,7 @@ dependencies = [
   "nshmdb>=2026.09.1",
   "oq_wrapper>=2026.05.2",
   "qcore-utils>=2025.12.2",
+  "hf-simulation @ git+https://github.com/ucgmsim/high-frequency",
   "site-calculation>=2026.7.1",
   "source_modelling>=2026.08.1",
   # Data Formats

tests/test_hf.py

@@ -1,60 +1,37 @@
-from pathlib import Path
 from types import SimpleNamespace
 
 import numpy as np
-import pytest
+from hf_simulation import PathDurationModel, Ray
 from hypothesis import given
 from hypothesis import strategies as st
 
 from workflow.realisations import (
     HFConfig,
-    Resolution,
     RuptureVelocity,
 )
 from workflow.scripts import hf_sim
 
 
-def test_build_hf_input_serialisation() -> None:
-    stoch_ffp = Path("/path/to/stoch")
-    velocity_model = Path("/path/to/vmodel")
+def test_build_config_mirrors_the_realisation() -> None:
+    """The realisation's `hf` section reaches `hf_simulation.HfConfig` unchanged.
 
+    The two structures mirror each other group for group, so `build_config` is a splat plus
+    the two values the realisation deliberately does not carry: the record duration, which
+    the domain computes, and the rupture-velocity multipliers, which live in their own
+    section because SRF generation reads them too. This pins both halves of that.
+    """
     hf_config = HFConfig(
-        sdrop=50.0,
-        rayset=[1, 2],
-        no_siteamp=False,
-        nbu=1,
-        ift=0,
-        flo=0.1,
-        fhi=10.0,
-        fmax=20.0,
-        kappa=0.045,
-        qfexp=0.6,
-        czero=2.5,
-        calpha=0.0,
-        mom=None,
-        rupv=1.2,
-        vs_moho=3.5,
-        nl_skip=0,
-        vp_sig=0.1,
-        vsh_sig=0.1,
-        rho_sig=0.1,
-        qs_sig=0.1,
-        ic_flag=True,
-        velocity_name="test_vel",
-        fa_sig1=0.2,
-        fa_sig2=0.2,
-        rv_sig1=0.1,
-        path_dur=11,
-        t_sec=0.0,
-        site_specific=False,
-        dpath_pert=0,
-        stress_parameter_adjustment_fault_area=None,
-        stress_parameter_adjustment_target_magnitude=None,
-        stress_parameter_adjustment_tect_type=0,
+        source={
+            "stress_drop_bars": 50.0,
+            "corner_frequency_constant": 2.5,
+            "corner_frequency_alpha": 0.1,
+            "rupture_velocity": {"sigma": 0.1},
+        },
+        path={"rayset": [1, 2], "q_frequency_exponent": 0.6, "path_duration_model": 11},
+        site={"kappa_s": 0.045, "fmax_hz": 20.0},
+        record={"dt": 0.005},
     )
-
-    res = Resolution(resolution=0.1)
-    rv = RuptureVelocity(
+    rupture_velocity = RuptureVelocity(
         rvfrac=0.8,
         rvfrac_shal=0.7,
         rvfrac_deep=0.9,
@@ -64,32 +41,26 @@ def test_build_hf_input_serialisation() -> None:
         deep_transition_range=1,
         rvfrac_slip_sig=None,
     )
-    # Rather than create DomainParameters with a bounding box, we simplify with a mock object
+    # A bounding box is not needed to read one field off it.
     domain = SimpleNamespace(duration=100.0)
 
-    result = hf_sim.build_hf_input(
-        stoch_ffp,
-        velocity_model,
-        res,
-        hf_config,
-        rv,
-        domain,  # ty: ignore[invalid-argument-type]
-    )
-
-    lines = result.split("\n")
-
-    assert lines[1] == "50.0"  # sdrop
-    assert lines[2] == "{station_input_file}"  # placeholder for station file
-    assert lines[3] == "{output_file}"  # placeholder for output file
-    assert lines[4] == "2 1 2"  # rayset count + rays
-    assert lines[5] == "1"  # int(not no_siteamp) -> int(not False) -> 1
-    assert lines[7] == "{seed}"  # seed placeholder
-    assert lines[9] == "100.0 0.005 20.0 0.045 0.6"  # Domain and resolution parameters
-    assert lines[10] == "0.8 0.7 0.9 2.5 0.0"  # rupture velocity + czero,alpha
-    assert lines[11] == "-1 1.2"  # mom (None -> -1) and rupv
-    assert lines[12] == str(stoch_ffp)  # Stoch file path
-    assert lines[15] == "0 0.1 0.1 0.1 0.1 1"  # Sigs and ic_flag (True -> 1)
-    assert lines[20] == "-1 -1 -1"  # Optional stress parameters
+    config = hf_sim.build_config(hf_config, rupture_velocity, domain)  # ty: ignore[invalid-argument-type]
+
+    # Splatted through unchanged.
+    assert config.source.stress_drop_bars == 50.0
+    assert config.source.corner_frequency_constant == 2.5
+    assert config.site.fmax_hz == 20.0
+    assert config.record.dt == 0.005
+    # Ints become the enums the simulation takes.
+    assert config.path.rayset == (Ray.DIRECT, Ray.MOHO_REFLECTION)
+    assert config.path.path_duration_model is PathDurationModel.BOORE_THOMPSON_2014
+    # Injected, because the `hf` section does not carry them.
+    assert config.record.duration_s == 100.0
+    assert config.source.rupture_velocity.fraction == 0.8
+    assert config.source.rupture_velocity.shallow == 0.7
+    assert config.source.rupture_velocity.deep == 0.9
+    # ... but the sigma does come from the `hf` section.
+    assert config.source.rupture_velocity.sigma == 0.1
 
 
 STATION_STRATEGY = st.text(
@@ -99,15 +70,16 @@ def test_build_hf_input_serialisation() -> None:
 
 def test_station_seeds() -> None:
     seed = hf_sim.station_seeds(0, ["station"])
-    assert seed.dtype == np.int32
+    assert seed.dtype == np.uint64
     assert seed.shape == (1,)
     # Seeds should be referentially transparent: i.e. depend only on the seed and station name
     seed_1 = hf_sim.station_seeds(0, ["station"])
     assert seed.item() == seed_1.item()
 
 
 @given(
-    seed=st.integers(min_value=-(1 << 31), max_value=(1 << 31) - 1),
+    # Non-negative: SeedSequence rejects negative entropy, which `station_seeds` says.
+    seed=st.integers(min_value=0, max_value=(1 << 31) - 1),
     stations=st.lists(STATION_STRATEGY, min_size=1, unique=True),
 )
 def test_station_seeds_on_name_only(seed: int, stations: list[str]) -> None:
@@ -123,48 +95,3 @@ def test_station_seeds_on_name_only(seed: int, stations: list[str]) -> None:
     # one.
     for station, expected_seed in zip(stations, station_seeds):
         assert hf_sim.station_seeds(seed, [station]).item() == expected_seed
-
-
-def test_create_hf_dataset_structure() -> None:
-    # 1. Setup Mock Data
-    n_stations = 2
-    n_components = 3  # Fixed by function logic
-    n_time = 100
-
-    names = ["station_a", "station_b"]
-    waveform = np.random.rand(n_components, n_stations, n_time).astype(np.float32)
-    lat = np.array([-43.5, -43.6])
-    lon = np.array([172.6, 172.7])
-    dist = np.array([10.5, 20.1])
-    seeds = np.array([123, 456])
-    vrefs = np.array([300.0, 350.0])
-    dt = 0.02
-    start_sec = 0.0
-
-    ds = hf_sim.create_hf_dataset(
-        waveform=waveform,
-        latitude=lat,
-        longitude=lon,
-        names=names,
-        epicentre_distance=dist,
-        seed=seeds,
-        vref=vrefs,
-        dt=dt,
-        start_sec=start_sec,
-    )
-
-    assert ds.sizes == {"component": 3, "station": 2, "time": 100}
-
-    np.testing.assert_array_equal(ds.station.values, names)
-    np.testing.assert_array_equal(ds.component.values, ["x", "y", "z"])
-    assert ds.time.values[1] == pytest.approx(0.02)
-    assert ds.lat.dims == ("station",)
-    assert ds.lon.dims == ("station",)
-
-    assert "waveform" in ds.data_vars
-    assert ds.waveform.dims == ("component", "station", "time")
-    np.testing.assert_allclose(ds.waveform.values, waveform)
-
-    assert ds.attrs["dt"] == dt
-    assert ds.attrs["nt"] == n_time
-    assert ds.attrs["units"] == "cm/s^2"

workflow/default_parameters/root/defaults.yaml

@@ -60,39 +60,31 @@ emod3d:
   yseis: 0
   zseis: 0
   pertbfile: none.pertb
+# Mirrors `hf_simulation.HfConfig` group for group and field for field, so this block can
+# be deserialised straight into it rather than translated. Two of that class's inputs are
+# absent on purpose: `record.duration_s` is computed from the domain, and the three
+# `source.rupture_velocity` multipliers live in `rupture_velocity:` below because SRF
+# generation reads the same physical values. hf-sim injects both.
 hf:
-  nbu: 4
-  ift: 0
-  flo: 0.02
-  fhi: 19.9
-  nl_skip: -99
-  vp_sig: 0.0
-  vsh_sig: 0.0
-  rho_sig: 0.0
-  qs_sig: 0.0
-  ic_flag: true
-  velocity_name: "-1"
-  t_sec: 0.0
-  sdrop: 50.0
-  rayset: [1]
-  no_siteamp: false
-  fmax: 10.0
-  kappa: 0.045
-  qfexp: 0.6
-  czero: 2.1
-  calpha: -99.0
-  mom: null
-  rupv: null
-  site_specific: false
-  vs_moho: 999.9
-  fa_sig1: 0.0
-  fa_sig2: 0.0
-  rv_sig1: 0.1
-  path_dur: 11
-  dpath_pert: 0.0
-  stress_parameter_adjustment_tect_type: 0
-  stress_parameter_adjustment_target_magnitude: null
-  stress_parameter_adjustment_fault_area: null
+  source:
+    stress_drop_bars: 50.0
+    corner_frequency_constant: 2.1
+    # 0.1, not -99.0. The Fortran read any value below -1.0 as "use the built-in default"
+    # (`if(Calpha.lt.-1.0) Calpha = Calpha_default`), so -99.0 here MEANT 0.1. That decode
+    # lived in the deck reader, which is gone -- hf-simulation takes the value literally
+    # now, and -99.0 makes alpha_T negative for anything but a vertical strike-slip fault.
+    corner_frequency_alpha: 0.1
+    rupture_velocity:
+      sigma: 0.1
+  path:
+    rayset: [1]
+    q_frequency_exponent: 0.6
+    path_duration_model: 11
+  site:
+    kappa_s: 0.045
+    fmax_hz: 10.0
+  record:
+    dt: 0.005
 stoch:
   stoch_dx: 2.0
   stoch_dy: 2.0
@@ -1530,6 +1522,8 @@ hf_velocity_model_1d:
       rho: 3.33
       Qp: 394.80
       Qs: 197.40
+  # Truncate the model where Vs reaches this. 999.9 means no layer does.
+  vs_moho: 999.9
 bb:
   fmin: 0.2
   fmidbot: 0.5

workflow/default_parameters/v24_2_2_1/defaults.yaml

@@ -3,3 +3,6 @@ resolution:
   resolution: 0.1
 bb:
   flo: 1.0
+hf:
+  record:
+    dt: 0.005

workflow/default_parameters/v24_2_2_2/defaults.yaml

@@ -3,3 +3,6 @@ resolution:
   resolution: 0.2
 bb:
   flo: 0.5
+hf:
+  record:
+    dt: 0.005

workflow/default_parameters/v24_2_2_4/defaults.yaml

@@ -3,3 +3,6 @@ resolution:
   resolution: 0.4
 bb:
   flo: 0.25
+hf:
+  record:
+    dt: 0.005

workflow/default_parameters/v26_7_1Hz/defaults.yaml

@@ -105,3 +105,6 @@ sw4:
     - name: "imagehdf5"
       parameters:
         { mode: "vmax", z: 0.0, file: "surf_vmax", precision: "float" }
+hf:
+  record:
+    dt: 0.005