truncation.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.io import read_h5file
from support.io import read_params

from support.iterations import iterations
from support.friction_stats import friction_stats

def build_time(h5file, final_time, interval = []):
    # read time
    time = np.array(h5file['relativeTime']) * FINAL_TIME
    time = np.delete(time, 0)
    
    if len(interval) == 0:
        interval = [0, FINAL_TIME]
    
    t = [i for i in range(len(time)) if time[i]>=interval[0] and time[i]<=interval[1]]
    return t, time[t]

def build_patch(coords, percentage = 1.0):
    x_coords = coords[:, 0]
    xmin = np.min(x_coords)
    xmax = np.max(x_coords)
    delta_x = (1 - percentage)*(xmax - xmin)/2

    xmin = xmin + delta_x
    xmax = xmax - delta_x

    return [i for i in range(len(x_coords)) if x_coords[i]>=xmin and x_coords[i]<=xmax]

# read h5 output file
h5file = read_h5file()
print(list(h5file.keys()))

# read problem parameters
params = read_params('foam.cfg')
print(params)

TANGENTIAL_COORDS = 1
FINAL_TIME = params['finalTime']
NBODIES = params['bodyCount']

t, time = build_time(h5file, FINAL_TIME, [0, FINAL_TIME])

for body_ID in range(NBODIES):
    body = 'body' + str(body_ID)

    if body not in h5file:
        continue

    coords = np.array(h5file[body + '/coordinates'])
    patch = build_patch(coords, 1.0)

    fig = plt.figure()

    # velocity and state data
    v = abs(np.array(h5file[body + '/velocity']))[:,:,TANGENTIAL_COORDS]
    alpha = np.array(h5file[body + '/state'])

    vmin = 0 #params['V0'] / np.exp( ( params['mu0'] + params['b'] * np.array(alpha)) / params['a'] )
    
    v_truncation = (v<=vmin)[t,:]
    v_truncated = np.sum(v_truncation, axis=1)

    v_comp = (np.abs(v-vmin)<1e-14)[t,:]
    v_small = np.sum(v_comp, axis=1)

    weighted_normal_stress = np.array(h5file[body + '/weightedNormalStress'])
    second_deriv_large = (v<0)[t,:] #((- weighted_normal_stress[None, ...] * params['a'] / v)>1e10)[0,t,:]
    second_deriv = np.sum(second_deriv_large, axis=1) 

    total_truncated = np.sum(v_truncation | v_comp | second_deriv_large, axis=1)

    # plot v_small
    ax_v_small = fig.add_subplot(4, 1, 1)
    ax_v_small.plot(time, v_small, color='black', linestyle='-')
    ax_v_small.set_ylabel('# v-vmin < 1e-14')
    #-------------------------

    # plot v_truncated
    ax_v_truncated = fig.add_subplot(4, 1, 2  )
    ax_v_truncated.plot(time, v_truncated, color='black', linestyle='-')
    ax_v_truncated.set_ylabel('# v <= vmin')
    #-------------------------

    # plot second_deriv_large
    ax_second_deriv = fig.add_subplot(4, 1, 3)
    ax_second_deriv.plot(time, second_deriv, color='black', linestyle='-')
    ax_second_deriv.set_ylabel('# second_deriv large')
    #-------------------------

    # plot total_truncated
    ax_truncated = fig.add_subplot(4, 1, 4)
    ax_truncated.plot(time, total_truncated, color='black', linestyle='-')
    ax_truncated.set_ylabel('# total truncated')
    #-------------------------
plt.show()

h5file.close()