.

data/tools/config.ini

 #/home/joscha/software/dune/build-release/dune-tectonic/src/foam/output/test/
 #/home/joscha/software/dune/build-debug/dune-tectonic/src/foam/output/test/
 #/home/joscha/Desktop/strikeslip/viscosity-1e3/
+ #/home/joscha/Downloads/
 
 [directories]
-simulation = /home/joscha/software/dune/build-release/dune-tectonic/src/foam/output/tresca/
+simulation = /home/joscha/software/dune/build-release/dune-tectonic/src/foam/output/tresca-0.6/
 experiment = ~/group/publications/2016-RosenauCorbiDominguezRudolfRitterPipping
 output     = generated

data/tools/elias.py

+import configparser as cp
+import os
+import numpy as np
+import csv
+import h5py
+import matplotlib.pyplot as plt
+
+from debug.outliers import outliers
+from debug.friction import truncated_friction
+from debug.state import aging_law
+from debug.diffplot import diffplot
+
+from support.maximum import maximum
+from support.norm import norm
+from support.find_quakes import find_quakes
+from support.slip_beginnings import slip_beginnings
+from support.slip_endings import slip_endings
+from support.max_distance import max_distance
+
+from support.iterations import iterations
+from support.friction_stats import friction_stats
+
+def build_patch(coords, percentage):
+    x_coords = coords[:, 0]
+    xmin = np.min(x_coords)
+    xmax = np.max(x_coords)
+    delta_x = (1 - percentage)*(xmax - xmin)/2
+
+    xmin = xmin + delta_x
+    xmax = xmax - delta_x
+
+    return [i for i in range(len(x_coords)) if x_coords[i]>=xmin and x_coords[i]<=xmax]
+
+FINAL_TIME = 1000  # s
+FINAL_VELOCITY = 5e-5  # m/s
+THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY  # 1000e-6 + FINAL_VELOCITY
+
+TANGENTIAL_COORDS = 0
+
+# friction params
+params = {
+    'L'  : 1e-5,
+    'V0' : 1e-6,
+    'mu0': 0.6,
+    'a'  : 0.010,
+    'b'  : 0.015
+}
+
+# read config ini
+config = cp.ConfigParser()
+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')
+
+# read hdf5 output file
+h5path = os.path.join(sim_path)
+h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')
+
+print(list(h5file.keys()))
+print(list(h5file['frictionalBoundary'].keys()))
+
+iterations(h5file, FINAL_TIME)
+
+coords = np.array(h5file['frictionalBoundary/coordinates'])
+patch = build_patch(coords, 0.05)
+
+friction_stats(h5file, 0, FINAL_TIME, [], [0, 50], TANGENTIAL_COORDS)
+    
+plt.show()
+
+h5file.close()
\ No newline at end of file

data/tools/main.py

@@ -17,6 +17,8 @@ from support.slip_beginnings import slip_beginnings
 from support.slip_endings import slip_endings
 from support.max_distance import max_distance
 
+from support.io import read_h5file
+
 from support.iterations import iterations
 from support.friction_stats import friction_stats
 
@@ -33,7 +35,7 @@ def build_patch(coords, percentage):
 
 
 NBODIES = 2
-FINAL_TIME = 15  # s
+FINAL_TIME = 150  # s
 FINAL_VELOCITY = 2e-4  # m/s
 THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY  # 1000e-6 + FINAL_VELOCITY
 
@@ -48,20 +50,8 @@ params = {
     'b'  : 0.015
 }
 
-# read config ini
-config = cp.ConfigParser()
-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')
-
-# read hdf5 output file
-h5path = os.path.join(sim_path)
-h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')
-
+h5file = read_h5file()
 print(list(h5file.keys()))
-print(list(h5file['body1'].keys()))
 
 iterations(h5file, FINAL_TIME)
 
@@ -74,7 +64,7 @@ for body_ID in range(NBODIES):
     coords = np.array(h5file[body + '/coordinates'])
     patch = build_patch(coords, 0.05)
 
-    friction_stats(h5file, body_ID, FINAL_TIME, [], [0, 3])
+    friction_stats(h5file, body_ID, FINAL_TIME, [], [0, FINAL_TIME])
     
 plt.show()
 

data/tools/support/friction_stats.py

 import numpy as np
 import matplotlib.pyplot as plt
 
-TANGENTIAL_COORDS = 1
-
-def friction_stats(h5file, body_ID, FINAL_TIME, patch=[], interval=[]):
-    body = 'body' + str(body_ID)
+def friction_stats(h5file, body_ID, FINAL_TIME, patch=[], interval=[], TANGENTIAL_COORDS=1):
+    body = 'body' + str(body_ID) # 'frictionalBoundary'     'body' + str(body_ID)
 
     coords = np.array(h5file[body + '/coordinates'])
     if len(patch) == 0:
@@ -30,12 +28,17 @@ def friction_stats(h5file, body_ID, FINAL_TIME, patch=[], interval=[]):
     max_v = np.max(v_tx, axis=1)
     # plot
     ax_slip = fig.add_subplot(3, 1, 1)
-    ax_slip.plot(time, min_v, color='gray', linestyle='--')
+    #ax_slip.plot(time, min_v, color='gray', linestyle='--')
     ax_slip.plot(time, avg_v, color='black', linestyle='-')
-    ax_slip.plot(time, max_v, color='gray', linestyle='--')
+    #ax_slip.plot(time, max_v, color='gray', linestyle='--')
     ax_slip.set_ylabel('slip rate')
+    #ax_slip.set_yscale('log')
     #-------------------------
 
+    print(np.min(min_v))
+    print(np.average(avg_v))
+    print(np.max(max_v))
+
     # state
     states = np.array(h5file[body + '/state'])
     states_t = states[t,:]

data/tools/support/io.py

+import configparser as cp
+import os
+import h5py
+
+def remove_comment(str, char = "#"):
+    split_str = str.split(char, 1)
+    return split_str[0]
+
+# tag = 'simulation', 'experiment', 'output'
+def read_h5file(tag = 'simulation'):
+    # read config ini
+    config = cp.ConfigParser()
+    config_path = os.path.join('tools/config.ini')
+    config.read(config_path)
+    path = config.get('directories', tag)
+
+    # read hdf5 output file
+    h5path = os.path.join(path)
+
+    return h5py.File(os.path.join(h5path, 'output.h5'), 'r')
+
+# tag = 'simulation', 'experiment', 'output'
+def read_params(file_name, tag = 'simulation'):
+    # read config ini
+    config = cp.ConfigParser()
+    config_path = os.path.join('tools/config.ini')
+    config.read(config_path)
+    path = config.get('directories', tag)
+
+    # read cfg parameter file
+    params = cp.ConfigParser()
+    params_path = os.path.join(path, file_name)
+    params.read(params_path)
+
+    #with open(params_path) as stream:
+    #    params.read_string("[top]\n" + stream.read()) 
+
+    out = {
+        'L'  : float(remove_comment(params.get('boundary.friction', 'L'))),
+        'V0' : float(remove_comment(params.get('boundary.friction', 'V0'))),
+        'mu0': float(remove_comment(params.get('boundary.friction', 'mu0'))), 
+        'a'  : float(remove_comment(params.get('boundary.friction.weakening', 'a'))), 
+        'b'  : float(remove_comment(params.get('boundary.friction.weakening', 'b'))),
+        'bodyCount' : int(remove_comment(params.get('problem', 'bodyCount'))),
+        'finalTime' : float(remove_comment(params.get('problem', 'finalTime'))),
+        'finalVelocity' : float(remove_comment(params.get('boundary.dirichlet', 'finalVelocity')))
+    }
+
+    return out
\ No newline at end of file

data/tools/truncation.py

+import configparser as cp
+import os
+import numpy as np
+import csv
+import h5py
+import matplotlib.pyplot as plt
+
+from debug.outliers import outliers
+from debug.friction import truncated_friction
+from debug.state import aging_law
+from debug.diffplot import diffplot
+
+from support.maximum import maximum
+from support.norm import norm
+from support.find_quakes import find_quakes
+from support.slip_beginnings import slip_beginnings
+from support.slip_endings import slip_endings
+from support.max_distance import max_distance
+
+from support.io import read_h5file
+from support.io import read_params
+
+from support.iterations import iterations
+from support.friction_stats import friction_stats
+
+def build_time(h5file, final_time, interval = []):
+    # 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]]
+    return t, time[t]
+
+def build_patch(coords, percentage = 1.0):
+    x_coords = coords[:, 0]
+    xmin = np.min(x_coords)
+    xmax = np.max(x_coords)
+    delta_x = (1 - percentage)*(xmax - xmin)/2
+
+    xmin = xmin + delta_x
+    xmax = xmax - delta_x
+
+    return [i for i in range(len(x_coords)) if x_coords[i]>=xmin and x_coords[i]<=xmax]
+
+# read h5 output file
+h5file = read_h5file()
+print(list(h5file.keys()))
+
+# read problem parameters
+params = read_params('foam.cfg')
+print(params)
+
+TANGENTIAL_COORDS = 1
+FINAL_TIME = params['finalTime']
+NBODIES = params['bodyCount']
+
+t, time = build_time(h5file, FINAL_TIME, [0, FINAL_TIME])
+
+for body_ID in range(NBODIES):
+    body = 'body' + str(body_ID)
+
+    if body not in h5file:
+        continue
+
+    coords = np.array(h5file[body + '/coordinates'])
+    patch = build_patch(coords, 1.0)
+
+    fig = plt.figure()
+
+    # velocity and state data
+    v = abs(np.array(h5file[body + '/velocity']))[:,:,TANGENTIAL_COORDS]
+    alpha = np.array(h5file[body + '/state'])
+
+    vmin = 0 #params['V0'] / np.exp( ( params['mu0'] + params['b'] * np.array(alpha)) / params['a'] )
+    
+    v_truncation = (v<=vmin)[t,:]
+    v_truncated = np.sum(v_truncation, axis=1)
+
+    v_comp = (np.abs(v-vmin)<1e-14)[t,:]
+    v_small = np.sum(v_comp, axis=1)
+
+    weighted_normal_stress = np.array(h5file[body + '/weightedNormalStress'])
+    second_deriv_large = (v<0)[t,:] #((- weighted_normal_stress[None, ...] * params['a'] / v)>1e10)[0,t,:]
+    second_deriv = np.sum(second_deriv_large, axis=1) 
+
+    total_truncated = np.sum(v_truncation | v_comp | second_deriv_large, axis=1)
+
+    # plot v_small
+    ax_v_small = fig.add_subplot(4, 1, 1)
+    ax_v_small.plot(time, v_small, color='black', linestyle='-')
+    ax_v_small.set_ylabel('# v-vmin < 1e-14')
+    #-------------------------
+
+    # plot v_truncated
+    ax_v_truncated = fig.add_subplot(4, 1, 2  )
+    ax_v_truncated.plot(time, v_truncated, color='black', linestyle='-')
+    ax_v_truncated.set_ylabel('# v <= vmin')
+    #-------------------------
+
+    # plot second_deriv_large
+    ax_second_deriv = fig.add_subplot(4, 1, 3)
+    ax_second_deriv.plot(time, second_deriv, color='black', linestyle='-')
+    ax_second_deriv.set_ylabel('# second_deriv large')
+    #-------------------------
+
+    # plot total_truncated
+    ax_truncated = fig.add_subplot(4, 1, 4)
+    ax_truncated.plot(time, total_truncated, color='black', linestyle='-')
+    ax_truncated.set_ylabel('# total truncated')
+    #-------------------------
+plt.show()
+
+h5file.close()
\ No newline at end of file