Implementation of GMM of Ambraseys 1996 and Sapetta_Pugliese 1996

doc/api-reference/openquake.hazardlib.gsim.rst

@@ -126,6 +126,14 @@ allen_2022
     :undoc-members:
     :show-inheritance:
 
+ambraseys_1996
+------------------------------------------
+
+.. automodule:: openquake.hazardlib.gsim.ambraseys_1996
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
 ambraseys_2005
 ------------------------------------------
 
@@ -1034,6 +1042,15 @@ rietbrock_edwards_2019
     :undoc-members:
     :show-inheritance:
 
+sabetta_pugliese_1996
+-------------------------------------------
+
+.. automodule:: openquake.hazardlib.gsim.sabetta_pugliese_1996
+
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
 sandikkaya_akkar_2017
 -------------------------------------------
 

openquake/hazardlib/gsim/ambraseys_1996.py

@@ -0,0 +1,179 @@
+# -*- coding: utf-8 -*-
+# vim: tabstop=4 shiftwidth=4 softtabstop=4
+#
+# Copyright (C) 2014-2023 GEM Foundation
+#
+# OpenQuake 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.
+#
+# OpenQuake 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 OpenQuake. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+Module exports :class:`AmbraseysEtAl1996`.
+"""
+import numpy as np
+from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
+from openquake.hazardlib.gsim import utils
+from openquake.hazardlib import const
+from openquake.hazardlib.imt import PGA, SA
+
+
+
+def _compute_magnitude(ctx, C):
+    return C['c1'] + (C['c2'] * ctx.mag)
+
+def _compute_distance(ctx, C):
+    #convertion from repi to rjb if M>6 (Montaldo et al 2005, eq. 3)
+
+    R = np.where(ctx.mag >= 6, -3.5525 + (0.8845 * ctx.repi), ctx.repi)
+    R = np.where(R < 0, 0, R)
+    return (C['c4']* np.log10(np.sqrt(R ** 2 + C['h'] ** 2)))
+
+def _site_amplification(ctx, C):
+    #site condition
+    ssa = np.zeros(len(ctx.vs30))
+    sss = np.zeros(len(ctx.vs30))
+
+    idx = (ctx.vs30 > 360.0) & (ctx.vs30 <= 750.0)
+    ssa[idx] = 1.0
+    idx = (ctx.vs30 <= 360.0)
+    sss[idx] = 1.0
+
+    return C['ca'] * ssa + C['cs'] * sss
+
+def _get_mechanism(ctx, C):
+    #style of faulting correction as in Montale et al 2005 (Table 2)
+    SS, NS, RS = utils.get_fault_type_dummy_variables(ctx)
+
+    idx = (ctx.mag < 6)
+    SS[idx] = 0
+    NS[idx] = 0
+    RS[idx] = 0
+
+    Fr = np.log10(1.13)
+    Fn = np.log10(0.88)
+    Fss = np.log10(0.93)
+
+    return Fr*RS + Fn*NS + Fss*SS
+
+
+class AmbraseysEtAl1996(GMPE):
+    """
+    Implements GMPE developed by Ambraseys, N. N., Simpson, K. U., & 
+    Bommer, J. J. (1996). Prediction of horizontal response spectra in Europe.
+    Earthquake Engineering & Structural Dynamics, 25(4), 371-400.
+    SA are given up to 2 s.
+    The regressions are developed considering the largest horizontal component
+
+    GMPE implemented as adopted in MPS04 (Montaldo et al 2005)
+    """
+    #: Supported tectonic region type is 'active shallow crust'
+
+    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST
+
+    #: Set of :mod:`intensity measure types <openquake.hazardlib.imt>`
+    #: this GSIM can calculate. A set should contain classes from module
+    #: :mod:`openquake.hazardlib.imt`.
+    DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, SA}
+
+    #: Supported intensity measure component is the largest
+    #: horizontal component
+    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.GREATER_OF_TWO_HORIZONTAL
+
+    #: Supported standard deviation type is total
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+
+    #: Required site parameter is only Vs30
+    REQUIRES_SITES_PARAMETERS = {'vs30'}
+
+    #: Required rupture parameters are magnitude and rake.
+    REQUIRES_RUPTURE_PARAMETERS = {'mag', 'rake'}
+
+    #: Required distance measure is Repi
+    REQUIRES_DISTANCES = {'repi'}
+
+
+    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.
+        """
+
+        for m, imt in enumerate(imts):
+            C = self.COEFFS[imt]
+            imean = _compute_magnitude(ctx, C) + \
+            _compute_distance(ctx, C) + \
+            _site_amplification(ctx, C) +\
+            _get_mechanism(ctx, C)
+
+            # convert from base 10 to base e (already in g)
+            mean[m] = np.log(10.0 ** imean)
+
+            # Return stddevs in terms of natural log scaling
+            sig[m] = np.log(10.0 ** C['sigma'])
+
+
+
+    #: Coefficients for PGA and SA
+
+    COEFFS = CoeffsTable(sa_damping=5, table="""
+         IMT       c1      c2     h     c4     ca      cs  sigma
+         pga    -1.48   0.266   3.5 -0.922  0.117   0.124   0.25
+        0.10    -0.84   0.219   4.5 -0.954  0.078   0.027   0.27
+        0.11    -0.86   0.221   4.5 -0.945  0.098   0.036   0.27
+        0.12    -0.87   0.231   4.7 -0.960  0.111   0.052   0.27
+        0.13    -0.87   0.238   5.3 -0.981  0.131   0.068   0.27
+        0.14    -0.94   0.244   4.9 -0.955  0.136   0.077   0.27
+        0.15    -0.98   0.247   4.7 -0.938  0.143   0.085   0.27
+        0.16    -1.05   0.252   4.4 -0.907  0.152   0.101   0.27
+        0.17    -1.08   0.258   4.3 -0.896  0.140   0.102   0.27
+        0.18    -1.13   0.268   4.0 -0.901  0.129   0.107   0.27
+        0.19    -1.19   0.278   3.9 -0.907  0.133   0.130   0.28
+        0.20    -1.21   0.284   4.2 -0.922  0.135   0.142   0.27
+        0.22    -1.28   0.295   4.1 -0.911  0.120   0.143   0.28
+        0.24    -1.37   0.308   3.9 -0.916  0.124   0.155   0.28
+        0.26    -1.40   0.318   4.3 -0.942  0.134   0.163   0.28
+        0.28    -1.46   0.326   4.4 -0.946  0.134   0.158   0.29
+        0.30    -1.55   0.338   4.2 -0.933  0.133   0.148   0.30
+        0.32    -1.63   0.349   4.2 -0.932  0.125   0.161   0.31
+        0.34    -1.65   0.351   4.4 -0.939  0.118   0.163   0.31
+        0.36    -1.69   0.354   4.5 -0.936  0.124   0.160   0.31
+        0.38    -1.82   0.364   3.9 -0.900  0.132   0.164   0.31
+        0.40    -1.94   0.377   3.6 -0.888  0.139   0.172   0.31
+        0.42    -1.99   0.384   3.7 -0.897  0.147   0.180   0.32
+        0.44    -2.05   0.393   3.9 -0.908  0.153   0.187   0.32
+        0.46    -2.11   0.401   3.7 -0.911  0.149   0.191   0.32
+        0.48    -2.17   0.410   3.5 -0.920  0.150   0.197   0.32
+        0.50    -2.25   0.420   3.3 -0.913  0.147   0.201   0.32
+        0.55    -2.38   0.434   3.1 -0.911  0.134   0.203   0.32
+        0.60    -2.49   0.438   2.5 -0.881  0.124   0.212   0.32
+        0.65    -2.58   0.451   2.8 -0.901  0.122   0.215   0.32
+        0.70    -2.67   0.463   3.1 -0.914  0.116   0.214   0.33
+        0.75    -2.75   0.477   3.5 -0.942  0.113   0.212   0.32
+        0.80    -2.86   0.485   3.7 -0.925  0.127   0.218   0.32
+        0.85    -2.93   0.492   3.9 -0.920  0.124   0.218   0.32
+        0.90    -3.03   0.502   4.0 -0.920  0.124   0.225   0.32
+        0.95    -3.10   0.503   4.0 -0.892  0.121   0.217   0.32
+        1.00    -3.17   0.508   4.3 -0.885  0.128   0.219   0.32
+        1.10    -3.30   0.513   4.0 -0.857  0.123   0.206   0.32
+        1.20    -3.38   0.513   3.6 -0.851  0.128   0.214   0.31
+        1.30    -3.43   0.514   3.6 -0.848  0.115   0.200   0.31
+        1.40    -3.52   0.522   3.4 -0.839  0.109   0.197   0.31
+        1.50    -3.61   0.524   3.0 -0.817  0.109   0.204   0.31
+        1.60    -3.68   0.520   2.5 -0.781  0.108   0.206   0.31
+        1.70    -3.74   0.517   2.5 -0.759  0.105   0.206   0.31
+        1.80    -3.79   0.514   2.4 -0.730  0.104   0.204   0.32
+        1.90    -3.80   0.508   2.8 -0.724  0.103   0.194   0.32
+        2.00    -3.79   0.503   3.2 -0.728  0.101   0.182   0.32
+    """)
+
+

openquake/hazardlib/gsim/sabetta_pugliese_1996.py

@@ -0,0 +1,150 @@
+# -*- coding: utf-8 -*-
+# vim: tabstop=4 shiftwidth=4 softtabstop=4
+#
+# Copyright (C) 2014-2023 GEM Foundation
+#
+# OpenQuake 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.
+#
+# OpenQuake 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 OpenQuake. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+Module exports :class:`AmbraseysEtAl1996`.
+"""
+import numpy as np
+from scipy.constants import g, pi
+
+from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
+from openquake.hazardlib.gsim import utils
+from openquake.hazardlib import const
+from openquake.hazardlib.imt import PGA, PGV, SA
+
+
+
+
+def _compute_magnitude(ctx, C):
+    return C['a'] + (C['b'] * ctx.mag)
+
+def _compute_distance(ctx, C):
+    return (C['c']* np.log10(np.sqrt(ctx.repi ** 2 + C['h'] ** 2)))
+
+def _site_amplification(ctx, C):
+    #site condition
+    sshallow = np.zeros(len(ctx.vs30))
+    sdeep = np.zeros(len(ctx.vs30))
+    idx = (ctx.vs30 >= 400.0) & (ctx.vs30 < 800.0)
+    sshallow[idx] = 1.0
+    idx = (ctx.vs30 < 400.0)
+    sdeep[idx] = 1.0
+    return C['e1'] * sshallow + C['e2'] * sdeep
+
+def _get_mechanism(ctx, C):
+    #style of faulting correction as in Montale et al 2005 (Table 2)
+    SS, NS, RS = utils.get_fault_type_dummy_variables(ctx)
+    idx = (ctx.mag < 6)
+    SS[idx] = 0
+    NS[idx] = 0
+    RS[idx] = 0
+    Fr = np.log10(1.15)
+    Fn = np.log10(0.89)
+    Fss = np.log10(0.94)
+    return Fr*RS + Fn*NS + Fss*SS
+
+
+
+
+class SabettaPugliese1996(GMPE):
+    """
+    Implements GMPE developed by Sabetta, F., & Pugliese, A. (1996). 
+    Estimation of response spectra and simulation of nonstationary earthquake 
+    ground motions. Bulletin of the Seismological Society of America, 86(2), 337-352.
+    SA are given up to 4 s.
+    The regressions are developed considering the largest horizontal component
+
+    GMPE implemented as adopted in MPS04 (Montaldo et al 2005)
+    """
+    #: Supported tectonic region type is 'active shallow crust'
+
+    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST
+
+    #: Set of :mod:`intensity measure types <openquake.hazardlib.imt>`
+    #: this GSIM can calculate. A set should contain classes from module
+    #: :mod:`openquake.hazardlib.imt`.
+    DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, PGV, SA}
+
+    #: Supported intensity measure component is the largest
+    #: horizontal component
+    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.GREATER_OF_TWO_HORIZONTAL
+
+    #: Supported standard deviation type is total
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+
+    #: Required site parameter is only Vs30
+    REQUIRES_SITES_PARAMETERS = {'vs30'}
+
+    #: Required rupture parameter is magnitude.
+    REQUIRES_RUPTURE_PARAMETERS = {'mag', 'rake'}
+
+    #: Required distance measure is Repi (eq. 1).
+    REQUIRES_DISTANCES = {'repi'}
+
+
+    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.
+        """
+
+        for m, imt in enumerate(imts):
+            C = self.COEFFS[imt]
+            imean = _compute_magnitude(ctx, C) + \
+            _compute_distance(ctx, C) + \
+            _site_amplification(ctx, C) + \
+            _get_mechanism(ctx, C)
+
+            # Convert units to g (only pga is already in g)
+            # Sa in PSV (cm/s)
+            if imt.string.startswith(('SA')):
+                mean[m] = np.log((10.0 ** (imean - 2.0)) * (2*pi/imt.period) / g)
+            elif imt == PGV():
+                #PGV - convert in g unit and from base 10 to base e
+                mean[m] = np.log((10.0 ** (imean - 2.0)) / g)
+            else:
+                # PGA - from base 10 to base e
+                mean[m] = np.log(10.0 ** imean)
+
+            # Return stddevs in terms of natural log scaling
+            sig[m] = np.log(10.0 ** C['sigma'])
+
+
+    #: Coefficients for PGA and SA
+
+    COEFFS = CoeffsTable(sa_damping=5, table="""
+        IMT      a      b       c       e1      e2      h       sigma
+        4.0000  -2.5000 0.7250  -1.0000 0.0000  0.1000  2.6000  0.3190
+        3.0303  -2.2500 0.7150  -1.0000 0.0000  0.1080  3.0000  0.3190
+        2.0000  -1.9000 0.6870  -1.0000 0.0000  0.1500  3.6000  0.3190
+        1.4925  -1.6470 0.6600  -1.0000 0.0100  0.1750  4.0000  0.3150
+        1.0000  -1.2800 0.6120  -1.0000 0.0500  0.2080  4.4000  0.3080
+        0.7519  -1.0000 0.5700  -1.0000 0.1200  0.1900  4.7000  0.3030
+        0.5000  -0.5950 0.5000  -1.0000 0.2300  0.1240  5.0000  0.2900
+        0.4000  -0.2810 0.4450  -1.0000 0.2220  0.0780  5.2000  0.2800
+        0.3003  0.1000  0.3770  -1.0000 0.1850  0.0200  5.4000  0.2600
+        0.2000  0.2960  0.3230  -1.0000 0.1610  0.0000  5.7000  0.2340
+        0.1499  0.2220  0.3100  -1.0000 0.1610  0.0000  5.9000  0.2200
+        0.1000  -0.0190 0.3040  -1.0000 0.1610  0.0000  6.2000  0.2080
+        0.0667  -0.3120 0.3040  -1.0000 0.1610  0.0000  6.3000  0.2000
+        0.0400  -0.8170 0.3300  -1.0000 0.1610  0.0000  4.7000  0.1950
+        pga     -1.845  0.363   -1.0000 0.1950  0.0000  5.0000  0.1900
+        pgv     -0.828  0.489   -1.0000 0.1160  0.1160  3.9000  0.2490
+    """)
+

openquake/hazardlib/tests/gsim/ambraseys_1996_test.py

@@ -0,0 +1,37 @@
+# -*- coding: utf-8 -*-
+# vim: tabstop=4 shiftwidth=4 softtabstop=4
+#
+# Copyright (C) 2013-2023 GEM Foundation
+#
+# OpenQuake 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.
+#
+# OpenQuake 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 OpenQuake. If not, see <http://www.gnu.org/licenses/>.
+
+from openquake.hazardlib.gsim.ambraseys_1996 import AmbraseysEtAl1996
+from openquake.hazardlib.tests.gsim.utils import BaseGSIMTestCase
+
+# Test data provided by INGV
+
+
+class AmbraseysEtAl1996TestCase(BaseGSIMTestCase):
+    GSIM_CLASS = AmbraseysEtAl1996
+
+    # Tables provided by INGV (Giovanni Lanzano- Vera D'Amico)
+    MEAN_FILE = "AMB96/AMB96_MPS04_MEAN.csv"
+    STD_FILE = "AMB96/AMB96_MPS04_STD_TOTAL.csv"
+
+    def test_all(self):
+        self.check(self.MEAN_FILE, self.STD_FILE,
+                   max_discrep_percentage=0.1)
+
+
+