From 2d5cf5ffcf36e77c0193a75825febe5b50b6c8ee Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 31 Aug 2026 20:28:54 +1200 Subject: [PATCH] Describe SW4 grids in the realisation, and add the v26.7.1Hz profile Introduces the realisation sections an SW4 run is configured by, the geometry derived from them, and the first shipped set of values. They are one change because the geometry exists only to read the config classes, and the defaults are the shipped instance of both -- the tests that pin the geometry assert against those numbers, not synthetic ones. **Configuration.** `Refinements` describes the stack of vertical grids SW4 solves on, in the abstract -- the layers a domain of any depth would be given. `refinements_for_depth` resolves that against a particular domain: it truncates at the domain floor, extends with `unbounded_refinement_resolution` if the listed layers do not reach that far, and guarantees the last layer is at least two cells deep so a domain ending just past a refinement boundary does not degenerate. `SW4Parameters` holds SW4's input file as a list of `SW4Command`s rather than a fixed set of typed fields. SW4 has a large and growing command vocabulary and we do not want a field per command; anything SW4 accepts can be written as `{"name": ..., "parameters": {...}}` without touching this code. `render()` handles the two places SW4 differs from Python's `str()`: booleans go out as 0/1, and `None` parameters are omitted rather than rendered as "None". **Geometry.** `workflow.sw4` is the single definition of the supergrid sponge, which three later stages need and must not each re-derive. `minimum_fault_buffer_m` is `sponge + 5h`, additive rather than a multiple: the two terms have different origins -- the sponge is where the equation changes, the five grid points are the source stencil plus the dissipation operator applied to its outermost point -- and a multiplicative margin would collapse as the grid refines while the stencil still spans five points. `supergrid_width` follows SW4's own precedence: `width=` beats `gp=`, and with neither, SW4's built-in 30 gridpoint default. `absorbed_period` is advisory and nothing enforces it. Clearing the sponge is not the same as the sponge working: the layer absorbs adiabatically only while `W cos(theta) / lambda >> 0.431`, so a 12 km sponge absorbs below ~8 s at normal incidence and ~4 s at 60 degrees. That is a property of the period band being asked for, not of the buffer. **Values.** The load-bearing one is `fault_buffer: 14.0 km`, against the root default of 2.0 km. Two criteria bear on it and only the first is enforced (by a later commit): HARD FLOOR, 14 km. A source must sit outside the sponge plus the stencil margin: (30 + 5) * 400 m, at the widest sponge any supported domain produces. Below this the run is not a ground motion prediction -- it is the standing-wave failure seen in validation_results_24-08. LONG-PERIOD CRITERION, ~80 km. Sound runs in that campaign sat 83-113 km from the domain edge, broken ones 4-5 km. Deliberately NOT enforced: enforcing it would reject every affordable domain. The 2.0 km root default stays. It is the EMOD3D value, kept for CyberShake reproducibility; v24.2.2.x have no absorbing layer to clear. The sponge is stated as `width: 12000.0` rather than `gp: 30`, because `gp` is measured on the coarsest grid and would vary with domain depth. `SKIP_PAIRS` records the two genuinely version-specific sections: `refinements` describes a grid the EMOD3D versions do not have, and `resolution` is a single uniform spacing an SW4 run does not have. No pipeline stage reads any of this yet. Co-Authored-By: Claude Opus 5 --- tests/test_realisation.py | 145 ++++++++- tests/test_sw4.py | 199 +++++++++++++ .../default_parameters/root/defaults.yaml | 4 + .../default_parameters/v26_7_1Hz/__init__.py | 0 .../v26_7_1Hz/defaults.yaml | 110 +++++++ workflow/defaults.py | 1 + workflow/domain.py | 19 ++ workflow/realisations.py | 164 ++++++++++ workflow/schemas.py | 46 +++ workflow/sw4.py | 279 ++++++++++++++++++ 10 files changed, 966 insertions(+), 1 deletion(-) create mode 100644 tests/test_sw4.py create mode 100644 workflow/default_parameters/v26_7_1Hz/__init__.py create mode 100644 workflow/default_parameters/v26_7_1Hz/defaults.yaml create mode 100644 workflow/domain.py create mode 100644 workflow/sw4.py diff --git a/tests/test_realisation.py b/tests/test_realisation.py index 4a2223d4..18da2867 100644 --- a/tests/test_realisation.py +++ b/tests/test_realisation.py @@ -16,7 +16,7 @@ from source_modelling import magnitude_scaling, rupture_propagation from velocity_modelling import bounding_box from workflow import defaults, realisations, schemas -from workflow.realisations import SourceConfig +from workflow.realisations import SourceConfig, SW4Command def test_bounding_box_example(tmp_path: Path) -> None: @@ -886,6 +886,133 @@ def test_resolution(tmp_path: Path) -> None: assert realisations.Resolution.read_from_realisation(realisation_path) == resolution +def test_refinements(tmp_path: Path) -> None: + refinements = realisations.Refinements( + refinements=[ + realisations.Refinement(resolution=50.0, bottom=2000.0), + realisations.Refinement(resolution=100.0, bottom=5000.0), + realisations.Refinement(resolution=200.0, bottom=25000.0), + ], + unbounded_refinement_resolution=400.0, + ) + + realisation_path = tmp_path / "realisation.json" + refinements.write_to_realisation(realisation_path) + with open(realisation_path, "r") as realisation_handle: + assert json.load(realisation_handle) == { + "refinements": { + "refinements": [ + {"resolution": 50.0, "bottom": 2000.0}, + {"resolution": 100.0, "bottom": 5000.0}, + {"resolution": 200.0, "bottom": 25000.0}, + ], + "unbounded_refinement_resolution": 400.0, + } + } + + assert ( + realisations.Refinements.read_from_realisation(realisation_path) == refinements + ) + + +def test_sw4_command() -> None: + command = SW4Command( + "imagehdf5", + { + "mode": "hmax", + "plane": "z", + "plane_value": 0, + "file": "surf_hmax", + "time": None, + }, + ) + assert ( + command.render() == "imagehdf5 mode=hmax plane=z plane_value=0 file=surf_hmax" + ) + + merged = command.merged(time=10.0, cycle=None) + assert merged.parameters["time"] == 10.0 + assert ( + merged.render() + == "imagehdf5 mode=hmax plane=z plane_value=0 file=surf_hmax time=10.0" + ) + # merged() must not mutate the original command. + assert command.parameters["time"] is None + + +def test_sw4_command_renders_bools_as_ints() -> None: + """SW4 expects 0/1 for boolean flags, not Python's True/False.""" + command = SW4Command("developer", {"reporttiming": True, "failonnan": False}) + assert command.render() == "developer reporttiming=1 failonnan=0" + + +def test_sw4_parameters(tmp_path: Path) -> None: + sw4 = realisations.SW4Parameters( + verbose=2, + printcycle=10, + nz_min=12, + commands=[ + SW4Command( + "grid", + { + "proj": "tmerc", + "ellps": "GRS80", + "lon_p": 173.0, + "lat_p": 0.0, + "scale": 0.9996, + }, + ), + SW4Command("developer", {"cfl": 0.9, "reporttiming": True}), + SW4Command( + "topography", + {"order": 3}, + ), + SW4Command( + "imagehdf5", + { + "mode": "hmax", + "plane": "z", + "plane_value": 0, + "file": "surf_hmax", + "precision": "float", + }, + ), + ], + ) + + realisation_path = tmp_path / "realisation.json" + sw4.write_to_realisation(realisation_path) + with open(realisation_path, "r") as realisation_handle: + written = json.load(realisation_handle) + assert written["sw4"]["verbose"] == 2 + assert len(written["sw4"]["commands"]) == 4 + assert written["sw4"]["commands"][1] == { + "name": "developer", + "parameters": {"cfl": 0.9, "reporttiming": True}, + } + + assert realisations.SW4Parameters.read_from_realisation(realisation_path) == sw4 + + +def test_sw4_parameters_defaults_loadable() -> None: + """SW4Parameters should load from v26_7_1Hz defaults and raise for older versions.""" + sw4 = realisations.SW4Parameters.read_from_defaults( + defaults.DefaultsVersion.v26_7_1Hz + ) + assert sw4.verbose == 2 + developer = realisations.find_command(sw4.commands, "developer") + assert developer is not None + assert developer.parameters["reporttiming"] is True + assert developer.parameters["cfl"] == 0.9 + assert sum(1 for command in sw4.commands if command.name == "imagehdf5") == 10 + + for version in defaults.DefaultsVersion: + if version == defaults.DefaultsVersion.v26_7_1Hz: + continue + with pytest.raises(realisations.RealisationParseError): + realisations.SW4Parameters.read_from_defaults(version) + + def test_sources(tmp_path: Path) -> None: realisation_ffp = tmp_path / "realisation.json" source_json = { @@ -934,6 +1061,17 @@ def test_sources(tmp_path: Path) -> None: assert json.load(f_old) == json.load(f_new) +SKIP_PAIRS = { + # `refinements` and `sw4` describe an SW4 grid; the EMOD3D-only versions + # have neither. `resolution` is the uniform EMOD3D grid spacing, which an + # SW4 run has no single value for. + (defaults.DefaultsVersion.v24_2_2_1, realisations.Refinements), + (defaults.DefaultsVersion.v24_2_2_2, realisations.Refinements), + (defaults.DefaultsVersion.v24_2_2_4, realisations.Refinements), + (defaults.DefaultsVersion.v26_7_1Hz, realisations.Resolution), +} + + @pytest.mark.parametrize( "realisation_config", [ @@ -947,6 +1085,7 @@ def test_sources(tmp_path: Path) -> None: realisations.HFVelocityModel1D, realisations.Resolution, realisations.RuptureVelocity, + realisations.Refinements, ], ) @pytest.mark.parametrize("defaults_version", list(defaults.DefaultsVersion)) @@ -955,4 +1094,8 @@ def test_defaults_are_loadable( realisation_config: realisations.RealisationConfiguration, defaults_version: defaults.DefaultsVersion, ) -> None: + if (defaults_version, realisation_config) in SKIP_PAIRS: + pytest.skip( + f"Configuration {realisation_config} unsupported for defaults {defaults_version}" + ) realisation_config.read_from_defaults(defaults_version) diff --git a/tests/test_sw4.py b/tests/test_sw4.py new file mode 100644 index 00000000..ef92b97e --- /dev/null +++ b/tests/test_sw4.py @@ -0,0 +1,199 @@ +"""Tests for `workflow.sw4`, the shared supergrid geometry. + +These are decision tests. The numbers here are the ones the rest of the SW4 +pipeline is built on, and the point of pinning them is that a silent change to +any of them puts a source back inside the absorbing layer. +""" + +import pytest + +from workflow import defaults, sw4 +from workflow.realisations import ( + Refinements, + SW4Command, + SW4Parameters, + VelocityModelParameters, +) + +DEEPEST_SUPPORTED_DOMAIN_KM = 350.0 +"""The deepest domain any realisation can ask for. + +One static `fault_buffer` default has to clear the widest sponge any run can +produce, and the sponge is widest on the coarsest grid, which is the one a very +deep domain gets. +""" + + +def sw4_parameters(**supergrid_parameters: float) -> SW4Parameters: + """Build minimal SW4 parameters carrying a single `supergrid` command. + + Parameters + ---------- + **supergrid_parameters : float + Parameters for the `supergrid` command. Pass none to omit the command + entirely. + + Returns + ------- + SW4Parameters + The parameters. + """ + commands = [SW4Command("grid", {"proj": "tmerc"})] + if supergrid_parameters: + commands.append(SW4Command("supergrid", dict(supergrid_parameters))) + return SW4Parameters(verbose=2, printcycle=10, nz_min=12, commands=commands) + + +def test_supergrid_width_from_gridpoints() -> None: + """`gp=` is a thickness on the coarsest grid, so it scales with resolution.""" + parameters = sw4_parameters(gp=30) + assert sw4.supergrid_width(parameters, 400.0) == 12000.0 + assert sw4.supergrid_width(parameters, 200.0) == 6000.0 + + +def test_supergrid_width_from_metres() -> None: + """`width=` is already metres and must not scale with resolution.""" + parameters = sw4_parameters(width=12000.0) + assert sw4.supergrid_width(parameters, 400.0) == 12000.0 + assert sw4.supergrid_width(parameters, 200.0) == 12000.0 + + +def test_supergrid_width_prefers_width_over_gridpoints() -> None: + """SW4 rejects both together, but if they appear, `width=` is what it uses.""" + parameters = sw4_parameters(gp=30, width=6000.0) + assert sw4.supergrid_width(parameters, 400.0) == 6000.0 + + +def test_supergrid_width_falls_back_to_the_sw4_default() -> None: + """With no `supergrid` command SW4 still builds a sponge, at its own default.""" + parameters = sw4_parameters() + assert ( + sw4.supergrid_width(parameters, 400.0) + == sw4.SW4_DEFAULT_SUPERGRID_GRIDPOINTS * 400.0 + ) + # An empty `supergrid` command is the same situation. + assert sw4.supergrid_width(sw4_parameters(dc=0.02), 400.0) == 12000.0 + + +def test_minimum_fault_buffer_is_additive() -> None: + """`sponge + 5h`, not `k * sponge`. + + The two terms have different physical origins, so they must not be folded + into a multiplier: a multiplicative margin collapses to nothing as the grid + refines, even though the stencil still spans five points. + """ + parameters = sw4_parameters(gp=30) + assert sw4.minimum_fault_buffer_m(parameters, 400.0) == 14000.0 + assert sw4.minimum_fault_buffer_m(parameters, 200.0) == 7000.0 + + # Stated explicitly, so the additive form cannot be refactored away. + for resolution in (100.0, 200.0, 400.0): + assert ( + sw4.minimum_fault_buffer_m(parameters, resolution) + == (sw4.SW4_DEFAULT_SUPERGRID_GRIDPOINTS + sw4.STENCIL_MARGIN_GRIDPOINTS) + * resolution + ) + + +def test_check_fault_buffer_boundary() -> None: + """14.0 km is exactly enough on a 400 m grid; 13.9 km is not.""" + parameters = sw4_parameters(gp=30) + sw4.check_fault_buffer(14.0, parameters, 400.0) + with pytest.raises(ValueError, match="supergrid absorbing layer"): + sw4.check_fault_buffer(13.9, parameters, 400.0) + + +def test_check_fault_buffer_message_names_the_remedy() -> None: + """The error has to say what to change, not just that something is wrong.""" + with pytest.raises(ValueError) as error: + sw4.check_fault_buffer(2.0, sw4_parameters(gp=30), 400.0) + message = str(error.value) + assert "fault_buffer" in message + assert "14.000 km" in message + + +def test_coarsest_resolution_is_the_bottom_refinement() -> None: + """The sponge is measured on SW4's `mGridSize[0]`, the deepest layer.""" + refinements = Refinements.read_from_defaults(defaults.DefaultsVersion.v26_7_1Hz) + assert sw4.coarsest_resolution(refinements, 3.0) == 100.0 + assert sw4.coarsest_resolution(refinements, 20.0) == 200.0 + assert sw4.coarsest_resolution(refinements, 60.0) == 400.0 + assert sw4.coarsest_resolution(refinements, DEEPEST_SUPPORTED_DOMAIN_KM) == 400.0 + + +def test_default_fault_buffer_is_the_derived_minimum() -> None: + """The YAML holds the number, Python holds the derivation, this pins them. + + A static YAML cannot hold a derived value, so the only thing keeping + `v26_7_1Hz`'s `fault_buffer` honest is this test. It is evaluated at the + deepest supported domain, because one default has to clear the widest + sponge any run can produce. + """ + version = defaults.DefaultsVersion.v26_7_1Hz + sw4_params = SW4Parameters.read_from_defaults(version) + refinements = Refinements.read_from_defaults(version) + velocity_model = VelocityModelParameters.read_from_defaults(version) + + coarsest = sw4.coarsest_resolution(refinements, DEEPEST_SUPPORTED_DOMAIN_KM) + minimum = sw4.minimum_fault_buffer_m(sw4_params, coarsest) + + assert minimum == 14000.0 + assert velocity_model.fault_buffer * 1000.0 == minimum + + # And it must actually pass its own gate at every supported depth. + for depth in (5.0, 25.0, 60.0, DEEPEST_SUPPORTED_DOMAIN_KM): + sw4.check_fault_buffer( + velocity_model.fault_buffer, + sw4_params, + sw4.coarsest_resolution(refinements, depth), + ) + + +def test_root_fault_buffer_is_left_alone() -> None: + """The EMOD3D-only versions keep 2.0 km for CyberShake reproducibility.""" + for version in defaults.DefaultsVersion: + if version == defaults.DefaultsVersion.v26_7_1Hz: + continue + assert VelocityModelParameters.read_from_defaults(version).fault_buffer == 2.0 + + +def test_check_lateral_gridpoints() -> None: + """A grid whose two sponges meet has no interior to simulate in.""" + parameters = sw4_parameters(width=12000.0) + # 100 km domain padded to 124 km: 100 km of interior. + sw4.check_lateral_gridpoints(124000.0, 124000.0, parameters, 400.0) + # Exactly 2 * 5 gridpoints of interior is the boundary case. + sw4.check_lateral_gridpoints(28000.0, 28000.0, parameters, 400.0) + with pytest.raises(ValueError, match="supergrid sponges"): + sw4.check_lateral_gridpoints(27999.0, 28000.0, parameters, 400.0) + with pytest.raises(ValueError, match="supergrid sponges"): + sw4.check_lateral_gridpoints(28000.0, 27999.0, parameters, 400.0) + + +def test_absorbed_period() -> None: + """`T_max = W cos(theta) / (0.431 c)`, pinned against edits to the constant.""" + parameters = sw4_parameters(width=12000.0) + assert sw4.absorbed_period(parameters, 400.0, 3.5) == pytest.approx(7.96, abs=5e-3) + assert sw4.absorbed_period(parameters, 400.0, 3.5, 60.0) == pytest.approx( + 3.98, abs=5e-3 + ) + # Grazing incidence destroys absorption; a nominally adequate sponge is not + # adequate for a wave running along it. + assert sw4.absorbed_period(parameters, 400.0, 3.5, 85.0) == pytest.approx( + 0.69, abs=5e-3 + ) + assert sw4.absorbed_period(parameters, 400.0, 3.5, 90.0) == pytest.approx(0.0) + + +def test_adiabatic_coefficient() -> None: + """`max|Psi0'| / 2pi` with `max|Psi0'| = 2772/1024`.""" + assert sw4.ADIABATIC_COEFFICIENT == pytest.approx(0.4308374, abs=1e-7) + + +def test_stencil_margin_matches_sw4s_own_margin() -> None: + """`src_reach(3) + sgd_reach(2)` at 4th order, the order SW4 defaults to. + + This constant is duplicated in SW4's own source-in-sponge check. If one + moves without the other, one of the two guards becomes wrong. + """ + assert sw4.STENCIL_MARGIN_GRIDPOINTS == 5 diff --git a/workflow/default_parameters/root/defaults.yaml b/workflow/default_parameters/root/defaults.yaml index d34d270b..0bc55587 100644 --- a/workflow/default_parameters/root/defaults.yaml +++ b/workflow/default_parameters/root/defaults.yaml @@ -232,6 +232,10 @@ velocity_model: topo_type: "SQUASHED_TAPERED" vs30: 500.0 s_wave_velocity: 3500.0 + # NOTE: This 2.0 km buffer is the EMOD3D value and is kept for CyberShake + # reproducibility (v24.2.2.x are EMOD3D-only and have no absorbing sponge). + # It is far too small for SW4, which needs the buffer to clear the supergrid + # layer; v26.7.1Hz overrides this. See `workflow.sw4.minimum_fault_buffer_m`. fault_buffer: 2.0 rrup_interpolants: [ diff --git a/workflow/default_parameters/v26_7_1Hz/__init__.py b/workflow/default_parameters/v26_7_1Hz/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/workflow/default_parameters/v26_7_1Hz/defaults.yaml b/workflow/default_parameters/v26_7_1Hz/defaults.yaml new file mode 100644 index 00000000..635c53d6 --- /dev/null +++ b/workflow/default_parameters/v26_7_1Hz/defaults.yaml @@ -0,0 +1,110 @@ +--- +bb: + flo: 1.0 + fmidbot: 0.5 + fmin: 0.2 + site_amp_version: "ba2018" +velocity_model: + # Overrides the root default of 2.0 km, which is an EMOD3D value and has no + # absorbing layer to clear. Two criteria bear on this number, and only the + # first is enforced (`workflow.sw4.check_fault_buffer`): + # + # 1. HARD FLOOR, 14 km. A source must sit outside the SW4 supergrid sponge, + # inside which SW4 solves a damped, coordinate-stretched equation rather + # than the wave equation, plus a discretisation margin for its own stencil: + # `(gp + STENCIL_MARGIN_GRIDPOINTS) * coarsest_resolution` + # = (30 + 5) * 400 m = 14 km at the widest sponge any supported domain + # produces. Below this the run is not a ground motion prediction; it is + # the standing-wave failure seen in validation_results_24-08. + # + # 2. LONG-PERIOD CRITERION, ~80 km. Clearing the sponge is not the same as + # the sponge working. The layer absorbs adiabatically only while + # `W cos(theta) / lambda >> max|Psi0'| / 2pi = 0.431`, so a 12 km sponge + # absorbs only below ~8 s at normal incidence and ~4 s at 60 degrees. + # Empirically, sound runs in that campaign were 83-113 km from the domain + # edge and broken ones 4-5 km. This is NOT enforced: it is a property of + # the period band being asked for, not of the buffer alone, and enforcing + # it would reject every affordable domain. + fault_buffer: 14.0 +refinements: + refinements: + [ + { resolution: 100.0, bottom: 5000.0 }, + { resolution: 200.0, bottom: 25000.0 }, + ] + unbounded_refinement_resolution: 400.0 +sw4: + verbose: 2 + printcycle: 10 + nz_min: 12 + commands: + - name: "grid" + parameters: + { + proj: "tmerc", + ellps: "GRS80", + lon_p: 173.0, + lat_p: 0.0, + scale: 0.9996, + } + # The sponge width is stated in metres rather than as `gp: 30`, which would + # be 30 gridpoints of the coarsest (400 m) grid and so vary with the + # domain depth. `gp` and `width` are mutually exclusive in SW4. + - name: "supergrid" + parameters: { width: 12000.0 } + - name: "attenuation" + parameters: { maxfreq: 10.0, phasefreq: 0.5, nmech: 3 } + - name: "developer" + parameters: { cfl: 0.9, reporttiming: true, failonnan: true } + - name: "prefilter" + parameters: { order: 2, passes: 2, fc1: null, fc2: 1.0, type: "lowpass" } + - name: "topography" + parameters: { order: 3 } + - name: "imagehdf5" + parameters: + { mode: "topo", z: 0.0, file: "topo", cycle: 0, precision: "float" } + - name: "imagehdf5" + parameters: + { mode: "grid", z: 0.0, file: "grid", cycle: 0, precision: "float" } + - name: "imagehdf5" + parameters: + { mode: "p", z: 0.0, file: "surf_vp", cycle: 0, precision: "float" } + - name: "imagehdf5" + parameters: + { mode: "s", z: 0.0, file: "surf_vs", cycle: 0, precision: "float" } + - name: "imagehdf5" + parameters: + { mode: "rho", z: 0.0, file: "surf_rho", cycle: 0, precision: "float" } + - name: "imagehdf5" + parameters: + { + mode: "mag", + z: 0.0, + file: "surf_mag", + timeInterval: 0.5, + precision: "float", + } + - name: "imagehdf5" + parameters: + { + mode: "velmag", + z: 0.0, + file: "surf_velmag", + timeInterval: 0.5, + precision: "float", + } + - name: "imagehdf5" + parameters: + { + mode: "uz", + z: 0.0, + file: "surf_uz", + timeInterval: 0.5, + precision: "float", + } + - name: "imagehdf5" + parameters: + { mode: "hmax", z: 0.0, file: "surf_hmax", precision: "float" } + - name: "imagehdf5" + parameters: + { mode: "vmax", z: 0.0, file: "surf_vmax", precision: "float" } diff --git a/workflow/defaults.py b/workflow/defaults.py index 023c0926..19242e9a 100644 --- a/workflow/defaults.py +++ b/workflow/defaults.py @@ -16,6 +16,7 @@ class DefaultsVersion(StrEnum): v24_2_2_1 = "24.2.2.1" v24_2_2_2 = "24.2.2.2" v24_2_2_4 = "24.2.2.4" + v26_7_1Hz = "26.7.1Hz" # noqa: N815 - mirrors the version string def load_defaults(version: DefaultsVersion) -> dict[str, int | float | str]: diff --git a/workflow/domain.py b/workflow/domain.py new file mode 100644 index 00000000..16a1ea08 --- /dev/null +++ b/workflow/domain.py @@ -0,0 +1,19 @@ +from workflow.realisations import ( + DomainParameters, + Refinements, +) + + +def gridpoints_from_domain( + domain_parameters: DomainParameters, refinements: Refinements +) -> int: + depth = domain_parameters.depth + area = domain_parameters.domain.area * (1000**2) + domain_refinements = refinements.refinements_for_depth(depth) + top = 0.0 + gridpoints = 0 + for refinement in domain_refinements: + volume = (refinement.bottom - top) * area + gridpoints += int(volume // (refinement.resolution) ** 3) + top = refinement.bottom + return gridpoints diff --git a/workflow/realisations.py b/workflow/realisations.py index 364b64f2..ccc1e9cf 100644 --- a/workflow/realisations.py +++ b/workflow/realisations.py @@ -922,6 +922,82 @@ def to_dict(self) -> dict: return param_dict +@dataclasses.dataclass +class Refinement: + """One vertical mesh refinement layer.""" + + resolution: float + """Grid spacing within this layer (metres).""" + bottom: float + """Depth of the bottom of this layer (metres).""" + + +@dataclasses.dataclass +class Refinements(RealisationConfiguration): + """The vertical mesh refinements, from the surface down. + + SW4 solves on a stack of grids that coarsen with depth. This describes + that stack in the abstract -- the layers a domain of *any* depth would + be given -- and `refinements_for_depth` resolves it against a + particular domain. + """ + + _config_key: ClassVar[str] = "refinements" + _schema: ClassVar[Schema] = schemas.REFINEMENTS_SCHEMA + refinements: list[Refinement] + unbounded_refinement_resolution: float + + def __post_init__(self) -> None: + """Coerce refinements read from JSON into `Refinement` instances.""" + if self.refinements and not isinstance(self.refinements[0], Refinement): + self.refinements = [Refinement(**item) for item in self.refinements] # type: ignore + + def refinements_for_depth(self, depth: float) -> list[Refinement]: + """Resolve the refinement stack against a domain of a given depth. + + Layers below `depth` are dropped and the last one is truncated to + it. If the stack does not reach `depth`, a final layer at + `unbounded_refinement_resolution` extends to the bottom. The last + layer is always given at least two cells, so a domain that ends + just past a refinement boundary does not produce a degenerate grid. + + Parameters + ---------- + depth : float + The domain depth, in kilometres. + + Returns + ------- + list of Refinement + The layers covering `depth`, from the surface down. + """ + depth_m = depth * 1000.0 + refinements = [] + for refinement in self.refinements: + refinements.append( + dataclasses.replace(refinement, bottom=min(refinement.bottom, depth_m)) + ) + if refinement.bottom > depth_m: + break + else: + # This block only runs when we finish the loop without breaking, i.e. we + # exhaust the refinement list. + refinements.append( + Refinement( + resolution=self.unbounded_refinement_resolution, bottom=depth_m + ) + ) + + match refinements: + case [*_, previous_layer, last_layer]: + # Ensure a minimum amount in the last layer. + last_layer.bottom = max( + previous_layer.bottom + last_layer.resolution * 2, last_layer.bottom + ) + + return refinements + + @dataclasses.dataclass class VelocityModelParameters(RealisationConfiguration): """Parameters defining the velocity model.""" @@ -1281,6 +1357,94 @@ class BroadbandParameters(RealisationConfiguration): site_amp_version: str +@dataclasses.dataclass +class SW4Command: + """A single SW4 input file command, e.g. `attenuation maxfreq=10 nmech=3`.""" + + name: str + """The SW4 command name (e.g. `attenuation`, `supergrid`, `imagehdf5`).""" + parameters: dict[str, str | int | float | bool | None] = dataclasses.field( + default_factory=dict + ) + """The command's key=value parameters. None values are omitted when rendered.""" + + def render(self) -> str: + """Render this command as a single SW4 input file line. + + Returns + ------- + str + The command name followed by its non-None `key=value` parameters. + Booleans render as `0`/`1`, SW4's expected format, rather than + Python's `False`/`True`. + """ + parts = [self.name] + for key, value in self.parameters.items(): + if value is None: + continue + if isinstance(value, bool): + value = int(value) + parts.append(f"{key}={value}") + return " ".join(parts) + + def merged(self, **overrides: str | int | float | bool | None) -> "SW4Command": + """Return a copy of this command with `overrides` merged into its parameters. + + Parameters + ---------- + **overrides : str | int | float | bool | None + Parameter values to overlay on top of the existing parameters. + + Returns + ------- + SW4Command + A new command with the merged parameters. + """ + return dataclasses.replace(self, parameters={**self.parameters, **overrides}) + + +def find_command(commands: list[SW4Command], name: str) -> SW4Command | None: + """Find the first command with the given name. + + Parameters + ---------- + commands : list[SW4Command] + The commands to search. + name : str + The command name to look for. + + Returns + ------- + SW4Command | None + The first matching command, or None if no command has this name. + """ + return next((command for command in commands if command.name == name), None) + + +@dataclasses.dataclass +class SW4Parameters(RealisationConfiguration): + """Parameters for SW4 simulation.""" + + _config_key: ClassVar[str] = "sw4" + _schema: ClassVar[Schema] = schemas.SW4_PARAMETERS_SCHEMA + + verbose: int + """Fileio verbosity level.""" + printcycle: int + """Output fileio print cycle.""" + nz_min: int + """Minimum vertical cells in each refinement layer.""" + commands: list[SW4Command] + """All other SW4 input file commands (grid projection, attenuation, supergrid, + developer, prefilter, topography, imagehdf5 outputs, and any other non-testing + SW4 command). Add any SW4 command here as `{"name": ..., "parameters": {...}}`.""" + + def __post_init__(self) -> None: + """Coerce commands read from JSON into `SW4Command` instances.""" + if self.commands and not isinstance(self.commands[0], SW4Command): + self.commands = [SW4Command(**item) for item in self.commands] # type: ignore + + @dataclasses.dataclass class IntensityMeasureCalculationParameters(RealisationConfiguration): """Intensity measure calculation parameters.""" diff --git a/workflow/schemas.py b/workflow/schemas.py index 81460cf0..fb98dbbd 100644 --- a/workflow/schemas.py +++ b/workflow/schemas.py @@ -1214,3 +1214,49 @@ def _corners_to_array(corners_spec: list[dict[str, float]]) -> np.ndarray: ) } ) + +REFINEMENT_SCHEMA = Schema( + { + Literal( + "resolution", description="Vertical mesh resolution in this layer." + ): And(NUMBER, _is_positive), + Literal( + "bottom", description="Bottom depth of this refinement layer (m)." + ): And(NUMBER, _is_positive), + } +) + +REFINEMENTS_SCHEMA = Schema( + { + Literal( + "refinements", + description="List of vertical mesh refinements from top to bottom.", + ): [REFINEMENT_SCHEMA], + Literal( + "unbounded_refinement_resolution", + description="Resolution below the last refinement layer.", + ): And(NUMBER, _is_positive), + } +) + +SW4_COMMAND_SCHEMA = Schema( + { + "name": str, + "parameters": {str: Or(str, int, float, bool, None)}, + } +) + +SW4_PARAMETERS_SCHEMA = Schema( + { + Literal("verbose", description="Fileio verbosity level."): int, + Literal("printcycle", description="Output fileio print cycle."): int, + Literal( + "nz_min", + description="Minimum vertical cells in each refinement layer.", + ): int, + Literal( + "commands", + description="List of SW4 input file commands (grid, attenuation, supergrid, developer, prefilter, topography, imagehdf5, or any other non-testing SW4 command).", + ): [SW4_COMMAND_SCHEMA], + } +) diff --git a/workflow/sw4.py b/workflow/sw4.py new file mode 100644 index 00000000..b7ea4d3d --- /dev/null +++ b/workflow/sw4.py @@ -0,0 +1,279 @@ +"""SW4 supergrid geometry. + +SW4 surrounds the simulation domain with a *supergrid* absorbing layer (a +"sponge"), inside which it deliberately solves a damped, coordinate-stretched +equation rather than the wave equation. Anything inside that layer — a source, a +receiver — is not a ground motion prediction. + +This module is the single shared definition of the sponge's geometry for the +workflow. It lives here rather than in `workflow.domain` because +`workflow.domain`'s one function is EMOD3D-specific, and because three separate +consumers need these numbers: + +- `workflow.scripts.sw4_template` pads the SW4 grid laterally and vertically by + the sponge width, so the requested domain becomes the grid's *interior*; +- `workflow.scripts.nzvm_input_template` pads the *velocity model* by at least + as much again, so SW4 never queries outside the sfile; +- `workflow.scripts.generate_domain` checks the fault buffer against the sponge + before a domain is ever written. + +Validation lives here and not in `workflow.schemas`: the module docstring of +`workflow.realisations` reserves the schemas for loose, field-level validation, +and these checks are cross-field and depend on the resolved refinements. +""" + +import math + +from workflow.realisations import Refinements, SW4Parameters, find_command + +SW4_DEFAULT_SUPERGRID_GRIDPOINTS = 30 +"""SW4's own default supergrid thickness, in grid points (`sw4/src/EW.C`). + +Used when the realisation's SW4 commands set neither `gp=` nor `width=` on the +`supergrid` command, because that is what SW4 itself would then use. +""" + +STENCIL_MARGIN_GRIDPOINTS = 5 +"""Grid points of clearance required between a source and the sponge. + +A source must be far enough from the sponge that neither its own stencil nor the +dissipation operator applied to its outermost point reaches into the region +where the stretching function is not the identity. At 4th order that is +`src_reach(3) + sgd_reach(2) = 5` grid points: the moment tensor spans +`ic-2..ic+3`, and `addsgd4` sweeps `+/-2` with `reach=1`. The two terms sum +rather than max, because they are consecutive stencils, not alternatives. + +This number must stay equal to `margin_pts` at 4th order in SW4's own source +check. +""" + +ADIABATIC_COEFFICIENT = (2772.0 / 1024.0) / (2.0 * math.pi) +"""The constant in the supergrid's adiabatic absorption criterion. + +SW4's stretching function has `Psi0'(xi) = 2772 * xi**5 * (1 - xi)**5`, so +`max|Psi0'| = 2772/1024` at `xi = 0.5`. The layer absorbs a wave adiabatically +while `W * cos(theta) / lambda >> max|Psi0'| / (2 pi)`, which gives a longest +absorbable period of `W * cos(theta) / (ADIABATIC_COEFFICIENT * c)`. +""" + + +def supergrid_width(sw4_params: SW4Parameters, coarsest_resolution: float) -> float: + """Compute the supergrid sponge width SW4 will use, in metres. + + SW4's `supergrid` command accepts either a thickness in grid points (`gp=`) + or a width in metres (`width=`), and `CHECK_INPUT` in `parseInputFile.C` + makes them mutually exclusive. When both are somehow present, `width=` wins, + matching SW4's own precedence. When neither is given (or there is no + `supergrid` command at all) SW4 falls back to its own default `gp`. + + The grid-point form is measured on SW4's *coarsest* grid (`mGridSize[0]`), + which is the bottom refinement layer, because a single scalar sponge width + serves every grid and every face. + + Parameters + ---------- + sw4_params : SW4Parameters + The SW4 parameters read from the realisation (or defaults). + coarsest_resolution : float + The coarsest grid spacing in the run, in metres. See + `coarsest_resolution`. + + Returns + ------- + float + The sponge width, in metres. + """ + command = find_command(sw4_params.commands, "supergrid") + parameters = command.parameters if command is not None else {} + + width = parameters.get("width") + if width is not None: + return float(width) + + gridpoints = parameters.get("gp") + if gridpoints is None: + gridpoints = SW4_DEFAULT_SUPERGRID_GRIDPOINTS + + return float(gridpoints) * coarsest_resolution + + +def coarsest_resolution(refinements: Refinements, depth_km: float) -> float: + """Find the coarsest grid spacing SW4 will use for a domain, in metres. + + This is SW4's `mGridSize[0]`: the grid spacing of the deepest (and so + coarsest) mesh refinement once the theoretical refinements have been + resolved against the domain depth. + + Parameters + ---------- + refinements : Refinements + The theoretical mesh refinements from the realisation (or defaults). + depth_km : float + The domain depth, in kilometres. + + Returns + ------- + float + The coarsest grid spacing, in metres. + """ + return max( + refinement.resolution + for refinement in refinements.refinements_for_depth(depth_km) + ) + + +def minimum_fault_buffer_m( + sw4_params: SW4Parameters, coarsest_resolution: float +) -> float: + """Compute the smallest fault buffer that clears the supergrid sponge. + + The buffer is **additive**, `sponge + STENCIL_MARGIN_GRIDPOINTS * h`, not a + multiple of the sponge width. The two terms have different physical origins: + the sponge is a geometric exclusion (SW4 solves a different PDE inside it), + while the stencil margin is a discretisation margin that scales with the + grid spacing. A multiplicative `1.2 * sponge` would collapse to a 1.2 km + margin on a 200 m grid even though the stencil still spans five points. + + At SW4's default `gp=30` this is `(30 + 5) * h`, i.e. exactly 14 km on a + 400 m grid and 7 km on a 200 m grid. + + Parameters + ---------- + sw4_params : SW4Parameters + The SW4 parameters read from the realisation (or defaults). + coarsest_resolution : float + The coarsest grid spacing in the run, in metres. + + Returns + ------- + float + The minimum fault buffer, in metres. + """ + return ( + supergrid_width(sw4_params, coarsest_resolution) + + STENCIL_MARGIN_GRIDPOINTS * coarsest_resolution + ) + + +def check_fault_buffer( + fault_buffer_km: float, sw4_params: SW4Parameters, coarsest_resolution: float +) -> None: + """Check that a fault buffer keeps every source clear of the supergrid sponge. + + Parameters + ---------- + fault_buffer_km : float + The `velocity_model.fault_buffer` value, in kilometres. The domain edge + is guaranteed to be at least this far from every source. + sw4_params : SW4Parameters + The SW4 parameters read from the realisation (or defaults). + coarsest_resolution : float + The coarsest grid spacing in the run, in metres. + + Raises + ------ + ValueError + If the buffer is narrower than the sponge plus the stencil margin, so a + source could sit inside the absorbing layer. + """ + minimum = minimum_fault_buffer_m(sw4_params, coarsest_resolution) + if fault_buffer_km * 1000.0 < minimum: + sponge = supergrid_width(sw4_params, coarsest_resolution) + raise ValueError( + f"The fault buffer of {fault_buffer_km:.3f} km is smaller than the " + f"{minimum / 1000.0:.3f} km needed to keep sources out of the SW4 " + f"supergrid absorbing layer. On a {coarsest_resolution:.0f} m " + f"coarsest grid the sponge is {sponge / 1000.0:.3f} km wide, and a " + f"source needs a further {STENCIL_MARGIN_GRIDPOINTS} grid points " + f"({STENCIL_MARGIN_GRIDPOINTS * coarsest_resolution / 1000.0:.3f} km) " + "of clearance for its own stencil and the dissipation operator. " + "Inside the layer SW4 solves a damped, coordinate-stretched " + "equation, so the result is not a ground motion. Raise " + f"velocity_model.fault_buffer to at least {minimum / 1000.0:.3f} km, " + "or narrow the supergrid." + ) + + +def check_lateral_gridpoints( + x_m: float, y_m: float, sw4_params: SW4Parameters, coarsest_resolution: float +) -> None: + """Check that a SW4 grid has a usable interior between its lateral sponges. + + Both lateral axes carry a sponge on each side, so the interior of an axis is + the extent less twice the sponge width. That interior has to be wide enough + for the stencil margin on each side, otherwise the two sponges effectively + meet and there is nowhere in the grid where SW4 solves the wave equation + undisturbed. + + Parameters + ---------- + x_m, y_m : float + The lateral extents of the SW4 grid, in metres, in SW4's own axis + convention (`x` is north). + sw4_params : SW4Parameters + The SW4 parameters read from the realisation (or defaults). + coarsest_resolution : float + The coarsest grid spacing in the run, in metres. + + Raises + ------ + ValueError + If either axis has too little interior left between its two sponges. + """ + sponge = supergrid_width(sw4_params, coarsest_resolution) + minimum_interior = 2 * STENCIL_MARGIN_GRIDPOINTS * coarsest_resolution + + for axis, extent in (("x", x_m), ("y", y_m)): + interior = extent - 2 * sponge + if interior < minimum_interior: + raise ValueError( + f"The SW4 grid's {axis} extent of {extent / 1000.0:.3f} km leaves " + f"only {interior / 1000.0:.3f} km between its two " + f"{sponge / 1000.0:.3f} km supergrid sponges, but at least " + f"{minimum_interior / 1000.0:.3f} km " + f"({2 * STENCIL_MARGIN_GRIDPOINTS} grid points on a " + f"{coarsest_resolution:.0f} m grid) is needed for a usable " + "interior. Widen the domain or narrow the supergrid." + ) + + +def absorbed_period( + sw4_params: SW4Parameters, + coarsest_resolution: float, + vs_km_s: float, + incidence_degrees: float = 0.0, +) -> float: + """Compute the longest period the supergrid sponge can absorb, in seconds. + + The supergrid absorbs adiabatically while the layer is long compared to a + wavelength measured along the layer normal, i.e. while + `W cos(theta) / lambda >> ADIABATIC_COEFFICIENT`. Setting that ratio to one + gives the longest period the layer still absorbs: + `W cos(theta) / (ADIABATIC_COEFFICIENT * c)`. Longer periods are reflected + rather than absorbed, and with a source in or near the layer they ring. + + Parameters + ---------- + sw4_params : SW4Parameters + The SW4 parameters read from the realisation (or defaults). + coarsest_resolution : float + The coarsest grid spacing in the run, in metres. + vs_km_s : float + The shear wave speed at the layer, in km/s. The slowest material + against the layer is the conservative choice. + incidence_degrees : float, default 0.0 + The angle between the ray and the layer normal, in degrees. Grazing + incidence degrades absorption by `cos(theta)`. + + Returns + ------- + float + The longest absorbable period, in seconds. + """ + width = supergrid_width(sw4_params, coarsest_resolution) + speed = vs_km_s * 1000.0 + return ( + width + * math.cos(math.radians(incidence_degrees)) + / (ADIABATIC_COEFFICIENT * speed) + )