Merge pull request #10025 from gem/CY14

Modernized the regionalization of Chiou Youngs (2014)

debian/changelog

@@ -1,4 +1,6 @@
   [Michele Simionato]
+  * Modernized the regionalization of Chiou Youngs (2014)
+  * Modernized the regionalization of Campbell Bozorgnia (2014)
   * Internal: made it possible to override CoeffsTable
   * Added a memory check in disaggregation calculations
   * Made `scientific_format` resilient against encoding errors

openquake/hazardlib/gsim/can20/can_shm6_active_crust.py

@@ -162,19 +162,19 @@ def shm6_site_correction(C, mean, ctx, imt):
     mean[vs30_ge1100] = factor + cy14_760 + mean[vs30_ge1100]
 
 
-def get_mean_stddevs_cy14(name, C, ctx, conf):
+def get_mean_stddevs_cy14(region, C, ctx, conf):
     """
     Return mean and standard deviation values
     """
     # Get ground motion on reference rock
-    ln_y_ref = CY14.get_ln_y_ref(name, C, ctx, conf)
+    ln_y_ref = CY14.get_ln_y_ref(region, C, ctx, conf)
     y_ref = np.exp(ln_y_ref)
 
     # Set basin depth to 0
     f_z1pt0 = 0.0
 
     # Get linear amplification term
-    f_lin = CY14.get_linear_site_term(name, C, ctx)
+    f_lin = CY14.get_linear_site_term(region, C, ctx)
 
     # Get nonlinear amplification term
     f_nl, f_nl_scaling = CY14.get_nonlinear_site_term(C, ctx, y_ref)
@@ -184,7 +184,7 @@ def get_mean_stddevs_cy14(name, C, ctx, conf):
 
     # Get standard deviations
     sig, tau, phi = CY14.get_stddevs(
-        name, C, ctx, ctx.mag, y_ref, f_nl_scaling)
+        conf['peer'], C, ctx, ctx.mag, y_ref, f_nl_scaling)
 
     return mean, sig, tau, phi
 
@@ -215,14 +215,13 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
         Notes:
             - Centered cdpp already defaulted to 0
         """
-        name = self.__class__.__name__
 
         # Checking the IMTs used to compute ground-motion
         _check_imts(imts)
 
         # Reference to page 1144, PSA might need PGA value
         pga_mean, pga_sig, pga_tau, pga_phi = get_mean_stddevs_cy14(
-            name, self.COEFFS[PGA()], ctx, self.conf)
+            self.region, self.COEFFS[PGA()], ctx, self.conf)
 
         # Processing IMTs
         for m, imt in enumerate(imts):
@@ -234,8 +233,8 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
 
                 sig[m], tau[m], phi[m] = pga_sig, pga_tau, pga_phi
             else:
-                imt_mean, imt_sig, imt_tau, imt_phi = \
-                    get_mean_stddevs_cy14(name, self.COEFFS[imt], ctx, self.conf)
+                imt_mean, imt_sig, imt_tau, imt_phi = get_mean_stddevs_cy14(
+                    self.region, self.COEFFS[imt], ctx, self.conf)
                 # reference to page 1144
                 # Predicted PSA value at T ≤ 0.3s should be set equal to the
                 # value of PGA when it falls below the predicted PGA
@@ -245,7 +244,6 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
 
                 # Site term correction fos SHM6
                 shm6_site_correction(self.COEFFS[imt], mean[m], ctx, imt)
-
                 sig[m], tau[m], phi[m] = imt_sig, imt_tau, imt_phi
 
 # =============================================================================

openquake/hazardlib/gsim/chiou_youngs_2014.py

@@ -18,9 +18,6 @@
 
 """
 Module exports :class:`ChiouYoungs2014`
-               :class:`ChiouYoungs2014Japan`
-               :class:`ChiouYoungs2014Italy`
-               :class:`ChiouYoungs2014Wenchuan`
                :class:`ChiouYoungs2014PEER`
                :class:`ChiouYoungs2014NearFaultEffect`
 """
@@ -29,29 +26,21 @@
 import numpy as np
 
 from openquake.baselib.general import CallableDict
-from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
+from openquake.hazardlib.gsim.base import GMPE, CoeffsTable, add_alias
 from openquake.hazardlib.gsim.abrahamson_2014 import get_epistemic_sigma
 from openquake.hazardlib import const
 from openquake.hazardlib.imt import PGA, PGV, SA
 
 CONSTANTS = {"c2": 1.06, "c4": -2.1, "c4a": -0.5, "crb": 50.0,
              "c8a": 0.2695, "c11": 0.0, "phi6": 300.0, "phi6jp": 800.0}
 
-def _get_centered_cdpp(clsname, ctx):
-    """
-    Returns the centred dpp term (zero by default)
-    """
-    if clsname.endswith("NearFaultEffect"):
-        return ctx.rcdpp
-    return np.zeros(ctx.rrup.shape)
-
 
-def _get_centered_z1pt0(clsname, ctx):
+def _get_centered_z1pt0(region, ctx):
     """
     Get z1pt0 centered on the Vs30- dependent average z1pt0(m)
     California and non-Japan regions
     """
-    if clsname.endswith("Japan"):
+    if region == "JPN":
         mean_z1pt0 = (-5.23 / 2.) * np.log(((ctx.vs30 ** 2.) + 412.39 ** 2.)
                                            / (1360 ** 2. + 412.39 ** 2.))
         return ctx.z1pt0 - np.exp(mean_z1pt0)
@@ -151,30 +140,31 @@ def _get_mean(ctx, C, ln_y_ref, exp1, exp2):
     return ln_y
 
 
-def get_basin_depth_term(clsname, C, centered_z1pt0):
+def get_basin_depth_term(region, C, centered_z1pt0):
     """
     Returns the basin depth scaling
     """
-    if clsname.endswith("Japan"):
+    if region == "JPN":
         return C["phi5jp"] * (1.0 - np.exp(-centered_z1pt0 /
                                            CONSTANTS["phi6jp"]))
     return C["phi5"] * (1.0 - np.exp(-centered_z1pt0 /
                                      CONSTANTS["phi6"]))
 
 
-def get_directivity(clsname, C, ctx):
+def get_directivity(C, ctx):
     """
-    Returns the directivity term.
+    Returns the directivity term, if any.
 
     The directivity prediction parameter is centered on the average
     directivity prediction parameter. Here we set the centered_dpp
     equal to zero, since the near fault directivity effect prediction is
     off by default in our calculation.
-    """
-    cdpp = _get_centered_cdpp(clsname, ctx)
-    if not np.any(cdpp > 0.0):
+    """    
+    try:
+        cdpp =  ctx.rcdpp
+    except AttributeError:
         # No directivity term
-        return 0.0
+        return 0.
     f_dir = np.exp(-C["c8a"] * ((ctx.mag - C["c8b"]) ** 2.)) * cdpp
     f_dir *= np.clip((ctx.mag - 5.5) / 0.8, 0., 1.)
     rrup_max = ctx.rrup - 40.
@@ -293,28 +283,14 @@ def get_hanging_wall_term(C, ctx):
     return fhw
 
 
-def get_linear_site_term(clsname, C, ctx):
+def get_linear_site_term(region, C, ctx):
     """
     Returns the linear site scaling term
     """
-    if clsname.endswith("Japan"):
+    if region == "JPN":
         return C["phi1jp"] * np.log(ctx.vs30 / 1130).clip(-np.inf, 0.0)
     return C["phi1"] * np.log(ctx.vs30 / 1130).clip(-np.inf, 0.0)
 
-
-def get_region(clsname):
-    """
-    Returns the region parameter
-    """
-    if clsname.endswith("Italy"):
-        return "ITA"
-    elif clsname.endswith("Japan"):
-        return "JPN"
-    elif clsname.endswith("Wenchuan"):
-        return "WEN"
-    else:
-        return "CAL"
-
 
 def get_delta_c1(rrup, imt, mag):
     """
@@ -393,7 +369,7 @@ def _get_delta_g(delta_gamma_tab, ctx, imt):
     return delta_g
 
 
-def get_ln_y_ref(clsname, C, ctx, conf):
+def get_ln_y_ref(region, C, ctx, conf):
     """
     Returns the ground motion on the reference rock, described fully by
     Equation 11 in CY14 (page 1131).
@@ -405,7 +381,6 @@ def get_ln_y_ref(clsname, C, ctx, conf):
     alpha_nm = conf.get('alpha_nm')
 
     # Get the region name from the name of the class
-    region = get_region(clsname)
     delta_ztor = _get_centered_ztor(ctx)
 
     # If stress correction required get delta cm
@@ -434,7 +409,7 @@ def get_ln_y_ref(clsname, C, ctx, conf):
            get_geometric_spreading(C, ctx.mag, ctx.rrup) +
            get_far_field_distance_scaling(
                region, C, ctx.mag, ctx.rrup, delta_g) +
-           get_directivity(clsname, C, ctx))
+           get_directivity(C, ctx))
 
     # Adjust ground-motion for the hanging wall effect
     if use_hw:
@@ -510,11 +485,11 @@ def get_source_scaling_terms(C, ctx, delta_ztor, alpha_nm):
     return f_src
 
 
-def get_stddevs(clsname, C, ctx, mag, y_ref, f_nl_scaling):
+def get_stddevs(peer, C, ctx, mag, y_ref, f_nl_scaling):
     """
     Returns the standard deviation model described in equation 13
     """
-    if clsname == 'ChiouYoungs2014PEER':
+    if peer:
         # the standard deviation, which is fixed at 0.65 for every site
         return [0.65 * np.ones_like(ctx.vs30), 0, 0]
 
@@ -545,23 +520,23 @@ def get_tau(C, mag):
     return C['tau1'] + (C['tau2'] - C['tau1']) / 1.5 * mag_test
 
 
-def get_mean_stddevs(name, C, ctx, imt, conf):
+def get_mean_stddevs(region, C, ctx, imt, conf):
     """
     Return mean and standard deviation values
     """
     # Get ground motion on reference rock
-    ln_y_ref = get_ln_y_ref(name, C, ctx, conf)
+    ln_y_ref = get_ln_y_ref(region, C, ctx, conf)
     y_ref = np.exp(ln_y_ref)
 
     # Get basin depth
-    dz1pt0 = _get_centered_z1pt0(name, ctx)
+    dz1pt0 = _get_centered_z1pt0(region, ctx)
 
     # for Z1.0 = 0.0 no deep soil correction is applied
     dz1pt0[ctx.z1pt0 <= 0.0] = 0.0
-    f_z1pt0 = get_basin_depth_term(name, C, dz1pt0)
+    f_z1pt0 = get_basin_depth_term(region, C, dz1pt0)
 
     # Get linear amplification term
-    f_lin = get_linear_site_term(name, C, ctx)
+    f_lin = get_linear_site_term(region, C, ctx)
 
     # Get nonlinear amplification term
     f_nl, f_nl_scaling = get_nonlinear_site_term(C, ctx, y_ref)
@@ -571,7 +546,7 @@ def get_mean_stddevs(name, C, ctx, imt, conf):
 
     # Get standard deviations
     sig, tau, phi = get_stddevs(
-        name, C, ctx, ctx.mag, y_ref, f_nl_scaling)
+        conf['peer'], C, ctx, ctx.mag, y_ref, f_nl_scaling)
 
     return mean, sig, tau, phi
 
@@ -637,15 +612,19 @@ class ChiouYoungs2014(GMPE):
     #: Reference shear wave velocity
     DEFINED_FOR_REFERENCE_VELOCITY = 1130
 
-    def __init__(self, sigma_mu_epsilon=0.0, use_hw=True, add_delta_c1=False,
-                 alpha_nm=1.0, stress_par_host=None, stress_par_target=None,
-                 delta_gamma_tab=None):
-
-        # Add sigma_mu_epsilon 
+    def __init__(self, region='CAL', sigma_mu_epsilon=0.0, use_hw=True,
+                 add_delta_c1=False, alpha_nm=1.0, stress_par_host=None,
+                 stress_par_target=None, delta_gamma_tab=None):
+
+        # set region
+        self.region = region
+
+        # set sigma_mu_epsilon 
         self.sigma_mu_epsilon = sigma_mu_epsilon
 
-        # Adding into the conf dictionary
+        # set the conf dictionary
         self.conf = {}
+        self.conf['peer'] = self.__class__.__name__.endswith('PEER')
         self.conf['use_hw'] = use_hw 
         self.conf['alpha_nm'] = alpha_nm
         self.conf['add_delta_c1'] = add_delta_c1
@@ -690,11 +669,10 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
         <.base.GroundShakingIntensityModel.compute>`
         for spec of input and result values.
         """
-        name = self.__class__.__name__
         # Reference to page 1144, PSA might need PGA value
         self.conf['imt'] = PGA()
         pga_mean, pga_sig, pga_tau, pga_phi = get_mean_stddevs(
-            name, self.COEFFS[PGA()], ctx, PGA(), self.conf)
+            self.region, self.COEFFS[PGA()], ctx, PGA(), self.conf)
         # compute
         for m, imt in enumerate(imts):
             self.conf['imt'] = imt
@@ -704,7 +682,7 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
                 sig[m], tau[m], phi[m] = pga_sig, pga_tau, pga_phi
             else:
                 imt_mean, imt_sig, imt_tau, imt_phi = get_mean_stddevs(
-                    name, self.COEFFS[imt], ctx, imt, self.conf)
+                    self.region, self.COEFFS[imt], ctx, imt, self.conf)
                 # Reference to page 1144
                 # Predicted PSA value at T ≤ 0.3s should be set equal to the
                 # value of PGA when it falls below the predicted PGA
@@ -747,28 +725,20 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
 """)
 
 
-class ChiouYoungs2014Japan(ChiouYoungs2014):
-    """
-    Regionalisation of the Chiou & Youngs (2014) GMPE for use with the
-    Japan far-field distance attuation scaling and site model
-    """
-
-
-class ChiouYoungs2014Italy(ChiouYoungs2014):
-    """
-    Adaption of the Chiou & Youngs (2014) GMPE for the the Italy far-field
-    attenuation scaling, but assuming the California site amplification model
-    """
+# Regionalisation of the Chiou & Youngs (2014) GMPE for use with the
+# Japan far-field distance attuation scaling and site model
+add_alias('ChiouYoungs2014Japan', ChiouYoungs2014, region='JPN')
 
+# Adaption of the Chiou & Youngs (2014) GMPE for the the Italy far-field
+# attenuation scaling, but assuming the California site amplification model
+add_alias('ChiouYoungs2014Italy', ChiouYoungs2014, region='ITA')
 
-class ChiouYoungs2014Wenchuan(ChiouYoungs2014):
-    """
-    Adaption of the Chiou & Youngs (2014) GMPE for the Wenchuan far-field
-    attenuation scaling, but assuming the California site amplification model.
-    It should be note that according to Chiou & Youngs (2014) this adjustment
-    is calibrated only for the M7.9 Wenchuan earthquake, so application to
-    other scenarios is at the user's own risk
-    """
+# Adaption of the Chiou & Youngs (2014) GMPE for the Wenchuan far-field
+# attenuation scaling, but assuming the California site amplification model.
+# It should be note that according to Chiou & Youngs (2014) this adjustment
+# is calibrated only for the M7.9 Wenchuan earthquake, so application to
+#  other scenarios is at the user's own risk
+add_alias('ChiouYoungs2014Wenchuan', ChiouYoungs2014, region='WEN')
 
 
 class ChiouYoungs2014PEER(ChiouYoungs2014):
@@ -811,9 +781,9 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
             C = self.COEFFS[imt]
             # intensity on a reference soil is used for both mean
             # and stddev calculations.
-            ln_y_ref = _get_ln_y_ref(self.__class__.__name__, ctx, C)
+            ln_y_ref = _get_ln_y_ref(self.region, ctx, C)
             # exp1 and exp2 are parts of eq. 12 and eq. 13,
             # calculate it once for both.
             exp1 = np.exp(C['phi3'] * (ctx.vs30.clip(-np.inf, 1130) - 360))
             exp2 = np.exp(C['phi3'] * (1130 - 360))
-            mean[m] = _get_mean(ctx, C, ln_y_ref, exp1, exp2)
+            mean[m] = _get_mean(ctx, C, ln_y_ref, exp1, exp2)

openquake/hazardlib/gsim/nz22/stafford_2022.py

@@ -347,7 +347,7 @@ def get_ln_y_ref(
         + get_far_field_distance_scaling(
             T, C, ctx, mu_branch, adjust_c1, adjust_chm, adjust_c7, adjust_cg1
         )
-        + get_directivity("Stafford2022", C, ctx)
+        + get_directivity(C, ctx)
     )