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")