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) + )