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 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()