diff --git a/data/generated/events1.csv b/data/generated/events1.csv
new file mode 100644
index 0000000000000000000000000000000000000000..39ddcc413264bddc2f7ba692707a059b8650245a
--- /dev/null
+++ b/data/generated/events1.csv
@@ -0,0 +1,10 @@
+start,end,peakSlipRate,ruptureWidth
+26.75999999999787,26.784999999997865,0.027083990074991897,0.7066604616124774
+28.98999999999762,29.049999999997617,0.03226818570782891,0.7066460434893268
+31.56499999999734,31.614999999997334,0.042600764506162535,0.7065953907642707
+34.18999999999705,34.24499999999704,0.039804029857929246,0.7065001704582403
+36.69999999999678,36.75499999999677,0.04534969591160434,0.7063800319864028
+39.32999999999648,39.39499999999648,0.038848015455120025,0.7062439809396165
+42.71999999999611,42.7699999999961,0.04713727822866609,0.7063968867127027
+45.77499999999578,45.83999999999577,0.04896737991925046,0.7064612510276674
+48.69499999999545,48.75499999999545,0.04528040064825157,0.7062545440454634
diff --git a/data/tools/2d_event_writer.py b/data/tools/2d_event_writer.py
index f79edf660786f0443b772ec6cdbbd2a52262e533..6d4e34f6498094abb511042b4de594ebd211eb2f 100644
--- a/data/tools/2d_event_writer.py
+++ b/data/tools/2d_event_writer.py
@@ -4,6 +4,9 @@ import numpy as np
import csv
import h5py
+from support.io import read_h5file
+from support.io import read_params
+
from support.maximum import maximum
from support.norm import norm
from support.find_quakes import find_quakes
@@ -11,20 +14,10 @@ from support.slip_beginnings import slip_beginnings
from support.slip_endings import slip_endings
from support.max_distance import max_distance
-
-NBODIES = 2
-FINAL_TIME = 1000 # s
-FINAL_VELOCITY = 1e-5 # m/s
-THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY # 1000e-6 + FINAL_VELOCITY
-
-TANGENTIAL_COORDS = 1
-
-# read config ini
+# read outpath from config ini
config = cp.ConfigParser()
-config_path = os.path.join('config.ini')
+config_path = os.path.join('tools/config.ini')
config.read(config_path)
-sim_path = config.get('directories', 'simulation')
-exp_path = config.get('directories', 'experiment')
out_path = config.get('directories', 'output')
# create output directory
@@ -32,12 +25,27 @@ out_dir = os.path.join(out_path)
if not os.path.exists(out_dir):
os.mkdir(out_dir)
-h5path = os.path.join(sim_path)
-h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')
+# read problem parameters
+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.5*FINAL_TIME, FINAL_TIME]
# read time
-relative_time = np.array(h5file['relativeTime'])
-real_time = relative_time * FINAL_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]
for body_ID in range(NBODIES):
body = 'body' + str(body_ID)
@@ -48,14 +56,19 @@ for body_ID in range(NBODIES):
# read data
coordinates = np.array(h5file[body + '/coordinates'])
displacement = np.array(h5file[body + '/displacement'])
+
velocity = np.array(h5file[body + '/velocity'])
+ #print(velocity[t,:,TANGENTIAL_COORDS:])
+ rate = np.linalg.norm(velocity[t,:,TANGENTIAL_COORDS:], axis=2)
+ #print(rate)
+ avg_rate = np.average(rate, axis=1)
+ max_rate = np.max(rate, axis=1)
+ print(avg_rate)
num_vertices = displacement.shape[1]
- velocity_norm = norm(velocity[:, :, TANGENTIAL_COORDS:])
- maximum_velocity = maximum(velocity_norm)
- [quake_starts, quake_ends] = find_quakes(THRESHOLD_VELOCITY,
- maximum_velocity)
+ [quake_starts, quake_ends] = find_quakes(avg_rate, THRESHOLD_VELOCITY)
+ print("Number of quakes: " + str(len(quake_starts)))
quakes = {}
for quake_ID in range(len(quake_starts)):
@@ -63,23 +76,23 @@ for body_ID in range(NBODIES):
quake_end = int(quake_ends[quake_ID])
quake = {}
- quake['start'] = real_time[quake_start]
- quake['end'] = real_time[quake_end]
+ quake['start'] = time[quake_start]
+ quake['end'] = time[quake_end]
- local_slipping_times = velocity_norm[quake_start:quake_end, :] \
+ local_slipping_times = rate[quake_start:quake_end, :] \
> THRESHOLD_VELOCITY
quake_displacement = displacement[quake_start:quake_end, :, :]
- slip = np.zeros(num_vertices)
+ #slip = np.zeros(num_vertices)
- for i in range(num_vertices):
- if any(local_slipping_times[:, i]):
- starts = slip_beginnings(local_slipping_times[:, i])
- ends = slip_endings(local_slipping_times[:, i])
- slip[i] = np.linalg.norm(quake_displacement[ends, i, :]
- - quake_displacement[starts, i, :])
+ # for i in range(num_vertices):
+ # if any(local_slipping_times[:, i]):
+ # starts = slip_beginnings(local_slipping_times[:, i])
+ # ends = slip_endings(local_slipping_times[:, i])
+ # slip[i] = np.linalg.norm(quake_displacement[ends, i, :]
+ # - quake_displacement[starts, i, :])
- quake['peakSlip'] = max(slip)
- quake['peakSlipRate'] = max(maximum_velocity[quake_start:quake_end])
+ # quake['peakSlip'] = max(slip)
+ quake['peakSlipRate'] = max(max_rate[quake_start:quake_end])
max_rupture_width = 0
for time_step in range(local_slipping_times.shape[0]):
diff --git a/data/tools/config.ini b/data/tools/config.ini
index 0d5e43204cf1209fe9eea0435f41b1e11a335396..053017c1bb8accdfe5279a84acde4dfb8eb5931e 100644
--- a/data/tools/config.ini
+++ b/data/tools/config.ini
@@ -4,6 +4,6 @@
#/home/joscha/Downloads/
[directories]
-simulation = /home/joscha/software/dune/build-release/dune-tectonic/src/multi-body-problem/output/newmark/
+simulation = /home/joscha/software/dune/build-release/dune-tectonic/src/strikeslip/output/newmark-uniform-1e4
experiment = ~/group/publications/2016-RosenauCorbiDominguezRudolfRitterPipping
output = generated
diff --git a/data/tools/main.py b/data/tools/main.py
index 332d2747e5e047d77d2ed747209bb5a142c9153d..7bb3eef822c981c53c2c45b9afac320a350f6c08 100644
--- a/data/tools/main.py
+++ b/data/tools/main.py
@@ -43,7 +43,7 @@ params = read_params('strikeslip.cfg')
TANGENTIAL_COORDS = 1
FINAL_TIME = params['finalTime']
NBODIES = params['bodyCount']
-THRESHOLD_VELOCITY = 1e-3
+THRESHOLD_VELOCITY = 1e-2
h5file = read_h5file()
print(list(h5file.keys()))
@@ -66,7 +66,7 @@ for body_ID in range(NBODIES):
friction_stats(h5file, body_ID, FINAL_TIME, patch, interval)
#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)
+ #[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 81e57ce259326255357f6b82cc35c02f83db661d..380dd6821f1ae60d56d7fb91ff63320f6d537af8 100644
--- a/data/tools/support/find_quakes.py
+++ b/data/tools/support/find_quakes.py
@@ -11,42 +11,21 @@ def peak_slip(quake_start, quake_end, velocities):
#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
+def find_quakes(rate, threshold_rate):
+ slipping_times = rate > threshold_rate
quake_starts = slip_beginnings(slipping_times)
quake_ends = slip_endings(slipping_times)
- #print(quake_starts)
- #print(quake_ends)
+ 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))
+ return [quake_starts, quake_ends]
peak_v = np.zeros(min_len)
quake_times = np.zeros(min_len)
@@ -59,4 +38,4 @@ def find_quakes(h5file, body_ID, FINAL_TIME, patch, interval, threshold_velocity
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]
+