TITLE 'steady newtonian Fluid'    
COORDINATES cartesian2

SELECT busymon   ! forces monitors to be displayed on each newton iteration

DEFINITIONS
    U0 = 0.1

VARIABLES       
  	u
	v(U0/100)
	p        

DEFINITIONS    
	L=5 	!length of the channel
	H=1    !height of the channel
	dens=10   !density of the fluid
	visc=0.1   !dyn. viscosity 

	speed=sqrt(u^2+v^2)

	penalty=1000

	Re=globalmax(speed)*H*dens/visc	

INITIAL VALUES
	u=0
	v=0
	p=0

EQUATIONS        
  	u:	visc*div(grad(u))-dx(p)=dens*(u*dx(u)+v*dy(u))
	v:	visc*div(grad(v))-dy(p)=dens*(u*dx(v)+v*dy(v))
	p: div(grad(p))=penalty*(dx(u)+dy(v))


BOUNDARIES       
  REGION 1       
    START(0,0)  
   		value(u)=0 value(v)=0 load(p)=0   
	LINE TO (L,0)
		load(u)=0 {load(v)=0}  value(p)=0 value(v)=0 !! nobc(p)
	line to (L,H)
		value(u)=0 value(v)=0 load(p)=0 
	Line to (0,H)
		value(u)=U0    !*y*(H-y)*(6/H^3)    ! parabolic input velocity 
		value(v)=0 nobc(p)
	line to close 

MONITORS     
	contour (u) 	report (Re)
	contour (v) 	report (Re)
	contour (p) 	report (Re)
	contour (speed) 	report (Re)
	elevation(p) from(0,0) to (0,H)
	elevation(p) from(0,H/2) to (L,H/2)
	elevation(u) from(L/2,0) to (L/2,H)
	elevation(visc*div(grad(u)), dx(p), dens*(u*dx(u)+v*dy(u))) from(0,H/2) to (L,H/2)
	elevation(visc*div(grad(v)), dy(p), dens*(u*dx(v)+v*dy(v))) from(0,H/2) to (L,H/2)
PLOTS            
  	grid(x,y)
	contour(u) report(Re)
	contour(v) report(Re)
	contour(speed) painted report(Re)
	vector(u,v)	as "flow"	report(Re)
	contour(p) as "pressure" painted
	contour(dx(u)+dy(v))as"continuity error"
	elevation(u) from(0,0) to (0,H)
	elevation(u) from(L,0) to (L,H)
	
	elevation(u) from(L/2,0) to (L/2,H)
	elevation(p) from(0,H/2) to (L,H/2)


END