{fpe_2D.PDE  - FPE PDESolns - RGN's advice f: FIRSTPARTS }
title 'time-dependent 2D diffusion' {D.Bakewell 16:45 we/10/10/2000}
SELECT
errlim=8e-5  regrid=on    firstparts
COORDINATES
cartesian (x, y)
VARIABLES 
p     {solute probability or concentration}
DEFINITIONS
a=0.1 b=50.1  {rectangle horiz and vert dimnsions, respctivly}
d=10  d12=d/2 d15=d/5 d10=d/10 dA=0.83 d2=2*d  Vo=1
xc=d/2 { x center for this boundary is +ve elctrde right edge}
ktz=298*1.38/(6.00*pi*9.00*1.08)  { = 2.245 = KbT*(10^12)/zeta }
Fcm=0.737 {Clusius-Msoti Fctr} k=1510*Fcm {nrmlizd valu - see Book II}
Cdiff=ktz   {pi already defined} th=5000 {field trucation threshold}
 xh=(pi*x/d2)+(pi/4)  yh=pi*y/d2 { h = 'hat' } 
Ex=(2.*Vo/(pi*d))*(ARCTAN(sin(xh)/sinh(yh))-ARCTAN(cos(xh)/sinh(yh)))
Ey=(Vo/(pi*d))*(LN(cosh(yh)+cos(xh))-LN(cosh(yh)-cos(xh))+
  LN(cosh(yh)+sin(xh))-LN(cosh(yh)-sin(xh)))
Dpos=cosh(2*yh)+cos(2*xh)  Dneg=cosh(2*yh)-cos(2*xh)
Exy=2.*Vo*cosh(yh)/(d^2)*(cos(xh)/Dpos-sin(xh)/Dneg)  Eyx=Exy
Exx=2.*Vo*sinh(yh)/(d^2)*(cos(xh)/Dneg+sin(xh)/Dpos)  Eyy=-Exx
dxE2=2*(Ex*Exx+Ey*Eyx) dyE2=2*(Ex*Exy+Ey*Eyy)
Ux=k*dxE2 Uy=k*dyE2
Vx=th*SIGN(Ux)*USTEP(ABS(Ux)-th)+Ux*USTEP(th-ABS(Ux))
Vy=-th*USTEP(-Uy-th)+Uy*USTEP(Uy+th) {Areap=pi*(dA^2)/2}
Areap=(2*dA)*1  pA=integral(p,2)/Areap
{- steady state p(x,y) from analytical solution and flux monitoring -}
E2=Ex^2+Ey^2        pss=EXP(k*E2/Cdiff) 
 Cint=integral(pss) 	{dble intgrl of p}
 IQx=integral(Vx) IQy=integral(Vy) IE2=integral(E2/2)
ps=pss/Cint  psA=integral(ps,2)/Areap  {pInit=initial p} 
pInit=1/(d)*1/(b-a)  pIA=integral(pInit,2)/Areap 
pAtau=(psA-pIA)*(1.-exp(-1.))+pIA 
DpARel=(psA-pIA)/pIA
Jx= -(Cdiff*dx(p)-p*Vx)   {for monitoring Flux in X direction} 
Jy= -(Cdiff*dy(p)-p*Vy)    {for monitoring Flux in Y direction}
J= vector(p*Vx-Cdiff*dx(p),p*Vy-Cdiff*dy(p)) 
INITIAL VALUES
p=1/(d)*1/(b-a) 
EQUATIONS
Cdiff*div[grad(p)]-dx[p*Vx]-dy[Vy*p]=dt(p)      { 2D FPE equation }
BOUNDARIES
Region 1   natural(p) = 0  {walking outer boundary}
 start (xc-d12,a) line to (xc+d12,a) to (xc+d12,b) to (xc-d12,b) to finish 
 {start (xc-d,a) line to (xc+d,a) to (xc+d,b) to (xc-d,b) to finish}
Region 2 start  'Aint1' (xc,a) line to (xc-dA,a) to (xc-dA,a+1)
 to (xc+dA,a+1) to (xc+dA,a) to finish 
{arc (center=xc,a)  to (xc+dA,a) line to finish}
FEATURE
start (xc-.1,a) line to (xc-.09,a) to (xc-.08,a) to (xc-.07,a) to (xc-.06,a)
 to (xc-.05,a) line to (xc-.04,a) to (xc-.03,a) to (xc-.02,a) to (xc-.01,a)
 to (xc,a) to (xc+.01,a) to (xc+.02,a) to (xc+.03,a) to (xc+.04,a) to
 (xc+.05,a) to (xc+.06,a) to (xc+.07,a) to (xc+.08,a) to (xc+.09,a) to (xc+.1,a)
TIME 0 to 50 by 0.00000001     { time range and initial timestep (secs)}
MONITORS for cycle=10  contour (p) contour (LOG10(p)) 
{elevation (p) from (xc,a) to (xc,b)  as "X-section"}
elevation (LOG10(p)) from (xc,a) to (xc,b)  as "X-section"
elevation (LOG10(ps)) from (xc,a) to (xc,b)  as "X-section" 
contour(dxE2) contour(dyE2)
 contour(Ux) contour(Uy) {contour(LOG10(k*E2/Cdiff)) }
contour (ps ) contour(LOG10(ps)) contour ((p-ps)/ps)
elevation (log10(dxE2)) from (xc,a) to (xc,b)  as "log10(dxE2) X-section"
elevation (dxE2) from (xc-d12,a) to (xc+d12,a)  as "dxE2 X-section"
elevation (log10(-dyE2)) from (xc,a) to (xc,b)  as "log10(-dyE2) X-section"
elevation (dyE2) from (xc-d12,a) to (xc+d12,a)  as "dyE2 X-section"
HISTORIES history (p) at (xc,a)
 history (pA) report(k) report(Vo) report(pIA) report(psA) report(pAtau) 
report(DpARel)
PLOTS     { plot prints a tag-delimited data list to the file "filename.TXT":}
 history(pA) at (a,0) export format "#t#b#" file="c:\windows\desktop\0.45V0.5M1.66.txt"