diff --git a/src/main/java/gov/usgs/earthquake/nshmp/calc/CalcConfig.java b/src/main/java/gov/usgs/earthquake/nshmp/calc/CalcConfig.java index 37fb6e6da6cbc2cdf5897b85ba58548f9b3333ec..126015be8a593b9e5d8305c3aeb0f00b3270a250 100644 --- a/src/main/java/gov/usgs/earthquake/nshmp/calc/CalcConfig.java +++ b/src/main/java/gov/usgs/earthquake/nshmp/calc/CalcConfig.java @@ -542,7 +542,7 @@ public final class CalcConfig { * The return period (in years) at which disaggregations will be performed * unless period-specific ground motion levels are included in a sites file. */ - public final Integer returnPeriod; + public final Double returnPeriod; /** * The distance, magnitude, and epsilon bins into which contributing sources @@ -611,7 +611,7 @@ public final class CalcConfig { private static final class Builder { - Integer returnPeriod; + Double returnPeriod; Bins bins; Double contributorLimit; @@ -793,7 +793,7 @@ public final class CalcConfig { * * <p><b>Default:</b> [475, 975, 2475] */ - public final List<Integer> returnPeriods; + public final List<Double> returnPeriods; private Output(Builder b) { this.directory = b.directory; @@ -805,7 +805,7 @@ public final class CalcConfig { Path directory; Set<DataType> dataTypes; - List<Integer> returnPeriods; + List<Double> returnPeriods; Output build() { checkNotNull(directory); diff --git a/src/main/java/gov/usgs/earthquake/nshmp/calc/Transforms.java b/src/main/java/gov/usgs/earthquake/nshmp/calc/Transforms.java index b8c59a40f7f02d8d68eac768f08203f0dae34fd8..07a8fa3891c09787fd22f33600d0f42e4eba8e0c 100644 --- a/src/main/java/gov/usgs/earthquake/nshmp/calc/Transforms.java +++ b/src/main/java/gov/usgs/earthquake/nshmp/calc/Transforms.java @@ -509,6 +509,10 @@ final class Transforms { * for sources in a cluster. Note that this is only to be used with cluster * sources as the weight of each magnitude variant is stored in the * HazardInput.rate field, which is kinda KLUDGY, but works. + * + * TODO: Investigate sensitivity to collapsing magnitude variants on each + * fault in a cluster [F1(M1,M2,M3), F2(M1 M2,M3), ...] versus doing each + * magnitude branch separately [(F1M1,F2M1,...), (F1M2,F2M2,...)] */ private static final class ClusterGroundMotionsToCurves implements Function<ClusterGroundMotions, ClusterCurves> { diff --git a/src/main/java/gov/usgs/earthquake/nshmp/gmm/MotazedianAtkinson_2005.java b/src/main/java/gov/usgs/earthquake/nshmp/gmm/MotazedianAtkinson_2005.java index 1024432d595a8884f15f83b1b8ddee51dc7b39e9..6cb0affe6c0ffe402c209c6eda88fcec17921de6 100644 --- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/MotazedianAtkinson_2005.java +++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/MotazedianAtkinson_2005.java @@ -242,14 +242,14 @@ public class MotazedianAtkinson_2005 implements GroundMotionModel { // Get μ from MA_05_BASE for current IMT (interpolate if necessary) double μ = GroundMotions.combine(super.calc(in)).mean(); - double μPga = calcMean(coeffsPGA, in.rRup, in.Mw); + double pgaRef = Math.exp(calcMean(coeffsPGA, in.rRup, in.Mw)); double site = siteAmp.isPresent() - ? siteAmp.get().siteAmp(μPga, in.vs30) - : calcInterpolatedSite(super.coeffs.imt, μPga, in.vs30); + ? siteAmp.get().siteAmp(pgaRef, in.vs30) + : calcInterpolatedSite(super.coeffs.imt, pgaRef, in.vs30); return GroundMotions.createTree(μ + site, SIGMA); } - private static double calcInterpolatedSite(Imt imt, double μPga, double vs30) { + private static double calcInterpolatedSite(Imt imt, double pgaRef, double vs30) { Range<Imt> imtRange = INTERPOLATED_SITE_IMTS.get(imt); Imt imtLo = imtRange.lowerEndpoint(); Imt imtHi = imtRange.upperEndpoint(); @@ -263,8 +263,8 @@ public class MotazedianAtkinson_2005 implements GroundMotionModel { double tTarget = imt.period(); BooreAtkinson_2008 ba08lo = (BooreAtkinson_2008) BA_08_BASE.instance(imtLo); BooreAtkinson_2008 ba08hi = (BooreAtkinson_2008) BA_08_BASE.instance(imtHi); - double siteLo = ba08lo.siteAmp(μPga, vs30); - double siteHi = ba08hi.siteAmp(μPga, vs30); + double siteLo = ba08lo.siteAmp(pgaRef, vs30); + double siteHi = ba08hi.siteAmp(pgaRef, vs30); return Interpolator.findY(tLo, siteLo, tHi, siteHi, tTarget); } } diff --git a/src/main/java/gov/usgs/earthquake/nshmp/gmm/UsgsPrviBackbone2025.java b/src/main/java/gov/usgs/earthquake/nshmp/gmm/UsgsPrviBackbone2025.java index 8cb0141e79559e5cb49614e5c9dee001baf9e688..5adb04891e979f3e10612875ddd87dd42cc242bb 100644 --- a/src/main/java/gov/usgs/earthquake/nshmp/gmm/UsgsPrviBackbone2025.java +++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/UsgsPrviBackbone2025.java @@ -82,7 +82,7 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel { final double τ1, τ2, τ3, τ4; final double a, b; - final double φs2sPrvi, φs2sWus; + final double φs2sPrvi, φs2sNga; final double τSub, φSub, φSubSS, φSubS2S; Coefficients( @@ -104,7 +104,7 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel { Map<String, Double> coeffsPhiS2S = ccPhiS2S.get(imt); φs2sPrvi = coeffsPhiS2S.get("phiS2S_PRVI"); - φs2sWus = coeffsPhiS2S.get("phiS2S_WUS"); + φs2sNga = coeffsPhiS2S.get("phiS2S_WUS"); Map<String, Double> coeffsSubSigma = ccSubSigma.get(imt); τSub = coeffsSubSigma.get("tau_mean"); @@ -244,9 +244,9 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel { Coefficients c = super.coeffs; double τ = calcTau(in.Mw, c); double φSS = calcPhiSS(in.Mw, c); - double σWus = sqrt(τ * τ + c.φs2sWus * c.φs2sWus + φSS * φSS); + double σNga = sqrt(τ * τ + c.φs2sNga * c.φs2sNga + φSS * φSS); double σPrvi = sqrt(τ * τ + c.φs2sPrvi * c.φs2sPrvi + φSS * φSS); - return new double[] { σWus, σPrvi }; + return new double[] { σNga, σPrvi }; } // active crust tau @@ -400,9 +400,9 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel { @Override double[] calcSigmas(GmmInput in) { Coefficients c = super.coeffs; - double σ = sqrt(c.τSub * c.τSub + c.φSub * c.φSub); + double σNga = sqrt(c.τSub * c.τSub + c.φSub * c.φSub); double σPrvi = sqrt(c.τSub * c.τSub + c.φSubS2S * c.φSubS2S + c.φSubSS * c.φSubSS); - return new double[] { σ, σPrvi }; + return new double[] { σNga, σPrvi }; } } @@ -535,9 +535,9 @@ public abstract class UsgsPrviBackbone2025 implements GroundMotionModel { @Override double[] calcSigmas(GmmInput in) { Coefficients c = super.coeffs; - double σ = sqrt(c.τSub * c.τSub + c.φSub * c.φSub); + double σNga = sqrt(c.τSub * c.τSub + c.φSub * c.φSub); double σPrvi = sqrt(c.τSub * c.τSub + c.φSubS2S * c.φSubS2S + c.φSubSS * c.φSubSS); - return new double[] { σ, σPrvi }; + return new double[] { σNga, σPrvi }; } } diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/ClusterRuptureSet.java b/src/main/java/gov/usgs/earthquake/nshmp/model/ClusterRuptureSet.java index b17de10cc833448ad910515919f457ebdb07be02..85c9e9ece965efe154344d0309ec936f65fe802c 100644 --- a/src/main/java/gov/usgs/earthquake/nshmp/model/ClusterRuptureSet.java +++ b/src/main/java/gov/usgs/earthquake/nshmp/model/ClusterRuptureSet.java @@ -8,12 +8,16 @@ import static java.util.stream.Collectors.toList; import java.util.Iterator; import java.util.List; +import java.util.Map; +import java.util.TreeMap; import java.util.function.Predicate; +import gov.usgs.earthquake.nshmp.data.XySequence; import gov.usgs.earthquake.nshmp.geo.Location; import gov.usgs.earthquake.nshmp.geo.json.Feature; import gov.usgs.earthquake.nshmp.mfd.Mfd; -import gov.usgs.earthquake.nshmp.tree.LogicGroup; +import gov.usgs.earthquake.nshmp.mfd.Mfd.Type; +import gov.usgs.earthquake.nshmp.tree.Branch; import gov.usgs.earthquake.nshmp.tree.LogicTree; /** @@ -36,13 +40,36 @@ import gov.usgs.earthquake.nshmp.tree.LogicTree; */ public class ClusterRuptureSet extends AbstractRuptureSet<ClusterSource> { + /* + * Developer notes: + * + * A ClusterRuptureSet is composed of two or more FaultRuptureSets (fault + * sections) that rupture together according to some rate or logic tree of + * rates. Each nested FaultRuptureSet itself contains a logic tree of SINGLE + * MFDs, usually labelled 'M1', 'M2', etc. The MFDs for each FaultRuptureSet + * are the same size and have identical weights. Technically, the M1 fault + * branches across each cluster should be considered together, the M2 fault + * branches together, etc. in the cluster hazard calculation. However, the + * magnitude branch exceedance curves are currently collapsed prior to doing + * the cluster calculation. + * + * For the purposes of reporting an MFD tree for a ClusterRuptureSet, we + * create MFDs for each magnitude branch spanning the individual sections in + * the cluster. In some cases there may be identical magnitudes in a cluster + * such that the cluster rate reflected in the MFD may be doubled or more. + * + */ + + @Deprecated final ModelData data; // holds cluster rate tree private final List<ClusterSource> source; + private final LogicTree<Mfd> mfdTree; private ClusterRuptureSet(Builder builder) { super(builder); this.data = builder.data; this.source = List.of(builder.source); + this.mfdTree = builder.mfdTree; } @Override @@ -65,23 +92,7 @@ public class ClusterRuptureSet extends AbstractRuptureSet<ClusterSource> { @Override public LogicTree<Mfd> mfdTree() { - /* - * A ClusterSource is composed of multiple FaultRuptureSets that rupture - * together at some rate. For the purposes of reporting an MFD tree for the - * rupture set, we combine the SINGLE MFDs associated with each - * FaultRuptureSet and set the rate of each magnitude to its weight. We then - * build a logic group where the weight (scale) of each branch is set to the - * cluster rate. - */ - ClusterSource clusterSource = source.get(0); - double rate = clusterSource.rate(); - - LogicGroup.Builder<Mfd> builder = LogicGroup.builder(); - for (RuptureSet<? extends Source> ruptureSet : clusterSource.ruptureSets()) { - Mfd mfd = ModelTrees.reduceMfdTree(ruptureSet.mfdTree()); - builder.addBranch(ruptureSet.name(), mfd, rate); - } - return builder.build(); + return mfdTree; } /* @@ -168,6 +179,7 @@ public class ClusterRuptureSet extends AbstractRuptureSet<ClusterSource> { private ModelData data; // holds cluster rate-tree private ClusterSource source; + private LogicTree<Mfd> mfdTree; private Builder() { super(SourceType.FAULT_CLUSTER); @@ -205,8 +217,37 @@ public class ClusterRuptureSet extends AbstractRuptureSet<ClusterSource> { checkNotNull(source, "%s cluster source", label); } + LogicTree<Mfd> createMfdTree() { + LogicTree.Builder<Mfd> builder = LogicTree.builder("cluster-mfd"); + LogicTree<Mfd> model = source.ruptureSets().get(0).mfdTree(); + /* loop magnitude branches */ + for (int i = 0; i < model.size(); i++) { + /* collect cluster magnitudes and rate; all M have cluster rate */ + Map<Double, Double> magMap = new TreeMap<>(); + String magId = model.get(i).id(); + double magWeight = model.get(i).weight(); + for (RuptureSet<? extends Source> ruptureSet : source.ruptureSets()) { + Branch<Mfd> mfd = ruptureSet.mfdTree().get(i); + /* check magnitude branch IDs and weights match */ + checkState(mfd.id().equals(magId)); + checkState(mfd.value().properties().type() == Type.SINGLE); + magMap.merge( + mfd.value().data().x(0), + source.rate(), + (v1, v2) -> v1 + v2); + } + XySequence mfdXy = XySequence.create( + magMap.keySet(), + magMap.values()); + Mfd clusterMfd = Mfd.create(mfdXy); + builder.addBranch(magId, clusterMfd, magWeight); + } + return builder.build(); + } + ClusterRuptureSet build() { validateAndInit("ClusterRuptureSet.Builder"); + mfdTree = createMfdTree(); return new ClusterRuptureSet(this); } } diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/InterfaceRuptureSet.java b/src/main/java/gov/usgs/earthquake/nshmp/model/InterfaceRuptureSet.java index 8b965edb7d8904623e340f37fd5928f284b133ac..1351e5239af8dd69410de2efa14bc3f49130aacf 100644 --- a/src/main/java/gov/usgs/earthquake/nshmp/model/InterfaceRuptureSet.java +++ b/src/main/java/gov/usgs/earthquake/nshmp/model/InterfaceRuptureSet.java @@ -14,8 +14,6 @@ import java.util.Map; import java.util.function.Predicate; import gov.usgs.earthquake.nshmp.Earthquakes; -import gov.usgs.earthquake.nshmp.fault.surface.ApproxGriddedSurface; -import gov.usgs.earthquake.nshmp.fault.surface.GriddedSurface; import gov.usgs.earthquake.nshmp.geo.Location; import gov.usgs.earthquake.nshmp.geo.LocationList; import gov.usgs.earthquake.nshmp.geo.Locations; @@ -40,8 +38,6 @@ public class InterfaceRuptureSet extends AbstractRuptureSet<InterfaceSource> { final LogicTree<Mfd> mfdTree; final Mfd mfdTotal; - final List<Integer> sectionIds;// reference: actually needed? - final GriddedSurface surface; private final List<InterfaceSource> source; InterfaceRuptureSet(Builder builder) { @@ -52,8 +48,6 @@ public class InterfaceRuptureSet extends AbstractRuptureSet<InterfaceSource> { this.mfdTree = builder.mfdTree; this.mfdTotal = builder.mfdTotal; - this.sectionIds = builder.sectionIds; - this.surface = builder.surface; this.source = List.of(createSource()); } @@ -193,7 +187,6 @@ public class InterfaceRuptureSet extends AbstractRuptureSet<InterfaceSource> { /* created on build */ private LogicTree<Mfd> mfdTree; private Mfd mfdTotal; - private GriddedSurface surface; private Builder() { super(SourceType.INTERFACE); @@ -254,12 +247,6 @@ public class InterfaceRuptureSet extends AbstractRuptureSet<InterfaceSource> { mfdTotal = ModelTrees.reduceMfdTree(mfdTree); feature = createFeature(); - - SourceConfig.Interface config = data.sourceConfig().asInterface(); - surface = new ApproxGriddedSurface( - feature.traces.get(0), - feature.traces.get(1), - config.surfaceSpacing); } private SourceFeature.Interface createFeature() { diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/Models.java b/src/main/java/gov/usgs/earthquake/nshmp/model/Models.java index aa70ac9930da148d1383415fd84cf6475edd2a95..b823cdb5a66ddf1978fe70c6a89b74cdf93f8046 100644 --- a/src/main/java/gov/usgs/earthquake/nshmp/model/Models.java +++ b/src/main/java/gov/usgs/earthquake/nshmp/model/Models.java @@ -17,9 +17,6 @@ import java.util.function.Predicate; import java.util.stream.Collectors; import java.util.stream.IntStream; -import com.google.gson.Gson; -import com.google.gson.GsonBuilder; - import gov.usgs.earthquake.nshmp.data.XySequence; import gov.usgs.earthquake.nshmp.geo.json.Feature; import gov.usgs.earthquake.nshmp.geo.json.FeatureCollection; @@ -37,10 +34,6 @@ import gov.usgs.earthquake.nshmp.tree.LogicTree; */ public class Models { - static final Gson GSON = new GsonBuilder() - .setPrettyPrinting() - .create(); - /** * Returns an object for JSON serialization with the name and ID of all source * logic tree groups in the supplied model organized by tectonic setting and @@ -226,7 +219,7 @@ public class Models { /* * Map source tree into a logic tree of MFDs. Note that the mfdTree - * created below is actually a LogicGroup becuase a SourceTree may contain + * created below is actually a LogicGroup because a SourceTree may contain * null (do-nothing) branches such that weights do not sum to one. */ LogicTree<Mfd> mfdTree = ModelTrees.toMfdTree(tree); diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/SystemRuptureSet.java b/src/main/java/gov/usgs/earthquake/nshmp/model/SystemRuptureSet.java index a9943527037e70d8e5a1bccffe03b125ab92a8e5..43591f0850ee9f7b20970001550a5c55124cfe7e 100644 --- a/src/main/java/gov/usgs/earthquake/nshmp/model/SystemRuptureSet.java +++ b/src/main/java/gov/usgs/earthquake/nshmp/model/SystemRuptureSet.java @@ -84,7 +84,6 @@ public class SystemRuptureSet extends AbstractRuptureSet<SystemRuptureSet.System final Map<Integer, SystemSection> sectionMap; private final Path ruptures; - // from SSSb final GriddedSurface[] sections; final String[] sectionNames; @@ -115,14 +114,16 @@ public class SystemRuptureSet extends AbstractRuptureSet<SystemRuptureSet.System this.sectionMap = builder.sectionMap; this.ruptures = builder.ruptures; - Function<SystemSection, GriddedSurface> surfaceFunction = (type() == SourceType.FAULT_SYSTEM) - ? SystemRuptureSet::featureToCrustalSurface - : SystemRuptureSet::featureToInterfaceSurface; + SourceConfig config = builder.data.sourceConfig(); this.sections = sectionMap.keySet().stream() .sorted() .map(sectionMap::get) - .map(surfaceFunction) + .map(f -> { + return (type() == SourceType.FAULT_SYSTEM) + ? featureToCrustalSurface(f, config.asFault().surfaceSpacing) + : featureToInterfaceSurface(f, config.asInterface().surfaceSpacing); + }) .toArray(GriddedSurface[]::new); this.sectionNames = sectionMap.keySet().stream() @@ -219,7 +220,10 @@ public class SystemRuptureSet extends AbstractRuptureSet<SystemRuptureSet.System return IntStream.rangeClosed(min, max); } - private static GriddedSurface featureToCrustalSurface(SourceFeature.SystemSection section) { + private static GriddedSurface featureToCrustalSurface( + SourceFeature.SystemSection section, + double spacing) { + return DefaultGriddedSurface.builder() .depth(section.upperDepth) .lowerDepth(section.lowerDepth) @@ -227,15 +231,18 @@ public class SystemRuptureSet extends AbstractRuptureSet<SystemRuptureSet.System .dip(section.dip) .dipDir(section.dipDirection) .trace(section.traces.get(0)) - .spacing(1.0) + .spacing(spacing) .build(); } - private static GriddedSurface featureToInterfaceSurface(SourceFeature.SystemSection section) { + private static GriddedSurface featureToInterfaceSurface( + SourceFeature.SystemSection section, + double spacing) { + return new ApproxGriddedSurface( section.traces.get(0), section.traces.get(1), - 5.0); + spacing); } @Override @@ -473,17 +480,8 @@ public class SystemRuptureSet extends AbstractRuptureSet<SystemRuptureSet.System private Map<Integer, SystemSection> sectionMap; private Path ruptures; - /* Unfiltered UCERF3: FM31 = 253,706 FM32 = 305,709 */ - static final int RUP_SET_SIZE = 306000; - static final String ID = "SystemRuptureSet.Builder"; - private List<GriddedSurface> sections; - private List<String> sectionNames; - - private double mMin = Double.POSITIVE_INFINITY; - private double mMax = Double.NEGATIVE_INFINITY; - private Builder(SourceType type) { super(type); } @@ -520,20 +518,6 @@ public class SystemRuptureSet extends AbstractRuptureSet<SystemRuptureSet.System return this; } - Builder sections(List<GriddedSurface> sections) { - checkNotNull(sections, "Section surface list is null"); - checkArgument(sections.size() > 0, "Section surface list is empty"); - this.sections = sections; - return this; - } - - Builder sectionNames(List<String> names) { - checkNotNull(names, "Section name list is null"); - checkArgument(names.size() > 0, "Section name list is empty"); - this.sectionNames = names; - return this; - } - private void validateAndInit(String label) { validateState(label); checkNotNull(data, "%s model data", label); diff --git a/src/test/java/gov/usgs/earthquake/nshmp/calc/CalcConfigTests.java b/src/test/java/gov/usgs/earthquake/nshmp/calc/CalcConfigTests.java index 5324603329921e1ae26f07d50756d12631cf0a5b..53f2b189af0f9dd6f957a543f134eb68b8677a6b 100644 --- a/src/test/java/gov/usgs/earthquake/nshmp/calc/CalcConfigTests.java +++ b/src/test/java/gov/usgs/earthquake/nshmp/calc/CalcConfigTests.java @@ -370,7 +370,7 @@ class CalcConfigTests { @Test void testOutputMember() { - List<Integer> defaultReturnPeroiods = List.of(475, 975, 2475, 10000); + List<Double> defaultReturnPeroiods = List.of(475.0, 975.0, 2475.0, 10000.0); CalcConfig.Output def = DEFAULTS.output; assertEquals(Path.of("hazout"), def.directory); @@ -381,7 +381,7 @@ class CalcConfigTests { assertEquals(Path.of("custom"), def.directory); // also tests automatic addition of TOTAL and MAP assertEquals(Set.of(DataType.TOTAL, DataType.MAP, DataType.GMM), def.dataTypes); - assertEquals(List.of(100, 200), def.returnPeriods); + assertEquals(List.of(100.0, 200.0), def.returnPeriods); def = EXTENDS_EMPTY.output; assertEquals(Path.of("hazout"), def.directory);