Fig_13.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 find_idx_nearest(x,v):
    return np.argmin(np.abs(x-v))

def smooth(c):
    w = signal.windows.get_window('hann',100)
    w = w/np.sum(w)
    return signal.convolve(c,w,mode='same')

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

mosaic = """
AABB
.CC.
"""
fig, axdict = plt.subplot_mosaic(mosaic,figsize=(4.4,3.0),layout="tight")

axes = (axdict['A'],axdict['B'],axdict['C'])

for ax in axes:
    ax.axhline(1,linewidth=0.5,color='k')

amp_list = [ f'a{a/10:.1f}cm' for a in [5,10,15,20] ]
linestyles = ['-','--','-.',':']
colorstyles = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']

for amp,linestyle,color in zip(amp_list,linestyles,colorstyles):
    # 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

    line_kwargs = {
        'color':color,
        'linestyle':'-',
        'linewidth':1,
    }
    
    # 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)
    space1 = 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_val1 = 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   = (space1*0.3)/h0
    t   = (time/125)/t0
    a0  = float(amp[len('a'):-len('cm')])*10/h0
    
    xoff_idx = find_idx_nearest(x,-12)
    aoff = c_val1[xoff_idx]
    
    eta = elev/aoff
    c   = c_val1/aoff # - 1
    
    ax.axhline(0,color='k',linewidth=0.5)
    ax.plot(x,c,**line_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)
    space2 = 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_val2 = 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   = (space2*0.3)/h0
    t   = (time/125)/t0
    a0  = float(amp[len('a'):-len('cm')])*10/h0
    
    xoff_idx = find_idx_nearest(x,-12)
    aoff = c_val2[xoff_idx]
    
    eta = elev/aoff
    c   = c_val2/aoff # - 1
    
    ax.axhline(0,color='k',linewidth=0.5)
    ax.plot(x,smooth(c),**line_kwargs)
    
    # choose the last axis
    ax = axes[2]
    
    xvalid, xidx1, xidx2 = np.intersect1d(space1,space2,return_indices=True)
    ratio = c_val2[xidx2]/c_val1[xidx1]
    
    ax.plot(xvalid*0.3/h0, smooth(ratio),**line_kwargs)


# combined paramaters
xlim = (-14,14)
ylim = (0.95, 1.30)

axes[0].set_ylabel(r'$a/a_i$')
axes[2].set_ylabel(r'$a_w/a_n$')

axes[0].set_ylim(*ylim)
axes[1].set_ylim(*ylim)
axes[2].set_ylim(ymin=1)

for ax in axes:
    ax.set_xlim(*xlim)
    ax.set_xticks([-10,-5,0,5,10])
    ax.set_xlabel(r'$x$')

lines  = [ Line2D([0],[0],color=color,linestyle='-') for color in colorstyles ]
labels = [ f'{a:.1f}' for a in [0.1,0.2,0.3,0.4] ]
axes[2].legend(lines,labels,title=r'$a_0$',bbox_to_anchor=(1, 0.5),loc='center left')

#for label, ax in zip(['a)', 'b)','c)'], axes):
#    trans = mtransforms.ScaledTranslation(-35/72, 0, 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.savefig(f'Fig_13.pdf')
plt.show()