From 8ea3315795b4ced755c1233540ffa69510cb3a09 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Fri, 19 Feb 2021 16:29:19 +0100
Subject: [PATCH] python based debug tool
---
data/tools/2d_event_writer.py | 2 +-
data/tools/config.ini | 6 +-
data/tools/debug.py | 129 ++++++++++++++++++++++++++++++++++
data/tools/debug/__init__.py | 0
data/tools/debug/friction.py | 8 +++
data/tools/debug/outliers.py | 6 ++
data/tools/debug/state.py | 25 +++++++
7 files changed, 174 insertions(+), 2 deletions(-)
create mode 100644 data/tools/debug.py
create mode 100644 data/tools/debug/__init__.py
create mode 100644 data/tools/debug/friction.py
create mode 100644 data/tools/debug/outliers.py
create mode 100644 data/tools/debug/state.py
diff --git a/data/tools/2d_event_writer.py b/data/tools/2d_event_writer.py
index d36e85be..f79edf66 100644
--- a/data/tools/2d_event_writer.py
+++ b/data/tools/2d_event_writer.py
@@ -1,4 +1,4 @@
-import ConfigParser as cp
+import configparser as cp
import os
import numpy as np
import csv
diff --git a/data/tools/config.ini b/data/tools/config.ini
index d8a6165b..40b64b13 100644
--- a/data/tools/config.ini
+++ b/data/tools/config.ini
@@ -1,4 +1,8 @@
+#/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
diff --git a/data/tools/debug.py b/data/tools/debug.py
new file mode 100644
index 00000000..cd37940a
--- /dev/null
+++ b/data/tools/debug.py
@@ -0,0 +1,129 @@
+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
diff --git a/data/tools/debug/__init__.py b/data/tools/debug/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/data/tools/debug/friction.py b/data/tools/debug/friction.py
new file mode 100644
index 00000000..47df0de2
--- /dev/null
+++ b/data/tools/debug/friction.py
@@ -0,0 +1,8 @@
+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
diff --git a/data/tools/debug/outliers.py b/data/tools/debug/outliers.py
new file mode 100644
index 00000000..e7e47a6a
--- /dev/null
+++ b/data/tools/debug/outliers.py
@@ -0,0 +1,6 @@
+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
diff --git a/data/tools/debug/state.py b/data/tools/debug/state.py
new file mode 100644
index 00000000..bb17d028
--- /dev/null
+++ b/data/tools/debug/state.py
@@ -0,0 +1,25 @@
+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
--
GitLab