{Fpe2dTh.PDE  - FPE PDESolns - RGN's advice f: FIRSTPARTS }
title 'Fpe2DTh: 2D FPE  Thesis' {D.Bakewell 18:30 Su/25/2/2001}
SELECT
errlim=8e-5  regrid=on    firstparts
COORDINATES
cartesian (x, y)
VARIABLES 
p     {solute probability or concentration}
DEFINITIONS
a=0.1 b=10.0  {rectangle horiz and vert dimnsions, respctivly}
d=10  d12=d/2 d15=d/5 d10=d/10 dA=0.83 d2=2*d  Vo=0.5
xc=d/2 { x center for this boundary is +ve elctrde right edge}
Fcm=0.04 {Clusius-Msoti Fctr} r=.108 {r=.0465}
Kdiff=0.2424/r  { = 2.245 = KbT*(10^12)/zeta }
Kfpe=(1.637e5)*(r^2)*(Vo^2)*Fcm/(d^3) {nrmlizd valu - see Book III, p.40}
Kss=(5.373e4)*(r^3)*(Vo^2)*Fcm/(d^2)
 xh=(pi*x/d2)+(pi/4)  yh=pi*y/d2 { h = 'hat' } 
ThT=(ARCTAN(sin(xh)/sinh(yh))-ARCTAN(cos(xh)/sinh(yh)))
ThL=(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)
XiS=sinh(yh)*(cos(xh)/Dneg+sin(xh)/Dpos) {Xi sum}
XiD=cosh(yh)*(cos(xh)/Dpos-sin(xh)/Dneg)  {Xi diff}
GamX=2.*ThT*XiS+ThL*XiD GamY=2.*ThT*XiD-ThL*XiS
Vx=Kfpe*GamX  Vy=Kfpe*GamY
Areap=(2*dA)*1  {Areap=pi*(dA^2)/2}  pA=integral(p,2)/Areap
{- steady state p(x,y) from analytical solution and flux monitoring -}
ThSq=4.*ThT^2+ThL^2        pss=EXP(Kss*ThSq) 
 Cint=integral(pss) 	{dble intgrl of p}
 IQx=integral(Vx) IQy=integral(Vy) IThSq=integral(ThSq)
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= -(Kdiff*dx(p)-p*Vx)   {for monitoring Flux in X direction} 
Jy= -(Kdiff*dy(p)-p*Vy)    {for monitoring Flux in Y direction}
J= vector(p*Vx-Kdiff*dx(p),p*Vy-Kdiff*dy(p)) 
INITIAL VALUES
p=1/(d)*1/(b-a) 
EQUATIONS
Kdiff*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 500 by 0.00000001     { time range and initial timestep (secs)}
MONITORS for cycle=1  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(GamX) contour(GamY)
! contour(Vx) contour(Vy) {contour(LOG10(k*E2/Cdiff)) }
 contour (ps ) contour(LOG10(ps)) contour ((p-ps)/ps)
!elevation (log10(GamX)) from (xc,a) to (xc,b)  as "log10(dxE2) X-section"
!elevation (Gamx) from (xc-d12,a) to (xc+d12,a)  as "dxE2 X-section"
!elevation (log10(-GamY)) from (xc,a) to (xc,b)  as "log10(-dyE2) X-section"
!elevation (GamY) from (xc-d12,a) to (xc+d12,a)  as "dyE2 X-section"
HISTORIES history (p) at (xc,a)
 history (pA) report(r) report(Vo) report (Fcm) 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.5V108r0.04F.txt"