New GSIM Idini2017

debian/changelog

@@ -1,3 +1,7 @@
+  [Pablo Heresi]
+  * Implemented Idini et al (2017) GSIM.
+  * Included 'soiltype' into site.py and adds it in site_test.py and oqvalidation.py
+
   [Claudia Mascandola]
   * Fix to LanzanoEtAl2016 in presence of a "bas" term in the site model
 

doc/sphinx/openquake.hazardlib.gsim.rst

@@ -650,6 +650,14 @@ hong_goda_2007
     :undoc-members:
     :show-inheritance:
 
+idini_2017
+---------------------------------------------
+
+.. automodule:: openquake.hazardlib.gsim.idini_2017
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
 idriss_2014
 -------------------------------------------
 

openquake/calculators/tests/classical_test.py

@@ -37,7 +37,7 @@
     case_34, case_35, case_36, case_37, case_38, case_39, case_40, case_41,
     case_42, case_43, case_44, case_45, case_46, case_47, case_48, case_49,
     case_50, case_51, case_52, case_53, case_54, case_55, case_56, case_57,
-    case_58, case_59, case_60, case_61, case_62, case_64)
+    case_58, case_59, case_60, case_61, case_62, case_63, case_64)
 
 aac = numpy.testing.assert_allclose
 
@@ -868,6 +868,13 @@ def test_case_62(self):
         # multisurface with kite faults
         self.run_calc(case_62.__file__, 'job.ini')
         [f] = export(('hcurves/mean', 'csv'), self.calc.datastore)
+        self.assertEqualFiles('expected/hcurve-mean.csv', f)    
+
+    def test_case_63(self):
+        # test soiltype
+        self.run_calc(case_63.__file__, 'job.ini')
+        [f] = export(('hcurves/mean', 'csv'), self.calc.datastore)
+        self.assertEqualFiles('expected/hcurve-mean.csv', f)    
         self.assertEqualFiles('expected/hcurve-mean.csv', f)
 
     def test_case_64(self):

openquake/commonlib/oqvalidation.py

@@ -474,6 +474,13 @@
   Example: *reference_siteclass = "A"*.
   Default: "D"
 
+reference_soiltype:
+  Used when there is no site model to specify a global soil type.
+  The soiltype is an integer between 1 (rock) and 6, and it is used in some 
+  GMPEs, like Idini (2017).
+  Example: *reference_soiltype = 3*.
+  Default: 1
+
 reference_vs30_type:
   Used when there is no site model to specify a global vs30 type.
   The choices are "inferred" or "measured"
@@ -706,7 +713,8 @@ class OqParam(valid.ParamSet):
         z1pt0='reference_depth_to_1pt0km_per_sec',
         z2pt5='reference_depth_to_2pt5km_per_sec',
         siteclass='reference_siteclass',
-        backarc='reference_backarc')
+        backarc='reference_backarc',
+        soiltype='reference_soiltype')
     aggregate_by = valid.Param(valid.namelist, [])
     amplification_method = valid.Param(
         valid.Choice('convolution', 'kernel'), None)
@@ -799,6 +807,7 @@ class OqParam(valid.ParamSet):
         valid.positivefloat, numpy.nan)
     reference_siteclass = valid.Param(valid.Choice('A', 'B', 'C', 'D'), 'D')
     reference_backarc = valid.Param(valid.boolean, False)
+    reference_soiltype = valid.Param(valid.positiveints, 1)
     region = valid.Param(valid.wkt_polygon, None)
     region_grid_spacing = valid.Param(valid.positivefloat, None)
     risk_imtls = valid.Param(valid.intensity_measure_types_and_levels, {})

openquake/hazardlib/gsim/idini_2017.py

@@ -0,0 +1,210 @@
+# -*- coding: utf-8 -*-
+# vim: tabstop=4 shiftwidth=4 softtabstop=4
+#
+# Copyright (C) 2015-2021 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:`IdiniEtAl2017SInter`
+               :class:`IdiniEtAl2017SSlab`
+"""
+import numpy as np
+
+from openquake.hazardlib import const
+from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
+from openquake.hazardlib.imt import PGA, SA
+
+class IdiniEtAl2017SInter(GMPE):
+    """
+    Implements the GMPE developed by Idini et al. (2017) for subduction 
+    interface earthquakes, publised as:
+
+    Idini, B., F. Rojas, S. Ruiz, and C. Pastén. 2017. “Ground motion 
+    prediction equations for the Chilean subduction zone.” Bull. Earthq. Eng. 
+    15(5): 1853–1880.
+    """
+
+    #: Supported tectonic region type is subduction interface
+    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.SUBDUCTION_INTERFACE
+
+    #: Supported intensity measure types are spectral acceleration,
+    #: and peak ground acceleration
+    DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, SA}
+
+    #: Supported intensity measure component is the geometric mean component
+    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.AVERAGE_HORIZONTAL
+
+    #: Supported standard deviation types are inter-event, intra-event
+    #: and total
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {
+        const.StdDev.TOTAL, const.StdDev.INTER_EVENT, const.StdDev.INTRA_EVENT}
+
+    #: Site amplification is dependent on the Site Class and Vs30
+    REQUIRES_SITES_PARAMETERS = {'soiltype','vs30'}
+
+    #: Required rupture parameters are only magnitude for the interface model
+    REQUIRES_RUPTURE_PARAMETERS = {'mag'}
+
+    #: Required distance measure is closest distance to rupture, for
+    #: interface events M>=7.7, and hypocentral distance for interface events 
+    #: M<7.7
+    REQUIRES_DISTANCES = {'rrup','rhypo'}
+
+    def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
+        """
+        See :meth:`superclass method
+        <.base.GroundShakingIntensityModel.get_mean_and_stddevs>`
+        for spec of input and result values.
+        """
+        # extracting dictionary of coefficients specific to required
+        # intensity measure type.
+        C = self.COEFFS[imt]
+
+        mean = (self._get_magnitude_term(C, rup) +
+                self._get_distance_term(C, rup, dists) +
+                self._get_site_term(C, sites))
+        # Convert from log10 to ln
+        mean = np.log(10 ** mean)
+
+        stddevs = self._get_stddevs(C, len(sites.vs30), stddev_types)
+
+        return mean, stddevs
+
+    def _get_magnitude_term(self, C, rup):
+        """
+        Returns the magnitude scaling term defined in Equation (3)
+        """
+        mag = rup.mag
+        return C['c1'] + C['c2'] * mag + C['c9'] * mag ** 2.0
+
+
+    def _get_distance_term(self, C, rup, dists):
+        """
+        Returns the magnitude dependent distance scaling term defined in 
+        Equation (5)
+        """
+        mag = rup.mag
+        if mag < 7.7:
+            R = dists.rhypo
+        else:
+            R = dists.rrup
+
+        g = C['c3'] + self.CONSTS['c4'] * (mag - self.CONSTS['Mr'])
+        Ro = self.CONSTS['c6'] * 10 ** (self.CONSTS['c7'] * (mag -\
+             self.CONSTS['Mr']))
+        return g * np.log10(R + Ro) + C['c5'] * R
+
+    def _get_site_term(self, C, sites):
+        """
+        Returns the site scaling term defined in Equation (18)
+        """
+        soiltype = sites.soiltype
+        vs30 = sites.vs30
+
+        # sT* depends on the soil type.
+        # If soiltype > 6, use soiltype = 1 (rock)
+        sT = [C['s%i' % st] if st > 1 and st <= 6 else 0 for st in soiltype]
+
+        return sT * np.log10(vs30 / self.CONSTS['Vref'])
+
+    def _get_stddevs(self, C, n_sites, stddev_types):
+        """
+        Returns the standard deviations
+        """
+        # sigma_e and sigma_r are in log10 base, so we need to transform them
+        tau = np.log(10 ** C['sigma_e']) + np.zeros(n_sites)
+        phi = np.log(10 ** C['sigma_r']) + np.zeros(n_sites)
+        stddevs = []
+        for stddev_type in stddev_types:
+            assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
+            if stddev_type == const.StdDev.TOTAL:
+                sigma = np.sqrt(tau ** 2.0 + phi ** 2.0)
+                stddevs.append(sigma)
+            elif stddev_type == const.StdDev.INTRA_EVENT:
+                stddevs.append(phi)
+            elif stddev_type == const.StdDev.INTER_EVENT:
+                stddevs.append(tau)
+        return stddevs
+
+    COEFFS = CoeffsTable(sa_damping=5, table="""\
+    imt	          c1	     c2	       c9	      c8	    dc1	      dc2	sigma_e	sigma_t	       c3	       c5	      dc3	sigma_r	     s2	     s3	     s4	     s5	     s6                         
+    pga	    -2.85480	0.77410	 -0.039580	0.005860	2.56990	 -0.47610	0.1720	0.2890	-0.975580	-0.001740	-0.527450	0.2320	-0.5840	-0.3220	-0.1090	-0.0950	-0.2120
+    0.01	-2.84240	0.80520	 -0.041350	0.005840	2.73700	 -0.51910	0.1730	0.2880	-1.029930	-0.001750	-0.504660	0.2310	-0.5230	-0.2620	-0.1000	-0.0920	-0.1930
+    0.02	-2.83370	0.83830	 -0.043250	0.005830	2.90870	 -0.56400	0.1760	0.2920	-1.085670	-0.001760	-0.480430	0.2330	-0.4590	-0.2080	-0.0920	-0.0890	-0.1770
+    0.03	-2.82350	0.88380	 -0.045950	0.005860	3.07350	 -0.62270	0.1780	0.2950	-1.159510	-0.001760	-0.424900	0.2350	-0.3900	-0.1600	-0.0850	-0.0880	-0.1640
+    0.05	-2.73580	0.95390	 -0.050330	0.006210	3.21470	 -0.70790	0.1900	0.3070	-1.286400	-0.001780	-0.312390	0.2410	-0.3060	-0.0880	-0.0750	-0.0900	-0.1460
+    0.07	-2.60040	0.98080	 -0.052250	0.006030	3.08510	 -0.74250	0.2130	0.3290	-1.346440	-0.001810	-0.179950	0.2510	-0.3510	-0.0560	-0.0690	-0.0960	-0.1410
+    0.10	-2.48910	0.95440	 -0.050600	0.005710	2.80910	 -0.70550	0.1950	0.3210	-1.323530	-0.001820	-0.132080	0.2550	-0.5240	-0.0870	-0.0700	-0.1130	-0.1560
+    0.15	-2.65050	0.92320	 -0.048790	0.005600	2.62600	 -0.62700	0.1600	0.3020	-1.176870	-0.001830	-0.264510	0.2550	-0.6910	-0.3360	-0.0950	-0.1660	-0.2450
+    0.20	-3.00960	0.94260	 -0.050340	0.005730	2.60630	 -0.59760	0.1570	0.3100	-1.045080	-0.001820	-0.391050	0.2680	-0.6710	-0.5470	-0.1270	-0.2090	-0.3590
+    0.25	-3.33210	0.95780	 -0.051430	0.005070	2.36540	 -0.58200	0.1420	0.2990	-0.943630	-0.001780	-0.343480	0.2640	-0.5840	-0.6740	-0.1780	-0.2350	-0.4440
+    0.30	-3.54220	0.94410	 -0.050520	0.004280	2.20170	 -0.54120	0.1410	0.2960	-0.848140	-0.001730	-0.366950	0.2600	-0.5060	-0.7300	-0.2580	-0.2340	-0.4910
+    0.40	-3.39850	0.77730	 -0.038850	0.003080	1.63670	 -0.34480	0.1570	0.3060	-0.692780	-0.001660	-0.463010	0.2630	-0.3860	-0.7180	-0.4230	-0.1640	-0.5350
+    0.50	-2.80410	0.50690	 -0.019730	0.002570	0.76210	 -0.06170	0.1520	0.3020	-0.578990	-0.001610	-0.540980	0.2610	-0.3000	-0.6350	-0.5370	-0.1100	-0.5570
+    0.75	-4.45880	0.86910	 -0.041790	0.001350	2.10030	 -0.43490	0.1460	0.2910	-0.568870	-0.001580	-0.462660	0.2520	-0.2760	-0.3950	-0.5750	-0.3580	-0.5990
+    1.00	-5.33910	1.01670	 -0.049990	0.000450	2.56100	 -0.56780	0.1530	0.2900	-0.532820	-0.001540	-0.423140	0.2470	-0.2750	-0.2540	-0.4620	-0.6700	-0.5840
+    1.50	-6.12040	1.10050	 -0.054260	0.000680	2.89230	 -0.58980	0.1520	0.2890	-0.462630	-0.001450	-0.585190	0.2460	-0.2490	-0.2380	-0.3000	-0.8010	-0.5220
+    2.00	-7.03340	1.25010	 -0.063560	0.000510	3.39410	 -0.70090	0.1570	0.2910	-0.405940	-0.001390	-0.659990	0.2450	-0.2180	-0.2310	-0.2200	-0.7460	-0.4790
+    3.00	-8.25070	1.46520	 -0.077970	0.000660	4.00330	 -0.84650	0.1550	0.2790	-0.339570	-0.001370	-0.790040	0.2310	-0.1800	-0.2190	-0.2100	-0.6280	-0.4610
+    4.00	-8.74330	1.48270	 -0.078630	0.000630	3.93370	 -0.81340	0.1600	0.2790	-0.264790	-0.001370	-0.865450	0.2280	-0.1710	-0.2180	-0.2120	-0.5310	-0.4480
+    5.00	-8.99270	1.46300	 -0.076380	0.000670	3.75760	 -0.76420	0.1670	0.2860	-0.223330	-0.001370	-0.887350	0.2320	-0.1680	-0.2180	-0.2030	-0.4380	-0.4390
+    7.50	-9.82450	1.63830	 -0.086200	0.001080	4.39480	 -0.93130	0.1640	0.2830	-0.303460	-0.001310	-0.912590	0.2310	-0.1680	-0.2180	-0.1530	-0.2560	-0.4350
+    10.00	-9.86710	1.58770	 -0.081680	0.000140	4.38750	 -0.88920	0.1760	0.2700	-0.337710	-0.001170	-0.963630	0.2040	-0.1680	-0.2180	-0.1250	-0.2310	-0.4350    
+    """)
+
+
+    CONSTS = {'Mr': 5,
+              'Vref': 1530,
+              'c4': 0.1,
+              'c6': 5,
+              'c7': 0.35,
+              'h0': 50}
+
+
+class IdiniEtAl2017SSlab(IdiniEtAl2017SInter):
+    """
+    Implements the GMPE developed by Idini et al. (2017) for subduction 
+    inslab (intraslab) earthquakes, publised as:
+
+    Idini, B., F. Rojas, S. Ruiz, and C. Pastén. 2017. “Ground motion 
+    prediction equations for the Chilean subduction zone.” Bull. Earthq. Eng. 
+    15(5): 1853–1880.
+    """
+
+    #: Supported tectonic region type is subduction in-slab
+    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.SUBDUCTION_INTRASLAB
+
+    #: Required rupture parameters for the in-slab model are magnitude and top
+    # of rupture depth
+    REQUIRES_RUPTURE_PARAMETERS = {'mag', 'hypo_depth'}
+
+    def _get_magnitude_term(self, C, rup):
+        """
+        Returns the magnitude scaling term defined in Equation (3)
+        """
+        mag = rup.mag
+        H = rup.hypo_depth
+        return C['c1'] + C['c2'] * mag + C['c8'] * (H - self.CONSTS['h0']) + \
+            C['dc1'] + C['dc2'] * mag
+
+    def _get_distance_term(self, C, rup, dists):
+        """
+        Returns the magnitude dependent distance scaling term defined in 
+        Equation (5)
+        """
+        mag = rup.mag
+        R = dists.rhypo
+        g = C['c3'] + self.CONSTS['c4'] * (mag - self.CONSTS['Mr']) + C['dc3']
+        return g * np.log10(R) + C['c5'] * R

openquake/hazardlib/site.py

@@ -125,6 +125,7 @@ def _extract(array_or_float, indices):
     'z1pt4': numpy.float64,
     'backarc': numpy.bool,
     'xvf': numpy.float64,
+    'soiltype': numpy.uint32,
     'bas': numpy.bool,
 
     # Parameters for site amplification
@@ -252,6 +253,7 @@ def from_points(cls, lons, lats, depths=None, sitemodel=None,
             self._set('z2pt5', sitemodel.reference_depth_to_2pt5km_per_sec)
             self._set('siteclass', sitemodel.reference_siteclass)
             self._set('backarc', sitemodel.reference_backarc)
+            self._set('soiltype', sitemodel.reference_soiltype)
         else:
             for name in sitemodel.dtype.names:
                 if name not in ('lon', 'lat'):