Dev jk

.github/workflows/send_update_signal.yml

@@ -0,0 +1,20 @@
+name: Dispatch Event on Master Update
+
+on:
+  push:
+    branches:
+      - master
+
+jobs:
+  dispatch:
+    runs-on: ubuntu-latest    
+    steps:
+      - name: Send repository dispatch event
+        run: |
+          curl -L \
+          -X POST \
+          -H "Accept: application/vnd.github+json" \
+          -H "Authorization: Bearer ${{ secrets.CICD_PAT }}" \
+          -H "X-GitHub-Api-Version: 2022-11-28" \
+          https://api.github.com/repos/UO-Geophysics/on_demand_fakequakes/dispatches \
+          -d '{"event_type":"repo-update","client_payload":{"unit":false,"integration":true}}'

src/fk/bessel.f

@@ -1,10 +1,51 @@
+# 0 "bessel.FF"
+# 0 "<built-in>"
+# 0 "<command-line>"
+# 1 "/usr/include/stdc-predef.h" 1 3 4
+
+# 17 "/usr/include/stdc-predef.h" 3 4
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+# 45 "/usr/include/stdc-predef.h" 3 4
+
+# 55 "/usr/include/stdc-predef.h" 3 4
+
+
+
+
+
+
+
+
+
+# 0 "<command-line>" 2
 # 1 "bessel.FF"
+<<<<<<< HEAD
 # 1 "<built-in>" 1
 # 1 "<built-in>" 3
-# 390 "<built-in>" 3
+# 423 "<built-in>" 3
 # 1 "<command line>" 1
 # 1 "<built-in>" 2
 # 1 "bessel.FF" 2
+=======
+>>>>>>> 3159b946cd1d6ea15965e9cff3ef38de094c9b59
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c bessel.FF:	Compute Bessel function Jn(z) for n=0,1,2
 c Reivsion History
@@ -13,19 +54,10 @@
       subroutine besselFn(z, aj0, aj1, aj2)
 	IMPLICIT NONE
         real z, aj0, aj1, aj2
-
-
-
-
-
-
-
-
+# 17 "bessel.FF"
 	aj0 = BesJ0(z)
 	aj1 = BesJ1(z)
 	aj2 = BesJN(2,z)
-
 # 30 "bessel.FF"
         return
       end
-

src/fk/fk.pl

@@ -53,8 +53,10 @@
     taper applies a low-pass cosine filter at fc=(1-taper)*f_Niquest ($taper).
 -P: specify the min. and max. slownesses in term of 1/vs_at_the_source ($pmin/$pmax)
     and optionally kmax at zero frequency in term of 1/hs ($kmax).
+-G: specify sigma, small imaginary frequency in 1/T ($sigma).
 -R: receiver depth ($r_depth).
 -S: 0=explosion; 1=single force; 2=double couple ($src).
+-T: specify the number of samples before the first arrival ($tb).
 -U: 1=down-going wave only; -1=up-going wave only ($updn).
 -X: dump the input to cmd for debug ($fk).
 Examples
@@ -73,6 +75,8 @@
    my @value = split(/\//,substr($_,2));
    if ($opt eq "D") {
      $deg2km = 6371*3.14159/180.;
+   } elsif ($opt eq "G") {
+     $sigma = $value[0];
    } elsif ($opt eq "H") {
      $f1 = $value[0]; $f2 = $value[1];
    } elsif ($opt eq "M") {
@@ -96,6 +100,8 @@
      $rdep = "_$r_depth";
    } elsif ($opt eq "S") {
      $src = $value[0];
+   } elsif ($opt eq "T") {
+     $tb = $value[0];
    } elsif ($opt eq "U") {
      $updn = $value[0];
    } elsif ($opt eq "X") {

src/fk/model.h

@@ -19,7 +19,7 @@ c		xi--- mu/bulk_modulus.
 c		si(3,6) ---- source coefs. of n=0,1,2.
       integer mb,stype,src,rcv,nlay,ndis,nt,updn
 c max. # of layers and receivers and time length
-      parameter(nlay=20, ndis=5000, nt=4096)
+      parameter(nlay=500, ndis=5000, nt=4096)
       real d(nlay),rho(nlay),mu(nlay),xi(nlay),si(3,6),epsilon
       complex ka(nlay),kb(nlay)
       PARAMETER (epsilon=0.0001)

src/python/mudpy/fakequakes.py

@@ -533,12 +533,11 @@ def select_faults(whole_fault,Dstrike,Ddip,target_Mw,num_modes,scaling_law,
 
 
     #Select random subfault as center of the locus of L and W
-    if subfault_hypocenter==None: 
+    if subfault_hypocenter is None: 
         hypo_fault=randint(0,len(whole_fault)-1)
     else: #get subfault closest to hypo
         hypo_fault=subfault_hypocenter
 
-
     if force_area==True and no_random==False: #Use entire fault model  nothing more to do here folks
 
         selected_faults =  arange(len(whole_fault))
@@ -627,10 +626,10 @@ def select_faults(whole_fault,Dstrike,Ddip,target_Mw,num_modes,scaling_law,
                       width=60+10**normal(0,width_std)
                       length=a/width
 
-    #so which subfault ended up being the middle?
-    center_subfault=hypo_fault
+    # #so which subfault ended up being the middle?
+    # center_subfault=hypo_fault  # I'm not sure why this is getting defined here?
 
-    hypo_fault=int(hypo_fault)
+    # hypo_fault=int(hypo_fault)
 
     #Get max/min distances from hypocenter to all faults
     dstrike_max=Dstrike[:,hypo_fault].max()
@@ -673,7 +672,7 @@ def select_faults(whole_fault,Dstrike,Ddip,target_Mw,num_modes,scaling_law,
     hypo_found=False
 
 
-    if len(selected_faults)>1: #there is mroe than 1 subfault
+    if len(selected_faults)>1: #there is more than 1 subfault
 
         i=randint(0,len(selected_faults)-1)
         hypo_fault=selected_faults[i]
@@ -701,7 +700,7 @@ def select_faults(whole_fault,Dstrike,Ddip,target_Mw,num_modes,scaling_law,
     if use_hypo_fraction == True:
 
         #Need to find "center" of the selected faults. To do this look for the subfault
-        #witht he lowest combined maximum along strike and along dip distance to all
+        #with the lowest combined maximum along strike and along dip distance to all
         #other faults
 
 
@@ -865,7 +864,7 @@ def get_rise_times(M0,slip,fault_array,rise_time_depths,stoc_rake,rise_time='MH2
 
 
 
-def get_rupture_onset(home,project_name,slip,fault_array,model_name,hypocenter,
+def get_rupture_onset(home,project_name,slip,fault_array,model_name,hypocenter,shypo,
         rise_time_depths,M0,velmod,sigma_rise_time=0.2,shear_wave_fraction_shallow=0.49,shear_wave_fraction_deep=0.8):
 
     home,project_name,slip,fault_array,model_name,hypocenter,rise_time_depths,M0,velmod,shear_wave_fraction_shallow,shear_wave_fraction_deep
@@ -883,7 +882,6 @@ def get_rupture_onset(home,project_name,slip,fault_array,model_name,hypocenter,
     # swf_r=randint(1,11)
     # shear_wave_fraction_shallow=1/60/60/24*swf_r
     # shear_wave_fraction_deep=1/60/60/24*swf_r
-    print(shear_wave_fraction_deep*86400)
 
     #I don't condone it but this cleans up the warnings
     warnings.filterwarnings("ignore")
@@ -914,6 +912,7 @@ def get_rupture_onset(home,project_name,slip,fault_array,model_name,hypocenter,
     i_same_as_hypo=where(fault_array[:,3]==hypocenter[2])[0]
     dist=((fault_array[:,1]-hypocenter[0])**2+(fault_array[:,2]-hypocenter[1])**2)**0.5
     i_hypo=argmin(dist)
+
     #Get faults at same depth that are NOT the hypo
     i_same_as_hypo=setxor1d(i_same_as_hypo,i_hypo)
     #perturb
@@ -923,7 +922,6 @@ def get_rupture_onset(home,project_name,slip,fault_array,model_name,hypocenter,
     R=rand(len(i_same_as_hypo))
     fault_array[i_same_as_hypo,3]=fault_array[i_same_as_hypo,3]+delta*R
 
-
     #Loop over all faults
     t_onset=zeros(len(slip))
     length2fault=zeros(len(slip))
@@ -1263,7 +1261,7 @@ def generate_ruptures(home,project_name,run_name,fault_name,slab_name,mesh_name,
 		Lstrike,num_modes,Nrealizations,rake,rise_time,rise_time_depths,time_epi,
 		max_slip,source_time_function,lognormal,slip_standard_deviation,scaling_law,ncpus,
 		force_magnitude=False,force_area=False,mean_slip_name=None,hypocenter=None,
-		slip_tol=1e-2,force_hypocenter=False,no_random=False,shypo=None,use_hypo_fraction=True,
+		slip_tol=1e-2,force_hypocenter=False,no_random=False,use_hypo_fraction=True,
 		shear_wave_fraction_shallow=0.49,shear_wave_fraction_deep=0.8,max_slip_rule=True):
     '''
     Set up rupture generation-- use ncpus if available
@@ -1290,7 +1288,7 @@ def generate_ruptures(home,project_name,run_name,fault_name,slab_name,mesh_name,
     Lstrike,num_modes,Nrealizations,rake,rise_time,rise_time_depths,time_epi,
     max_slip,source_time_function,lognormal,slip_standard_deviation,scaling_law,ncpus,
     force_magnitude,force_area,mean_slip_name,hypocenter,slip_tol,force_hypocenter,
-    no_random,shypo,use_hypo_fraction,shear_wave_fraction_shallow,shear_wave_fraction_deep,max_slip_rule)
+    no_random,use_hypo_fraction,shear_wave_fraction_shallow,shear_wave_fraction_deep,max_slip_rule)
 
 
 
@@ -1301,7 +1299,7 @@ def run_generate_ruptures_parallel(home,project_name,run_name,fault_name,slab_na
         Lstrike,num_modes,Nrealizations,rake,rise_time,rise_time_depths,time_epi,
         max_slip,source_time_function,lognormal,slip_standard_deviation,scaling_law,ncpus,
         force_magnitude,force_area,mean_slip_name,hypocenter,slip_tol,force_hypocenter,
-        no_random,shypo,use_hypo_fraction,shear_wave_fraction_shallow,shear_wave_fraction_deep,max_slip_rule):
+        no_random,use_hypo_fraction,shear_wave_fraction_shallow,shear_wave_fraction_deep,max_slip_rule):
 
     from numpy import ceil
     from os import environ
@@ -1322,7 +1320,7 @@ def run_generate_ruptures_parallel(home,project_name,run_name,fault_name,slab_na
     print("MPI: Starting " + str(Nrealizations_parallel*ncpus) + " FakeQuakes Rupture Generations on ", ncpus, "CPUs")
     mud_source=environ['MUD']+'/src/python/mudpy/'
 
-    mpi='mpiexec -n '+str(ncpus)+' python '+mud_source+'generate_ruptures_parallel.py run_parallel_generate_ruptures '+home+' '+project_name+' '+run_name+' '+fault_name+' '+str(slab_name)+' '+str(mesh_name)+' '+str(load_distances)+' '+distances_name+' '+UTM_zone+' '+str(tMw)+' '+model_name+' '+str(hurst)+' '+Ldip+' '+Lstrike+' '+str(num_modes)+' '+str(Nrealizations_parallel)+' '+str(rake)+' '+str(rise_time)+' '+str(rise_time_depths0)+' '+str(rise_time_depths1)+' '+str(time_epi)+' '+str(max_slip)+' '+source_time_function+' '+str(lognormal)+' '+str(slip_standard_deviation)+' '+scaling_law+' '+str(ncpus)+' '+str(force_magnitude)+' '+str(force_area)+' '+str(mean_slip_name)+' "'+str(hypocenter)+'" '+str(slip_tol)+' '+str(force_hypocenter)+' '+str(no_random)+' '+str(shypo)+' '+str(use_hypo_fraction)+' '+str(shear_wave_fraction_shallow)+' '+str(shear_wave_fraction_deep)+' '+str(max_slip_rule)
+    mpi='mpiexec -n '+str(ncpus)+' python '+mud_source+'generate_ruptures_parallel.py run_parallel_generate_ruptures '+home+' '+project_name+' '+run_name+' '+fault_name+' '+str(slab_name)+' '+str(mesh_name)+' '+str(load_distances)+' '+distances_name+' '+UTM_zone+' '+str(tMw)+' '+model_name+' '+str(hurst)+' '+Ldip+' '+Lstrike+' '+str(num_modes)+' '+str(Nrealizations_parallel)+' '+str(rake)+' '+str(rise_time)+' '+str(rise_time_depths0)+' '+str(rise_time_depths1)+' '+str(time_epi)+' '+str(max_slip)+' '+source_time_function+' '+str(lognormal)+' '+str(slip_standard_deviation)+' '+scaling_law+' '+str(ncpus)+' '+str(force_magnitude)+' '+str(force_area)+' '+str(mean_slip_name)+' "'+str(hypocenter)+'" '+str(slip_tol)+' '+str(force_hypocenter)+' '+str(no_random)+' '+str(use_hypo_fraction)+' '+str(shear_wave_fraction_shallow)+' '+str(shear_wave_fraction_deep)+' '+str(max_slip_rule)
     mpi=split(mpi)
     p=subprocess.Popen(mpi)
     p.communicate()
@@ -1522,7 +1520,7 @@ def run_generate_ruptures(home,project_name,run_name,fault_name,slab_name,mesh_n
             #Calculate rupture onset times
             if force_hypocenter==False: #Use random hypo, otehrwise force hypo to user specified
                 hypocenter=whole_fault[hypo_fault,1:4]
-
+            
             t_onset=get_rupture_onset(home,project_name,slip,fault_array,model_name,hypocenter,
                                       rise_time_depths,M0,velmod,shear_wave_fraction)
             fault_out[:,12]=0
@@ -1563,4 +1561,4 @@ def run_generate_ruptures(home,project_name,run_name,fault_name,slab_name,mesh_n
 
             realization+=1
 
-
+

src/python/mudpy/forward.py

@@ -299,6 +299,8 @@ def waveforms_fakequakes(home,project_name,fault_name,rupture_list,GF_list,
         all_sources=array([rupture_name])  
     else:
         all_sources=genfromtxt(home+project_name+'/data/'+rupture_list,dtype='U')
+        all_sources = array(all_sources, ndmin=1)  # in case only 1 entry
+
 
 
 
@@ -806,7 +808,8 @@ def make_parallel_hfsims(home,project_name,rupture_name,ncpus,sta,sta_lon,sta_la
     for k in range(ncpus):
         i=arange(k,len(fault),ncpus)
         mpi_source=fault[i,:]
-        fmt='%d\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6E'
+        #fmt='%d\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6E' #old
+        fmt='%d\t%10.6f\t%10.6f\t%10.4f\t%10.2f\t%10.2f\t%10.6f\t%10.6E\t%10.6E\t%10.6E\t%10.6f\t%10.6f\t%10.6E\t%10.6E\t%10.6E' #new
         savetxt(home+project_name+'/output/ruptures/mpi_rupt.'+str(k)+'.'+rupture_name,mpi_source,fmt=fmt)
     #Make mpi system call
     print("MPI: Starting Stochastic High Frequency Simulation on ", ncpus, "CPUs")
@@ -2516,7 +2519,7 @@ def get_mu_and_area(home,project_name,fault_name,model_name):
     return mu,area
 
 
-def get_source_time_function(mu,area,rise_time,t0,slip):
+def get_source_time_function(mu,area,rise_time,t0,slip,single_force=0):
     '''
     Compute source time function for a given rise time, right now it assumes 1m of slip
     and a triangle STF
@@ -2528,7 +2531,10 @@ def get_source_time_function(mu,area,rise_time,t0,slip):
     t=linspace(t0,t0+rise_time,1000)
     Mdot=zeros(t.shape)
     #Triangle gradient
-    m=4*mu*area/(rise_time**2)
+    if single_force == 1:
+        m = 4/(rise_time**2)
+    else:
+        m=4*mu*area/(rise_time**2)
     #Upwards intercept
     b1=-m*t0
     #Downwards intercept
@@ -3032,12 +3038,26 @@ def ssds2rake(ss,ds):
     return rake
 
 def makefault(fout,strike,dip,nstrike,dx_dip,dx_strike,epicenter,num_updip,num_downdip,rise_time):
-    '''
-    Make a planar fault
-
-    strike - Strike angle (degs)
-    dip - Dip angle (degs)200/5
-    '''
+    """
+    Create a planar fault and write its properties to a file.
+
+    Parameters:
+
+    fout (str): Output file path to save the fault properties.
+    strike (float): Strike angle in degrees.
+    dip (float): Dip angle in degrees.
+    nstrike (int): Number of subfaults along the strike direction.
+    dx_dip (float): Distance between subfaults along the dip direction.
+    dx_strike (float): Distance between subfaults along the strike direction.
+    epicenter (tuple): Coordinates of the epicenter (longitude, latitude, depth).
+    num_updip (int): Number of subfaults updip from the epicenter.
+    num_downdip (int): Number of subfaults downdip from the epicenter.
+    rise_time (float): Rise time for the fault slip.
+
+    Returns:
+    None
+    """
+
     from numpy import arange,sin,cos,deg2rad,r_,ones,arctan,rad2deg,zeros,isnan,unique,where,argsort
     import pyproj
 
@@ -3125,6 +3145,11 @@ def makefault(fout,strike,dip,nstrike,dx_dip,dx_strike,epicenter,num_updip,num_d
     L=ones(loout.shape)*dx_strike*1000
     W=ones(loout.shape)*dx_dip*1000
     f=open(fout,'w')
+    # Write header line
+    header = "# No\tLongitude\tLatitude\tDepth(km)\tStrike\tDip\ttype\tRise Time\tLength(m)\tWidth(m)\n"
+    print('Hi')
+    f.write(header)
+
     for k in range(len(x)):   
         out='%i\t%.6f\t%.6f\t%.3f\t%.2f\t%.2f\t%.1f\t%.1f\t%.2f\t%.2f\n' % (k+1,loout[k],laout[k],zout[k],strike[k],dip[k],tw[k],rise[k],L[k],W[k])
         f.write(out)