import scipy.io as sio
import numpy as np
import sys
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 26})
plt.interactive(True)

# Channel flow
nu=1/5200
data=np.loadtxt("y_u_k_eps_nut_5200_ymax4.txt")

y=data[:,0]   # cell centers
u=data[:,1]   # velocity
k=data[:,2]   # turbulent kinetic energy
eps=data[:,3] # dissipation
nut=data[:,4] # turbulent viscosity

# Q2/1
# compute the friction velocity for the upper wall.
ymax=4
dy=ymax-y[-1]
dudy=u[-1]/dy
ustar=(nu*dudy)**0.5

kplus=k/ustar**2
yplus=(ymax-y)*ustar/nu

#************
fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.30,bottom=0.20)
plt.plot(kplus,yplus,'b-')
plt.axis([0,5,0,5200])
plt.xlabel('$k^+$')
plt.ylabel('$y^+_{max} - y^+$')
plt.savefig('kplus_yplus.eps')


# Q2/2
#************
fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.30,bottom=0.20)
dudy=np.gradient(u,y)
uv_visc=nu*dudy
uv_turb=-nut*dudy
plt.plot(uv_visc,y,'b-')
plt.plot(-uv_turb,y,'r--')
plt.axis([0,0.3,0,0.01])
plt.xlabel('shear stress')
plt.ylabel('$y$')
plt.savefig('shear-stress.eps')

# they cross between y[7] adnd y[8], i.e. between 0.00363279 and  0.00454428



# Q2/3
#************
fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.30,bottom=0.20)
cmu=0.09
L=cmu**0.75*k**1.5/eps  #see https://www.cfd-online.com/Wiki/Turbulence_length_scale
# L=cmu*k**1.5/eps  # is also correct
plt.plot(L,y,'b-')
plt.axis([0,0.25,0,1])
plt.xlabel('$L$')
plt.ylabel('$y$')
plt.savefig('L.eps')