Source code for openquake.hazardlib.gsim.parker_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:`ParkerEtAl2020SInter` :class:`ParkerEtAl2020SInterB` :class:`ParkerEtAl2020SSlab` :class:`ParkerEtAl2020SSlabB`"""importmathimportnumpyasnpfromscipy.specialimporterffromopenquake.baselib.generalimportCallableDictfromopenquake.hazardlibimportconstfromopenquake.hazardlib.gsim.baseimportGMPE,CoeffsTable,add_aliasfromopenquake.hazardlib.imtimportPGA,SA,PGVCONSTANTS={"b4":0.1,"f3":0.05,"Vb":200,"vref_fnl":760,"V1":270,"vref":760}_a0=CallableDict()@_a0.add(const.TRT.SUBDUCTION_INTERFACE)def_a0_1(trt,region,basin,C,C_PGA):""" Regional anelastic coefficient, a0 """ifregionisNoneorregion=="Cascadia":a0=C["a0"]a0_pga=C_PGA["a0"]else:a0=C[region+"_a0"]a0_pga=C_PGA[region+"_a0"]returna0,a0_pga@_a0.add(const.TRT.SUBDUCTION_INTRASLAB)def_a0_2(trt,region,basin,C,C_PGA):""" Regional anelastic coefficient for subduction slab, a0 """ifregionisNone:returnC["a0slab"],C_PGA["a0slab"]returnC[region+"_a0slab"],C_PGA[region+"_a0slab"]def_basin_term(region,basin,C,ctx=None):""" Basin term main handler. """ifnothasattr(ctx,'z2pt5'):return0ifregion=="JP":return_get_basin_term_factors(3.05,-0.8,500,0.33,C["J_e1"],C["J_e2"],C["J_e3"],ctx)ifregion=="Cascadia":ifbasinisNone:return_get_basin_term_factors(3.94,-0.42,200,0.2,C["C_e1"],C["C_e2"],C["C_e3"],ctx)ifbasin=="out":dn=C["del_None"]return_get_basin_term_factors(3.94,-0.42,200,0.2,C["C_e1"],C["C_e2"]+dn,C["C_e3"]+dn,ctx)ifbasin=="Seattle":ds=C["del_Seattle"]return_get_basin_term_factors(3.94,-0.42,200,0.2,C["C_e1"],C["C_e2"]+ds,C["C_e3"]+ds,ctx)return0_c0=CallableDict()@_c0.add(const.TRT.SUBDUCTION_INTERFACE)def_c0_1(trt,region,saturation_region,C,C_PGA):""" c0 factor. """ifsaturation_regionisNone:c0_col="c0"else:c0_col=saturation_region+"_c0"returnC[c0_col],C_PGA[c0_col]@_c0.add(const.TRT.SUBDUCTION_INTRASLAB)def_c0_2(trt,region,saturation_region,C,C_PGA):""" c0 factor. """ifsaturation_regionisNone:c0_col="c0slab"elifregionin["AK","SA"]:c0_col=saturation_region+"_c0slab"else:# no more specific region availablec0_col=region+"_c0slab"returnC[c0_col],C_PGA[c0_col]_depth_scaling=CallableDict()@_depth_scaling.add(const.TRT.SUBDUCTION_INTERFACE)def_depth_scaling_1(trt,C,ctx):""" Depth scaling is for slab. """return0@_depth_scaling.add(const.TRT.SUBDUCTION_INTRASLAB)def_depth_scaling_2(trt,C,ctx):res=C["m"]*(ctx.hypo_depth-C["db"])+C["d"]res[ctx.hypo_depth>=C["db"]]=C["d"]res[ctx.hypo_depth<=20]=C["m"]*(20-C["db"])+C["d"]returnresdef_get_basin_term_factors(theta0,theta1,vmu,vsig,e1,e2,e3,ctx):""" Basin term for given factors. """btf=np.zeros_like(ctx.vs30)select=ctx.z2pt5!=0iflen(select)==0:returnbtfvs30=ctx.vs30[select]z2pt5=ctx.z2pt5[select]z2pt5_pred=10**(theta0+theta1*(1+erf((np.log10(vs30)-math.log10(vmu))/(vsig*math.sqrt(2)))))del_z2pt5=np.log(z2pt5)-np.log(z2pt5_pred)btf[select]=np.where(del_z2pt5<=(e1/e3),e1,np.where(del_z2pt5>=(e2/e3),e2,e3*del_z2pt5))returnbtfdef_linear_amplification(region,C,vs30):""" Linear site term. """# site coefficientsv1=CONSTANTS["V1"]vref=CONSTANTS["vref"]ifregionisNoneorregion=="CAM":s2=C["s2"]s1=s2elifregion=="TW"orregion=="JP":s2=C[region+"_s2"]s1=C[region+"_s1"]else:s2=C[region+"_s2"]s1=s2# linear site termfnl=np.where(vs30<=v1,s1*np.log(vs30/v1)+s2*math.log(v1/vref),0)fnl=np.where((v1<vs30)&(vs30<=C["V2"]),s2*np.log(vs30/vref),fnl)fnl=np.where(vs30>C["V2"],s2*math.log(C["V2"]/vref),fnl)returnfnldef_magnitude_scaling(sfx,C,C_PGA,mag,m_b):""" Magnitude scaling factor. """m_diff=mag-m_bfm=np.where(m_diff>0,C["c6"+sfx]*m_diff,C["c4"+sfx]*m_diff+C["c5"+sfx]*m_diff**2)fm_pga=np.where(m_diff>0,C_PGA["c6"+sfx]*m_diff,C_PGA["c4"+sfx]*m_diff+C_PGA["c5"+sfx]*m_diff**2)returnfm,fm_pgadef_non_linear_term(C,imt,vs30,fp,fm,c0,fd=0):""" Non-linear site term. """# fd for slab onlypgar=np.exp(fp+fm+c0+fd)ifhasattr(imt,"period")andimt.period>=3:fnl=0else:fnl=C["f4"]*(np.exp(C["f5"]*(np.minimum(vs30,CONSTANTS["vref_fnl"])-CONSTANTS["Vb"]))-math.exp(C["f5"]*(CONSTANTS["vref_fnl"]-CONSTANTS["Vb"])))fnl*=np.log((pgar+CONSTANTS["f3"])/CONSTANTS["f3"])returnfnldef_path_term(trt,region,basin,suffix,C,C_PGA,mag,rrup,m_b):""" Path term. """h=_path_term_h(trt,mag,m_b)r=np.sqrt(rrup**2+h**2)# log(R / Rref)r_rref=np.log(r/np.sqrt(1+h**2))a0,a0_pga=_a0(trt,region,basin,C,C_PGA)c1n="c1"+suffixfp=C[c1n]*np.log(r)+(CONSTANTS["b4"]*mag)*r_rref+a0*rfp_pga=C_PGA[c1n]*np.log(r)+(CONSTANTS["b4"]*mag)*r_rref+a0_pga*rreturnfp,fp_pga_path_term_h=CallableDict()@_path_term_h.add(const.TRT.SUBDUCTION_INTERFACE)def_path_term_h_1(trt,mag,m_b=None):""" H factor for path term. """return10**(-0.82+0.252*mag)@_path_term_h.add(const.TRT.SUBDUCTION_INTRASLAB)def_path_term_h_2(trt,mag,m_b=None):""" H factor for path term, subduction slab. """m=(math.log10(35)-math.log10(3.12))/(m_b-4)returnnp.where(mag<=m_b,10**(m*(mag-m_b)+math.log10(35)),35)
[docs]defget_stddevs(C,rrup,vs30):""" Returns the standard deviations. Generate tau, phi, and total sigma computed from both total and partitioned phi models. """# define period-independent coefficients for phi modelsv1=200v2=500r1=200r2=500# total Phiphi_rv=np.zeros(len(vs30))fori,vs30iinenumerate(vs30):ifrrup[i]<=r1:phi_rv[i]=C["phi21"]elifrrup[i]>=r2:phi_rv[i]=C["phi22"]else:phi_rv[i]=((C["phi22"]-C["phi21"])/(math.log(r2)-math.log(r1))) \
*(math.log(rrup[i])-math.log(r1))+C["phi21"]ifvs30i<=v1:phi_rv[i]+=C["phi2V"]* \
(math.log(r2/max(r1,min(r2,rrup[i])))/math.log(r2/r1))elifvs30i<v2:phi_rv[i]+=C["phi2V"]* \
((math.log(v2/min(v2,vs30i)))/math.log(v2/v1))* \
(math.log(r2/max(r1,min(r2,rrup[i])))/math.log(r2/r1))phi_tot=np.sqrt(phi_rv)return[np.sqrt(C["Tau"]**2+phi_tot**2),C["Tau"],phi_tot]
[docs]classParkerEtAl2020SInter(GMPE):""" Implements Parker et al. (2020) for subduction interface. """DEFINED_FOR_TECTONIC_REGION_TYPE=const.TRT.SUBDUCTION_INTERFACE#: Supported intensity measure types are spectral acceleration,#: peak ground acceleration and peak ground velocityDEFINED_FOR_INTENSITY_MEASURE_TYPES={PGV,PGA,SA}#: Supported intensity measure component is the geometric mean componentDEFINED_FOR_INTENSITY_MEASURE_COMPONENT=const.IMC.GEOMETRIC_MEANDEFINED_FOR_STANDARD_DEVIATION_TYPES={const.StdDev.TOTAL,const.StdDev.INTER_EVENT,const.StdDev.INTRA_EVENT}#: Site amplification is dependent only upon Vs30REQUIRES_SITES_PARAMETERS={'vs30'}#: Required rupture parameters are only magnitude for the interface modelREQUIRES_RUPTURE_PARAMETERS={'mag'}#: Required distance measure is closest distance to rupture, for#: interface eventsREQUIRES_DISTANCES={'rrup'}REQUIRES_ATTRIBUTES={'region','saturation_region','basin'}def__init__(self,region=None,saturation_region=None,basin=None):""" Enable setting regions to prevent messy overriding and code duplication. """self.region=regionifsaturation_regionisNone:self.saturation_region=regionelse:self.saturation_region=saturation_regionself.basin=basin
[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. """trt=self.DEFINED_FOR_TECTONIC_REGION_TYPEC_PGA=self.COEFFS[PGA()]form,imtinenumerate(imts):C=self.COEFFS[imt]# Regional Mb factorifself.saturation_regioninself.MB_REGIONS:m_b=self.MB_REGIONS[self.saturation_region]else:m_b=self.MB_REGIONS["default"]c0,c0_pga=_c0(trt,self.region,self.saturation_region,C,C_PGA)fm,fm_pga=_magnitude_scaling(self.SUFFIX,C,C_PGA,ctx.mag,m_b)fp,fp_pga=_path_term(trt,self.region,self.basin,self.SUFFIX,C,C_PGA,ctx.mag,ctx.rrup,m_b)fd=_depth_scaling(trt,C,ctx)fd_pga=_depth_scaling(trt,C_PGA,ctx)fb=_basin_term(self.region,self.basin,C,ctx)flin=_linear_amplification(self.region,C,ctx.vs30)fnl=_non_linear_term(C,imt,ctx.vs30,fp_pga,fm_pga,c0_pga,fd_pga)# The output is the desired median model prediction in LN units# Take the exponential to get PGA, PSA in g or the PGV in cm/smean[m]=fp+fnl+fb+flin+fm+c0+fdsig[m],tau[m],phi[m]=get_stddevs(C,ctx.rrup,ctx.vs30)
[docs]classParkerEtAl2020SInterB(ParkerEtAl2020SInter):""" For Cascadia and Japan where basins are defined (also require z2pt5). """REQUIRES_SITES_PARAMETERS={'vs30','z2pt5'}
[docs]classParkerEtAl2020SSlab(ParkerEtAl2020SInter):""" Modifications for subduction slab. """DEFINED_FOR_TECTONIC_REGION_TYPE=const.TRT.SUBDUCTION_INTRASLAB# slab also requires hypo_depthREQUIRES_RUPTURE_PARAMETERS={'mag','hypo_depth'}# constant table suffixSUFFIX="slab"MB_REGIONS={"Aleutian":7.98,"AK":7.2,"Cascadia":7.2,"CAM_S":7.6,"CAM_N":7.4,"JP_Pac":7.65,"JP_Phi":7.55,"SA_N":7.3,"SA_S":7.25,"TW_W":7.7,"TW_E":7.7,"default":7.6}
[docs]classParkerEtAl2020SSlabB(ParkerEtAl2020SSlab):""" For Cascadia and Japan where basins are defined (also require z2pt5). """REQUIRES_SITES_PARAMETERS={'vs30','z2pt5'}
