update data evaluation

data/generated/events1.csv

+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

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]):

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

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

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]
+