Gmm service updates 44

Refactored and cleaned up gmm service classes

Resolves #48 (closed)

src/main/java/gov/usgs/earthquake/nshmp/www/gmm/GmmCalc.java

 package gov.usgs.earthquake.nshmp.www.gmm;
 
+import static java.lang.Math.cos;
+import static java.lang.Math.hypot;
+import static java.lang.Math.sin;
+import static java.lang.Math.tan;
+import static java.lang.Math.toRadians;
+
 import java.util.ArrayList;
+import java.util.Arrays;
 import java.util.EnumMap;
 import java.util.List;
 import java.util.Map;
@@ -8,41 +15,34 @@ import java.util.Set;
 import java.util.stream.Collectors;
 
 import com.google.common.collect.ImmutableList;
+import com.google.common.primitives.Doubles;
 
+import gov.usgs.earthquake.nshmp.Maths;
 import gov.usgs.earthquake.nshmp.gmm.Gmm;
 import gov.usgs.earthquake.nshmp.gmm.GmmInput;
 import gov.usgs.earthquake.nshmp.gmm.GroundMotion;
+import gov.usgs.earthquake.nshmp.gmm.GroundMotionModel;
 import gov.usgs.earthquake.nshmp.gmm.GroundMotions;
 import gov.usgs.earthquake.nshmp.gmm.Imt;
 import gov.usgs.earthquake.nshmp.tree.Branch;
 import gov.usgs.earthquake.nshmp.tree.LogicTree;
+import gov.usgs.earthquake.nshmp.www.gmm.GmmService.Distance;
+import gov.usgs.earthquake.nshmp.www.gmm.GmmService.Magnitude;
+import gov.usgs.earthquake.nshmp.www.gmm.GmmService.Request;
+import jakarta.inject.Singleton;
 
-/**
- * Factory class for computing deterministic response spectra.
+/*
+ * GMM service calculators.
  *
  * @author U.S. Geological Survey
  */
-class ResponseSpectra {
+@Singleton
+class GmmCalc {
 
   private static final double PGA_PERIOD = 0.001;
 
-  /**
-   * Compute the spectra of ground motions and their standard deviations for
-   * multiple models and a source. This method provides the option to compute
-   * ground motion values either for the set of common spectral accelerations
-   * supported by the {@link Gmm}s specified, or for every spectral acceleration
-   * supported by each {@link Gmm}. PGA is included in the results with a
-   * spectral period of 0.001s.
-   *
-   * @param gmms {@code GroundMotionModel}s to use
-   * @param input source and site parameterization
-   * @param commonImts {@code true} if only ground motions corresponding to the
-   *        spectral accelerations common to all {@code gmms} should be
-   *        computed; {@code false} if all spectral accelerations supported by
-   *        each gmm should be used.
-   * @return a {@link Result} data container
-   */
-  static Map<Gmm, GmmData> spectra(Set<Gmm> gmms, GmmInput input, boolean commonImts) {
+  /* Compute ground motion response spectra. */
+  static Map<Gmm, GmmData> spectra(Request request, boolean commonImts) {
 
     /*
      * NOTE: At present, program assumes that all supplied Gmms support PGA.
@@ -51,13 +51,13 @@ class ResponseSpectra {
      */
 
     /* Common imts and periods; may not be used. */
-    Set<Imt> saImts = Gmm.responseSpectrumImts(gmms);
+    Set<Imt> saImts = Gmm.responseSpectrumImts(request.gmms);
     List<Double> periods = new ArrayList<>();
     periods.add(PGA_PERIOD);
     periods.addAll(Imt.periods(saImts));
 
     Map<Gmm, GmmData> gmmSpectra = new EnumMap<>(Gmm.class);
-    for (Gmm gmm : gmms) {
+    for (Gmm gmm : request.gmms) {
       if (!commonImts) {
         saImts = gmm.responseSpectrumImts();
         periods = ImmutableList.<Double> builder()
@@ -67,9 +67,9 @@ class ResponseSpectra {
       }
 
       List<LogicTree<GroundMotion>> imtTrees = new ArrayList<>();
-      imtTrees.add(gmm.instance(Imt.PGA).calc(input));
+      imtTrees.add(gmm.instance(Imt.PGA).calc(request.input));
       for (Imt imt : saImts) {
-        imtTrees.add(gmm.instance(imt).calc(input));
+        imtTrees.add(gmm.instance(imt).calc(request.input));
       }
 
       GmmData dataGroup = treesToDataGroup(periods, imtTrees);
@@ -78,8 +78,8 @@ class ResponseSpectra {
     return gmmSpectra;
   }
 
-  static GmmData treesToDataGroup(
-      List<Double> periods,
+  private static GmmData treesToDataGroup(
+      List<Double> xValues,
       List<LogicTree<GroundMotion>> trees) {
 
     // Can't use Trees.transpose() because some
@@ -90,9 +90,9 @@ class ResponseSpectra {
         .map(Branch::id)
         .collect(Collectors.toList());
 
-    double[] Ts = periods.stream().mapToDouble(Double::valueOf).toArray();
-    double[] μs = new double[periods.size()];
-    double[] σs = new double[periods.size()];
+    double[] xs = xValues.stream().mapToDouble(Double::valueOf).toArray();
+    double[] μs = new double[xValues.size()];
+    double[] σs = new double[xValues.size()];
 
     // build combined μ and σ
     for (int i = 0; i < trees.size(); i++) {
@@ -111,9 +111,9 @@ class ResponseSpectra {
       List<double[]> σBranches = new ArrayList<>();
       List<double[]> weights = new ArrayList<>();
       for (int i = 0; i < modelTree.size(); i++) {
-        μBranches.add(new double[periods.size()]);
-        σBranches.add(new double[periods.size()]);
-        weights.add(new double[periods.size()]);
+        μBranches.add(new double[xValues.size()]);
+        σBranches.add(new double[xValues.size()]);
+        weights.add(new double[xValues.size()]);
       }
 
       // imt indexing
@@ -138,24 +138,122 @@ class ResponseSpectra {
       }
     }
 
-    GmmData gmmData = new GmmData(Ts, μs, σs, epiBranches);
+    GmmData gmmData = new GmmData(xs, μs, σs, epiBranches);
     return gmmData;
   }
 
+  /* Compute ground motions over a range of distances. */
+  static Map<Gmm, GmmData> distance(Distance.Request request, double[] distances) {
+    var inputList = hangingWallDistances(request.input, distances);
+    return calculateGroundMotions(request.gmms, inputList, request.imt, distances);
+  }
+
+  /* Compute ground motions over a range of magnitudes. */
+  static Map<Gmm, GmmData> magnitude(
+      Magnitude.Request request,
+      double[] magnitudes,
+      double distance) {
+
+    var gmmInputs = Arrays.stream(magnitudes)
+        .mapToObj(Mw -> GmmInput.builder().fromCopy(request.input).mag(Mw).build())
+        .map(gmmInput -> hangingWallDistance(gmmInput, distance))
+        .collect(Collectors.toList());
+    return calculateGroundMotions(request.gmms, gmmInputs, request.imt, magnitudes);
+  }
+
+  private static List<GmmInput> hangingWallDistances(
+      GmmInput inputModel,
+      double[] rValues) {
+
+    return Arrays.stream(rValues)
+        .mapToObj(r -> hangingWallDistance(inputModel, r))
+        .collect(Collectors.toList());
+  }
+
+  /* Compute distance metrics for a fault. */
+  private static GmmInput hangingWallDistance(GmmInput in, double r) {
+    /* Dip in radians */
+    double δ = toRadians(in.dip);
+
+    /* Horizontal and vertical widths of fault */
+    double h = cos(δ) * in.width;
+    double v = sin(δ) * in.width;
+
+    /* Depth to bottom of rupture */
+    double zBot = in.zTor + v;
+
+    /* Distance range over which site is normal to fault plane */
+    double rCutLo = tan(δ) * in.zTor;
+    double rCutHi = tan(δ) * zBot + h;
+
+    /* rRup values corresponding to cutoffs above */
+    double rRupLo = Maths.hypot(in.zTor, rCutLo);
+    double rRupHi = Maths.hypot(zBot, rCutHi - h);
+
+    double rJB = (r < 0) ? -r : (r < h) ? 0.0 : r - h;
+    double rRup = (r < rCutLo)
+        ? hypot(r, in.zTor)
+        : (r > rCutHi)
+            ? hypot(r - h, zBot)
+            : rRupScaled(
+                r, rCutLo, rCutHi, rRupLo, rRupHi);
+
+    return GmmInput.builder()
+        .fromCopy(in)
+        .distances(rJB, rRup, r)
+        .build();
+  }
+
+  private static Map<Gmm, GmmData> calculateGroundMotions(
+      Set<Gmm> gmms,
+      List<GmmInput> gmmInputs,
+      Imt imt,
+      double[] distances) {
+
+    Map<Gmm, GmmData> gmValues = new EnumMap<>(Gmm.class);
+
+    for (Gmm gmm : gmms) {
+      List<LogicTree<GroundMotion>> trees = new ArrayList<>();
+      GroundMotionModel model = gmm.instance(imt);
+      for (GmmInput gmmInput : gmmInputs) {
+        trees.add(model.calc(gmmInput));
+      }
+      GmmData dataGroup = GmmCalc.treesToDataGroup(Doubles.asList(distances), trees);
+      gmValues.put(gmm, dataGroup);
+    }
+    return gmValues;
+  }
+
+  /*
+   * Computes rRup for a surface distance r. The range [rCutLo, rCutHi] must
+   * contain r; rRupLo and rRupHi are rRup at rCutLo and rCutHi, respectively.
+   */
+  private static double rRupScaled(
+      double r,
+      double rCutLo,
+      double rCutHi,
+      double rRupLo,
+      double rRupHi) {
+
+    double rRupΔ = rRupHi - rRupLo;
+    double rCutΔ = rCutHi - rCutLo;
+    return rRupLo + (r - rCutLo) / rCutΔ * rRupΔ;
+  }
+
   static class GmmData {
 
-    final double[] periods; // should be any x-value
+    final double[] xs;
     final double[] μs;
     final double[] σs;
     final List<EpiBranch> tree;
 
     GmmData(
-        double[] periods,
+        double[] xs,
         double[] μs,
         double[] σs,
         List<EpiBranch> tree) {
 
-      this.periods = periods;
+      this.xs = xs;
       this.μs = μs;
       this.σs = σs;
       this.tree = tree;