# -*- coding: utf-8 -*-
"""
Created on Mon Aug 26 16:33:21 2019

@author: ddesantis
"""

from fipy import *

L = 0.05 #m
H = 0.007 #m
Po = 20*1000 # kg/m/s^2
PL = 10*1000 # kg/m/s^2 
dP = Po-PL

nx = 50
ny=nx

mesh = Grid2D(Lx=L,nx=nx,ny=ny,Ly=H)

rho = 1
mu = 0.013 #kg/m/s
Top_plate_vel = .5 #m/s

#%% - Cell Variables
vx = CellVariable(mesh=mesh,hasOld=True,name = "x - vel")
P = CellVariable(mesh=mesh,hasOld=True,name = "pressure")

#%% Constraints
P.constrain(Po,mesh.facesLeft)
P.constrain(PL,mesh.facesRight)

vx.constrain(Top_plate_vel,mesh.facesTop)
vx.constrain(0,mesh.facesBottom)

#%% Equations

# Continuity Equation
"""
Continuity equation: Mass flow in x direction only, so \nabla(\rho*v) = -\frac{\partial \rho}{\partial t} becomes
\frac{\partial \rho v_x}{\partial x} = -\frac{\partial \rho}{\partial t}
At steady state:
    \frac{\partial \rho v_x}{\partial x} = 0

"""
Eq_C = PowerLawConvectionTerm(coeff = (rho,rho),var = vx) == 0

# Momentum equation 
"""
Set up the diffusion term in the x direction only
"""
Eq_M = PowerLawConvectionTerm(coeff=(rho,rho),var=vx) == -PowerLawConvectionTerm((rho,0),var=P) + DiffusionTerm(coeff=mu,var = vx)

CpldEq = Eq_C & Eq_M

#%% Set up solver

tsd = .1
sweeps = 25
viewer = Viewer(vars=vx,datamin = 0,datamax=1)

for sweep in range(sweeps):
    vx.updateOld()
    P.updateOld()
    CpldEq.sweep(dt=tsd)
    viewer.plot()