From 791500eb692f69c53be61da4a5aa018730faff37 Mon Sep 17 00:00:00 2001 From: yararossi Date: Tue, 25 Nov 2025 09:29:38 +0100 Subject: [PATCH 1/3] making code run for a single fault. --- src/python/mudpy/fakequakes.py | 5 ++++- src/python/mudpy/forward.py | 12 ++++++++++-- src/python/mudpy/generate_ruptures_parallel.py | 4 +++- src/python/mudpy/hfsims.py | 3 ++- src/python/mudpy/hfsims_parallel.py | 3 ++- src/python/mudpy/parallel.py | 11 +++++++---- 6 files changed, 28 insertions(+), 10 deletions(-) diff --git a/src/python/mudpy/fakequakes.py b/src/python/mudpy/fakequakes.py index bcbf284..ae0238a 100644 --- a/src/python/mudpy/fakequakes.py +++ b/src/python/mudpy/fakequakes.py @@ -114,7 +114,7 @@ def subfault_distances_3D(home,project_name,fault_name,slab_name,projection_zone """ - from numpy import sqrt,sin,cos,deg2rad,zeros,meshgrid,linspace,where,c_,unravel_index,sort,diff,genfromtxt,sign,argmin + from numpy import sqrt,sin,cos,deg2rad,zeros,meshgrid,linspace,where,c_,unravel_index,sort,diff,genfromtxt,sign,argmin,atleast_2d from scipy.interpolate import griddata from matplotlib import pyplot as plt from scipy.spatial.distance import cdist @@ -124,6 +124,9 @@ def subfault_distances_3D(home,project_name,fault_name,slab_name,projection_zone if slab_name==None: #Read fault geometry data fault=genfromtxt(home+project_name+'/data/model_info/'+fault_name) + + # make fault be at least 2d + fault = atleast_2d(fault) #Initalize distance output arrays nsubfaults = len(fault) diff --git a/src/python/mudpy/forward.py b/src/python/mudpy/forward.py index 85153e6..df06c2c 100644 --- a/src/python/mudpy/forward.py +++ b/src/python/mudpy/forward.py @@ -764,6 +764,10 @@ def hf_waveforms(home,project_name,fault_name,rupture_list,GF_list,model_name,ru comp=['N','E','Z'] #components to loop through #Now loop over stations + # TODO: don't loop but run these in parallel + # if there are more stations than faults, go for station ncpu and + # if there are more fault than stations go for fault ncpu + # or if there are more components, then go for that. for ksta in range(len(sta)): #Now loop over components N,E,Z @@ -790,7 +794,7 @@ def make_parallel_hfsims(home,project_name,rupture_name,ncpus,sta,sta_lon,sta_la ''' Set up for MPI calculation of HF stochastics ''' - from numpy import savetxt,arange,genfromtxt,where + from numpy import savetxt,arange,genfromtxt,where,atleast_2d from os import environ import subprocess from shlex import split @@ -798,6 +802,7 @@ def make_parallel_hfsims(home,project_name,rupture_name,ncpus,sta,sta_lon,sta_la #Calculate the necessary full-fault parameters before splitting up the faults over your ncpus rupture=home+project_name+'/output/ruptures/'+rupture_name fault=genfromtxt(rupture) + fault = atleast_2d(fault) slip=(fault[:,8]**2+fault[:,9]**2)**0.5 subfault_M0=slip*fault[:,10]*fault[:,11]*fault[:,13] subfault_M0=subfault_M0*1e7 #to dyne-cm @@ -1267,7 +1272,7 @@ def get_fakequakes_G_and_m(Gimpulse,home,project_name,rupture_name,time_epi,GF_l G: Fully assembled GF matrix ''' - from numpy import genfromtxt,convolve,where,zeros,arange,unique,r_,sort,ones + from numpy import genfromtxt,convolve,where,zeros,arange,unique,r_,sort,ones,atleast_2d from numpy import expand_dims,squeeze,roll,float64,array import dask.array as da @@ -1275,6 +1280,9 @@ def get_fakequakes_G_and_m(Gimpulse,home,project_name,rupture_name,time_epi,GF_l source=genfromtxt(home+project_name+'/forward_models/'+rupture_name) else: source=genfromtxt(home+project_name+'/output/ruptures/'+rupture_name) + + + source = atleast_2d(source) rise_times=source[:,7] fraction=source[:,6] rupture_onset=source[:,12] diff --git a/src/python/mudpy/generate_ruptures_parallel.py b/src/python/mudpy/generate_ruptures_parallel.py index dfa1a02..3f6cef8 100755 --- a/src/python/mudpy/generate_ruptures_parallel.py +++ b/src/python/mudpy/generate_ruptures_parallel.py @@ -18,7 +18,7 @@ def run_parallel_generate_ruptures(home,project_name,run_name,fault_name,slab_na 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,argmin + from numpy import load,save,genfromtxt,log10,cos,sin,deg2rad,savetxt,zeros,where,argmin,atleast_2d from time import gmtime, strftime from numpy.random import shuffle from mudpy import fakequakes @@ -56,6 +56,8 @@ def run_parallel_generate_ruptures(home,project_name,run_name,fault_name,slab_na #Read fault and prepare output variable whole_fault=genfromtxt(home+project_name+'/data/model_info/'+fault_name) + whole_fault = atleast_2d(whole_fault) + #Get structure model vel_mod_file=home+project_name+'/structure/'+model_name diff --git a/src/python/mudpy/hfsims.py b/src/python/mudpy/hfsims.py index 358adba..77c840e 100644 --- a/src/python/mudpy/hfsims.py +++ b/src/python/mudpy/hfsims.py @@ -7,7 +7,7 @@ def stochastic_simulation(home,project_name,rupture_name,sta,sta_lon,sta_lat,com stress parameter is in bars ''' - from numpy import genfromtxt,pi,logspace,log10,mean,where,exp,arange,zeros,argmin,rad2deg,arctan2,real + from numpy import genfromtxt,pi,logspace,log10,mean,where,exp,arange,zeros,argmin,rad2deg,arctan2,real,atleast_2d from pyproj import Geod from obspy.geodetics import kilometer2degrees from obspy.taup import TauPyModel @@ -55,6 +55,7 @@ def stochastic_simulation(home,project_name,rupture_name,sta,sta_lon,sta_lat,com fault=genfromtxt(home+project_name+'/output/ruptures/'+rupture_name) #Onset times for each subfault + fault = atleast_2d(fault) onset_times=fault[:,12] #load velocity structure diff --git a/src/python/mudpy/hfsims_parallel.py b/src/python/mudpy/hfsims_parallel.py index 3a7fbcd..50e7b81 100644 --- a/src/python/mudpy/hfsims_parallel.py +++ b/src/python/mudpy/hfsims_parallel.py @@ -12,7 +12,7 @@ def run_parallel_hfsims(home,project_name,rupture_name,N,M0,sta,sta_lon,sta_lat, stress parameter is in bars ''' - from numpy import genfromtxt,pi,logspace,log10,mean,where,exp,arange,zeros,argmin,rad2deg,arctan2,real,savetxt,c_ + from numpy import genfromtxt,pi,logspace,log10,mean,where,exp,arange,zeros,argmin,rad2deg,arctan2,real,savetxt,c_,atleast_2d from pyproj import Geod from obspy.geodetics import kilometer2degrees from obspy.taup import TauPyModel @@ -81,6 +81,7 @@ def run_parallel_hfsims(home,project_name,rupture_name,N,M0,sta,sta_lon,sta_lat, fault=genfromtxt(mpi_rupt) #Onset times for each subfault + fault = atleast_2d(fault) onset_times=fault[:,12] #load velocity structure diff --git a/src/python/mudpy/parallel.py b/src/python/mudpy/parallel.py index 271d19c..48060ab 100644 --- a/src/python/mudpy/parallel.py +++ b/src/python/mudpy/parallel.py @@ -21,7 +21,7 @@ def run_parallel_green(home,project_name,station_file,model_name,dt,NFFT,static, import subprocess from os import chdir from shutil import copy,rmtree - from numpy import genfromtxt,zeros + from numpy import genfromtxt,zeros,atleast_2d from shlex import split from shutil import copy from glob import glob @@ -48,7 +48,9 @@ def run_parallel_green(home,project_name,station_file,model_name,dt,NFFT,static, print(out) #read your corresponding source file source=genfromtxt(home+project_name+'/data/model_info/mpi_source.'+str(rank)+'.fault') - for ksource in range(len(source)): + source = atleast_2d(source) + + for ksource in range(len(source[:,0])): #Where should I be working boss? depth='%.4f' % source[ksource,3] subfault=str(int(source[ksource,0])).rjust(4,'0') @@ -132,7 +134,7 @@ def run_parallel_synthetics(home,project_name,station_file,model_name,integrate, import subprocess from pandas import DataFrame as df from mudpy.forward import get_mu - from numpy import array,genfromtxt,loadtxt,savetxt,log10,zeros,sin,cos,ones,deg2rad + from numpy import array,genfromtxt,loadtxt,savetxt,log10,zeros,sin,cos,ones,deg2rad,atleast_2d from obspy import read,Stream,Trace from shlex import split from mudpy.green import src2sta,rt2ne,origin_time,okada_synthetics @@ -163,7 +165,8 @@ def run_parallel_synthetics(home,project_name,station_file,model_name,integrate, #Read your corresponding source file mpi_source=genfromtxt(home+project_name+'/data/model_info/mpi_source.'+str(rank)+'.fault') - + mpi_source = atleast_2d(mpi_source) + #Constant parameters rakeDS=90+beta #90 is thrust, -90 is normal rakeSS=0+beta #0 is left lateral, 180 is right lateral From 56973a40bfd218f71a79d1d2a485573a436e094b Mon Sep 17 00:00:00 2001 From: yararossi Date: Tue, 25 Nov 2025 09:30:08 +0100 Subject: [PATCH 2/3] making code run for a single fault. And adjust the ncpus to one if only a single fault. --- src/python/mudpy/runslip.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/python/mudpy/runslip.py b/src/python/mudpy/runslip.py index c168b66..5b0dd6b 100644 --- a/src/python/mudpy/runslip.py +++ b/src/python/mudpy/runslip.py @@ -197,7 +197,7 @@ def make_parallel_green(home,project_name,station_file,fault_name,model_name,dt, OUT: Nothing ''' - from numpy import arange,savetxt,genfromtxt + from numpy import arange,savetxt,genfromtxt,atleast_2d from os import path,makedirs,environ from shlex import split import subprocess @@ -208,6 +208,8 @@ def make_parallel_green(home,project_name,station_file,fault_name,model_name,dt, fault_file=home+project_name+'/data/model_info/'+fault_name #Load source model for station-event distance computations source=genfromtxt(fault_file) + source = atleast_2d(source) + #Create all output folders for k in range(len(source)): strdepth='%.4f' % source[k,3] @@ -228,6 +230,10 @@ def make_parallel_green(home,project_name,station_file,fault_name,model_name,dt, #It doesn't, make it, don't be lazy makedirs(subfault_folder) #Create individual source files + num_faults = len(source) + if num_faults < ncpus: + ncpus=len(source) + print(print('Cutting back to ' + str(ncpus) + ' cpus for ' + str(num_faults) + ' subfaults')) for k in range(ncpus): i=arange(k+hot_start,len(source),ncpus) mpi_source=source[i,:] @@ -401,6 +407,10 @@ def make_parallel_synthetics(home,project_name,station_file,fault_name,model_nam #First read fault model file source=loadtxt(fault_file,ndmin=2) #Create individual source files + num_faults = len(source) + if num_faults < ncpus: + ncpus=len(source) + print(print('Cutting back to ' + str(ncpus) + ' cpus for ' + str(num_faults) + ' subfaults')) for k in range(ncpus): i=arange(k+hot_start,len(source),ncpus) mpi_source=source[i,:] From 64d6c550c1f06cfca4cad96163f7acba825801c7 Mon Sep 17 00:00:00 2001 From: yararossi Date: Tue, 25 Nov 2025 09:30:41 +0100 Subject: [PATCH 3/3] I don't know why this changed, I guess this happenend when I initiated the fortran files. --- src/fk/bessel.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fk/bessel.f b/src/fk/bessel.f index 862bd25..23dd89b 100644 --- a/src/fk/bessel.f +++ b/src/fk/bessel.f @@ -1,7 +1,7 @@ # 1 "bessel.FF" # 1 "" 1 # 1 "" 3 -# 423 "" 3 +# 417 "" 3 # 1 "" 1 # 1 "" 2 # 1 "bessel.FF" 2