Source code for openquake.hazardlib.gsim.phung_2020
# -*- coding: utf-8 -*-# vim: tabstop=4 shiftwidth=4 softtabstop=4## Copyright (C) 2015-2018 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:`PhungEtAl2020SInter` :class:`PhungEtAl2020SSlab` :class:`PhungEtAl2020Asc`"""importmathimportnumpyasnpfromopenquake.hazardlibimportconstfromopenquake.hazardlib.gsim.baseimportGMPE,CoeffsTablefromopenquake.hazardlib.imtimportPGA,SA
[docs]defget_stddevs(C):""" Return standard deviations. """phi_tot=math.sqrt(C['phiss']**2+C['phis2s']**2)return[np.sqrt(C["tau"]**2+phi_tot**2),C["tau"],phi_tot]
def_get_basin_term(C,ctx,region):""" Basin term [16]. """vs30=ctx.vs30ifregion=='glb':return0ifregion=='tw':ez_1=np.exp(-3.73/2*np.log((vs30**2+290.53**2)/(1750**2+290.53**2)))phi6=300elifregionin['ca','jp']:ez_1=np.exp(-5.23/2*np.log((vs30**2+412.39**2)/(1360**2+412.39**2)))ifregion=='ca':phi6=300else:phi6=800d_z1=ctx.z1pt0-ez_1returnnp.where(ctx.z1pt0<0,0,C['phi5'+region]*(1-np.exp(-d_z1/phi6)))def_distance_attenuation(s,region,aftershocks,C,mag,rrup,ztor):""" Distance scaling and attenuation term [8, 9, 10]. """del_c5=C['dp']*np.maximum(ztor/50-20/50,0)* \
(ztor>20)*(mag<7.0)cns=(C['c5']+del_c5)*np.cosh(C['c6']*np.maximum(mag-C['c_hm'],0))f8=s['c4']*np.log(rrup+cns)f9=(s['c4_a']-s['c4'])*np.log(np.sqrt(rrup**2+s['c_rb']**2))f10=(C['c_g1'+region]+aftershocks*C['dc_g1as']+C['c_g2']/(np.cosh(np.maximum(mag-C['c_g3'],0))))*rrupreturnf8+f9+f10def_fsof_ztor(C,mag,rake,ztor):""" Factors requiring type of fault. """f234=np.zeros_like(rake)e_ztor=ztor.copy()pos=(30<=rake)&(rake<=150)e_ztor[pos]=np.maximum(3.5384-2.60*np.maximum(mag[pos]-5.8530,0),0)**2neg=(-120<=rake)&(rake<=-60)e_ztor[neg]=np.maximum(2.7482-1.7639*np.maximum(mag[neg]-5.5210,0),0)**2# reverse fault [2]f234[pos]+=C['c1_a']+C['c1_c']/np.cosh(2*np.maximum(mag[pos]-4.5,0))# normal fault [3]f234[neg]+=C['c1_b']+C['c1_d']/np.cosh(2*np.maximum(mag[neg]-4.5,0))ztor=np.where(ztor<0,e_ztor,ztor)delta_ztor=ztor-e_ztor# [4]f234+=(C['c7']+C['c7_b']/np.cosh(2*np.maximum(mag-4.5,0)))*delta_ztorreturnf234,ztor
[docs]classPhungEtAl2020Asc(GMPE):""" Implements Phung et al. (2020) for crustal. """DEFINED_FOR_TECTONIC_REGION_TYPE=const.TRT.ACTIVE_SHALLOW_CRUSTDEFINED_FOR_INTENSITY_MEASURE_TYPES=set([PGA,SA])DEFINED_FOR_INTENSITY_MEASURE_COMPONENT=const.IMC.GEOMETRIC_MEANDEFINED_FOR_REFERENCE_VELOCITY=1130DEFINED_FOR_STANDARD_DEVIATION_TYPES={const.StdDev.TOTAL,const.StdDev.INTER_EVENT,const.StdDev.INTRA_EVENT}#: rx for hanging wallREQUIRES_DISTANCES={'rjb','rrup','rx'}REQUIRES_RUPTURE_PARAMETERS={'dip','mag','rake','ztor'}REQUIRES_SITES_PARAMETERS={'vs30','z1pt0'}def__init__(self,region='glb',aftershocks=False,d_dpp=0):# region options:# 'glb', 'tw', 'ca', 'jp' (global, Taiwan, California, Japan)self.region=region# only for Taiwan regionself.aftershocks=aftershocksandregion=='tw'# direct point parameter for directivity effectself.d_dpp=d_dpp
[docs]defcompute(self,ctx:np.recarray,imts,mean,sig,tau,phi):""" See :meth:`superclass method <.base.GroundShakingIntensityModel.compute>` for spec of input and result values. """s=self.CONSTANTSform,imtinenumerate(imts):C=self.COEFFS[imt]lnmed,ztor=_fsof_ztor(C,ctx.mag,ctx.rake,ctx.ztor)# main shock [1]lnmed+=C['c1']# dip term [5]lnmed+=(C['c11']+C['c11_b']/np.cosh(2*np.maximum(ctx.mag-4.5,0)))*np.cos(np.radians(ctx.dip))**2# hanging wall term [12]lnmed+=(ctx.rx>=0)*C['c9']*np.cos(np.radians(ctx.dip))\
*(C['c9_a']+(1-C['c9_a'])*np.tanh(ctx.rx/C['c9_b']))\
*(1-np.sqrt(ctx.rjb**2+ztor**2)/(ctx.rrup+1))# directivity [11]lnmed+=C['c8'] \
*np.maximum(1-np.maximum(ctx.rrup-40,0)/30,0) \
*np.clip((ctx.mag-5.5)/0.8,0.,1.) \
*np.exp(s['c8_a']*(ctx.mag-C['c8_b'])**2) \
*self.d_dpp# fmag [6, 7]lnmed+=s['c2']*(ctx.mag-6)lnmed+=(s['c2']-C['c3'])/C['c_n'] \
*np.log(1+np.exp(C['c_n']*(C['c_m']-ctx.mag)))lnmed+=_distance_attenuation(s,self.region,self.aftershocks,C,ctx.mag,ctx.rrup,ztor)sa1130=np.exp(lnmed)# site response [14, 15]lnmed+=C['phi1'+self.region]*np.minimum(np.log(ctx.vs30/1130),0)lnmed+=C['phi2']*(np.exp(C['phi3']*(np.minimum(ctx.vs30,1130)-360))-np.exp(C['phi3']*(1130-360)))*np.log((sa1130+C['phi4'])/C['phi4'])# basin term [16]lnmed+=_get_basin_term(C,ctx,self.region)mean[m]=lnmedsig[m],tau[m],phi[m]=get_stddevs(C)
def_get_stddevs(region,C):""" Return standard deviations. """phi_tot=np.sqrt(C['phiss4'+region]**2+C['phis2s4'+region]**2)return[np.sqrt(C["tau4"+region]**2+phi_tot**2),C["tau4"+region],phi_tot]def_fmag(region,C,mag,pga=False):""" Magnitude term. """a4=C['a4']+(C['a4_del']ifregion=='jptw'else0)returnnp.where(mag<=C['mref'],a4*(mag-C['mref'])+C['a13']*(10-mag)**2,C['a5']*(mag-C['mref'])+C['a13']*(10-mag)**2)def_fp(trt,region,s,C,mag,rrup,pga=False):""" Path term for subduction interface. """a7,a14=(C['a7'],C['a14'])if(trt==const.TRT.SUBDUCTION_INTRASLAB)else(0,0)a1=C['a1']ifregion=='jptw':a1+=C['a1_del']r='tw'ifpgaelseregionreturna1+a7+(C['a2']+a14+s['a3']*(mag-7.8)) \
*np.log(rrup+s['c4']*np.exp(s['a9']*(mag-6))) \
+C['a6'+r]*rrupdef_fsite(s,region,C,vs30,pga_1000):""" Non-linear component. """result=C['a12'+region]*np.log(vs30/C['vlin'])idx=vs30<C['vlin']result[idx]+= \
-C['b']*np.log(pga_1000[idx]+s['c']) \
+C['b']*np.log(pga_1000[idx]+s['c']*(vs30[idx]/C['vlin'])**s['n'])result[~idx]+=C['b']*s['n']*np.log(vs30[~idx]/C['vlin'])returnresultdef_fz1pt0(region,C,vs30,z1pt0):""" Basin depth term. """result=np.zeros_like(z1pt0)idx=np.where(z1pt0>=0)ifregion=='tw':ez_1=np.exp(-3.96/2*np.log((vs30**2+352.7**2)/(1750**2+352.7**2)))elifregion=='jptw':ez_1=np.exp(-5.23/2*np.log((vs30**2+412.39**2)/(1360**2+412.39**2)))result[idx]=C['a8'+region]*np.minimum(np.log(z1pt0[idx]/ez_1),1)returnresultdef_fztor(trt,C,ztor):""" Ztor factor for subduction interface. """iftrt==const.TRT.SUBDUCTION_INTRASLAB:returnC['a11']*(np.minimum(ztor,80)-40)returnC['a10']*(np.minimum(ztor,40)-20)
[docs]defpga_rock(trt,region,s,C_PGA,mag,rrup,ztor,vs30):""" PGA at Vs30 (as Taiwan region, C_PGA) """f_mag=_fmag(region,C_PGA,mag,pga=True)f_ztor=_fztor(trt,C_PGA,ztor)f_p=_fp(trt,region,s,C_PGA,mag,rrup,pga=True)# site functionf_site=C_PGA['a12tw']*np.log(vs30/C_PGA['vlin']) \
+C_PGA['b']*s['n']*np.log(vs30/C_PGA['vlin'])returnnp.exp(f_mag+f_p+f_ztor+f_site)
[docs]classPhungEtAl2020SInter(GMPE):""" Implements Phung et al. (2020) for Subduction Interface. """DEFINED_FOR_TECTONIC_REGION_TYPE=const.TRT.SUBDUCTION_INTERFACEDEFINED_FOR_INTENSITY_MEASURE_TYPES={PGA,SA}DEFINED_FOR_INTENSITY_MEASURE_COMPONENT=const.IMC.GEOMETRIC_MEANDEFINED_FOR_REFERENCE_VELOCITY=1000DEFINED_FOR_STANDARD_DEVIATION_TYPES={const.StdDev.TOTAL,const.StdDev.INTER_EVENT,const.StdDev.INTRA_EVENT}#: rjb and rx not required for subductionREQUIRES_DISTANCES={'rrup'}#: dip and rake not required for subductionREQUIRES_RUPTURE_PARAMETERS={'mag','ztor'}REQUIRES_SITES_PARAMETERS={'vs30','z1pt0'}def__init__(self,region='jptw'):# region options:# 'jptw', 'tw' (Japan/Taiwan joined, Taiwan)self.region=region
[docs]defcompute(self,ctx:np.recarray,imts,mean,sig,tau,phi):""" See :meth:`superclass method <.base.GroundShakingIntensityModel.compute>` for spec of input and result values. """# extract dictionaries of coefficients specific to IM typetrt=self.DEFINED_FOR_TECTONIC_REGION_TYPEC_PGA=self.COEFFS[PGA()]s=self.CONSTANTSform,imtinenumerate(imts):C=self.COEFFS[imt]f_mag=_fmag(self.region,C,ctx.mag)f_ztor=_fztor(trt,C,ctx.ztor)# path termf_p=_fp(trt,self.region,s,C,ctx.mag,ctx.rrup)# PGA at rock with Vs30 = 1000 m/spga_1000=pga_rock(trt,self.region,self.CONSTANTS,C_PGA,ctx.mag,ctx.rrup,ctx.ztor,1000)# site effectvs30=np.minimum(ctx.vs30,1000)# non-linear componentf_site=_fsite(s,self.region,C,vs30,pga_1000)# basin depth termf_z1pt0=_fz1pt0(self.region,C,vs30,ctx.z1pt0)# median total and stddevmean[m]=f_mag+f_p+f_ztor+f_site+f_z1pt0sig[m],tau[m],phi[m]=_get_stddevs(self.region,C)
[docs]classPhungEtAl2020SSlab(PhungEtAl2020SInter):""" Implements Phung et al. (2020) for Subduction Intraslab. """DEFINED_FOR_TECTONIC_REGION_TYPE=const.TRT.SUBDUCTION_INTRASLAB
