diff --git a/debian/changelog b/debian/changelog index ffb15c987c8f..f2f0ceca11f3 100644 --- a/debian/changelog +++ b/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 diff --git a/openquake/hazardlib/gsim/can20/can_shm6_active_crust.py b/openquake/hazardlib/gsim/can20/can_shm6_active_crust.py index 9dff51498af5..530645fc1dfc 100644 --- a/openquake/hazardlib/gsim/can20/can_shm6_active_crust.py +++ b/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 # ============================================================================= diff --git a/openquake/hazardlib/gsim/chiou_youngs_2014.py b/openquake/hazardlib/gsim/chiou_youngs_2014.py index eed9cd43ca8c..883f0f59d6d5 100644 --- a/openquake/hazardlib/gsim/chiou_youngs_2014.py +++ b/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,7 +26,7 @@ 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 @@ -37,21 +34,13 @@ 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) \ No newline at end of file + mean[m] = _get_mean(ctx, C, ln_y_ref, exp1, exp2) diff --git a/openquake/hazardlib/gsim/nz22/stafford_2022.py b/openquake/hazardlib/gsim/nz22/stafford_2022.py index c29a393b9fec..ff351ffbbc02 100644 --- a/openquake/hazardlib/gsim/nz22/stafford_2022.py +++ b/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) ) diff --git a/openquake/hazardlib/tests/gsim/chiou_youngs_2014_test.py b/openquake/hazardlib/tests/gsim/chiou_youngs_2014_test.py index 9940f9251d47..1c4eceecdd7d 100644 --- a/openquake/hazardlib/tests/gsim/chiou_youngs_2014_test.py +++ b/openquake/hazardlib/tests/gsim/chiou_youngs_2014_test.py @@ -20,8 +20,7 @@ import numpy as np import copy from openquake.hazardlib.gsim.chiou_youngs_2014 import ( - ChiouYoungs2014, ChiouYoungs2014PEER, ChiouYoungs2014NearFaultEffect, - ChiouYoungs2014Japan, ChiouYoungs2014Italy, ChiouYoungs2014Wenchuan) + ChiouYoungs2014, ChiouYoungs2014PEER, ChiouYoungs2014NearFaultEffect) from openquake.hazardlib.gsim.chiou_youngs_2014 import ( _get_delta_cm, get_magnitude_scaling, _get_delta_g) @@ -90,34 +89,18 @@ def test_total_stddev_refactor(self): # Note that in the regionalisation cases the discrepancy percentage is raised # to 1 % to allow for a different interpretation of the deltaZ1.0 when Z1.0 = 0 # when compared to the verification code -class ChiouYoungs2014JapanTestCase(BaseGSIMTestCase): - GSIM_CLASS = ChiouYoungs2014Japan +# Data generated from implementation from Yue Hua +# https://web.stanford.edu/~bakerjw/GMPEs/CY_2014_nga.m +class ChiouYoungs2014RegionTestCase(BaseGSIMTestCase): + GSIM_CLASS = ChiouYoungs2014 - def test_mean_japan(self): - # Data generated from implementation from Yue Hua - # https://web.stanford.edu/~bakerjw/GMPEs/CY_2014_nga.m + def test_mean(self): self.check('NGA/CY14/CY14_Japan_MEAN.csv', - max_discrep_percentage=1.0) - - -class ChiouYoungs2014ItalyTestCase(BaseGSIMTestCase): - GSIM_CLASS = ChiouYoungs2014Italy - - def test_mean_italy(self): - # Data generated from implementation from Yue Hua - # https://web.stanford.edu/~bakerjw/GMPEs/CY_2014_nga.m + max_discrep_percentage=1.0, region='JPN') self.check('NGA/CY14/CY14_Italy_MEAN.csv', - max_discrep_percentage=1.0) - - -class ChiouYoungs2014WenchuanTestCase(BaseGSIMTestCase): - GSIM_CLASS = ChiouYoungs2014Wenchuan - - def test_mean_wenchuan(self): - # Data generated from implementation from Yue Hua - # https://web.stanford.edu/~bakerjw/GMPEs/CY_2014_nga.m + max_discrep_percentage=1.0, region='ITA') self.check('NGA/CY14/CY14_Wenchuan_MEAN.csv', - max_discrep_percentage=1.0) + max_discrep_percentage=1.0, region='WEN') class ChiouYoungs2014PEERTestCase(BaseGSIMTestCase):