TITLE     'Viscous Flow past a Circular Cylinder'	         { sim221.pde }
SELECT     errlim=1e-3     ngrid=1     spectral_colors     stages=2
VARIABLES     vx      vy       p
DEFINITIONS
   Lx=2.0     Ly=1.0     a=0.2     visc=1e4
   delp=100								{ Driving pressure }
   dens=1e3     Re=dens*globalmax( vx)*2*Lx/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 }
   natp=
        if stage=2 then visc*[ nx*div( grad( vx))+ ny*div( grad( vy))]  else 0
EQUATIONS
   vx:		dx( p)- visc*div( grad( vx))=0
   vy:		dy( p)- visc*div( grad( vy))=0
   p:		div( grad( p))- 1e4*visc/Ly^02*div( v)=0
BOUNDARIES
region 'domain'  start 'outer' (-Lx,Ly)
   natural( vx)=0   natural( vy)=0    value(p)=delp			{ In }
   line to (-Lx,-Ly)   value( vx)=0   value( vy)=0   natural(p)=natp { Wall}
   line to (Lx,-Ly)  natural( vx)=0   natural( vy)=0   value(p)=0  	{ Out }
   line to (Lx,Ly)   value( vx)=0   value( vy)=0   natural(p)=natp	{ Wall }
   line to close
   start 'outline' (a,0)						{ Exclude cylinder }
   value( vx)=0  value( vy)=0   natural(p)=natp
   arc( center=0,0) angle=360 close
PLOTS
   contour( vx) report( Re)     contour( vy)
   contour( vm) painted     contour( p)
   vector( v) norm     vector(v) norm  zoom(-2*a,-2*a, 4*a,4*a)
   contour( div( v))     contour( curl( v)) painted
   elevation( vx, vy) from (-Lx,0) to (Lx,0)
END