diff --git a/data/tools/main.py b/data/tools/main.py
index 31aabba6a8716c47ef042cf2e834bc59ea7024b6..332d2747e5e047d77d2ed747209bb5a142c9153d 100644
--- a/data/tools/main.py
+++ b/data/tools/main.py
@@ -24,6 +24,8 @@ from support.iterations import iterations
from support.friction_stats import friction_stats
from support.slip_rates import slip_rates
+from support.find_quakes import find_quakes
+
def build_patch(coords, percentage):
x_coords = coords[:, 0]
xmin = np.min(x_coords)
@@ -36,16 +38,17 @@ def build_patch(coords, percentage):
return [i for i in range(len(x_coords)) if x_coords[i]>=xmin and x_coords[i]<=xmax]
# read problem parameters
-params = read_params('multi-body-problem.cfg')
+params = read_params('strikeslip.cfg')
TANGENTIAL_COORDS = 1
FINAL_TIME = params['finalTime']
NBODIES = params['bodyCount']
+THRESHOLD_VELOCITY = 1e-3
h5file = read_h5file()
print(list(h5file.keys()))
-interval = [0, FINAL_TIME] #[0.75*FINAL_TIME, FINAL_TIME]
+interval = [0.5*FINAL_TIME, FINAL_TIME]
iterations(h5file, FINAL_TIME, interval)
@@ -58,9 +61,12 @@ for body_ID in range(NBODIES):
coords = np.array(h5file[body + '/coordinates'])
patch = build_patch(coords, 1.0)
+ print("patch length: " + str(len(patch)))
+
friction_stats(h5file, body_ID, FINAL_TIME, patch, interval)
- slip_rates(h5file, body_ID, FINAL_TIME, patch, interval, TANGENTIAL_COORDS)
+ #slip_rates(h5file, body_ID, FINAL_TIME, patch, interval, TANGENTIAL_COORDS)
+ [quake_starts, quake_ends] = find_quakes(h5file, body_ID, FINAL_TIME, patch, interval, THRESHOLD_VELOCITY, TANGENTIAL_COORDS)
plt.show()
h5file.close()
\ No newline at end of file
diff --git a/data/tools/support/find_quakes.py b/data/tools/support/find_quakes.py
index aa4653fd4384a89db8be9a258c882dc7583c1495..81e57ce259326255357f6b82cc35c02f83db661d 100644
--- a/data/tools/support/find_quakes.py
+++ b/data/tools/support/find_quakes.py
@@ -1,16 +1,62 @@
-from support.slip_beginnings import slip_beginnings
-from support.slip_endings import slip_endings
+import array as ar
+import numpy as np
+from .slip_beginnings import slip_beginnings
+from .slip_endings import slip_endings
-def find_quakes(threshold_velocity, maximum_velocities):
- slipping_times = maximum_velocities > threshold_velocity
+def peak_slip(quake_start, quake_end, velocities):
+ quake_v = velocities[quake_start-1:quake_end]
+ max_v = np.max(quake_v, axis=1)
+ return [quake_start+np.argmax(max_v), np.max(max_v)]
+
+#def rupture_width
+
+def find_quakes(h5file, body_ID, FINAL_TIME, patch, interval, threshold_velocity, tangential_coords):
+ body = 'body' + str(body_ID) # 'frictionalBoundary' 'body' + str(body_ID)
+
+ coords = np.array(h5file[body + '/coordinates'])
+ if len(patch) == 0:
+ patch = np.linspace(0, len(coords)-1, len(coords), endpoint=True, dtype='int8')
+
+ # read time
+ time = np.array(h5file['relativeTime']) * FINAL_TIME
+ time = np.delete(time, 0)
+ if len(interval) == 0:
+ interval = [0, FINAL_TIME]
+ t = [i for i in range(len(time)) if time[i]>=interval[0] and time[i]<=interval[1]]
+ #t = t[::5]
+ time = time[t]
+
+ # velocity data
+ v = abs(np.array(h5file[body + '/velocity']))
+ v_t = v[t,:,tangential_coords]
+ v_tx = v_t[:,patch]
+ avg_v = np.average(v_tx, axis=1)
+
+ slipping_times = avg_v > threshold_velocity
quake_starts = slip_beginnings(slipping_times)
quake_ends = slip_endings(slipping_times)
+ #print(quake_starts)
+ #print(quake_ends)
+
# remove incomplete quakes
min_len = min(len(quake_starts), len(quake_ends))
quake_ends = quake_ends[0:min_len]
quake_starts = quake_starts[0:min_len]
+ print("Number of quakes: " + str(min_len))
+
+ peak_v = np.zeros(min_len)
+ quake_times = np.zeros(min_len)
+ for i in range(min_len):
+ [quake_times[i], peak_v[i]] = peak_slip(quake_starts[i], quake_ends[i], v_tx)
+ print("peak slip velocity: " + str(peak_v.mean()) + " +- " + str(peak_v.std()))
+
+ recurrence_time = np.zeros(min_len-1)
+ for i in range(min_len-1):
+ recurrence_time[i] = time[int(quake_times[i+1])] - time[int(quake_times[i])]
+ print("recurrence time: " + str(recurrence_time.mean()) + " +- " + str(recurrence_time.std()))
+
return [quake_starts, quake_ends]
diff --git a/data/tools/support/friction_stats.py b/data/tools/support/friction_stats.py
index 7f0a13c48a69df130df7731b3f20957d12fd914d..1bff3d341bc249b9581239a05713d9c8eee10e65 100644
--- a/data/tools/support/friction_stats.py
+++ b/data/tools/support/friction_stats.py
@@ -14,6 +14,7 @@ def friction_stats(h5file, body_ID, FINAL_TIME, patch=[], interval=[], TANGENTIA
if len(interval) == 0:
interval = [0, FINAL_TIME]
t = [i for i in range(len(time)) if time[i]>=interval[0] and time[i]<=interval[1]]
+ #t = t[::5]
time = time[t]
fig = plt.figure()
@@ -32,6 +33,7 @@ def friction_stats(h5file, body_ID, FINAL_TIME, patch=[], interval=[], TANGENTIA
ax_slip.plot(time, avg_v, color='black', linestyle='-')
#ax_slip.plot(time, max_v, color='gray', linestyle='--')
ax_slip.set_ylabel('slip rate V [m/s]')
+ ax_slip.set_xlabel('time t [s]')
ax_slip.set_yscale('log')
#ax_slip.set_ylim([1e-6,1e-2])
#-------------------------
diff --git a/data/tools/support/iterations.py b/data/tools/support/iterations.py
index adf0550ad00ab743cd3d0683bf72cc05531c4c87..4df9adc656778b3d12721071ae626eff5b55b5aa 100644
--- a/data/tools/support/iterations.py
+++ b/data/tools/support/iterations.py
@@ -19,6 +19,9 @@ def iterations(h5file, FINAL_TIME, interval = []):
fpi_final = np.delete(fpi_final, 0)
#fpi_total = np.array(h5file['iterations/fixedPoint/total'])
+ print("FPI average: " + str(np.average(fpi_final)))
+ print("FPI max: " + str(np.max(fpi_final)))
+
# read multigrid iterations
multigrid_final = np.array(h5file['iterations/multiGrid/final'])
multigrid_final = np.delete(multigrid_final, 0)
@@ -38,9 +41,13 @@ def iterations(h5file, FINAL_TIME, interval = []):
ax_mg.set_xlabel('time [s]')
#-------------------------
+ tau_t = tau[t]
ax_tau = fig.add_subplot(3, 1, 3)
- ax_tau.plot(time, tau[t], color='black', linestyle='-')
+ ax_tau.plot(time, tau_t, color='black', linestyle='-')
ax_tau.set_ylabel('tau')
+ ax_tau.set_yscale('log')
+
+ print("tau_max / tau_min: " + str(np.max(tau_t)/np.min(tau_t)))
#-------------------------
fig.canvas.draw()
\ No newline at end of file
diff --git a/data/tools/support/slip_beginnings.py b/data/tools/support/slip_beginnings.py
index 255adf44d900595def438864124830b7134306c9..dff7ae015a0093cde12ba450178321f22ac65a34 100644
--- a/data/tools/support/slip_beginnings.py
+++ b/data/tools/support/slip_beginnings.py
@@ -1,6 +1,5 @@
import array as ar
-
def slip_beginnings(x):
# returns indicies i for which x[i-1]=0 and x[i]>0
starts = ar.array('i', (0 for i in range(len(x))))
diff --git a/data/tools/support/slip_endings.py b/data/tools/support/slip_endings.py
index 9fdcaf88a6f8f7942a580489a33bf60f6ce1271a..fe42e5373b58e0f64272d277d44812530349f7c5 100644
--- a/data/tools/support/slip_endings.py
+++ b/data/tools/support/slip_endings.py
@@ -1,6 +1,5 @@
import array as ar
-
def slip_endings(x):
# returns indicies i for which x[i-1]>0 and x[i]=0
ends = ar.array('i', (0 for i in range(len(x))))
@@ -12,5 +11,4 @@ def slip_endings(x):
ends[length] = i
length += 1
prev = x[i]
-
- return ends[:length]
+ return ends[:length]
\ No newline at end of file