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

plt.interactive(True)

nu=1/200

y_u=np.loadtxt("y_u.dat")
x_2=y_u[:,0]
v_1=y_u[:,1]

################################## Q1/1
#************
# velocity profile plot
fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.25,bottom=0.20)
# plot the velocity profile 
plt.plot(v_1,x_2,'b-')
plt.xlabel("$v_1$")
plt.ylabel("$x_2$")

################################## Q1/2
# bulk velocity
# integrate
u_int=np.trapz(v_1,x_2)
# divide by x_2_max
u_bulk=u_int/x_2[-1]
# Reynolds number:
Re_H=u_bulk*x_2[-1]/nu
print('ubulk',u_bulk,' Re_bulk',Re_H)
# Here I create the Reynolds number with channel width. You may also create it
# with the half  channel widthm i.e. Re_h = u_bulk*x_2[-1]/2/nu/

################################## Q1/3
# wall shear stress 
# 
dudy=np.gradient(v_1,x_2)

# south wall
tauw_s=nu*dudy[0]

# north wall
tauw_n=nu*dudy[-1]

print('wall shear south wall',tauw_s)
print('\nwall shear north wall',tauw_n)

################################## Q1/4
#
# incompressible flow => diss = 2*mu*s_ij*s_ij 
# Fully developed channel flow: dudx=v=0 => s11=s22=0
# => diss = 2*mu*(s_12*s_12 + s_21*s_21)  
# rho=1 => mu=nu
#
s12=0.5*dudy
s21=s12
diss=2*nu*(s12**2+s21**2)

# plot dissipation
fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.25,bottom=0.20)
plt.plot(diss,x_2,'b-')
plt.xlabel(r"$\varepsilon$")
plt.ylabel("$x_2$")