Fig_12.ipynb

Kernel: Python 3 (system-wide)

import os

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

import scipy.integrate as integrate
from scipy import signal

from scipy.signal import find_peaks

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

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

def moving_average(x, w):
    return np.convolve(x, np.ones(w), 'same') / w

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

h0, g = 50, 9810
c0 = np.sqrt(g*h0)
t0 = h0/c0

fig, axes = plt.subplots(nrows=1,ncols=2,sharey='all',figsize=(5.0,1.5),layout='constrained')

window = 'bartlett'
window_length = 190

#amp_list = [ f'a{a/10:.1f}cm' for a in [5,10,15,20] ]
amp_list = [ 'a1.5cm' ]
linestyles = ['-','--','-.',':']

for amp,linestyle in zip(amp_list,linestyles):
    # choose only the same amplitude cases
    data_dirs = sorted([d for d in data_directories if amp in d])
    # lexical sort will choose no_flow first and with_flow second

    line_kwargs = {
        'color':'k',
        #'linestyle':linestyle,
        'linewidth':1,
    }

    other_kwargs = {
        'color':'k',
        'linestyle':'-',
        'linewidth':0.5,
    }
    
    # choose the first axis
    ax = axes[0]
    data_dir = data_dirs[0]

    # load data
    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=',')
    t_val = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(0,),dtype=int)
    c_val = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(1,))
    t_idx = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(2,),dtype=int)

    x   = (space*0.3)/h0
    t   = (time/125)/t0
    tv  = (t_val/125)/t0
    eta = elev/h0
    c   = c_val/h0
    a0  = float(amp[len('a'):-len('cm')])*10/h0

    adjust = np.empty(x.shape)
    for i,j in enumerate(t_idx):
        a,b,_ = np.polyfit([(-1/125)/t0,0,(1/125)/t0],eta[j-1:j+2,i],2).tolist()
        adjust[i] = -b/(2*a)

    dt = np.diff(tv+adjust)
    dx = x[1] - x[0]

    F = 1 + (1/2)*a0 - (3/20)*(a0**2) + (3/56)*(a0**3)

    w = signal.get_window(window,window_length)
    w /= np.sum(w)
    speed = np.convolve(dx/dt,w,'same')
    
    ax.plot(x[:-1],speed/F,'k-')
    ax.axhline(1,**other_kwargs)

    # choose the second axis
    ax = axes[1]
    data_dir = data_dirs[1]

    # load data
    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=',')
    t_val = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(0,),dtype=int)
    c_val = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(1,))
    t_idx = np.loadtxt(f"{data_dir}/crests.csv",delimiter=',',usecols=(2,),dtype=int)

    x   = (space*0.3)/h0
    t   = (time/125)/t0
    tv  = (t_val/125)/t0
    eta = elev/h0
    c   = c_val/h0
    a0  = float(amp[len('a'):-len('cm')])*10/h0

    adjust = np.empty(x.shape)
    for i,j in enumerate(t_idx):
        a,b,_ = np.polyfit([(-1/125)/t0,0,(1/125)/t0],eta[j-1:j+2,i],2).tolist()
        adjust[i] = -b/(2*a)

    dt = np.diff(tv+adjust)
    dx = x[1] - x[0]

    F = 1 + (1/2)*a0 - (3/20)*(a0**2) + (3/56)*(a0**3)

    w = signal.get_window(window,window_length)
    w /= np.sum(w)
    speed = np.convolve(dx/dt,w,'same')

    ax.plot(x[:-1],speed/F,'k-')
    ax.axhline(1,**other_kwargs)

for ax in axes:
    ax.set_xlim(-20,20)
    ax.set_ylim(0.7,1.3)
    ax.set_xlabel(r'$x$')

axes[0].set_ylabel(r'$c/F$')

#for label, ax in zip(['a)', 'b)'], axes):
#    ax.set_title(label,loc='left',fontfamily='Times New Roman',fontsize='9.0')

plt.savefig('Fig_12.pdf')