diff --git a/data/tools/2d_event_writer.py b/data/tools/2d_event_writer.py index d36e85beaa0ab513ca90d46d73fc2f89347fa490..f79edf660786f0443b772ec6cdbbd2a52262e533 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 d8a6165b28fd6ef91665bfedf4653f55827fc05f..40b64b13c922d79bb6319d5c1907f75384e49378 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 0000000000000000000000000000000000000000..cd37940adb0fb7b092f0c31a48e92c7ac07ba3b5 --- /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 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/data/tools/debug/friction.py b/data/tools/debug/friction.py new file mode 100644 index 0000000000000000000000000000000000000000..47df0de21313c984ee8bb0bc9307d614e64f006f --- /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 0000000000000000000000000000000000000000..e7e47a6a25463e2246feeda3ba46d6df3a8eb937 --- /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 0000000000000000000000000000000000000000..bb17d02832d6a65ec0b4ac2b65f4febf6e2e63ec --- /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