Fig_04.ipynb

Kernel: Python 3 (system-wide)

import os

import numpy as np
import matplotlib.pyplot as plt

from scipy import signal
from scipy.signal import convolve2d
from scipy.stats import linregress

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")

h, g, b0 = 50, 9810, 229
c = np.sqrt(g*h)

x = (x_int/b0)
y = (y_int/b0)

u /= c
v /= c
m /= c

u_std /=c
v_std /= c
m_std /= c

N = 5
k = np.ones((N,N))/(N**2)

kwargs = {'mode':'same','boundary':'fill','fillvalue':np.nan}

u_smooth = convolve2d(u,k,**kwargs)
v_smooth = convolve2d(v,k,**kwargs)
m_smooth = convolve2d(m,k,**kwargs)

u_std_smooth = convolve2d(u_std,k,**kwargs)
v_std_smooth = convolve2d(v_std,k,**kwargs)
m_std_smooth = convolve2d(m_std,k,**kwargs)

def find_nearest(arr,val):
    idx = np.abs(arr - val).argmin()
    return idx, arr[idx]

x_vals = np.arange(-2200,-15,5,dtype=int)
b   = np.empty(x_vals.shape)
u_m = np.empty(x_vals.shape)

y_val_1 = 0
y_val_2 = 700

for idx, x_val in enumerate(x_vals):
    mask1 = np.logical_and(x_int==x_val,y_int>y_val_1)
    mask2 = y_int<y_val_2
    mask = np.logical_and(mask1,mask2)
    
    u_m[idx] = np.abs(u_smooth[mask]).max()
    b_idx, _ = find_nearest(u_smooth[mask],-u_m[idx]/2)
    b[idx] = y[mask][b_idx]
    
mask = x_vals < -1000

slope, intercept, r, p, se = linregress(x_vals[mask], b[mask])
fit = slope*x_vals + intercept
#print(slope, intercept)

mm_to_inch = 1/25.4

b0 = 229

w = signal.get_window('boxcar',15)
w = w/np.sum(w)
x_vals_smooth = signal.convolve(x_vals,w,'valid')
b_smooth = signal.convolve(b,w,'valid')

fig, ax = plt.subplots(figsize=(5.0,1.85))

y_val_1 = 0
y_val_2 = 700

x_vals = [-20,-200,-1000,-2000]
linestyles = ["-","--","-.",":"]

# plot theory
y_fit = np.linspace(0,3)
u_fit = np.exp(-np.log(2)*(y_fit)**2)
ax.plot(y_fit,u_fit,'r-')

# plot axis
ax.axhline(0,color='k',linewidth=0.5)

for x_val, linestyle in zip(x_vals,linestyles):
    # find B value from precalculated
    x_val_b = np.arange(-2200,-15,5,dtype=int)
    b_idx = np.argwhere(x_val==x_val_b)
    B = b[b_idx].flat
    #B = 230
    
    # construct a mask
    mask1 = np.logical_and(x_int==x_val,y_int>y_val_1)
    mask2 = y_int<y_val_2
    mask = np.logical_and(mask1,mask2)

    # find centerline velocity
    u_m_idx = np.argmax(np.abs(u_smooth[mask]))
    u_m = u_smooth[mask][u_m_idx]
    
    #print(x_val/b0,u_m)
    
    # plot normalized
    ax.plot(y[mask]/B,u_smooth[mask].ravel()/u_m,'k',linestyle=linestyle,label=f"{x_val/h:.2g}")

ax.set_xlim(0,3)
ax.set_ylim(-0.05,1.05)
ax.set(xlabel=r'$\beta$',ylabel=r'$u/u_m$')

plt.legend(title=r'$x$')

plt.tight_layout()

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