Merge pull request #55 from SCECcode/gp-mrf

GP method updates: Jbsim-3.0.0 and user-selectable merging frequency

Author: fabiolsilva

bbp/comps/anderson_gof.py

@@ -2,7 +2,7 @@
 """
 BSD 3-Clause License
 
-Copyright (c) 2023, University of Southern California
+Copyright (c) 2024, University of Southern California
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
@@ -51,6 +51,16 @@
 from scipy import signal, stats
 from scipy.signal import butter, filtfilt
 
+# Compatible with SciPy 1.4
+try:
+    cumulative_trapezoid = integrate.cumulative_trapezoid
+except AttributeError:
+    cumulative_trapezoid = integrate.cumtrapz
+try:
+    simpson = integrate.simpson
+except AttributeError:
+    simpson = integrate.simps
+
 # Import BBP modules
 import bband_utils
 import plot_config
@@ -135,15 +145,15 @@ def eval_func2(p1, p2):
 
     @staticmethod
     def integral(x, idt):
-        return integrate.simps(x, x=None, dx=idt)
+        return simpson(x, x=None, dx=idt)
 
     @staticmethod
     def integral2(x, idt):
-        return integrate.simps(x * x, x=None, dx=idt)
+        return simpson(x * x, x=None, dx=idt)
 
     @staticmethod
     def integ(array_in, idt):
-        return integrate.cumtrapz(array_in, dx=idt)
+        return cumulative_trapezoid(array_in, dx=idt)
 
     @staticmethod
     def shift_bit_length(x):

bbp/comps/arias_duration.py

@@ -1,7 +1,8 @@
+#!/usr/bin/env python3
 """
 BSD 3-Clause License
 
-Copyright (c) 2021, University of Southern California
+Copyright (c) 2024, University of Southern California
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
@@ -40,6 +41,12 @@
 import math
 import scipy.integrate
 
+# Compatible with SciPy 1.4
+try:
+    cumulative_trapezoid = scipy.integrate.cumulative_trapezoid
+except AttributeError:
+    cumulative_trapezoid = scipy.integrate.cumtrapz
+
 # Converting to cm units. Use approximation to g
 G_TO_CMS = 981.0 # %(cm/s)
 
@@ -97,9 +104,9 @@ def ad_from_acc(a_in_peer_file, a_out_ad):
     # Integrate and calculate vel and disp arrays
     #
     acc_cms = [value * G_TO_CMS for value in acc]
-    velo = list(dt * scipy.integrate.cumtrapz(acc_cms))
+    velo = list(dt * cumulative_trapezoid(acc_cms))
     velo.insert(0, 0)
-    disp = list(dt * scipy.integrate.cumtrapz(velo))
+    disp = list(dt * cumulative_trapezoid(velo))
     disp.insert(0, 0)
 
     #
@@ -138,7 +145,7 @@ def ad_from_acc(a_in_peer_file, a_out_ad):
     # Arias Intensities
     # Using the trapezoidal integration
     arias_intensity = [pow((value * G_TO_CMS), 2) for value in acc]
-    tsum = list(scipy.integrate.cumtrapz(arias_intensity) * dt)
+    tsum = list(cumulative_trapezoid(arias_intensity) * dt)
     tsum.insert(0, 0)
 
     arias_intensity = [num * math.pi / (2 * G_TO_CMS) for num in tsum]

bbp/comps/jbsim.py

@@ -2,7 +2,7 @@
 """
 BSD 3-Clause License
 
-Copyright (c) 2023, University of Southern California
+Copyright (c) 2024, University of Southern California
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
@@ -78,9 +78,18 @@ def run(self):
         a_statfile = os.path.join(install.A_IN_DATA_DIR,
                                   str(sim_id),
                                   self.r_stations)
-        a_srffile = os.path.join(install.A_IN_DATA_DIR,
-                                 str(sim_id),
-                                 self.r_srffile)
+        a_rupture_file = os.path.join(install.A_IN_DATA_DIR,
+                                      str(sim_id),
+                                      self.r_srffile)
+        # RWG 20241025: jbsim-v3.0.0 is backward compatible with jbsim-v2.0.0
+        # it can use rupture files in both SRF and MRF formats
+        if os.path.splitext(a_rupture_file)[-1].lower() == ".srf":
+            rupmod_type = "SRF"
+        elif os.path.splitext(a_rupture_file)[-1].lower() == ".mrf":
+            rupmod_type = "MRF"
+        else:
+            print("[ERROR]: Unknown rupture format in file %s" % (a_rupture_file))
+            sys.exit(-1)
 
         # Set directories, and make sure they exist
         a_indir = os.path.join(install.A_IN_DATA_DIR, str(sim_id))
@@ -106,9 +115,11 @@ def run(self):
             #
             progstring = ("%s latloncoords=1 slon=%f slat=%f " %
                           (os.path.join(install.A_GP_BIN_DIR,
-                                        "jbsim-v2.0.0"), slon, slat) +
+                                        "jbsim-v3.0.0"), slon, slat) +
                           "tshift_timedomain=1 use_closest_gf=1 " +
-                          "rupmodtype=SRF rupmodfile=%s " % a_srffile +
+                          "cgf_flag=%d " % (config.GF_CGFFLAG) +
+                          "rupmodtype=%s " % (rupmod_type) +
+                          "rupmodfile=%s " % (a_rupture_file) +
                           "moment=-1 outdir=%s stat=%s " % (a_veldir, site) +
                           "min_taper_range=0.0 max_taper_range=0.0 " +
                           "gftype=fk gflocs=%s gftimes=%s gf_swap_bytes=%d " %

bbp/comps/jbsim_cfg.py

@@ -1,8 +1,8 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 """
 BSD 3-Clause License
 
-Copyright (c) 2021, University of Southern California
+Copyright (c) 2024, University of Southern California
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
@@ -93,6 +93,11 @@ def __init__(self, vmodel_name):
             self.NTOUT = int(vmodel_params['NTOUT'])
         else:
             self.NTOUT = 4096
+        # RWG 20241025: Added parameter GF_CGFFLAG
+        if 'GF_CGFFLAG' in vmodel_params:
+            self.GF_CGFFLAG = int(vmodel_params['GF_CGFFLAG'])
+        else:
+            self.GF_CGFFLAG = 0
 
         # Now, configure other paramteres
         self.GF_SWAP_BYTES = 0

bbp/comps/match.py

@@ -2,7 +2,7 @@
 """
 BSD 3-Clause License
 
-Copyright (c) 2023, University of Southern California
+Copyright (c) 2024, University of Southern California
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
@@ -87,11 +87,25 @@ def run(self):
         dirs = [a_tmpdir]
         bband_utils.mkdirs(dirs, print_cmd=False)
 
+        # Get pointer to the velocity model object
+        vel_obj = velocity_models.get_velocity_model_by_name(self.vmodel_name)
+        if vel_obj is None:
+            raise bband_utils.ParameterError("Cannot find velocity model: %s" %
+                                             (self.vmodel_name))
+
+        # Check for velocity model-specific parameters
+        vmodel_params = vel_obj.get_codebase_params('gp')
+
         # Start with defaults
         self.phase = config.PHASE
         self.hf_fhi = config.HF_FHI
         self.lf_flo = config.LF_FLO
 
+        # Check if we have a different merging frequency
+        if 'MATCH_MERGING_FREQUENCY' in vmodel_params:
+            self.hf_fhi = float(vmodel_params['MATCH_MERGING_FREQUENCY'])
+            self.lf_flo = float(vmodel_params['MATCH_MERGING_FREQUENCY'])
+
         # Set match method
         if config.MATCH_METHOD == 1:
             self.phase = 1
@@ -109,15 +123,6 @@ def run(self):
         slo = StationList(a_statfile)
         site_list = slo.get_station_list()
 
-        # Get pointer to the velocity model object
-        vel_obj = velocity_models.get_velocity_model_by_name(self.vmodel_name)
-        if vel_obj is None:
-            raise bband_utils.ParameterError("Cannot find velocity model: %s" %
-                                             (self.vmodel_name))
-
-        # Check for velocity model-specific parameters
-        vmodel_params = vel_obj.get_codebase_params('gp')
-
         # Figure out what DT we should use when resampling
 
         # Figure out the LF DT value

bbp/comps/match_cfg.py

@@ -1,8 +1,8 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 """
 BSD 3-Clause License
 
-Copyright (c) 2021, University of Southern California
+Copyright (c) 2024, University of Southern California
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
@@ -45,13 +45,16 @@ class MatchCfg(object):
     Define the configuration parameters for the Match program
     """
     def __init__(self):
-        self.HF_FHI = 1.0
+        # Default merging frequency set
+        self.DEFAULT_MERGING_FREQUENCY = 1.0
+
+        self.HF_FHI = self.DEFAULT_MERGING_FREQUENCY
         self.HF_FLO = "1.0e+15"
         self.HF_ORD = 4
         self.HF_TSTART = 0.0
 
         self.LF_FHI = 0.0
-        self.LF_FLO = 1.0
+        self.LF_FLO = self.DEFAULT_MERGING_FREQUENCY
         self.LF_ORD = 4
         self.LF_TSTART = 0.0
 

bbp/comps/rzz2015.py

@@ -2,7 +2,7 @@
 """
 BSD 3-Clause License
 
-Copyright (c) 2023, University of Southern California
+Copyright (c) 2024, University of Southern California
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
@@ -45,6 +45,12 @@
     mpl.use('Agg') # Disable use of Tk/X11
 import pylab
 
+# Compatible with SciPy 1.4
+try:
+    cumulative_trapezoid = integrate.cumulative_trapezoid
+except AttributeError:
+    cumulative_trapezoid = integrate.cumtrapz
+
 # Import BBP modules
 import bband_utils
 import plot_config
@@ -106,7 +112,7 @@ def arias(self, F, dt, percent):
         n = len(F)
 
         a_i = [pow(value, 2) for value in F]
-        I = integrate.cumtrapz(a_i) * dt
+        I = cumulative_trapezoid(a_i) * dt
         # Arias Intensity
         Ia = (F[0]**2) * dt / 2.0 + I[n-2] + (F[n-1]**2) * dt / 2.0
         It = (percent / 100.0) * Ia

bbp/src/gp/JordanBailey/.gitignore

@@ -2,4 +2,5 @@
 jbsim
 jbsim3d
 ray_stimes
-jbsim-v2.0.0
+jbsim-v2.0.0
+jbsim-v3.0.0

bbp/src/gp/JordanBailey/defs.h

@@ -11,6 +11,7 @@
 #define   U1FLAG   1
 #define   U2FLAG   2
 #define   U3FLAG   3
+#define   MOMTEN_FLAG   4
 
 #define RDONLY_FLAGS    O_RDONLY
 #define RDWR_FLAGS      O_RDWR

bbp/src/gp/JordanBailey/function.h

@@ -36,6 +36,8 @@ void sum_4gf_v2(float *,int,float *,struct gfheader *,int,int,float *,float *,fl
 void timeshift_gf_v2(float *,int,float *,struct gfheader *,int,float *,float *,int);
 void mech_4gf_v2(float *,float *,struct gfheader *,struct gfparam gfp,int,struct mechparam,float *,float *);
 
+void momten_4gf(struct standrupformat *,int,int,float *,float *,struct gfheader *,struct gfparam,int,float *);
+
 void read_beroza(struct beroza *,char *rfile,float *);
 void get_brmpars(struct beroza *,int,int,float *,float *,float *,float *);
 void beroza_stf(struct beroza *,int,int,float *,float *,float *,int,float *,float *);
@@ -57,6 +59,9 @@ void get_srfpars(struct standrupformat *,int,int,float *,float *,float *,float *
 void get_srfpars_v2(struct standrupformat *,int,int,float *,float *,struct mechparam *);
 void get_srfparsOLD(struct standrupformat *,int,int,float *,float *,float *,int,int,struct mechparam *);
 void srf_stf(struct standrupformat *,int,int,float *,float *,float *,int,float *,struct mechparam,float *);
+void srf_stf_v2(struct standrupformat *,int,int,float *,float *,float *,int,float *,struct mechparam,float *);
+
+void get_indx_pow(float *,float *,float *,struct sgtindex *,struct geoprojection *,int);
 
 char *skipval(int,char *);
 void do_cnvlv(float *,float *,int,float *,int);