Fig_11.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
from matplotlib.lines import Line2D

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 solitary_amp(a0,h0,x,t,xoff=0,toff=0,g=9.81):
    '''Returns the amplitude of a 3rd order solitary wave on a meshgrid or as scalar'''
    assert isinstance(a0, (int, float)) and  isinstance(h0, (int, float)) \
        and isinstance(x, (int, float, np.ndarray)) and isinstance(t, (int, float, np.ndarray)) \
        and isinstance(toff, (int, float)) and isinstance(g, (int, float))
    
    from numpy import sqrt
    from numpy import cosh
    
    a = a0/h0
    c0 = sqrt(g*h0)
    F = 1 + (1/2)*a - (3/20)*(a**2) + (3/56)*(a**3)
    c = c0*F
    alpha = sqrt(3*a/4) * (1 - (5/8)*a + (71/128)*(a**2))
    xi = x/h0
    xioff = xoff/h0
    tau = c0*t/h0
    tauoff = c0*toff/h0
    X, T = np.meshgrid(xi,tau,indexing='ij')
    if X.size == 1 and T.size == 1:
        X = X.ravel()[0]
        T = T.ravel()[0]
    s = 1/(cosh(alpha*((X-xioff) - F*(T - tauoff))))
    return h0*(a*(s**2) - (3/4)*(a**2)*(s**2 - s**4) + (a**3)*( (5/8)*(s**2) - (151/80)*(s**4) + (101/80)*(s**6) ))


def solitary_vel(a0,h0,x,t,xoff=0,toff=0,g=9.81):
    '''Returns the flow velocity of a 3rd order solitary wave on a meshgrid or as scalar'''
    assert isinstance(a0, (int, float)) and  isinstance(h0, (int, float)) \
        and isinstance(x, (int, float, np.ndarray)) and isinstance(t, (int, float, np.ndarray)) \
        and isinstance(toff, (int, float)) and isinstance(g, (int, float))
    from math import sqrt
    eta = solitary_amp(a0,h0,x,t,xoff,toff,g)
    a = a0/h0
    F = 1 + (1/2)*a - (3/20)*(a**2) + (3/56)*(a**3)
    c = sqrt(g*h0)*F
    return (c*eta)/(h0+eta)

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

amp = 'a1.5cm'

# choose only the same amplitude cases
data_dirs = [d for d in data_directories if amp in d]
# lexical sort will choose no_flow first and with_flow second

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

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

step = 50
N_step = 3

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

print("Step between lines is: ",(step/125)/t0)

# 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
eta = elev/h0

x_zero_idx = np.argwhere(space==0).ravel()
t_zero_idx = t_idx[x_zero_idx][0]

start = t_zero_idx - step*N_step
stop  = t_zero_idx + step*(N_step+1)

ax.axhline(0,**line_kwargs)

for idx in range(start,stop,step):
    ax.plot(x,eta[idx,:],**line_kwargs)

compare_idx = t_zero_idx - step*(N_step-1)
# need to compute xoff
# which is the location of the crest at compare_idx
x_idx = np.argmax(eta[compare_idx,:])
xoff = x[x_idx]
# compute theoretical shape at xoff
theory = solitary_amp(c_val[compare_idx],h0,x*h0,0,xoff=xoff*h0,toff=0,g=9.81)/h0
ax.plot(x,theory,'k--',linewidth=1)
    
# 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
eta = elev/h0

x_zero_idx = np.argwhere(space==0).ravel()
t_zero_idx = t_idx[x_zero_idx][0]

start = t_zero_idx - step*N_step
stop  = t_zero_idx + step*(N_step+1)

ax.axhline(0,**line_kwargs)

for idx in range(start,stop,step):
    ax.plot(x,eta[idx,:],**line_kwargs)

compare_idx = t_zero_idx - step*(N_step-1)
# need to compute xoff
# which is the location of the crest at compare_idx
x_idx = np.argmax(eta[compare_idx,:])
xoff = x[x_idx]
# compute theoretical shape at xoff
# c_val[compare_idx]
theory = solitary_amp(eta[compare_idx,x_idx]*h0,h0,x*h0,0,xoff=xoff*h0,toff=0,g=9.81)/h0
ax.plot(x,theory,'k--',linewidth=1)

# combined paramaters
xlim = (-20,20)
ylim = (-0.07,0.47)

#plt.xlabel(r'$x$')
axes[0].set_ylabel(r'$\eta$')

for ax in axes:
    ax.set_xlim(*xlim)
    ax.set_ylim(*ylim)
    ax.set_yticks([0,0.1,0.2,0.3,0.4])
    ax.set_xlabel(r'$x$')
    #ax.set_ylabel(r'$\eta$')

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

Step between lines is:  5.602856414365801