diff --git a/debian/changelog b/debian/changelog index a4b5ba4e30ff..491574525f4c 100644 --- a/debian/changelog +++ b/debian/changelog @@ -1,3 +1,9 @@ + [Christopher Brooks] + * Added M9 Basin amplification term for ModifiableGMPE + and assoc. unit tests + * Added the CB14 basin amplification term for + ModifiableGMPE and assoc. unit tests + [Michele Simionato] * Added parameter `with_betw_ratio` * Added command `oq info peril` diff --git a/openquake/hazardlib/contexts.py b/openquake/hazardlib/contexts.py index 98355eb7de88..991aa2d67f4b 100644 --- a/openquake/hazardlib/contexts.py +++ b/openquake/hazardlib/contexts.py @@ -1686,7 +1686,7 @@ def __eq__(self, other): return False -# mock of a site collection used in the tests and in the SMTK +# mock of a site collection used in the tests and in the SMT class SitesContext(BaseContext): """ Sites calculation context for ground shaking intensity models. @@ -1699,7 +1699,7 @@ class SitesContext(BaseContext): Only those required parameters are made available in a result context object. """ - # _slots_ is used in hazardlib check_gsim and in the SMTK + # _slots_ is used in hazardlib check_gsim and in the SMT def __init__(self, slots='vs30 vs30measured z1pt0 z2pt5'.split(), sitecol=None): self._slots_ = slots @@ -1708,7 +1708,7 @@ def __init__(self, slots='vs30 vs30measured z1pt0 z2pt5'.split(), for slot in slots: setattr(self, slot, getattr(sitecol, slot)) - # used in the SMTK + # used in the SMT def __len__(self): return len(self.sids) @@ -1745,7 +1745,7 @@ def get_dists(ctx): # used to produce a RuptureContext suitable for legacy code, i.e. for calls -# to .get_mean_and_stddevs, like for instance in the SMTK +# to .get_mean_and_stddevs, like for instance in the SMT def full_context(sites, rup, dctx=None): """ :returns: a full RuptureContext with all the relevant attributes @@ -1785,7 +1785,7 @@ def get_mean_stds(gsim, ctx, imts, **kw): return out[:, 0] if single else out -# mock of a rupture used in the tests and in the SMTK +# mock of a rupture used in the tests and in the SMT class RuptureContext(BaseContext): """ Rupture calculation context for ground shaking intensity models. diff --git a/openquake/hazardlib/gsim/mgmpe/cb14_basin_term.py b/openquake/hazardlib/gsim/mgmpe/cb14_basin_term.py new file mode 100644 index 000000000000..acad2d59e72f --- /dev/null +++ b/openquake/hazardlib/gsim/mgmpe/cb14_basin_term.py @@ -0,0 +1,85 @@ +# The Hazard Library +# Copyright (C) 2012-2023 GEM Foundation +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Affero General Public License as +# published by the Free Software Foundation, either version 3 of the +# License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with this program. If not, see . + +""" +Module :mod:`openquake.hazardlib.mgmpe.cb14_basin_term` implements +:class:`~openquake.hazardlib.mgmpe.CB14BasinTerm` +""" +import numpy as np + +from openquake.hazardlib import const +from openquake.hazardlib.gsim.base import GMPE, registry +from openquake.hazardlib.gsim.campbell_bozorgnia_2014 import CampbellBozorgnia2014 + + +def _get_cb14_basin_term(ctx, C, jpn_flag=False): + """ + Get the basin response term defined in equation 20 of the Campbell and + Bozorgnia (2014) GMM paper. + + Currently the global basin term is provided (i.e. the Japan-regionalised + basin term is for now turned off). + """ + z2pt5 = ctx.z2pt5 + fb = np.zeros(len(z2pt5)) + idx = z2pt5 < 1.0 + fb[idx] = (C["c14"] + C["c15"] * jpn_flag) * (z2pt5[idx] - 1.0) + idx = z2pt5 > 3.0 + fb[idx] = C["c16"] * C["k3"] * np.exp(-0.75) * ( + 1. - np.exp(-0.25 * (z2pt5[idx] - 3.))) + + return fb + + +class CB14BasinTerm(GMPE): + """ + Implements a modified GMPE class that can be used to implement the Campbell + and Bozorgnia (2014) GMM's basin term + + :param gmpe_name: + The name of a GMPE class + """ + # Req Params + REQUIRES_SITES_PARAMETERS = {'z2pt5'} + + # Others are set from underlying GMM + REQUIRES_DISTANCES = set() + REQUIRES_RUPTURE_PARAMETERS = set() + DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = "" + DEFINED_FOR_INTENSITY_MEASURE_TYPES = set() + DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL} + DEFINED_FOR_TECTONIC_REGION_TYPE = "" + DEFINED_FOR_REFERENCE_VELOCITY = None + + def __init__(self, gmpe_name, **kwargs): + self.gmpe = registry[gmpe_name]() + self.set_parameters() + + # Need z2pt5 in req site params to ensure in ctx site col + if 'z2pt5' not in self.gmpe.REQUIRES_SITES_PARAMETERS: + self.REQUIRES_SITES_PARAMETERS = frozenset( + self.gmpe.REQUIRES_SITES_PARAMETERS | {'z2pt5'}) + + def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi): + """ + See :meth:`superclass method + <.base.GroundShakingIntensityModel.compute>` + for spec of input and result values. + """ + self.gmpe.compute(ctx, imts, mean, sig, tau, phi) + for m, imt in enumerate(imts): + C = CampbellBozorgnia2014.COEFFS[imt] + mean[m] += _get_cb14_basin_term(ctx, C) \ No newline at end of file diff --git a/openquake/hazardlib/gsim/mgmpe/cy14_site_term.py b/openquake/hazardlib/gsim/mgmpe/cy14_site_term.py index c13f4bd6e455..cd73fbff459b 100644 --- a/openquake/hazardlib/gsim/mgmpe/cy14_site_term.py +++ b/openquake/hazardlib/gsim/mgmpe/cy14_site_term.py @@ -26,7 +26,7 @@ from openquake.hazardlib.gsim.chiou_youngs_2014 import ChiouYoungs2014 -def _get_site_term(C, vs30, ln_y_ref): +def _get_cy14_site_term(C, vs30, ln_y_ref): """ Applies the linear and nonlinear site amplification term of Chiou & Youngs (2014) (excluding the basin amplification term) @@ -104,4 +104,4 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi): vs30 = ctx.vs30.copy() for m, imt in enumerate(imts): C = ChiouYoungs2014.COEFFS[imt] - mean[m] += _get_site_term(C, vs30, mean[m]) + mean[m] += _get_cy14_site_term(C, vs30, mean[m]) diff --git a/openquake/hazardlib/gsim/mgmpe/m9_basin_term.py b/openquake/hazardlib/gsim/mgmpe/m9_basin_term.py new file mode 100644 index 000000000000..e2df2291e25a --- /dev/null +++ b/openquake/hazardlib/gsim/mgmpe/m9_basin_term.py @@ -0,0 +1,77 @@ +# The Hazard Library +# Copyright (C) 2012-2023 GEM Foundation +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Affero General Public License as +# published by the Free Software Foundation, either version 3 of the +# License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with this program. If not, see . + +""" +Module :mod:`openquake.hazardlib.mgmpe.m9_basin_term` implements +:class:`~openquake.hazardlib.mgmpe.M9BasinTerm` +""" +import numpy as np + +from openquake.hazardlib import const +from openquake.hazardlib.gsim.base import GMPE, registry + + +def _apply_m9_basin_term(ctx, imt, mean): + if imt.period > 1.9: # Only apply to long-period SA + fb_m9 = np.log(2.0) + idx = ctx.z2pt5 >= 6.0 # Apply only to sites with z2pt5 >= 6 + mean[idx] += fb_m9 + + return mean + + +class M9BasinTerm(GMPE): + """ + Implements a modified GMPE class that can be used to implement the "M9" + US 2023 NSHM basin amplification adjustment. + + This implementation is based on the description of the M9 adjustment + within the Moschetti et al. (2024) EQ Spectra article on the conterminous + US 2023 NSHM GMC (pp. 1178). + + :param gmpe_name: + The name of a GMPE class + """ + # Req Params + REQUIRES_SITES_PARAMETERS = {'z2pt5'} + + # Others are set from underlying GMM + REQUIRES_DISTANCES = set() + REQUIRES_RUPTURE_PARAMETERS = set() + DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = "" + DEFINED_FOR_INTENSITY_MEASURE_TYPES = set() + DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL} + DEFINED_FOR_TECTONIC_REGION_TYPE = "" + DEFINED_FOR_REFERENCE_VELOCITY = None + + def __init__(self, gmpe_name, **kwargs): + self.gmpe = registry[gmpe_name]() + self.set_parameters() + + # Need z2pt5 in req site params to ensure in ctx site col + if 'z2pt5' not in self.gmpe.REQUIRES_SITES_PARAMETERS: + self.REQUIRES_SITES_PARAMETERS = frozenset( + self.gmpe.REQUIRES_SITES_PARAMETERS | {'z2pt5'}) + + def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi): + """ + See :meth:`superclass method + <.base.GroundShakingIntensityModel.compute>` + for spec of input and result values. + """ + self.gmpe.compute(ctx, imts, mean, sig, tau, phi) + for m, imt in enumerate(imts): + mean[m] = _apply_m9_basin_term(ctx, imt, mean[m]) \ No newline at end of file diff --git a/openquake/hazardlib/gsim/mgmpe/modifiable_gmpe.py b/openquake/hazardlib/gsim/mgmpe/modifiable_gmpe.py index c3f431ce3b55..9d5cf525be6f 100644 --- a/openquake/hazardlib/gsim/mgmpe/modifiable_gmpe.py +++ b/openquake/hazardlib/gsim/mgmpe/modifiable_gmpe.py @@ -26,8 +26,11 @@ from openquake.hazardlib.imt import from_string from openquake.hazardlib.gsim.mgmpe.nrcan15_site_term import ( NRCan15SiteTerm, BA08_AB06) -from openquake.hazardlib.gsim.mgmpe.cy14_site_term import _get_site_term +from openquake.hazardlib.gsim.mgmpe.cy14_site_term import _get_cy14_site_term from openquake.hazardlib.gsim.chiou_youngs_2014 import ChiouYoungs2014 +from openquake.hazardlib.gsim.mgmpe.cb14_basin_term import _get_cb14_basin_term +from openquake.hazardlib.gsim.campbell_bozorgnia_2014 import CampbellBozorgnia2014 +from openquake.hazardlib.gsim.mgmpe.m9_basin_term import _apply_m9_basin_term from openquake.hazardlib.gsim.nga_east import ( TAU_EXECUTION, get_phi_ss, TAU_SETUP, PHI_SETUP, get_tau_at_quantile, @@ -73,10 +76,25 @@ def cy14_site_term(ctx, imt, me, si, ta, phi): This function adds the CY14 site term to GMMs requiring it """ C = ChiouYoungs2014.COEFFS[imt] - fa = _get_site_term(C, ctx.vs30, me) # ref mean must be in natural log + fa = _get_cy14_site_term(C, ctx.vs30, me) # Ref mean must be in natural log me[:] += fa +def cb14_basin_term(ctx, imt, me, si, ta, phi): + """ + This function adds the CB14 basin term to GMMs requiring it. + """ + C = CampbellBozorgnia2014.COEFFS[imt] + me[:] += _get_cb14_basin_term(ctx, C) + + +def m9_basin_term(ctx, imt, me, si, ta, phi): + """ + This function applies the M9 basin adjustment + """ + me = _apply_m9_basin_term(ctx, imt, me) + + def add_between_within_stds(ctx, imt, me, si, ta, ph, with_betw_ratio): """ This adds the between and within standard deviations to a model which has @@ -248,6 +266,12 @@ def __init__(self, **kwargs): setattr(self, 'DEFINED_FOR_STANDARD_DEVIATION_TYPES', {StdDev.TOTAL, StdDev.INTRA_EVENT, StdDev.INTER_EVENT}) + if ('cb14_basin_term' in self.params or 'm9_basin_term' in self.params + ) and ( 'z2pt5' not in self.gmpe.REQUIRES_SITES_PARAMETERS): + tmp = list(self.gmpe.REQUIRES_SITES_PARAMETERS) + tmp.append('z2pt5') + self.gmpe.REQUIRES_SITES_PARAMETERS = frozenset(tmp) + # This is required by the `sigma_model_alatik2015` function key = 'sigma_model_alatik2015' if key in self.params: @@ -301,6 +325,7 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi): <.base.GroundShakingIntensityModel.compute>` for spec of input and result values. """ + # Set reference Vs30 if required if ('nrcan15_site_term' in self.params or 'cy14_site_term' in self.params): ctx_copy = ctx.copy() @@ -308,7 +333,7 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi): rock_vs30 = 760. elif 'cy14_site_term' in self.params: rock_vs30 = 1130. - ctx_copy.vs30 = np.full_like(ctx.vs30, rock_vs30) # rock + ctx_copy.vs30 = np.full_like(ctx.vs30, rock_vs30) # rock else: ctx_copy = ctx g = globals() diff --git a/openquake/hazardlib/tests/gsim/mgmpe/cb14_basin_term_test.py b/openquake/hazardlib/tests/gsim/mgmpe/cb14_basin_term_test.py new file mode 100644 index 000000000000..717343816af2 --- /dev/null +++ b/openquake/hazardlib/tests/gsim/mgmpe/cb14_basin_term_test.py @@ -0,0 +1,89 @@ +# The Hazard Library +# Copyright (C) 2012-2023 GEM Foundation +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Affero General Public License as +# published by the Free Software Foundation, either version 3 of the +# License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with this program. If not, see . + +import numpy as np +import unittest + +from openquake.hazardlib.tests.gsim.mgmpe.dummy import new_ctx +from openquake.hazardlib.contexts import simple_cmaker +from openquake.hazardlib import valid +from openquake.hazardlib.imt import PGA, SA +from openquake.hazardlib.const import TRT, IMC +from openquake.hazardlib.gsim.mgmpe.cb14_basin_term import CB14BasinTerm + +ae = np.testing.assert_equal +aae = np.testing.assert_almost_equal + +exp_gmm_origin = np.array([[-1.55067862, -3.8245714 , -5.48623229], + [-2.25723424, -3.72056054, -5.28051045], + [-3.45764676, -4.67077352, -5.84789307]]) + +exp_with_basin = np.array([[-1.35815033, -3.78220564, -5.293704 ], + [-1.86104115, -3.63337845, -4.88431736], + [-3.04727152, -4.58047065, -5.43751783]]) + +class M9BasinTermTestCase(unittest.TestCase): + + def test_instantiation(self): + mgmpe = CB14BasinTerm(gmpe_name='AtkinsonMacias2009') + + # Check the assigned IMTs + expected = set([PGA, SA]) + ae(mgmpe.DEFINED_FOR_INTENSITY_MEASURE_TYPES, expected) + # Check the TR + expected = TRT.SUBDUCTION_INTERFACE + ae(mgmpe.DEFINED_FOR_TECTONIC_REGION_TYPE, expected) + # Check the IM component + expected = IMC.RANDOM_HORIZONTAL + ae(mgmpe.DEFINED_FOR_INTENSITY_MEASURE_COMPONENT, expected) + # Check the required distances + expected = {'rrup'} + ae(mgmpe.REQUIRES_DISTANCES, expected) + + def test_all(self): + """ + Test that the M9 basin term applied to Kuehn et al. (2020) + provides the expected values (using sites with z2pt5 above + and below 6 km threshold and considering SAs with periods + above and below 1.9 s) + """ + am09 = valid.gsim('AtkinsonMacias2009') + + # Make original GMM + gmpe = valid.modified_gsim(am09) + + # Make GMM with basin term using the mgmpe class directly + mgmpe_cls = CB14BasinTerm(gmpe_name='AtkinsonMacias2009') + + # Make GMM with basin term using ModifiableGMPE and kwargs + mgmpe_val = valid.modified_gsim(am09, cb14_basin_term={}) + # Make the ctx + imts = ['PGA', 'SA(1.0)', 'SA(2.0)'] + cmaker = simple_cmaker([gmpe, mgmpe_cls, mgmpe_val], imts) + ctx = new_ctx(cmaker, 3) + ctx.dip = 60. + ctx.rake = 90. + ctx.z1pt0 = np.array([522.32, 516.98, 522.32]) + ctx.z2pt5 = np.array([6.32, 3.53, 6.32]) + ctx.rrup = np.array([50., 200., 500.]) + ctx.vs30 = 1100. + ctx.vs30measured = 1 + mea, _, _, _ = cmaker.get_mean_stds([ctx]) + + # Check expected is observed (original vs modified with basin term) + aae(mea[0], exp_gmm_origin) + aae(mea[1], exp_with_basin) + aae(mea[2], exp_with_basin) \ No newline at end of file diff --git a/openquake/hazardlib/tests/gsim/mgmpe/m9_basin_term_test.py b/openquake/hazardlib/tests/gsim/mgmpe/m9_basin_term_test.py new file mode 100644 index 000000000000..449c24a4708c --- /dev/null +++ b/openquake/hazardlib/tests/gsim/mgmpe/m9_basin_term_test.py @@ -0,0 +1,105 @@ +# The Hazard Library +# Copyright (C) 2012-2023 GEM Foundation +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Affero General Public License as +# published by the Free Software Foundation, either version 3 of the +# License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with this program. If not, see . + +import numpy as np +import unittest + +from openquake.hazardlib.tests.gsim.mgmpe.dummy import new_ctx +from openquake.hazardlib.contexts import simple_cmaker +from openquake.hazardlib import valid +from openquake.hazardlib.imt import PGA, PGV, SA +from openquake.hazardlib.const import TRT, IMC +from openquake.hazardlib.gsim.mgmpe.m9_basin_term import M9BasinTerm + +ae = np.testing.assert_equal +aae = np.testing.assert_almost_equal + + +class M9BasinTermTestCase(unittest.TestCase): + + def test_instantiation(self): + mgmpe = M9BasinTerm(gmpe_name='KuehnEtAl2020SInter') + + # Check the assigned IMTs + expected = set([PGA, SA, PGV]) + ae(mgmpe.DEFINED_FOR_INTENSITY_MEASURE_TYPES, expected) + # Check the TR + expected = TRT.SUBDUCTION_INTERFACE + ae(mgmpe.DEFINED_FOR_TECTONIC_REGION_TYPE, expected) + # Check the IM component + expected = IMC.RotD50 + ae(mgmpe.DEFINED_FOR_INTENSITY_MEASURE_COMPONENT, expected) + # Check the required distances + expected = {'rrup'} + ae(mgmpe.REQUIRES_DISTANCES, expected) + + + def test_all(self): + """ + Test that the M9 basin term applied to Kuehn et al. (2020) + provides the expected values (using sites with z2pt5 above + and below 6 km threshold and considering SAs with periods + above and below 1.9 s) + """ + k20 = valid.gsim('KuehnEtAl2020SInter') + + # Make original GMM + gmpe = valid.modified_gsim(k20) + + # Make GMM with basin term using the mgmpe class directly + mgmpe_cls = M9BasinTerm(gmpe_name='KuehnEtAl2020SInter') + + # Make GMM with basin term using ModifiableGMPE and kwargs + mgmpe_val = valid.modified_gsim(k20, m9_basin_term={}) + + # Make the ctx + imts = ['PGA', 'SA(1.0)', 'SA(2.0)'] + cmaker = simple_cmaker([gmpe, mgmpe_cls, mgmpe_val], imts) + ctx = new_ctx(cmaker, 3) + ctx.dip = 60. + ctx.rake = 90. + ctx.z1pt0 = np.array([522.32, 516.98, 522.32]) + ctx.z2pt5 = np.array([6.32, 3.53, 6.32]) + ctx.rrup = np.array([50., 200., 500.]) + ctx.vs30 = 1100. + ctx.vs30measured = 1 + mea, _, _, _ = cmaker.get_mean_stds([ctx]) + + # For SA with period less than 2 all means should be + # the same (amp term only applied for long periods) + ori_mea_pga = mea[0][0] + cls_mea_pga = mea[1][0] + val_mea_pga = mea[2][0] + ae(ori_mea_pga, cls_mea_pga, val_mea_pga) + ori_mea_sa1pt0 = mea[0][1] + cls_mea_sa1pt0 = mea[1][1] + val_mea_sa1pt0 = mea[2][1] + ae(ori_mea_sa1pt0, cls_mea_sa1pt0, val_mea_sa1pt0) + + # Check means using m9 term from mgmpe cls and + # valid.modifiable_gmpe are the same + ae(mea[1], mea[2]) + + # For SA(2.0) basin amplification should be added + ori_sa2pt0 = mea[0][2] + cls_sa2pt0 = mea[1][2] + val_sa2pt0 = mea[2][2] + exp_diff = np.array([np.log(2.0), 0., np.log(2.0)]) # Site 1 intensities + for idx_v, v in enumerate(exp_diff): # note be changed + diff_cls = cls_sa2pt0[idx_v] - ori_sa2pt0[idx_v] + diff_val = val_sa2pt0[idx_v] - ori_sa2pt0[idx_v] + aae(diff_cls, exp_diff[idx_v]) + aae(diff_val, exp_diff[idx_v]) \ No newline at end of file