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 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 = abs(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') #------------------------- diffplot(friction_coeff[1], friction_coeff_calc[1], 'diff friction coeff') iterations(h5file, FINAL_TIME) plt.show() h5file.close()