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/5000
data=np.loadtxt("y_u_k_eps_8000.dat")

j=data[:,0]
y=data[:,1] # cell centers
u=data[:,2]
k=data[:,3]
eps=data[:,4]

################################## Q2/1
# compute the friction velocity
dudy=np.gradient(u,y)
u_tau=(nu*dudy[0])**0.5

uplus=u/u_tau
yplus=y*u_tau/nu

#************
fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.30,bottom=0.20)
plt.plot(uplus,yplus,'b-')
plt.xlabel('$U^+$')
plt.ylabel('$y^+$')
plt.savefig('uplus_yplus.eps')


################################## Q2/2
# u'v' and u'u'. Use Boussinesq (Eq. 11.33 in the eBook). u'v'=-nu_t*dudy; u'u'=-2*nu_t*dudx+2/3*k.
# the flow is fully developed => dudx=0 =>  u'u'=2/3*k
# the AKN model is used. Compute fmu, see Eq. F.3 in eBOOK

ueps=(eps*nu)**0.25
# compute distance to closest wall
dist=np.minimum(y,y[-1]-y)
ystar=ueps*dist/nu
rt=k**2/eps/nu
fmu=((1.-np.exp(-ystar/14.))**2)*(1.+5./rt**0.75*np.exp(-(rt/200.)**2))

# compute nu_t
nu_t=fmu*0.09*k**2/eps

# compute u'v'
dudy=np.gradient(u,y)
uv=-nu_t*dudy

# compute u'u'
uu=2/3*k

fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.20,bottom=0.20)
plt.plot(uv,y,'b-')
plt.plot(uu,y,'r--')
plt.xlabel("$\overline{u'v'}, \quad \overline{u'u'}$")
plt.ylabel('$y$')
plt.savefig('uv-uu.eps')


################################# Q2/3. The two largest terms in the k eq are the production term and the dissipation term
prod=nu_t*dudy**2
fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.20,bottom=0.20)
plt.plot(prod,yplus,'b-')
plt.plot(-eps,yplus,'r--')
plt.xlabel(r"$P^k, \quad \varepsilon$")
plt.axis([-2000,2000,0,200])
plt.ylabel('$y$')
plt.savefig('prod-diss.eps')