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

# Channel flow
data=np.genfromtxt("channel_flow_data.dat", comments="%")

ni=199  # number of grid nodes in x_1 direction
nj=28  # number of grid nodes in x_2 direction
x1=data[:,0] #don't use this array
x2=data[:,1] #don't use this array
v1=data[:,2] #don't use this array
v2=data[:,3] #don't use this array
p=data[:,4]  #don't use this array

# transform the arrays from 1D fields x(n) to 2D fields x(i,j)
# the first index 'i', correponds to the x-direction
# the second index 'j', correponds to the y-direction

x1_2d=np.reshape(x1,(nj,ni)) #this is x_1 (streamwise coordinate)
x2_2d=np.reshape(x2,(nj,ni)) #this is x_2 (wall-normal coordinate)
v1_2d=np.reshape(v1,(nj,ni)) #this is v_1 (streamwise velocity component)
v2_2d=np.reshape(v2,(nj,ni)) #this is v_1 (wall-normal velocity component)
p_2d=np.reshape(p,(ni,nj))   #this is p   (pressure)

x1_2d=np.transpose(x1_2d)
x2_2d=np.transpose(x2_2d)
v1_2d=np.transpose(v1_2d)
v2_2d=np.transpose(v2_2d)
p_2d=np.transpose(p_2d)


########################## Q1/1.

# compute the bulk velocity. We know that due to contiuity it should be constant and equal to the inlet velocity = v1_2d[0,:]=0.9
# anbyway , let;s compute it. See Eq. F.1
ymax=x2_2d[1,-1]
ubulk=np.trapz(v1_2d,x2_2d,axis=1)/ymax

#************
fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.25,bottom=0.20)
# I plot is although I know is should be constant due to continuity. I find that it varies slightly due to numerical errors
plt.plot(x1_2d[:,1],ubulk,'b-')
plt.xlabel('$x_1$') 
plt.ylabel(r'$u_{bulk}}$') 
plt.savefig('ubulk.eps')


########################## Q1/2.
##  wall shear stress on the upper wall
dudy=np.gradient(v1_2d,x2_2d[1,:],axis=1)
# air of 20 degrees C: 
mu=18.2e-6
tauw=mu*dudy[:,-1]   # index -1 is the last element


#************
fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.25,bottom=0.20)
plt.plot(x1_2d[:,1],abs(tauw),'b-')
plt.xlabel('$x_1$') 
plt.ylabel(r'$\tau_w$') 
plt.savefig('tau.eps')

############################  Q1/3. vorticity vector on the upper wall
# it is defined as (eq. 1.13 in eBook): omega_3 = dv/dx-dudx. omega_1=omega_2=0 because the flow is 2D
dvdx=np.gradient(v2_2d,x1_2d[:,1],axis=0) # we know that at the wall, this is zero (but I compute it anyway)
omega_3=dvdx-dudy

fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.25,bottom=0.20)
plt.plot(x1_2d[:,:],omega_3[:,-1],'b-')
plt.xlabel('$x_1$') 
plt.ylabel(r'$\omega_3$') 
plt.savefig('omega_3.eps')