From 999dba56fb936f628c60f5724a1acc4514b12bf9 Mon Sep 17 00:00:00 2001
From: Jason Altekruse <jaltekruse@contractor.usgs.gov>
Date: Wed, 13 Jan 2021 16:48:56 -0700
Subject: [PATCH] fixed issues with median gm calc, used expected results from
official python implementation (modified to use updated coeffs). uncertainty
not yet implemented
---
.../earthquake/nshmp/gmm/ParkerEtAl_2020.java | 99 ++++++++++++++-----
.../resources/gmm/coeffs/PSHAB20_slab.csv | 2 +-
.../usgs/earthquake/nshmp/gmm/PSHAB20.java | 3 +-
src/test/resources/gmm/PSHAB20_results.csv | 46 ++++-----
4 files changed, 98 insertions(+), 52 deletions(-)
diff --git a/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java b/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java
index b00ea58c..f328e240 100644
--- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java
+++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java
@@ -21,6 +21,7 @@ import gov.usgs.earthquake.nshmp.gmm.GmmInput.Constraints;
/**
* Implementation of the PEER NGA-Subduction GMM, Parker et al., 2020 GMM
+ * (PSHAB20)
*
* Notes:
*
@@ -28,6 +29,11 @@ import gov.usgs.earthquake.nshmp.gmm.GmmInput.Constraints;
* analysis is planned to determine if there are attenuation differences between
* forearc and backarc zones.
*
+ * PSHAB20 uses zHyp (randomized source location on fault), but this can be
+ * replaced this event-specific zHyp with the mean depth expected for the fault
+ * plane, zHypBar, in Eq. 4.6 without modification to coefficients or
+ * uncertainty.
+ *
* <p><b>Reference:</b> Parker, G.A, Stewart, J.P., Boore, D.M., Atkinson, G.M.,
* Hassani, B. (2020) NGA-Subduction global ground-motion models with regional
* adjustment factors. Report no. 2020/03. Berkeley, CA: PEER, 131 p.
@@ -37,14 +43,6 @@ import gov.usgs.earthquake.nshmp.gmm.GmmInput.Constraints;
* @author U.S. Geological Survey
*/
-/*
- * Regions: global, Alaska, Cascadia. Unused Regions: CAM (Central America),
- * Japan, SA (South America), Taiwan
- *
- * Saturation Regions: global, Aleutian, Alaska, Cascadia. Unused Regions:
- * Central America S, Central America N, Japan Pacific, Japan Phillipines, South
- * America N, South America S, Taiwan W, Taiwan E
- */
public abstract class ParkerEtAl_2020 implements GroundMotionModel {
static final String NAME = "Parker et al. (2020) NGA-Subduction";
@@ -62,6 +60,14 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
.set(VS30, Range.closed(150.0, 2000.0))
.build();
+ /*
+ * Regions: global, Alaska, Cascadia, [CAM (Central America), Japan, SA (South
+ * America), Taiwan]
+ *
+ * Saturation Regions: global, Aleutian, Alaska, Cascadia, [Central America S,
+ * Central America N, Japan Pacific, Japan Phillipines, South America N, South
+ * America S, Taiwan W, Taiwan E]
+ */
static final class Region {
static final String GLOBAL = "Global";
static final String ALASKA = "Alaska";
@@ -253,11 +259,21 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
final double Mw,
final double rRup,
final double zHyp) {
- double h = nearSourceSaturationTerm(Mw, c.mc, isSlab);
- return c.c0 +
- getPathTerm(c, rRup, Mw, h) +
- getMagnitudeTerm(c, Mw) +
- getDepthTerm(c, isSlab, zHyp);
+ double h = nearSourceSaturation(Mw, c.mc, isSlab);
+ double Fp = getPathTermPy(c, rRup, Mw, h);
+ double Fm = getMagnitudeTerm(c, Mw);
+ double Fd = getDepthTerm(c, isSlab, zHyp);
+
+ // System.out.println("calcMean Reference for T = " + c.imt.toString());
+ // System.out.printf(" mc = %.2f\n", c.mc);
+ // System.out.printf(" h = %.6f\n", h);
+ // System.out.printf(" Fp = %.6f\n", Fp);
+ System.out.printf(" Fm(ref) = %.6f\n", Fm);
+ // System.out.printf(" Fd = %.6f\n", Fd);
+ // System.out.printf(" Fb = %.6f (reference calc)\n", 0.0);
+ // System.out.printf(" mu = %.6f\n", c.c0 + Fp + Fm + Fd);
+
+ return c.c0 + Fp + Fm + Fd;
}
/* Equation 4.1 */
@@ -272,28 +288,57 @@ public abstract class ParkerEtAl_2020 implements GroundMotionModel {
final double vs30,
final double z2p5) {
- /* Equation 4.7: Fs = Flin + Fnl + Fb */
- double siteTerm = getSiteTermLinear(c, vs30) + getSiteTermNonLinear(c, pgaRef, vs30) +
- getSiteBasinTerm(c, basin, vs30, z2p5);
-
- double h = nearSourceSaturationTerm(Mw, c.mc, isSlab);
-
- double μ = c.c0 +
- getPathTerm(c, rRup, Mw, h) +
- getMagnitudeTerm(c, Mw) +
- getDepthTerm(c, isSlab, zHyp) +
- siteTerm;
+ double h = nearSourceSaturation(Mw, c.mc, isSlab);
+ double Fp = getPathTermPy(c, rRup, Mw, h);
+ double Fm = getMagnitudeTerm(c, Mw);
+ double Fd = getDepthTerm(c, isSlab, zHyp);
- return μ;
+ /* Equation 4.7: Fs = Flin + Fnl + Fb */
+ double Fslin = getSiteTermLinear(c, vs30);
+ double Fsnl = getSiteTermNonLinear(c, pgaRef, vs30);
+ double Fsb = getSiteBasinTerm(c, basin, vs30, z2p5);
+ double Fs = Fslin + Fsnl + Fsb;
+
+ System.out.println("calcMean for T = " + c.imt.toString());
+ System.out.printf(" Vs30 = %.2f\n", vs30);
+ System.out.printf(" h = %.6f\n", h);
+ System.out.printf(" Fp = %.6f\n", Fp);
+ System.out.printf(" Fm = %.6f\n", Fm);
+ System.out.printf(" Fd = %.6f\n", Fd);
+ System.out.printf(" PGAref = %.6f\n", pgaRef);
+ System.out.printf(" Fslin = %.6e\n", Fslin);
+ System.out.printf(" Fsnl = %.6e\n", Fsnl);
+ System.out.printf(" Fssb = %.6e\n", Fsb);
+ System.out.printf(" Fs = %.6e\n", Fs);
+ System.out.printf(" mu = %.6f\n", c.c0 + Fp + Fm + Fd + Fs);
+ System.out.printf(" y = %.6f\n", exp(c.c0 + Fp + Fm + Fd + Fs));
+
+ return c.c0 + Fp + Fm + Fd + Fs;
}
/* Equations 4.2 and 4.3 */
private static final double getPathTerm(final Coefficients c, double rRup, double Mw, double h) {
double R = Maths.hypot(rRup, h);
- return c.c1 * log(R) + c.b4 * Mw * log(R) + c.a0 * R;
+ return c.c1 * log(R) + c.b4 * Mw * log(rRup / R) + c.a0 * R;
+ }
+
+ /* Equations 4.2 and 4.3, following official Python implementation */
+ private static final double getPathTermPy(final Coefficients c, double rRup, double Mw,
+ double h) {
+ double Rref = Maths.hypot(1, h);
+ double R = Maths.hypot(rRup, h);
+ double R_Rref = R / Rref;
+ // System.out.printf(" c1 = %.6e\n", c.c1);
+ // System.out.printf(" b4 = %.6e\n", c.b4);
+ // System.out.printf(" a0 = %.6e\n", c.a0);
+ // System.out.printf(" Rref = %.6e\n", R_Rref);
+ // System.out.printf(" R = %.6e\n", R);
+ // System.out.printf(" Mw = %.6e\n", Mw);
+
+ return c.c1 * log(R) + c.b4 * Mw * log(R_Rref) + c.a0 * R;
}
- private static final double nearSourceSaturationTerm(double Mw, double mc, boolean slab) {
+ private static final double nearSourceSaturation(double Mw, double mc, boolean slab) {
if (slab) {
/* Equation 4.4b */
return (Mw <= mc) ? pow(10, (1.050 / (mc - 4.0)) * (Mw - mc) + 1.544) : 35;
diff --git a/src/main/resources/gmm/coeffs/PSHAB20_slab.csv b/src/main/resources/gmm/coeffs/PSHAB20_slab.csv
index 78a98065..ab56cdbb 100644
--- a/src/main/resources/gmm/coeffs/PSHAB20_slab.csv
+++ b/src/main/resources/gmm/coeffs/PSHAB20_slab.csv
@@ -1,6 +1,6 @@
T,Global_c0,Alaska_c0,Aleutian_c0,Cascadia_c0,c1,b4,Global_a0,Alaska_a0,Cascadia_a0,c4,c5,c6,d,m,db,V1,V2,Vref,Global_s2,Alaska_s2,Cascadia_s2,f4,f5,C_e1,C_e2,C_e3,del_None,del_Seattle,Tau,phi21,phi22,phi2V,VM,phi2S2S0,a1,phi2SS1,phi2SS2,a2
PGV,13.194,12.79,13.6,12.874,-2.422,0.1,-1.90E-03,-2.38E-03,-1.09E-03,1.84,-0.05,0.8,0.2693,0.0252,67,270,850,760,-0.601,-1.031,-0.671,-0.31763,-0.0052,0,0.115,0.068,-0.115,0,0.477,0.348,0.288,-0.179,423,0.142,0.047,0.153,0.166,0.011
-PGV,9.907,9.404,9.912,9.6,-2.543,0.1,-2.55E-03,-2.27E-03,-3.54E-03,1.84,-0.05,0.4,0.3004,0.0314,67,270,1350,760,-0.498,-0.785,-0.572,-0.44169,-0.0052,0,0,1,0,0,0.48,0.396,0.565,-0.18,423,0.221,0.093,0.149,0.327,0.068
+PGA,9.907,9.404,9.912,9.6,-2.543,0.1,-2.55E-03,-2.27E-03,-3.54E-03,1.84,-0.05,0.4,0.3004,0.0314,67,270,1350,760,-0.498,-0.785,-0.572,-0.44169,-0.0052,0,0,1,0,0,0.48,0.396,0.565,-0.18,423,0.221,0.093,0.149,0.327,0.068
0.01,9.962,9.451,9.954,9.802,-2.554,0.1,-2.55E-03,-2.19E-03,-4.01E-03,1.84,-0.05,0.4,0.2839,0.0296,67,270,1300,760,-0.498,-0.803,-0.571,-0.4859,-0.0052,0,0,1,0,0,0.476,0.397,0.56,-0.18,423,0.223,0.098,0.148,0.294,0.071
0.02,10.099,9.587,10.086,9.933,-2.566,0.1,-2.55E-03,-2.19E-03,-4.01E-03,1.884,-0.05,0.415,0.2854,0.0298,67,270,1225,760,-0.478,-0.785,-0.575,-0.4859,-0.00518,0,0,1,0,0,0.482,0.401,0.563,-0.181,423,0.227,0.105,0.149,0.294,0.073
0.025,10.181,9.667,10.172,10.009,-2.578,0.1,-2.55E-03,-2.19E-03,-4.01E-03,1.884,-0.05,0.43,0.2891,0.0302,67,270,1200,760,-0.464,-0.745,-0.573,-0.4859,-0.00515,0,0,1,0,0,0.49,0.405,0.575,-0.183,423,0.231,0.12,0.15,0.31,0.076
diff --git a/src/test/java/gov/usgs/earthquake/nshmp/gmm/PSHAB20.java b/src/test/java/gov/usgs/earthquake/nshmp/gmm/PSHAB20.java
index 20816998..7b20051b 100644
--- a/src/test/java/gov/usgs/earthquake/nshmp/gmm/PSHAB20.java
+++ b/src/test/java/gov/usgs/earthquake/nshmp/gmm/PSHAB20.java
@@ -1,6 +1,7 @@
package gov.usgs.earthquake.nshmp.gmm;
import static gov.usgs.earthquake.nshmp.gmm.Gmm.PSHAB20_INTERFACE;
+import static gov.usgs.earthquake.nshmp.gmm.Imt.PGA;
import java.io.IOException;
import java.util.EnumSet;
@@ -33,7 +34,7 @@ class PSHAB20 implements ArgumentsProvider {
/* Result generation sets */
private static Set<Gmm> gmms = EnumSet.of(PSHAB20_INTERFACE);
- private static Set<Imt> imts = PSHAB20_INTERFACE.supportedIMTs();
+ private static Set<Imt> imts = EnumSet.of(PGA);// PSHAB20_INTERFACE.supportedIMTs();
public static void main(String[] args) throws IOException {
GmmTest.generateResults(gmms, imts, GMM_INPUTS, GMM_RESULTS);
diff --git a/src/test/resources/gmm/PSHAB20_results.csv b/src/test/resources/gmm/PSHAB20_results.csv
index d1c97d3b..90c40b89 100644
--- a/src/test/resources/gmm/PSHAB20_results.csv
+++ b/src/test/resources/gmm/PSHAB20_results.csv
@@ -1,23 +1,23 @@
-0-PSHAB20_INTERFACE-SA0P01,0.020897,0.885198
-0-PSHAB20_INTERFACE-SA0P02,0.021238,0.890687
-0-PSHAB20_INTERFACE-SA0P03,0.022458,0.907193
-0-PSHAB20_INTERFACE-SA0P05,0.025828,0.954874
-0-PSHAB20_INTERFACE-SA0P075,0.029636,0.984835
-0-PSHAB20_INTERFACE-SA0P1,0.032537,0.975488
-0-PSHAB20_INTERFACE-SA0P15,0.035565,0.936536
-0-PSHAB20_INTERFACE-SA0P2,0.036517,0.918151
-0-PSHAB20_INTERFACE-SA0P25,0.037127,0.923073
-0-PSHAB20_INTERFACE-SA0P3,0.037353,0.915458
-0-PSHAB20_INTERFACE-SA0P4,0.034904,0.913818
-0-PSHAB20_INTERFACE-SA0P5,0.031459,0.911078
-0-PSHAB20_INTERFACE-SA0P75,0.025667,0.921989
-0-PSHAB20_INTERFACE-SA1P0,0.020612,0.912175
-0-PSHAB20_INTERFACE-SA1P5,0.014567,0.898924
-0-PSHAB20_INTERFACE-SA2P0,0.010897,0.891664
-0-PSHAB20_INTERFACE-SA3P0,0.006269,0.876963
-0-PSHAB20_INTERFACE-SA4P0,0.004365,0.855607
-0-PSHAB20_INTERFACE-SA5P0,0.003441,0.836698
-0-PSHAB20_INTERFACE-SA7P5,0.001955,0.809978
-0-PSHAB20_INTERFACE-SA10P0,0.001252,0.795653
-0-PSHAB20_INTERFACE-PGA,0.020287,0.886792
-0-PSHAB20_INTERFACE-PGV,2.053897,0.857630
+0-PSHAB20_INTERFACE-SA0P01,2.08965604936e-02,0.0
+0-PSHAB20_INTERFACE-SA0P02,2.12376645069e-02,0.0
+0-PSHAB20_INTERFACE-SA0P03,2.24575034172e-02,0.0
+0-PSHAB20_INTERFACE-SA0P05,2.58276475982e-02,0.0
+0-PSHAB20_INTERFACE-SA0P075,2.96364966819e-02,0.0
+0-PSHAB20_INTERFACE-SA0P1,3.25369383278e-02,0.0
+0-PSHAB20_INTERFACE-SA0P15,3.55651022874e-02,0.0
+0-PSHAB20_INTERFACE-SA0P2,3.64296484198e-02,0.0
+0-PSHAB20_INTERFACE-SA0P25,3.69487943052e-02,0.0
+0-PSHAB20_INTERFACE-SA0P3,3.71144622014e-02,0.0
+0-PSHAB20_INTERFACE-SA0P4,3.46117876964e-02,0.0
+0-PSHAB20_INTERFACE-SA0P5,3.11738351662e-02,0.0
+0-PSHAB20_INTERFACE-SA0P75,2.56085017514e-02,0.0
+0-PSHAB20_INTERFACE-SA1P0,2.06123859520e-02,0.0
+0-PSHAB20_INTERFACE-SA1P5,1.45673139070e-02,0.0
+0-PSHAB20_INTERFACE-SA2P0,1.08967064936e-02,0.0
+0-PSHAB20_INTERFACE-SA3P0,6.26946254992e-03,0.0
+0-PSHAB20_INTERFACE-SA4P0,4.36498415871e-03,0.0
+0-PSHAB20_INTERFACE-SA5P0,3.44133448941e-03,0.0
+0-PSHAB20_INTERFACE-SA7P5,1.95537566087e-03,0.0
+0-PSHAB20_INTERFACE-SA10P0,1.25157789513e-03,0.0
+0-PSHAB20_INTERFACE-PGA,2.02870001012e-02,0.0
+0-PSHAB20_INTERFACE-PGV,1.95216532021e+00,0.0
--
GitLab