# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2012-2017 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:`ToroEtAl2002`, class:`ToroEtAl2002SHARE`.
"""
from __future__ import division
import numpy as np
# standard acceleration of gravity in m/s**2
from scipy.constants import g
from openquake.hazardlib.gsim.campbell_2003 import _compute_faulting_style_term
from openquake.hazardlib.gsim.base import CoeffsTable, GMPE
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, SA
[docs]class ToroEtAl2002(GMPE):
"""
Implements GMPE developed by G. R. Toro, N. A. Abrahamson, J. F. Sneider
and published in "Model of Strong Ground Motions from Earthquakes in
Central and Eastern North America: Best Estimates and Uncertainties"
(Seismological Research Letters, Volume 68, Number 1, 1997) and
"Modification of the Toro et al. 1997 Attenuation Equations for Large
Magnitudes and Short Distances" (available at:
http://www.riskeng.com/downloads/attenuation_equations)
The class implements equations for Midcontinent, based on moment magnitude.
SA at 3 and 4 s (not supported by the original equations) have been added
in the context of the SHARE project and they are obtained from SA at 2 s
scaled by specific factors for 3 and 4 s.
"""
#: Supported tectonic region type is stable continental crust,
#: given that the equations have been derived for central and eastern
#: north America
DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.STABLE_CONTINENTAL
#: Supported intensity measure types are spectral acceleration,
#: and peak ground acceleration, see table 2 page 47.
DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([
PGA,
SA
])
#: Supported intensity measure component is the geometric mean of
#two : horizontal components
#:attr:`~openquake.hazardlib.const.IMC.AVERAGE_HORIZONTAL`,
DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.AVERAGE_HORIZONTAL
#: Supported standard deviation type is only total.
DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([
const.StdDev.TOTAL
])
#: No site parameters required
REQUIRES_SITES_PARAMETERS = set()
#: Required rupture parameter is only magnitude.
REQUIRES_RUPTURE_PARAMETERS = set(('mag', ))
#: Required distance measure is rjb, see equation 4, page 46.
REQUIRES_DISTANCES = set(('rjb', ))
[docs] 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.
"""
assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
for stddev_type in stddev_types)
C = self.COEFFS[imt]
mean = self._compute_mean(C, rup.mag, dists.rjb)
stddevs = self._compute_stddevs(C, rup.mag, dists.rjb, imt,
stddev_types)
# apply decay factor for 3 and 4 seconds (not originally supported
# by the equations)
if isinstance(imt, SA):
if imt.period == 3.0:
mean /= 0.612
if imt.period == 4.0:
mean /= 0.559
return mean, stddevs
def _compute_term1(self, C, mag):
"""
Compute magnitude dependent terms (2nd and 3rd) in equation 3
page 46.
"""
mag_diff = mag - 6
return C['c2'] * mag_diff + C['c3'] * mag_diff ** 2
def _compute_term2(self, C, mag, rjb):
"""
Compute distance dependent terms (4th, 5th and 6th) in equation 3
page 46. The factor 'RM' is computed according to the 2002 model
(equation 4-3).
"""
RM = np.sqrt(rjb ** 2 + (C['c7'] ** 2) *
np.exp(-1.25 + 0.227 * mag) ** 2)
return (-C['c4'] * np.log(RM) -
(C['c5'] - C['c4']) *
np.maximum(np.log(RM / 100), 0) - C['c6'] * RM)
def _compute_mean(self, C, mag, rjb):
"""
Compute mean value according to equation 3, page 46.
"""
mean = (C['c1'] +
self._compute_term1(C, mag) +
self._compute_term2(C, mag, rjb))
return mean
def _compute_stddevs(self, C, mag, rjb, imt, stddev_types):
"""
Compute total standard deviation, equations 5 and 6, page 48.
"""
# aleatory uncertainty
sigma_ale_m = np.interp(mag, [5.0, 5.5, 8.0],
[C['m50'], C['m55'], C['m80']])
sigma_ale_rjb = np.interp(rjb, [5.0, 20.0], [C['r5'], C['r20']])
sigma_ale = np.sqrt(sigma_ale_m ** 2 + sigma_ale_rjb ** 2)
# epistemic uncertainty
if isinstance(imt, PGA) or (isinstance(imt, SA) and imt.period < 1):
sigma_epi = 0.36 + 0.07 * (mag - 6)
else:
sigma_epi = 0.34 + 0.06 * (mag - 6)
sigma_total = np.sqrt(sigma_ale ** 2 + sigma_epi ** 2)
stddevs = []
for _ in stddev_types:
stddevs.append(sigma_total)
return stddevs
#: Coefficient tables obtained by joining tables 2, 3, and 4, pages 47,
#: 50, 51.
COEFFS = CoeffsTable(sa_damping=5, table="""\
IMT c1 c2 c3 c4 c5 c6 c7 m50 m55 m80 r5 r20
pga 2.20 0.81 0.00 1.27 1.16 0.0021 9.3 0.55 0.59 0.50 0.54 0.20
0.03 4.00 0.79 0.00 1.57 1.83 0.0008 11.1 0.62 0.63 0.50 0.62 0.35
0.04 3.68 0.80 0.00 1.46 1.77 0.0013 10.5 0.62 0.63 0.50 0.57 0.29
0.10 2.37 0.81 0.00 1.10 1.02 0.0040 8.3 0.59 0.61 0.50 0.50 0.17
0.20 1.73 0.84 0.00 0.98 0.66 0.0042 7.5 0.60 0.64 0.56 0.45 0.12
0.40 1.07 1.05 -0.10 0.93 0.56 0.0033 7.1 0.63 0.68 0.64 0.45 0.12
1.00 0.09 1.42 -0.20 0.90 0.49 0.0023 6.8 0.63 0.64 0.67 0.45 0.12
2.00 -0.74 1.86 -0.31 0.92 0.46 0.0017 6.9 0.61 0.62 0.66 0.45 0.12
3.00 -0.74 1.86 -0.31 0.92 0.46 0.0017 6.9 0.61 0.62 0.66 0.45 0.12
4.00 -0.74 1.86 -0.31 0.92 0.46 0.0017 6.9 0.61 0.62 0.66 0.45 0.12
""")
[docs]class ToroEtAl2002SHARE(ToroEtAl2002):
#: Required rupture parameters are magnitude and rake
REQUIRES_RUPTURE_PARAMETERS = set(('mag', 'rake'))
[docs] 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.
"""
# extract faulting style and rock adjustment coefficients for the
# given imt
C_ADJ = self.COEFFS_FS_ROCK[imt]
mean, stddevs = super(ToroEtAl2002SHARE, self).\
get_mean_and_stddevs(sites, rup, dists, imt, stddev_types)
# apply faulting style and rock adjustment factor for mean and std
mean = np.log(np.exp(mean) *
_compute_faulting_style_term(C_ADJ['Frss'],
self.CONSTS_FS['pR'],
self.CONSTS_FS['Fnss'],
self.CONSTS_FS['pN'],
rup.rake) * C_ADJ['AFrock'])
stddevs = np.array(stddevs)
return mean, stddevs
#: Coefficients for faulting style and rock adjustment
COEFFS_FS_ROCK = CoeffsTable(sa_damping=5, table="""\
IMT Frss AFrock
pga 1.220000 0.735106
0.03 1.179400 0.423049
0.04 1.164000 0.477379
0.10 1.080000 0.888509
0.20 1.190000 1.197291
0.40 1.230000 1.308267
1.00 1.196667 1.265762
2.00 1.140000 1.215779
3.00 1.140000 1.215779
4.00 1.140000 1.215779
""")
#: Constants for faulting style adjustment
CONSTS_FS = {'Fnss': 0.95, 'pN': 0.01, 'pR': 0.81}