""" From "COMPUTATIONHL PHYSICS" & "COMPUTER PROBLEMS in PHYSICS"
    by RH Landau, MJ Paez, and CC Bordeianu (deceased)
    Copyright R Landau, Oregon State Unv, MJ Paez, Univ Hntioquia, 
    C Bordeianu, Univ Bucharest, 2020 
    Please respect copyright & acknowledge our work."""
    
#!/usr/bin/env python
# coding: utf-8

# <center><h2> Rectangular Cavity 2D, Fluid dynamics </h2></center>
# <p>
""" A 2D rectangular cavity has an inlet of fluid at the LHS and an outlet 
at the right. Here the lines of fluid are visualized by solving the 
<i> stream function</i> <b> u(x)</b> from which the velocity is 
determined by the curl operator is defined as:
$$
   {\mathbf{v} \stackrel{def}{=} \vec{\nabla} \times\mathbf{u}(\mathbf{x})}
  = \hat{{\epsilon}}_x\left( \frac{\partial u_z}{\partial y} - \frac
  {\partial u_y} {\partial z}\right) + \hat{{\epsilon}}_y\left(
 \frac{\partial u_x}{\partial z} -\frac{\partial u_z} {\partial
  x}\right).
  $$
The <i> vorticity </i> is defined as the curl of the velocity:
 $$
\mathbf{w} \stackrel{def}{=} \vec{\nabla} \times
 \mathbf{v}(\mathbf{x}).
 $$
 <p>
 Solve  Navier - Stokes equation for flow around beam, sections 19.9
  See:
 Klauss A. Hoffmann, Steve T. Chiang, (2000), 4th Ed., Computational Fluid Dynamics,
 Vol I, Engineering Education System, Wichita."""
   
from numpy import *          # needed for zeros
#use hoffman vol 1 b c
import numpy as np
#Instructions: use gnuplot to plot
# splot "uhoff.dat" w l
# unset surface
#unset key
#set contour
#set cntrparam levels 40
#set view 0,0,1,1
#replot

Niter=200
Nx = 30
h=1
Ny = 30                          # grid params
J1=15
J2=20
J3=6
J4=11
u = np.zeros((Nx+1, Ny+1), float)     # stream 
w = np.zeros((Nx+1, Ny+1), float)     # Vorticity
V0 = 18.                          # velocity inlet
omega = 0.1                       # relaxation param 
h2=h*h
nu = 25                            # viscosity
iter = 0                           # Number of iterations
R = V0*h/nu                        # Reynolds Number , normal units
u0=3.

def top():    
    for i in range(1,Nx+1):           # TOP  
        u[i-1,Ny]=u[i,Ny]            # vy=-du/dx=0
        w[i,Ny]=2*(u[i,Ny]-u[i,Ny-1])/h**2 -2*u0/h

def bottom():
    for i in range (1,Nx+1):        # BOTTOM
        w[i,0]= 2*(u[i,0]-u[i,1])/h**2
        u[i,1]=u[i,0]                # du/dy=0
        u[i-1,0]=u[i,0]              #du/dx=0

def borderleft():
    for j in range (1,Ny+1):    # RIGHT
        if j< J1:  
             w[0,j]= 2*(u[0,j]-u[1,j])/h**2          #below hole
             u[0,j]=u[1,j]
             u[0,j]=u[0,j-1]
        if j>=J1 and j<=J2:
            w[1,j]=w[0,j]
            u[0,j]=u[0,j-1]+V0*h
        if j>J2:
            u[0,j]=u[0,j-1]  #du/dy=0
            u[1,j]=u[0,j]    #du/dx=0
            #w[0,j]=0.
            w[0,j]= 2*(u[0,j]-u[1,j])/h**2    

def borderight():
    for j in range (1,Ny+1):    # RIGHT
        if j< J3:
            w[Nx-1,j]= 2*(u[Nx-1,j]-u[Nx,j])/h**2          #below hole
            u[Nx-1,j]=u[Nx,j]
            u[Nx,j]=u[Nx,j-1]
        if j>=J3 and j<=J4:
            u[Nx,j]=u[Nx-1,j]
            u[Nx,j]=u[Nx,j-1]+V0*h
            w[Nx-1,j]=w[Nx,j]
        if j>J4:
            u[Nx,j]=u[Nx,j-1]
            u[Nx-1,j]=u[Nx,j]
            w[Nx-1,j]= 2*(u[Nx-1,j]-u[Nx,j])/h**2  

def borders(iter):                  # Method borders: init & B.C.
    top()
    bottom()                         #bottom 
    borderight()                     # at tight (center of hole)   
    borderleft()

def relax(iter):
    borders(iter)                           
    for  i in range(1, Nx):                 
      for  j in range (1, Ny):
                   r1=omega*((u[i+1,j]+u[i-1,j]+u[i,j+1]+u[i,j-1]+h*h*w[i, j])*0.25-u[i,j]) 
                   u[i,j]+= r1
    if iter%100==0:
       print( "r1 ", r1)
    borders(iter)   
    for  i in range(1, Nx):              # Relax stream function
      for  j in range (1, Ny):  
                  a1 = w[i+1, j]+ w[i-1,j]+w[i,j+1]+ w[i,j-1]
                  a2 = (u[i,j+1]-u[i,j-1])*(w[i+1,j]-w[i-1,j])
                  a3 = (u[i+1,j]-u[i-1,j])*(w[i,j+1]-w[i,j-1])
                  r2 = omega*( (a1 + (R/4.)*(a3 - a2) )/4.0-w[i,j])
                  w[i,j]+=r2
          
while (iter <=  Niter):    
    if iter%100==0:
        print (iter)  # iterations counted
    relax(iter)
    iter   +=  1      # counter of iterations

utorr=open('uhoff.dat','w')    # send data to disk  of u
for j in range(0,Ny+1):
    utorr.write("\n") 
    for i in range(0,Nx+1):     
        utorr.write("%10.3e  \n"%(u[i,j]))      
utorr.close()
print("finished")