TITLE     'Constricted Channel with Divergence Term'      { sim213a.pde }
SELECT     errlim=1e-4     spectral_colors 	
VARIABLES     vx     vy     p
DEFINITIONS
   Lx=1.0     Ly=1.0     coef=0.5     visc=1e4
   delp=100							{ Driving pressure }
   dens=1e3     Re=dens*globalmax( vx)*2*Ly/visc	
   v=vector( vx, vy)     vm=magnitude( v)
   unit_x=vector(1,0)     unit_y=vector(0,1)	{ Unit vector fields }
   nx=normal( unit_x)     ny=normal( unit_y)	{Direction cosines }
{ Natural boundary condition for p: }
   natp=visc*[ nx*div( grad( vx))+ ny*div( grad( vy))]
EQUATIONS
   vx:		dx( p)- visc*div( grad( vx))=0
   vy:		dy( p)- visc*div( grad( vy))=0
   p:		div( grad( p))- 1e4*visc/Ly^2*div(v)=0			
BOUNDARIES
region 'domain' start 'outer' (0,Ly)
   natural( vx)=0   value( vy)=0   value( p)=delp     		{ In }   	
   line to (0,-Ly)  value( vx)=0   value( vy)=0   natural( p)=natp
   line to (Lx,-Ly) to (2*Lx,-Ly*coef) to (3*Lx,-Ly*coef)
   natural( vx)=0   value( vy)=0   value( p)=0			{ Out }
   line  to (3*Lx,Ly*coef)  value( vx)=0   value( vy)=0   natural( p)=natp
   line to (2*Lx,Ly*coef) to (Lx,Ly) close
PLOTS
   elevation( nx, ny) on 'outer'	 as 'direction cosines'			
   contour( vx) report(Re)     contour( vm)
   vector( v) norm     contour( p)
   contour( div( v))     contour( curl( v))
   elevation( vx) from (0.5*Lx,-Ly) to (0.5*Lx,Ly)
   elevation( vx) from (2.5*Lx,-Ly*coef) to (2.5*Lx,Ly*coef)
   elevation( natp) on 'outer'
END