Source code for openquake.hazardlib.gsim.climent_1994

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2014-2025 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:'ClimentEtAl1994'.
"""
import numpy as np
# standard acceleration of gravity in m/s**2
from scipy.constants import g


from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, SA


def _compute_term_1_2(ctx, C):
    """
    Compute terms 1 and 2 in equation 1 page 1.
    """
    return C['c1'] + C['c2'] * ctx.mag


def _compute_term_3_4(ctx, C):
    """
    Compute term 3 and 4 in equation 1 page 1.
    """
    cutoff = 6.056877878
    rhypo = ctx.rhypo.copy()
    rhypo[rhypo <= cutoff] = cutoff
    return C['c3'] * np.log(rhypo) + C['c4'] * rhypo


def _get_site_amplification(ctx, imt, C):
    """
    Compute the fith term of the equation (1), p. 1:
    ``c5 * S``
    """
    S = _get_site_type_dummy_variables(ctx)
    return (C['c5'] * S)


def _get_site_type_dummy_variables(ctx):
    """
    Get site type dummy variables, ``S`` (for rock and soil sites)
    """
    S = np.zeros_like(ctx.vs30)
    # S=0 for rock sites, S=1 otherwise pag 1.
    idxS = (ctx.vs30 < 760.0)
    S[idxS] = 1
    return S


def _compute_mean(C, ctx, imt):
    """
    Compute mean value for PGA and pseudo-velocity response spectrum,
    as given in equation 1. Converts also pseudo-velocity response
    spectrum values to SA, using:

      SA = (PSV * W)/ratio(SA_larger/SA_geo_mean)
       W = (2 * pi / T)
       T = period (sec)
    """
    mean = (_compute_term_1_2(ctx, C) +
            _compute_term_3_4(ctx, C) +
            _get_site_amplification(ctx, imt, C))

    # convert from m/s**2 to g for PGA and from m/s to g for PSV
    # and divided this value for the ratio(SA_larger/SA_geo_mean)
    if imt.string == "PGA":
        mean = (np.exp(mean) / g) / C['r_SA']
    else:
        W = (2. * np.pi) / imt.period
        mean = np.exp(mean) * W / g / C['r_SA']
    return np.log(mean)


[docs]class ClimentEtAl1994(GMPE):
    """
    Implements GMPE developed by Climent, A, W. Taylor, M. Ciudad Real,
    W. Strauch, M. Villagran, A. Dahle, and H. Bungum. Published as a
    NORSAR report: "Spectral strong motion attenuation in Central Ame-
    rica", NORSAR Technical Report No. 2-17, 46 pp.
    The original formulation predict PGA (m/s*s) and 5% damped PSV (m/s)
    for the largest component of horizontal ground motion.
    In this implementation:
    Spectral acceleration (SA) values are obtained from PSV ones using
    the following formula :

        SA = [PSV * (2 * pi/ T)]/ratio(SA_larger/SA_geo_mean)
        StdDev.TOTAL=StdDev.TOTAL/sd_ratio(SA_larger/SA_geo_mean)

    The ratio() and sd_ratio() from Beyer and Bommer(2006)
    """
    #: Supported tectonic region type is active shallow crust and/or
    #: interface subduction the authors did not distinction between shallow
    #: and sudbdution events (see topic 5.3 "Shallow crustal vs.subduction
    #: events, pag. 32).
    #: Any factor/parameter is used in the formulation to discriminate between
    #: shallow or interface tectonic regime, here this GMPE is implemented
    #: for active_shallow_crust only
    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST

    #: Supported intensity measure types are spectral acceleration,
    #: and peak ground acceleration. See Table 2 in page 1865
    DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, SA}

    #: Supported intensity measure component is the largest component of
    #: two horizontal components
    #: :attr:`openquake.hazardlib.const.IMC.GREATER_OF_TWO_HORIZONTAL`,
    #: see paragraph before table on Summary, page 1.
    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = (
        const.IMC.GREATER_OF_TWO_HORIZONTAL)

    #: Supported standard deviation types is total.
    #: See equation 1 on the Summary and Table 4.1, page 22.
    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}

    #: Required site parameters. The GMPE was developed for rock and soil
    #: site conditions. The parameter S in eq. 1 (see Summary) define the
    #: soil condition: S=0 for rock, S=1 for soil.
    #: Here we use the Vs30=760 as limit between the two soil conditions
    REQUIRES_SITES_PARAMETERS = {'vs30'}

    #: Required rupture parameters are magnitude.
    REQUIRES_RUPTURE_PARAMETERS = {'mag'}

    #: Required distance measure is Rhypo, explained in page 1(eq. 1)
    REQUIRES_DISTANCES = {'rhypo'}

[docs]    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]
            mean[m] = _compute_mean(C, ctx, imt)
            sig[m] = C['SigmaB'] / C['r_std']

    #: Equation coefficients, described in Table 4.1 on pp. 22
    #: the original imt values are defined as frequencies values
    #: the sigma_ls was excluded
    COEFFS = CoeffsTable(sa_damping=5, table="""\
     IMT     c1      c2       c3       c4       c5      SigmaB r_SA    r_std
     pga    -1.6870  0.5530  -0.5370  -0.00302  0.3270  0.750  1.1000  1.0200
     0.025  -7.2140  0.5530  -0.5370  -0.00302  0.3270  0.750  1.1000  1.0200
     0.050  -5.4870  0.4470  -0.5500  -0.00246  0.3090  0.780  1.1000  1.0200
     0.100  -4.7260  0.4830  -0.5810  -0.00199  0.3810  0.800  1.2020  1.0200
     0.200  -4.8760  0.6420  -0.6420  -0.00156  0.4700  0.820  1.2040  1.0200
     0.500  -5.8620  0.9170  -0.7260  -0.00107  0.5660  0.820  1.2100  1.0200
     1.000  -6.7440  1.0810  -0.7560  -0.00077  0.5880  0.820  1.2200  1.0200
     2.000  -7.3480  1.1280  -0.7280  -0.00053  0.5360  0.790  1.2400  1.0200
     4.000  -7.4410  1.0070  -0.6010  -0.00040  0.4960  0.730  1.2800  1.0200
    """)