Fig_07.ipynb

Kernel: Python 3 (system-wide)

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.transforms as mtransforms
import matplotlib.font_manager

from scipy import signal

plt.style.use('../jfm.mplstyle')

data_dir = '../Experimental_Dataset/with_flow/5cm_water_depth/particle_image_velocimitry'

x_int = np.loadtxt(f'{data_dir}/velocity_x_interp.txt',dtype=int)
y_int = np.loadtxt(f'{data_dir}/velocity_y_interp.txt',dtype=int)
u = np.loadtxt(f'{data_dir}/velocity_u_interp.txt',dtype="float32")
v = np.loadtxt(f'{data_dir}/velocity_v_interp.txt',dtype="float32")
m = np.loadtxt(f'{data_dir}/velocity_m_interp.txt',dtype="float32")
u_std = np.loadtxt(f'{data_dir}/velocity_u_std_interp.txt',dtype="float32")
v_std = np.loadtxt(f'{data_dir}/velocity_v_std_interp.txt',dtype="float32")
m_std = np.loadtxt(f'{data_dir}/velocity_m_std_interp.txt',dtype="float32")

h0, g, b0 = 50, 9810, 229
c0 = np.sqrt(g*h0)

mask = x_int==-500

#x_flow = (x_int/h0)
y_flow = (y_int[mask]/h0)

u_flow = u[mask]/c0

w = signal.get_window('boxcar',15)
w = w/np.sum(w)
u_flow_smooth = signal.convolve(u_flow,w,'valid')
y_flow_smooth = signal.convolve(y_flow,w,'valid')

!find ../Experimental_Dataset -type f -name "elevation.csv.gz" | xargs dirname | sort | grep "sol3rd" | tee "data_directories.txt" > /dev/null

with open("data_directories.txt","r") as data_directories_file:
    data_directories = data_directories_file.read().splitlines()

x_loc = "-50.00cm"
#x_loc = "-100.0cm"

amps = [5,10,15,20]
amp_strings = [ f"a{a/10:.1f}cm" for a in amps ]

fig, axes = plt.subplots(nrows=3,ncols=1,figsize=(5.3,4.0),sharex='all')

plt.xlim(0,12.5)

ax = axes[0]

ax.axhline(0,color='k',linewidth=0.5)
ax.plot(y_flow_smooth,u_flow_smooth,color='k')
ax.set_ylim(-0.1,0.005)
ax.set_ylabel(r'$u$')

ax = axes[1]

for amp_str, amp in zip(amp_strings,amps):
    data_dir = sorted([d for d in data_directories if x_loc in d and amp_str in d])[0]
    
    #time  = np.loadtxt(f"{data_dir}/time.csv",delimiter=',',dtype=int)
    space = np.loadtxt(f"{data_dir}/space.csv",delimiter=',',dtype=int)
    #elev  = np.loadtxt(f"{data_dir}/elevation.csv",delimiter=',')
    #c_t   = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',dtype=int,usecols=(0,))
    c_v   = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(1,))
    
    y_plif = space*0.3/h0
    a_plif = c_v/amp

    w = signal.get_window('boxcar',50)
    w = w/np.sum(w)
    a_plif_smooth = signal.convolve(a_plif,w,'valid')
    y_plif_smooth = signal.convolve(y_plif,w,'valid')
    
    ax.plot(y_plif_smooth,a_plif_smooth,label=f"{amp/h0:.1f}")
    

ax.axhline(1,color='k',linewidth=0.5)
ax.set_ylim(0.8,1.3)
ax.set_ylabel(r'$a/a_0$')

ax = axes[2]

for amp_str, amp in zip(amp_strings,amps):
    data_dir = sorted([d for d in data_directories if x_loc in d and amp_str in d])[0]
    
    time  = np.loadtxt(f"{data_dir}/time.csv",delimiter=',',dtype=int)
    space = np.loadtxt(f"{data_dir}/space.csv",delimiter=',',dtype=int)
    elev  = np.loadtxt(f"{data_dir}/elevation.csv",delimiter=',')
    c_t   = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',dtype=int,usecols=(0,))
    #c_v   = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(1,))
    
    y_plif = space*0.3/h0
    t0 = h0/c0
    a0 = amp
    F = 1 + (1/2)*(a0/h0) - (3/12)*((a0/h0)**2) + (3/56)*((a0/h0)**3)
    alpha = np.sqrt(3*a0/(4*h0))*(1 - (5/8)*(a0/h0) + (71/128)*((a0/h0)**2))

    # quatratic interpolate the arrival time to reduce stairstep artifact
    rng = range(len(space))
    c_idx = c_t - time[0]
    u = elev[c_idx-1,rng]
    c = elev[c_idx  ,rng]
    l = elev[c_idx+1,rng]
    n = u-l
    d = 2*u + 2*l - 4*c
    s = n/d
    s[np.where(np.abs(s)>0.5)] = 0.0 # ignore bad quadratic interps

    t_plif = (c_t/125 + s/125)/(t0/(alpha*F))
    t_plif -= t_plif[-1]

    #print(np.amax((c_t/125 + s/125)/t0) - np.amin((c_t/125 + s/125)/t0))

    w = signal.get_window('boxcar',50)
    w = w/np.sum(w)
    t_plif_smooth = signal.convolve(t_plif,w,'valid')
    y_plif_smooth = signal.convolve(y_plif,w,'valid')
    
    ax.plot(y_plif_smooth, t_plif_smooth, label=f"{a0/h0:.1f}")
    

ax.set_ylim(0,0.7)
ax.set_ylabel(r'$\Delta \tau$')
ax.set_xlabel(r'$y$')

fig.align_ylabels(axes)

#plt.legend(title=r'$a_0$')
axes[1].legend(bbox_to_anchor=(1, 0.5),loc='center left',title=r'$a_0$')

#for label, ax in zip(['a)', 'b)', 'c)'], axes):
#    trans = mtransforms.ScaledTranslation(-35/72, 7/72, fig.dpi_scale_trans)
#    ax.text(0.0, 1.0, label, transform=ax.transAxes + trans,
#            fontsize='9.0', verticalalignment='top', fontfamily='Times New Roman',
#            bbox=dict(facecolor='1.0', edgecolor='none', pad=2.0))

plt.tight_layout()

plt.savefig('Fig_07.pdf')
plt.show()