Aims

  • On-the-fly calculation
  • Multi-layered logic tree


Logic Tree layout


The weight of a branch is given by W_{ERF} * W_{AS} * W_{V} * W_{SS} * W_{SI}
where W_{ERF} is the weight of the ERF used and W_{AS}, W_{V}, ... are the weights of the GMMs choosen

This Logic tree can either used exhaustively or using a Monte Carlo sampling approach


Computation

To be able to compute results on the fly, evaluation of a single branch has to be very fast (i.e. fully vectorised).

This requires:

  • Fully vectorised GMMs across faults
  • Fully vectorised computation of hazard, disagg and UHS using multiple GMMs


Scribbles and code for vectorised hazard calculation below


import time

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

from empirical.util import empirical_factory
from empirical.util import classdef
from empirical.GMM_models import Bradley_2010_Sa as br10
from empirical.util import openquake_wrapper

import gmhazard_calc as gc

fault = classdef.Fault(
    Mw=6.0,
    hdepth=10.0,
    ztor=10.0,
    dip=90,
    rake=60,
    tect_type=classdef.TectType.ACTIVE_SHALLOW,
)

site = classdef.Site(rrup=120, rjb=120, rx=90, hw=True, rtvz=0.0, vs30=450, vs30measured=False, z1p0=classdef.estimate_z1p0(450))

faults_1 = []
for cur_mw in np.linspace(4, 8, 10000):
    faults_1.append(classdef.Fault(
        Mw=cur_mw,
        hdepth=10.0,
        ztor=10.0,
        dip=90,
        rake=60,
        tect_type=classdef.TectType.ACTIVE_SHALLOW,
    ))

faults_2 = []
for cur_mw in np.linspace(4, 8, 10000):
    faults_2.append(classdef.Fault(
        Mw=cur_mw,
        hdepth=10.0,
        ztor=10.0,
        dip=50,
        rake=0,
        tect_type=classdef.TectType.ACTIVE_SHALLOW,
    ))


sites = [site for ix in range(10000)]

rup_rates = np.random.uniform(0, 0.1, 10_000)

# N = Number of faults
# M = Number of IM values
# K = Number of GMMs
N = 10_000
M = 1000
K = 2

# Get GMM results, (N, K)
start_time = time.perf_counter()
gmm_results_1 = br10.calculate_Bradley_multi(sites, faults_1, np.asarray([0]))
gmm_results_2 = br10.calculate_Bradley_multi(sites, faults_2, np.asarray([0]))
mu = np.stack((gmm_results_1[0].ravel(), gmm_results_2[0].ravel()), axis=1)
sigma = np.stack((gmm_results_1[1][0].ravel(), gmm_results_2[1][0].ravel()), axis=1)
print(f"Took {time.perf_counter() - start_time} to get GMM results")

# Get IM values,  (M)
start_time = time.perf_counter()
im_values = gc.utils.get_im_values(gc.im.IM.from_str("PGA"), n_values=1000)
print(f"Took {time.perf_counter() - start_time} to get IM values")

# Compute GMEP
start_time = time.perf_counter()
gmep = stats.norm.sf(im_values, mu.reshape(-1, 1), sigma.reshape(-1, 1)).reshape(10_000, 2, 1000).transpose(0, 2, 1)
print(f"Took {time.perf_counter() - start_time} to compute GMEP")

# Compute exceedance probability
start_time = time.perf_counter()
excd_prob_rup = gmep * rup_rates.reshape((-1, 1, 1))
excd_prob = np.sum(excd_prob_rup, axis=0)
print(f"Took {time.perf_counter() - start_time} to compute annual exceedance prob")



  • No labels