Fixed bugs where a user-defined hypocenter was not being used to sele…

…ct the ruptue faults

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
@@ -929,6 +928,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
@@ -938,7 +938,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))
@@ -1278,7 +1277,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
@@ -1305,7 +1304,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)
 
 
 
@@ -1316,7 +1315,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
@@ -1337,7 +1336,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()
@@ -1537,7 +1536,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

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:
@@ -241,9 +250,7 @@ def run_parallel_generate_ruptures(home,project_name,run_name,fault_name,slab_na
                 shear_wave_fraction_shallow=1/60/60/24*2
                 shear_wave_fraction_deep=   1/60/60/24*2
 
-
-
-            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
@@ -395,19 +402,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':
@@ -418,7 +420,7 @@ 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")