diff --git a/tests/test_realisation.py b/tests/test_realisation.py index b7d198e..5f45398 100644 --- a/tests/test_realisation.py +++ b/tests/test_realisation.py @@ -114,7 +114,6 @@ def test_srf_config_example(tmp_path: Path) -> None: duration=60.0, ) srf_config = realisations.SRFConfig( - genslip_version="5.4.2", resolution=0.1, point_source_params=schemas.PointSourceParams( stype=schemas.Stype.cos, @@ -123,6 +122,19 @@ def test_srf_config_example(tmp_path: Path) -> None: risetimedep=0.0, inittime=0.0, ), + side_taper=0.02, + bot_taper=0.02, + top_taper=0.0, + alpha_rough=0.0, + gwid=[], + rvfac_seg=[], + seg_delay=False, + slip_sigma=1.0, + risetime_coef=1.6, + ymag_exp=None, + xmag_exp=None, + kx_corner=None, + ky_corner=None, ) realisation_ffp = tmp_path / "realisation.json" @@ -141,7 +153,6 @@ def test_srf_config_example(tmp_path: Path) -> None: "duration": 60.0, }, "srf": { - "genslip_version": "5.4.2", "resolution": 0.1, "point_source_params": { "stype": "cos", @@ -150,6 +161,19 @@ def test_srf_config_example(tmp_path: Path) -> None: "risetimedep": 0.0, "inittime": 0.0, }, + "side_taper": 0.02, + "bot_taper": 0.02, + "top_taper": 0.0, + "alpha_rough": 0.0, + "slip_sigma": 1.0, + "risetime_coef": 1.6, + "gwid": [], + "rvfac_seg": [], + "seg_delay": False, + "ymag_exp": None, + "xmag_exp": None, + "kx_corner": None, + "ky_corner": None, }, } diff --git a/workflow/default_parameters/develop/defaults.yaml b/workflow/default_parameters/develop/defaults.yaml index 1113b6d..7bb2a73 100644 --- a/workflow/default_parameters/develop/defaults.yaml +++ b/workflow/default_parameters/develop/defaults.yaml @@ -1,3 +1,4 @@ +--- resolution: resolution: 0.4 emod3d: @@ -100,8 +101,20 @@ hf: stoch_dx: 2.0 stoch_dy: 2.0 srf: - genslip_version: 5.4.2 + side_taper: 0.02 + bot_taper: 0.02 + top_taper: 0.0 + alpha_rough: 0.0 + gwid: [] + rvfac_seg: [] + seg_delay: false + ymag_exp: null + xmag_exp: null + kx_corner: null + ky_corner: null resolution: 0.1 + slip_sigma: 0.75 + risetime_coef: 1.6 # The following parameters are only used for the point source approximation. point_source_params: # Valid options for the slip time function (stype) are: diff --git a/workflow/default_parameters/v24_2_2_1/defaults.yaml b/workflow/default_parameters/v24_2_2_1/defaults.yaml index e6c38a4..8e73bf5 100644 --- a/workflow/default_parameters/v24_2_2_1/defaults.yaml +++ b/workflow/default_parameters/v24_2_2_1/defaults.yaml @@ -1,3 +1,4 @@ +--- resolution: resolution: 0.1 emod3d: @@ -100,8 +101,20 @@ hf: stoch_dx: 2.0 stoch_dy: 2.0 srf: - genslip_version: 5.4.2 + side_taper: 0.02 + bot_taper: 0.02 + top_taper: 0.0 + alpha_rough: 0.0 + gwid: [] + rvfac_seg: [] + seg_delay: false + ymag_exp: null + xmag_exp: null + kx_corner: null + ky_corner: null resolution: 0.1 + slip_sigma: 0.75 + risetime_coef: 1.6 # The following parameters are only used for the point source approximation. point_source_params: # Valid options for the slip time function (stype) are: diff --git a/workflow/default_parameters/v24_2_2_2/defaults.yaml b/workflow/default_parameters/v24_2_2_2/defaults.yaml index c71beb5..75289db 100644 --- a/workflow/default_parameters/v24_2_2_2/defaults.yaml +++ b/workflow/default_parameters/v24_2_2_2/defaults.yaml @@ -1,3 +1,4 @@ +--- resolution: resolution: 0.2 emod3d: @@ -100,8 +101,20 @@ hf: stoch_dx: 2.0 stoch_dy: 2.0 srf: - genslip_version: 5.4.2 + side_taper: 0.02 + bot_taper: 0.02 + top_taper: 0.0 + alpha_rough: 0.0 + gwid: [] + rvfac_seg: [] + seg_delay: false + ymag_exp: null + xmag_exp: null + kx_corner: null + ky_corner: null resolution: 0.1 + slip_sigma: 0.75 + risetime_coef: 1.6 # The following parameters are only used for the point source approximation. point_source_params: # Valid options for the slip time function (stype) are: diff --git a/workflow/default_parameters/v24_2_2_4/defaults.yaml b/workflow/default_parameters/v24_2_2_4/defaults.yaml index 1113b6d..7bb2a73 100644 --- a/workflow/default_parameters/v24_2_2_4/defaults.yaml +++ b/workflow/default_parameters/v24_2_2_4/defaults.yaml @@ -1,3 +1,4 @@ +--- resolution: resolution: 0.4 emod3d: @@ -100,8 +101,20 @@ hf: stoch_dx: 2.0 stoch_dy: 2.0 srf: - genslip_version: 5.4.2 + side_taper: 0.02 + bot_taper: 0.02 + top_taper: 0.0 + alpha_rough: 0.0 + gwid: [] + rvfac_seg: [] + seg_delay: false + ymag_exp: null + xmag_exp: null + kx_corner: null + ky_corner: null resolution: 0.1 + slip_sigma: 0.75 + risetime_coef: 1.6 # The following parameters are only used for the point source approximation. point_source_params: # Valid options for the slip time function (stype) are: diff --git a/workflow/realisations.py b/workflow/realisations.py index 3a5603b..24b3d83 100644 --- a/workflow/realisations.py +++ b/workflow/realisations.py @@ -370,15 +370,44 @@ class SRFConfig(RealisationConfiguration): _config_key: ClassVar[str] = "srf" _schema: ClassVar[Schema] = schemas.SRF_SCHEMA - genslip_version: str - """The version of genslip to use (currently supports "5.4.2").""" - resolution: float """The resolution of the SRF discretisation (different, in general, from the simulation resolution).""" point_source_params: schemas.PointSourceParams | None """Parameters for point source approximation, if applicable.""" + side_taper: float + """Side slip tapering, proportion of along-strike 0-1.""" + bot_taper: float + """Bottom slip tapering proportion of down-dip 0-1.""" + top_taper: float + """Top slip tapering proportion of down-dip 0-1.""" + + alpha_rough: float + """Roughness (0.0 = smooth)""" + + gwid: list[float] + """Width of delay zone around segment boundaries""" + rvfac_seg: list[float] + """Rupture speed reduction (proportion) at segment boundaries.""" + seg_delay: bool + """If true, delay rupture across slip boundaries according to specifications of rvfac_seg and gwid.""" + + slip_sigma: float + """The stddev of slip distribution.""" + + risetime_coef: float + """Risetime scaling coefficient.""" + + ymag_exp: float | None = None + """Corner magnitude exponent for along-strike slip correlation length. See genslip_v5.6.2c:1385""" + xmag_exp: float | None = None + """Corner magnitude exponent for down-dip slip correlation length. See genslip_v5.6.2c:1385""" + kx_corner: float | None = None + """Corner wavenumber for along-strike slip correlation length. See genslip_v5.6.2c:1385""" + ky_corner: float | None = None + """Corner wavenumber for down-dip slip correlation length. See genslip_v5.6.2c:1385""" + def to_dict(self) -> dict[str, Any]: """ Convert the object to a dictionary representation. diff --git a/workflow/schemas.py b/workflow/schemas.py index 54fa0ab..3538feb 100644 --- a/workflow/schemas.py +++ b/workflow/schemas.py @@ -18,28 +18,6 @@ from velocity_modelling.bounding_box import BoundingBox from workflow.defaults import DefaultsVersion -# NOTE: These functions seem silly and short, however there is a good -# reason for the choice to create functions like this. The reason is -# because when the schema library reports an error (such as the input -# file having a negative depth value) it prints the name of the -# function. So -# -# And(float, lambda x: x > 0).validate(-12) -# -# Would report the error -# -# schema.SchemaError: (-12) should evaluate to True -# -# But using the function `is_positive` we instead have -# -# And(float, is_positive).validate(-12) -# schema.SchemaError: is_positive(-12) should evaluate to True -# -# So using short functions with names improves the error reporting -# from the library. -# -# Accordingly, the most trivial of these functions lack docstrings. - class Stype(StrEnum): """Options for slip time function (stype) in generic_slip2srf.""" @@ -94,6 +72,29 @@ def from_dict(cls, params_dict: dict) -> "PointSourceParams": return cls(**params_dict) +# NOTE: These functions seem silly and short, however there is a good +# reason for the choice to create functions like this. The reason is +# because when the schema library reports an error (such as the input +# file having a negative depth value) it prints the name of the +# function. So +# +# And(float, lambda x: x > 0).validate(-12) +# +# Would report the error +# +# schema.SchemaError: (-12) should evaluate to True +# +# But using the function `is_positive` we instead have +# +# And(float, is_positive).validate(-12) +# schema.SchemaError: is_positive(-12) should evaluate to True +# +# So using short functions with names improves the error reporting +# from the library. +# +# Accordingly, the most trivial of these functions lack docstrings. + + def _is_positive(x: float) -> bool: # noqa: D103 # numpydoc ignore=GL08 return x > 0 @@ -130,6 +131,10 @@ def _is_valid_bearing(bearing: float) -> bool: # noqa: D103 # numpydoc ignore=G return 0 <= bearing <= 360 +def _is_proportion(x: float | int) -> bool: # noqa: D103 # numpydoc ignore=GL08 + return 0 <= x <= 1 + + def _is_correct_corner_shape(corners: np.ndarray) -> bool: """Check if the corner shape matches the corner shape for plane sources. @@ -373,21 +378,69 @@ def _corners_to_array(corners_spec: list[dict[str, float]]) -> np.ndarray: ) ) - SRF_SCHEMA = Schema( { - Literal("genslip_version", description="The version of genslip to use"): Or( - "5.4.2" - ), - Literal("resolution", description="The resolution of the SRF discretisation."): And( - NUMBER, _is_positive - ), + Literal( + "resolution", description="The resolution of the SRF discretisation." + ): And(NUMBER, _is_positive), Optional( Literal( "point_source_params", description="Parameters for point source approximation, if applicable", ) ): Or(None, POINT_SOURCE_PARAMS_SCHEMA), + Literal( + "side_taper", + description="Side slip tapering, proportion of along-strike 0-1.", + ): And(NUMBER, _is_proportion), + Literal( + "bot_taper", + description="Bottom slip tapering proportion of down-dip 0-1.", + ): And(NUMBER, _is_proportion), + Literal( + "top_taper", + description="Top slip tapering proportion of down-dip 0-1.", + ): And(NUMBER, _is_proportion), + Literal( + "alpha_rough", + description="Roughness (0.0 = smooth)", + ): And(NUMBER, _is_proportion), + Literal( + "gwid", + description="Width of delay zone around segment boundaries", + ): [And(NUMBER, _is_positive)], + Literal( + "rvfac_seg", + description="Rupture speed reduction (proportion) at segment boundaries.", + ): [And(NUMBER, _is_proportion)], + Literal( + "seg_delay", + description="If true, delay rupture across slip boundaries according to specifications of rvfac_seg and gwid.", + ): bool, + Literal( + "ymag_exp", + description="Corner magnitude exponent for along-strike slip correlation length. See genslip_v5.6.2c:1385", + ): Or(NUMBER, None), + Literal( + "xmag_exp", + description="Corner magnitude exponent for down-dip slip correlation length. See genslip_v5.6.2c:1385", + ): Or(NUMBER, None), + Literal( + "kx_corner", + description="Corner wavenumber for along-strike slip correlation length. See genslip_v5.6.2c:1385", + ): Or(NUMBER, None), + Literal( + "ky_corner", + description="Corner wavenumber for down-dip slip correlation length. See genslip_v5.6.2c:1385", + ): Or(NUMBER, None), + Literal( + "slip_sigma", + description="The stddev of slip distribution.", + ): And(NUMBER, _is_non_negative), + Literal( + "risetime_coef", + description="Risetime scaling coefficient.", + ): And(NUMBER, _is_positive), } ) diff --git a/workflow/scripts/realisation_to_srf.py b/workflow/scripts/realisation_to_srf.py index 8ddc932..c4569ee 100644 --- a/workflow/scripts/realisation_to_srf.py +++ b/workflow/scripts/realisation_to_srf.py @@ -399,6 +399,95 @@ def velocity_model_path(self) -> Path: # numpydoc ignore=RT01 return self.work_directory / "velocity_model" +def _build_genslip_command( + genslip_path: Path, + gsf_file_path: Path, + nx: int, + ny: int, + seed: int, + velocity_model_path: Path, + shypo: float, + dhypo: float, + magnitude: float, + dt: float, + srf_config: SRFConfig, +) -> list[str]: + """Build a genslip command to run. + + Parameters + ---------- + genslip_path : Path + Path to genslip executable. + gsf_file_path : Path + GSF file to read. + nx : int + Number of points along-strike. + ny : int + Number of points down-dip + seed : int + Genslip seed. + velocity_model_path : Path + Path to 1D velocity model. + shypo : float + Along-strike hypocentre. + dhypo : float + Down-dip hypocentre. + magnitude : float + SRF magnitude. + dt : float + Time resolution for SVFs. + srf_config : SRFConfig + Parameters for SRF generation. + + + Returns + ------- + list[str] + An argument list suitable to execute genslip with + `subprocess.run`. First element is the genslip path, + subsequent arguments are key=value pairs. + """ + cmd = [ + str(genslip_path), + "plane_header=1", + "srf_version=1.0", + "read_erf=0", + "write_srf=1", + "read_gsf=1", + "write_gsf=0", + "ns=1", + "nh=1", + f"infile={gsf_file_path}", + f"nstk={nx}", + f"ndip={ny}", + f"seed={seed}", + f"velfile={velocity_model_path}", + f"shypo={shypo}", + f"dhypo={dhypo}", + f"mag={magnitude}", + f"dt={dt}", + ] + skipped_fields = {"point_source_params"} + for field in dataclasses.fields(srf_config): + key = field.name + value = getattr(srf_config, key) + + if value is None or value == [] or key in skipped_fields: + continue + if isinstance(value, list): + serialised_value = ",".join([str(v) for v in value]) + elif isinstance(value, bool): + serialised_value = "1" if value else "0" + elif dataclasses.is_dataclass(value): + continue + else: + serialised_value = str(value) + + cmd.append(f"{key}={serialised_value}") + + return cmd + + def generate_fault_srf( name: str, params: SRFRealisationContext, environment: SRFEnvironmentContext ) -> None: @@ -439,34 +528,20 @@ def generate_fault_srf( genslip_hypocentre_coords = np.array([fault.length, fault.width]) * ( params.rupture_propagation_config.hypocentres[name] - np.array([1 / 2, 0]) ) - genslip_cmd = [ - str(environment.genslip_path), - "read_erf=0", - "write_srf=1", - "read_gsf=1", - "write_gsf=0", - f"infile={gsf_file_path}", - f"mag={params.magnitudes.magnitudes[name]}", - f"nstk={nx}", - f"ndip={ny}", - "ns=1", - "nh=1", - f"seed={environment.seeds.genslip_seed}", - f"velfile={environment.velocity_model_path}", - f"shypo={genslip_hypocentre_coords[0]}", - f"dhypo={genslip_hypocentre_coords[1]}", - f"dt={params.resolution.dt}", - "plane_header=1", - "srf_version=1.0", - "seg_delay={0}", - "rvfac_seg=-1", - "gwid=-1", - "side_taper=0.02", - "bot_taper=0.02", - "top_taper=0.0", - "rup_delay=0", - "alpha_rough=0.0", - ] + + genslip_cmd = _build_genslip_command( + genslip_path=environment.genslip_path, + gsf_file_path=gsf_file_path, + nx=nx, + ny=ny, + seed=environment.seeds.genslip_seed, + velocity_model_path=environment.velocity_model_path, + shypo=genslip_hypocentre_coords[0], + dhypo=genslip_hypocentre_coords[1], + magnitude=params.magnitudes.magnitudes[name], + dt=params.resolution.dt, + srf_config=params.srf_config, + ) srf_file_path = environment.srf_directory / (normalise_name(name) + ".srf") with open(srf_file_path, "w", encoding="utf-8") as srf_file_handle: