Refuse domains whose fault buffer sits inside the SW4 sponge

tests/test_generate_domain.py

@@ -7,10 +7,13 @@
 from hypothesis import strategies as st
 
 from source_modelling import magnitude_scaling, sources
+from workflow import defaults
 from workflow.realisations import (
     Magnitudes,
     Rakes,
+    Refinements,
     SourceConfig,
+    SW4Parameters,
     VelocityModelParameters,
 )
 from workflow.scripts import generate_domain
@@ -67,7 +70,7 @@
         source_config=source_config,  # ty: ignore[invalid-argument-type]
         rrups=rrups,
         nz_outline=nz_outline,
-        fault_buffer=2000.0,
+        fault_buffer=14.0,
     )
 
     assert result_domain.polygon.contains(fault_geom), (
@@ -99,7 +102,7 @@
         topo_type="BULLDOZED",
         ds_multiplier=1.2,
         vs30=500.0,
-        fault_buffer=2000,
+        fault_buffer=14.0,
         s_wave_velocity=3500,
         rrup_interpolants=np.array([[5.0, 8.0], [50.0, 50.0]]),
     )
@@ -110,3 +113,59 @@
         velocity_model_parameters,
     )
     assert shapely.contains(domain_parameters.domain.polygon, source.geometry)
+
+
+# Slow because of openquake import
+@pytest.mark.slow
+@pytest.mark.parametrize(
+    ("fault_buffer", "raises"), [(14.0, False), (13.9, True), (2.0, True)]
+)
+def test_generate_domain_gates_the_fault_buffer(
+    fault_buffer: float, raises: bool
+) -> None:
+    """Gate A: refuse to write a domain whose buffer sits inside the SW4 sponge.
+
+    The check is passed the SW4 parameters and refinements explicitly and both
+    default to None, because the EMOD3D-only defaults versions have neither.
+    """
+    version = defaults.DefaultsVersion.v26_7_1Hz
+    sw4_params = SW4Parameters.read_from_defaults(version)
+    refinements = Refinements.read_from_defaults(version)
+
+    source = sources.Point(
+        np.array([-43.0, -172.0, 10000.0]),
+        length_m=1000,
+        width_m=1000,
+        strike=90.0,
+        dip=45.0,
+        dip_dir=180.0,
+    )
+    arguments = (
+        SourceConfig(dict(source=source)),
+        Magnitudes(dict(source=magnitude_scaling.BoldM(6.0))),
+        Rakes(dict(source=180.0)),
+        VelocityModelParameters(
+            min_vs=500.0,
+            version="2.09",
+            topo_type="BULLDOZED",
+            ds_multiplier=1.2,
+            vs30=500.0,
+            fault_buffer=fault_buffer,
+            s_wave_velocity=3500,
+            rrup_interpolants=np.array([[5.0, 8.0], [50.0, 50.0]]),
+        ),
+    )
+
+    if raises:
+        with pytest.raises(ValueError, match="supergrid absorbing layer"):
+            generate_domain.generate_domain(
+                *arguments, sw4_params=sw4_params, refinements=refinements
+            )
+    else:
+        generate_domain.generate_domain(
+            *arguments, sw4_params=sw4_params, refinements=refinements
+        )
+
+    # Without the SW4 parameters there is nothing to check against, so even a
+    # 2 km buffer is accepted: that is the EMOD3D path, and it is unchanged.
+    generate_domain.generate_domain(*arguments)

workflow/scripts/generate_domain.py

@@ -40,13 +40,16 @@
 from source_modelling import magnitude_scaling, moment, sources
 from velocity_modelling import bounding_box
 from velocity_modelling.bounding_box import BoundingBox
-from workflow import log_utils, realisations, utils
+from workflow import log_utils, realisations, sw4, utils
 from workflow.realisations import (
     DomainParameters,
     Magnitudes,
     Rakes,
     RealisationMetadata,
+    RealisationParseError,
+    Refinements,
     SourceConfig,
+    SW4Parameters,
     VelocityModelParameters,
 )
 
@@ -519,6 +522,8 @@ def generate_domain(
     magnitudes: Magnitudes,
     rakes: Rakes,
     velocity_model_parameters: VelocityModelParameters,
+    sw4_params: SW4Parameters | None = None,
+    refinements: Refinements | None = None,
 ) -> DomainParameters:
     """
     Computes simulation domain spatial extent and temporal duration.
@@ -535,27 +540,54 @@ def generate_domain(
     velocity_model_parameters : VelocityModelParameters
         Parameters defining the velocity model, including Vs30, S-wave
         velocity, and the duration scaling multiplier.
+    sw4_params : SW4Parameters, optional
+        The SW4 parameters, when the realisation is destined for SW4. Given
+        together with `refinements`, the fault buffer is checked against the
+        supergrid sponge before a domain is written. Both default to None,
+        which skips the check, because the EMOD3D-only defaults versions have
+        no `sw4` or `refinements` section at all.
+    refinements : Refinements, optional
+        The theoretical mesh refinements, needed to know the coarsest grid
+        spacing the sponge will be measured on. See `sw4_params`.
 
     Returns
     -------
     DomainParameters
         An object containing the computed model domain (bounding box),
         the maximum simulation depth, and the estimated simulation duration.
+
+    Raises
+    ------
+    ValueError
+        If `sw4_params` and `refinements` are both given and the fault buffer
+        is too small to keep sources out of the SW4 supergrid sponge.
     """
 
     rupture_context = rupture_context_from(magnitudes, rakes, velocity_model_parameters)
 
     rrups = find_r_surfaces(
         source_config, magnitudes, velocity_model_parameters.rrup_interpolants
     )
+
+    # The depth does not depend on the lateral domain, so it can be computed
+    # first and used to resolve the refinements. Gate A: refuse to write a
+    # domain whose fault buffer cannot keep sources out of the SW4 supergrid
+    # sponge, rather than discovering it at the end of a run.
+    depth = domain_max_depth(source_config, magnitudes)
+    if sw4_params is not None and refinements is not None:
+        sw4.check_fault_buffer(
+            velocity_model_parameters.fault_buffer,
+            sw4_params,
+            sw4.coarsest_resolution(refinements, depth),
+        )
+
     nz_outline = utils.get_nz_outline_polygon()
     model_domain = estimate_domain(
         source_config, rrups, nz_outline, velocity_model_parameters.fault_buffer
     )
     sim_duration = estimate_simulation_duration(
         rupture_context, model_domain, source_config.source_geometries.values()
     )
-    depth = domain_max_depth(source_config, magnitudes)
 
     domain_parameters = DomainParameters(
         domain=model_domain,
@@ -602,8 +634,27 @@ def generate_domain_from_realisation(
 
     magnitudes = Magnitudes.read_from_realisation(realisation_ffp)
     rakes = Rakes.read_from_realisation(realisation_ffp)
+
+    # The EMOD3D-only defaults versions have no `sw4` or `refinements` section,
+    # so these are genuinely optional and their absence is not an error.
+    try:
+        sw4_params = SW4Parameters.read_from_realisation_or_defaults(
+            realisation_ffp, metadata.defaults_version
+        )
+        refinements = Refinements.read_from_realisation_or_defaults(
+            realisation_ffp, metadata.defaults_version
+        )
+    except RealisationParseError:
+        sw4_params = None
+        refinements = None
+
     domain_parameters = generate_domain(
-        source_config, magnitudes, rakes, velocity_model_parameters
+        source_config,
+        magnitudes,
+        rakes,
+        velocity_model_parameters,
+        sw4_params=sw4_params,
+        refinements=refinements,
     )
     domain_parameters.write_to_realisation(realisation_ffp)
     realisations.append_log_entry(realisation_ffp)