Merge pull request #90 from taranye96/master

Fixed a couple formatting errors plus a weird taupy bug

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

@@ -808,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")

src/python/mudpy/generate_ruptures_parallel.py

@@ -11,14 +11,14 @@ def run_parallel_generate_ruptures(home,project_name,run_name,fault_name,slab_na
         num_modes,Nrealizations,rake,rise_time,rise_time_depths0,rise_time_depths1,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,
+        no_random,use_hypo_fraction,shear_wave_fraction_shallow,shear_wave_fraction_deep,
         max_slip_rule,rank,size):
 
     '''
     Depending on user selected flags parse the work out to different functions
     '''
 
-    from numpy import load,save,genfromtxt,log10,cos,sin,deg2rad,savetxt,zeros,where
+    from numpy import load,save,genfromtxt,log10,cos,sin,deg2rad,savetxt,zeros,where,argmin
     from time import gmtime, strftime
     from numpy.random import shuffle
     from mudpy import fakequakes
@@ -38,12 +38,11 @@ def run_parallel_generate_ruptures(home,project_name,run_name,fault_name,slab_na
     else:
         time_epi=UTCDateTime(time_epi)
     rise_time_depths=[rise_time_depths0,rise_time_depths1]
-    #hypocenter=[hypocenter_lon,hypocenter_lat,hypocenter_dep]
     tMw=tMw.split(',')
     target_Mw=zeros(len(tMw))
     for rMw in range(len(tMw)):
         target_Mw[rMw]=float(tMw[rMw])
-
+        
 
     #Should I calculate or load the distances?
     if load_distances==1:  
@@ -63,6 +62,18 @@ def run_parallel_generate_ruptures(home,project_name,run_name,fault_name,slab_na
 
     #Get TauPyModel
     velmod = TauPyModel(model=home+project_name+'/structure/'+model_name.split('.')[0])
+
+    # Define the subfault hypocenter (if hypocenter is prescribed)
+    if hypocenter is None:
+        shypo=None
+    else:
+        dist=((whole_fault[:,1]-hypocenter[0])**2+(whole_fault[:,2]-hypocenter[1])**2)**0.5
+        shypo=argmin(dist)
+
+        # Re-define hypocenter as the coordinates of the hypocenter subfault in
+        # case the original hypocenter did not perfectly align with a subfault
+        hypocenter = whole_fault[shypo,1:4]
+
 
     #Now loop over the number of realizations
     realization=0
@@ -88,7 +99,6 @@ def run_parallel_generate_ruptures(home,project_name,run_name,fault_name,slab_na
                 current_target_Mw=target_Mw[kmag]
                 ifaults,hypo_fault,Lmax,Wmax,Leff,Weff,option,Lmean,Wmean=fakequakes.select_faults(whole_fault,Dstrike,Ddip,current_target_Mw,num_modes,scaling_law,
                                     force_area,no_shallow_epi=False,no_random=no_random,subfault_hypocenter=shypo,use_hypo_fraction=use_hypo_fraction)
-                print(option)
                 fault_array=whole_fault[ifaults,:]
                 Dstrike_selected=Dstrike[ifaults,:][:,ifaults]
                 Ddip_selected=Ddip[ifaults,:][:,ifaults]
@@ -224,8 +234,7 @@ def run_parallel_generate_ruptures(home,project_name,run_name,fault_name,slab_na
             #Calculate rupture onset times
             if force_hypocenter==False: #Use random hypo, otehrwise force hypo to user specified
                 hypocenter=whole_fault[hypo_fault,1:4]
-            else:
-                hypocenter=whole_fault[shypo,1:4]
+
 
             # edit ...
             # if rise_time==2:
@@ -243,9 +252,7 @@ def run_parallel_generate_ruptures(home,project_name,run_name,fault_name,slab_na
             else: #regular EQs, do nothing
                 pass
 
-
-
-            t_onset,length2fault=fakequakes.get_rupture_onset(home,project_name,slip,fault_array,model_name,hypocenter,rise_time_depths,
+            t_onset,length2fault=fakequakes.get_rupture_onset(home,project_name,slip,fault_array,model_name,hypocenter,shypo,rise_time_depths,
                                                  M0,velmod,shear_wave_fraction_shallow=shear_wave_fraction_shallow,
                                                  shear_wave_fraction_deep=shear_wave_fraction_deep)
             fault_out[:,12]=0
@@ -260,7 +267,7 @@ def run_parallel_generate_ruptures(home,project_name,run_name,fault_name,slab_na
 
             #Calculate average risetime
             rise = fault_out[:,7]
-            avg_rise = np.mean(rise)
+            avg_rise = np.mean(rise[np.where(rise>0)[0]])
 
             # Calculate average rupture velocity
             lon_array = fault_out[:,1]
@@ -397,19 +404,14 @@ def run_parallel_generate_ruptures(home,project_name,run_name,fault_name,slab_na
             no_random=True
         elif no_random=='False':
             no_random=False
-        shypo=sys.argv[36]
-        if shypo=='None':
-            shypo=None
-        else:
-            shypo=int(shypo)
-        use_hypo_fraction=sys.argv[37]
+        use_hypo_fraction=sys.argv[36]
         if use_hypo_fraction=='True':
             use_hypo_fraction=True
         if use_hypo_fraction=='False':
             use_hypo_fraction=False
-        shear_wave_fraction_shallow=float(sys.argv[38])
-        shear_wave_fraction_deep=float(sys.argv[39])
-        max_slip_rule=sys.argv[40]
+        shear_wave_fraction_shallow=float(sys.argv[37])
+        shear_wave_fraction_deep=float(sys.argv[38])
+        max_slip_rule=sys.argv[39]
         if max_slip_rule=='True':
             max_slip_rule=True
         if max_slip_rule=='False':
@@ -420,8 +422,8 @@ def run_parallel_generate_ruptures(home,project_name,run_name,fault_name,slab_na
         num_modes,Nrealizations,rake,rise_time,rise_time_depths0,rise_time_depths1,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,
+        no_random,use_hypo_fraction,shear_wave_fraction_shallow,shear_wave_fraction_deep,
         max_slip_rule,rank,size)
     else:
         print("ERROR: You're not allowed to run "+sys.argv[1]+" from the shell or it does not exist")
-
+

src/python/mudpy/hfsims_parallel.py

@@ -269,7 +269,11 @@ def run_parallel_hfsims(home,project_name,rupture_name,N,M0,sta,sta_lon,sta_lat,
             rake=rad2deg(arctan2(ds,ss))
 
             #Get ray paths for all direct P arrivals
-            Ppaths=velmod.get_ray_paths(zs,dist_in_degs,phase_list=['P','p'])
+            try:
+                Ppaths=velmod.get_ray_paths(zs,dist_in_degs,phase_list=['P','p'])
+            except:
+                zs = zs+0.0001
+                Ppaths=velmod.get_ray_paths(zs,dist_in_degs,phase_list=['P','p'])
 
             #Get ray paths for all direct S arrivals
             try:
@@ -647,7 +651,7 @@ def run_parallel_hfsims(home,project_name,rupture_name,N,M0,sta,sta_lon,sta_lat,
         Swave=sys.argv[21]
         if Swave=='True':
             Swave=True
-        high_stress_depth=int(sys.argv[22])
+        high_stress_depth=float(sys.argv[22])
         Qmethod=sys.argv[23]
         scattering=sys.argv[24]
         Qc_exp=float(sys.argv[25])