data/tools/generate-2d.bash

+#!/usr/bin/env bash
+
+set -e
+
+rr() {
+    echo "$(date)" Running: Rscript $@
+    Rscript --vanilla --default-packages=grDevices,methods,stats,utils tools/$@
+}
+
+# contours
+rr 2d-velocity-contours.R
+
+# dip
+rr 2d-dip-contours.R
+
+# iterations / adaptivity
+rr 2d-performance.R
+
+# box plot (one half)
+rr 2d-event-writer.R
+
+# fpi data
+rr 2d-fpi-tolerance.R
+
+date

data/tools/generate-3d.bash

+#!/usr/bin/env bash
+
+set -e
+
+rr() {
+    echo "$(date)" Running: Rscript $@
+    Rscript --vanilla --default-packages=grDevices,methods,stats,utils tools/$@
+}
+
+# performance
+rr 3d-performance.R
+
+# velocity contours
+rr 3d-velocity-contours.R
+
+date

data/tools/generate-others.bash

+#!/usr/bin/env bash
+
+set -e
+
+rr() {
+    echo "$(date)" Running: Rscript $@
+    Rscript --vanilla --default-packages=grDevices,methods,stats,utils tools/$@
+}
+
+# box plot (one half)
+rr comparison:lab-sim.R
+
+date

data/tools/generated/events1.csv

+start,peakSlipRate,end,ruptureWidth,peakSlip
+20.086059570312422,5.7045853762396602e-06,20.086669921874922,0,0.0
+20.090332031249918,5.6219949973608878e-06,20.092163085937418,0,0.0
+20.094604492187415,5.8412336903395239e-06,20.095825195312415,0,0.0
+20.099487304687411,5.3352579051826403e-06,20.100097656249908,0,0.0
+20.143737792968619,5.4285740142335498e-06,20.144958496093619,0,0.0
+20.148620605468615,5.6370561043864526e-06,20.149536132812365,0,0.0
+20.353851318359059,6.0525078256123743e-06,20.354461669921559,0,0.0
+20.358734130859055,5.1806538408711459e-06,20.359344482421555,0,0.0
+20.371856689452795,5.6547173619088325e-06,20.373077392577795,0,0.0
+20.425262451171495,6.4545930360194781e-06,20.428314208983991,0.088385894124722503,1.4298690605481336e-08
+20.479583740233945,5.8310594282664788e-06,20.481414794921445,0,0.0
+20.483245849608942,5.2537138171449348e-06,20.485076904296442,0,6.2000101383992609e-09
+20.531158447265149,5.2316818926703734e-06,20.531768798827649,0,0.0
+20.544586181640138,5.4092184353134343e-06,20.545196533202635,0,0.0
+20.593109130858842,6.190190600109647e-06,20.595550537108842,0,0.0
+20.608673095702581,5.5618674305029011e-06,20.609283447265081,0,0.0
+20.625152587890067,5.4920875934131527e-06,20.625762939452567,0,0.0
+20.630035400390064,5.1740512343599518e-06,20.630645751952564,0,0.0
+20.63369750976506,5.7392806837644885e-06,20.63491821289006,0,0.0

data/tools/main.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
+from support.slip_rates import slip_rates
+
+from support.find_quakes import find_quakes
+
+def build_patch(coords, percentage):
+    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 problem parameters
+params = read_params('strikeslip.cfg')
+
+TANGENTIAL_COORDS = 1
+FINAL_TIME = params['finalTime']
+NBODIES = params['bodyCount']
+THRESHOLD_VELOCITY = 1e-2 
+
+h5file = read_h5file()
+print(list(h5file.keys()))
+
+interval = [0.5*FINAL_TIME, FINAL_TIME]
+
+iterations(h5file, FINAL_TIME, interval)
+
+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)
+
+    print("patch length: " + str(len(patch)))
+
+    friction_stats(h5file, body_ID, FINAL_TIME, patch, interval)
+    #slip_rates(h5file, body_ID, FINAL_TIME, patch, interval, TANGENTIAL_COORDS)
+
+    #[quake_starts, quake_ends] = find_quakes(h5file, body_ID, FINAL_TIME, patch, interval, THRESHOLD_VELOCITY, TANGENTIAL_COORDS)
+plt.show()
+
+h5file.close()
\ No newline at end of file

data/tools/support/__init__.py

+ 

data/tools/support/find_quakes.py

+import array as ar
+import numpy as np
+
+from .slip_beginnings import slip_beginnings
+from .slip_endings import slip_endings
+
+def peak_slip(quake_start, quake_end, velocities):
+    quake_v = velocities[quake_start-1:quake_end]
+    max_v = np.max(quake_v, axis=1)
+    return [quake_start+np.argmax(max_v), np.max(max_v)]
+
+#def rupture_width
+
+def find_quakes(rate, threshold_rate):
+    slipping_times = rate > threshold_rate
+
+    quake_starts = slip_beginnings(slipping_times)
+    quake_ends = slip_endings(slipping_times)
+
+    print(quake_starts)
+    print(quake_ends)
+
+    # remove incomplete quakes
+    min_len = min(len(quake_starts), len(quake_ends))
+    quake_ends = quake_ends[0:min_len]
+    quake_starts = quake_starts[0:min_len]
+
+    return [quake_starts, quake_ends]
+
+    peak_v = np.zeros(min_len)
+    quake_times = np.zeros(min_len)
+    for i in range(min_len):
+        [quake_times[i], peak_v[i]] = peak_slip(quake_starts[i], quake_ends[i], v_tx)
+    print("peak slip velocity: " + str(peak_v.mean()) + " +- " + str(peak_v.std()))
+
+    recurrence_time = np.zeros(min_len-1)
+    for i in range(min_len-1):
+        recurrence_time[i] = time[int(quake_times[i+1])] - time[int(quake_times[i])]
+    print("recurrence time: " + str(recurrence_time.mean()) + " +- " + str(recurrence_time.std()))
+
+   

data/tools/support/find_quakes.py~

+from support.slip_beginnings import slip_beginnings
+from support.slip_endings import slip_endings
+
+
+def find_quakes(threshold_velocity, maximum_velocities):
+    slipping_times = maximum_velocities > threshold_velocity
+    print(slipping_times)
+
+    quake_starts = slip_beginnings(slipping_times)
+    quake_ends = slip_endings(slipping_times)
+
+    # remove incomplete quakes
+    min_len = min(len(quake_starts), len(quake_ends))
+    quake_ends = quake_ends[0:min_len]
+    quake_starts = quake_starts[0:min_len]
+
+    return [quake_starts, quake_ends]

data/tools/support/friction_stats.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+def friction_stats(h5file, body_ID, FINAL_TIME, patch=[], interval=[], TANGENTIAL_COORDS=1):
+    body = 'body' + str(body_ID)   # 'frictionalBoundary'     'body' + str(body_ID)
+
+    coords = np.array(h5file[body + '/coordinates'])
+    if len(patch) == 0:
+        patch = np.linspace(0, len(coords)-1, len(coords), endpoint=True, dtype='int8')
+
+    # 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]]
+    #t = t[::5]
+    time = time[t]
+
+    fig = plt.figure()
+
+    # velocity data
+    v = abs(np.array(h5file[body + '/velocity']))
+    v_t = v[t,:,TANGENTIAL_COORDS]
+    v_tx = v_t[:,patch]
+    # statistics
+    avg_v = np.average(v_tx, axis=1)
+    min_v = np.min(v_tx, axis=1)
+    max_v = np.max(v_tx, axis=1)
+    # plot
+    ax_slip = fig.add_subplot(1, 1, 1)
+    #ax_slip.plot(time, min_v, color='gray', linestyle='--')
+    ax_slip.plot(time, avg_v, color='black', linestyle='-')
+    #ax_slip.plot(time, max_v, color='gray', linestyle='--')
+    ax_slip.set_ylabel('slip rate V [m/s]')
+    ax_slip.set_xlabel('time t [s]')
+    ax_slip.set_yscale('log')
+    #ax_slip.set_ylim([1e-6,1e-2])
+    #-------------------------
+
+    print(np.min(min_v))
+    print(np.average(avg_v))
+    print(np.max(max_v))
+
+    # state
+    #states = np.array(h5file[body + '/state'])
+    #states_t = states[t,:]
+    #states_tx = states_t[:,patch]
+    # statistics
+    #avg_states = np.average(states_tx, axis=1)
+    #min_states = np.min(states_tx, axis=1)
+    #max_states = np.max(states_tx, axis=1)
+    # plot
+    #ax_state = fig.add_subplot(2, 1, 2)
+    #ax_state.plot(time, min_states, color='gray', linestyle='--')
+    #ax_state.plot(time, avg_states, color='black', linestyle='-')
+    #ax_state.plot(time, max_states, color='gray', linestyle='--')
+    #ax_state.set_ylabel('state')
+    #ax_state.set_xlabel('time [s]')
+    #-------------------------
+
+    # friction coefficient
+    # friction_coeff = np.array(h5file[body + '/coefficient'])
+    # friction_coeff_t = friction_coeff[t,:]
+    # friction_coeff_tx = friction_coeff_t[:,patch]
+    # # statistics
+    # avg_friction_coeff = np.average(friction_coeff_tx, axis=1)
+    # min_friction_coeff = np.min(friction_coeff_tx, axis=1)
+    # max_friction_coeff = np.max(friction_coeff_tx, axis=1)
+    # # plot
+    # ax_friction = fig.add_subplot(3, 1, 3)
+    # ax_friction.plot(time, min_friction_coeff, color='gray', linestyle='--')
+    # ax_friction.plot(time, avg_friction_coeff, color='black', linestyle='-')
+    # ax_friction.plot(time, max_friction_coeff, color='gray', linestyle='--')
+    # ax_friction.set_ylabel('friction coefficient')
+    #-------------------------
+
+    fig.canvas.draw()
\ No newline at end of file

data/tools/support/io.py

+import configparser as cp
+import os
+import h5py
+
+def remove_comment(str, char = "#"):
+    split_str = str.split(char, 1)
+    return split_str[0]
+
+# tag = 'simulation', 'experiment', 'output'
+def read_h5file(tag = 'simulation'):
+    # read config ini
+    config = cp.ConfigParser()
+    config_path = os.path.join('tools/config.ini')
+    config.read(config_path)
+    path = config.get('directories', tag)
+
+    # read hdf5 output file
+    h5path = os.path.join(path)
+
+    return h5py.File(os.path.join(h5path, 'output.h5'), 'r')
+
+# tag = 'simulation', 'experiment', 'output'
+def read_params(file_name, tag = 'simulation'):
+    # read config ini
+    config = cp.ConfigParser()
+    config_path = os.path.join('tools/config.ini')
+    config.read(config_path)
+    path = config.get('directories', tag)
+
+    # read cfg parameter file
+    params = cp.ConfigParser()
+    params_path = os.path.join(path, file_name)
+    params.read(params_path)
+
+    #with open(params_path) as stream:
+    #    params.read_string("[top]\n" + stream.read()) 
+
+    out = {
+        'L'  : float(remove_comment(params.get('boundary.friction', 'L'))),
+        'V0' : float(remove_comment(params.get('boundary.friction', 'V0'))),
+        'mu0': float(remove_comment(params.get('boundary.friction', 'mu0'))), 
+        'a'  : float(remove_comment(params.get('boundary.friction.weakening', 'a'))), 
+        'b'  : float(remove_comment(params.get('boundary.friction.weakening', 'b'))),
+        'bodyCount' : int(remove_comment(params.get('problem', 'bodyCount'))),
+        'finalTime' : float(remove_comment(params.get('problem', 'finalTime'))),
+        'finalVelocity' : float(remove_comment(params.get('boundary.dirichlet', 'finalVelocity')))
+    }
+
+    return out
\ No newline at end of file

data/tools/support/iterations.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+def iterations(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]]
+    time = time[t]
+
+    tau = np.array(h5file['relativeTimeIncrement']) * FINAL_TIME
+    tau = np.delete(tau, 0)
+
+    # read fpi iterations
+    fpi_final = np.array(h5file['iterations/fixedPoint/final'])
+    fpi_final = np.delete(fpi_final, 0)
+    #fpi_total = np.array(h5file['iterations/fixedPoint/total'])
+
+    print("FPI average: " + str(np.average(fpi_final)))
+    print("FPI max: " + str(np.max(fpi_final)))
+
+    # read multigrid iterations
+    multigrid_final = np.array(h5file['iterations/multiGrid/final'])
+    multigrid_final = np.delete(multigrid_final, 0)
+    #multigrid_total = np.array(h5file['iterations/multiGrid/total'])
+
+     # plot
+    fig = plt.figure()
+
+    ax_fpi = fig.add_subplot(3, 1, 1)
+    ax_fpi.plot(time, fpi_final[t], color='black', linestyle='-')
+    ax_fpi.set_ylabel('fpi')
+    #-------------------------
+
+    ax_mg = fig.add_subplot(3, 1, 2)
+    ax_mg.plot(time, multigrid_final[t], color='black', linestyle='-')
+    ax_mg.set_ylabel('multigrid iter.')
+    ax_mg.set_xlabel('time [s]')
+    #-------------------------
+
+    tau_t = tau[t]
+    ax_tau = fig.add_subplot(3, 1, 3)
+    ax_tau.plot(time, tau_t, color='black', linestyle='-')
+    ax_tau.set_ylabel('tau')
+    ax_tau.set_yscale('log')
+
+    print("tau_max / tau_min: " + str(np.max(tau_t)/np.min(tau_t)))
+    #-------------------------
+
+    fig.canvas.draw()
\ No newline at end of file

data/tools/support/max_distance.py

+from itertools import combinations
+
+
+def square_distance(x, y):
+    return sum([(xi-yi)**2 for xi, yi in zip(x, y)])
+
+
+def max_distance(points):
+    max_square_distance = 0
+    max_pair = (0, 0)
+    for pair in combinations(points, 2):
+        sq_dist = square_distance(*pair)
+
+        if sq_dist > max_square_distance:
+            max_square_distance = sq_dist
+            max_pair = pair
+    return max_pair

data/tools/support/max_distance.py~

+from itertools import combinations
+
+
+def square_distance(x, y):
+    return sum([(xi-yi)**2 for xi, yi in zip(x, y)])
+
+
+def max_distance(points):
+    max_square_distance = 0
+    max_pair = (0, 0)
+    for pair in combinations(points, 2):
+        sq_dist = square_distance(*pair)
+        print(sq_dist)
+        if sq_dist > max_square_distance:
+            max_square_distance = sq_dist
+            max_pair = pair
+    return max_pair

data/tools/support/max_velocity.py

+import numpy as np
+
+def max_velocity(x) :
+    return maximum(norm(x))
\ No newline at end of file

data/tools/support/max_velocity.py~

+import numpy as np
+
+def max_velocity(x) :
+    size = x.shape 
+    res = np.zeros(size[0])
+
+    for time_step in range(size[0]) :
+        for vertex in range(size[1]) :
+            res[time_step] = max(res[time_step], np.linalg.norm(x[time_step, vertex, :]))
+
+    return res
\ No newline at end of file

data/tools/support/maximum.py

+import numpy as np
+
+
+def maximum(x):
+    size = x.shape
+    res = np.zeros(size[0])
+
+    for time_step in range(size[0]):
+        res[time_step] = max(res[time_step], max(x[time_step, :]))
+
+    return res