python based debug tool

data/tools/2d_event_writer.py

-import ConfigParser as cp
+import configparser as cp
 import os
 import numpy as np
 import csv

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/
+
 [directories]
-simulation = /storage/mi/podlesny/software/dune/build-debug/dune-tectonic/src/strikeslip
+simulation = /home/joscha/software/dune/build-release/dune-tectonic/src/foam/output/test/
 experiment = ~/group/publications/2016-RosenauCorbiDominguezRudolfRitterPipping
 output     = generated

data/tools/debug.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 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
+
+
+NBODIES = 2
+FINAL_TIME = 15  # s
+FINAL_VELOCITY = 1e-5  # m/s
+THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY  # 1000e-6 + FINAL_VELOCITY
+
+TANGENTIAL_COORDS = 1
+
+# 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['body1'].keys()))
+
+# read time
+relative_time = np.array(h5file['relativeTime'])
+relative_tau = np.array(h5file['relativeTimeIncrement'])
+relative_tau = np.delete(relative_tau, 0)
+
+real_time = relative_time * FINAL_TIME
+real_tau = relative_tau * FINAL_TIME
+
+print(len(relative_time))
+
+for body_ID in range(NBODIES):
+    body = 'body' + str(body_ID)
+
+    if body not in h5file:
+        continue
+
+    # velocity data
+    v = np.array(h5file[body + '/velocity'])
+
+    # statistics
+    avg_v = np.average(v[:,:,TANGENTIAL_COORDS], axis=1)
+    min_v = np.min(v[:,:,TANGENTIAL_COORDS], axis=1)
+    max_v = np.max(v[:,:,TANGENTIAL_COORDS], axis=1)
+
+    # plot
+    plt.figure(1)
+    plt.subplot(311)
+    plt.plot(min_v, color='gray', linestyle='--')
+    plt.plot(avg_v, color='black', linestyle='-')
+    plt.plot(max_v, color='gray', linestyle='--')
+    plt.ylabel('slip rate')
+    #-------------------------
+
+    # state
+    states = np.array(h5file[body + '/state'])
+    states_calc = aging_law(params, states[0], v[:,:,TANGENTIAL_COORDS], real_tau)
+
+    # statistics
+    avg_states = np.average(states, axis=1)
+    avg_states_calc  = np.average(states_calc, axis=1)
+    min_states = np.min(states, axis=1)
+    max_states = np.max(states, axis=1)
+
+    # plot
+    plt.subplot(312)
+    plt.plot(min_states, color='gray', linestyle='--')
+    plt.plot(avg_states, color='black', linestyle='-')
+    plt.plot(max_states, color='gray', linestyle='--')
+    plt.plot(avg_states_calc, color='red', linestyle='-')
+    plt.ylabel('state')
+    #-------------------------
+
+    # friction coefficient
+    friction_coeff = np.array(h5file[body + '/coefficient'])
+
+    #weighted_normal_stress = np.array(h5file[body + '/weightedNormalStress'])
+    friction_coeff_calc = truncated_friction(params, v[:,:,TANGENTIAL_COORDS], states_calc)
+
+    # statistics
+    avg_friction_coeff = np.average(friction_coeff, axis=1)
+    avg_friction_coeff_calc = np.average(friction_coeff_calc, axis=1)
+    min_friction_coeff = np.min(friction_coeff, axis=1)
+    max_friction_coeff = np.max(friction_coeff, axis=1)
+
+    outliers_friction_coeff = outliers(avg_friction_coeff)
+
+    # plot
+    plt.subplot(313)
+    plt.plot(min_friction_coeff, color='gray', linestyle='--')
+    plt.plot(avg_friction_coeff, color='black', linestyle='-')
+    plt.plot(max_friction_coeff, color='gray', linestyle='--')
+    plt.plot(avg_friction_coeff_calc, color='red', linestyle='-')
+    plt.plot(outliers_friction_coeff[0], outliers_friction_coeff[1], color='red', marker='+')
+    plt.ylabel('friction coefficient')
+    #-------------------------
+
+    plt.show()
+
+h5file.close()
\ No newline at end of file

data/tools/debug/friction.py

+import numpy as np
+
+def truncated_friction(params, v, alpha):
+   vmin = params['V0'] / np.exp( ( params['mu0'] + params['b'] * np.array(alpha)) / params['a'] )
+   clipped_v = (v/vmin).clip(1)
+   return params['a'] * np.log(clipped_v)
+
+   #return (params['a'] * np.log(v/params['V0']) - params['mu0'] - (params['b']/params['a'])*np.array(alpha) ).clip(0)
\ No newline at end of file

data/tools/debug/outliers.py

+import numpy as np
+
+def outliers(data):
+    index = [i for i,x in enumerate(data) if np.abs(x)==np.inf]
+    val = np.zeros(len(index))
+    return [index, val]
\ No newline at end of file

data/tools/debug/state.py

+import numpy as np
+
+def lift_singularity(c, x):
+    def local_lift(y): 
+        if (y <= 0):
+            return -c
+        else:
+            return -np.expm1(c * y) / y
+    
+    return np.array([local_lift(y) for y in x])
+
+
+def aging_law(params, initial_alpha, v, time_steps):
+        #auto tangentVelocity = velocity_field[localToGlobal_[i]];
+        #tangentVelocity[0] = 0.0;
+        #double const V = tangentVelocity.two_norm();
+        
+    alpha = [initial_alpha]
+
+    for i,tau in enumerate(time_steps):
+        mtoL = -tau / params['L']
+        next_alpha = np.log(np.exp(alpha[-1] + v[i] * mtoL) + params['V0'] * lift_singularity(mtoL, v[i]))
+        alpha.append(next_alpha)
+
+    return alpha
\ No newline at end of file