{fpe_1D_Re.PDE   FPE Relaxation 1-D model PDE Solutions Inc.}
title 'FPE_1D_Re' {lastest update D.Bakewell 16:45 Sa/2/12/2000}
SELECT
errlim=8e-5  regrid=on  firstparts  {RGN's advice f: FIRSTPARTS}
COORDINATES
cartesian (x, y)
VARIABLES
p     {solute probability or concentration}
DEFINITIONS
a=0 b=80 {rectangle horiz & vert dimensions} 
c=0.5   {rectangle 1/2 width} 
ktz=298*1.38/(6.00*pi*9.00*1.08) { = 2.245 = KbT*(10^12)/zeta = Cdiff}
pas=1./(b-a)  ps=1./(b-a)
 { ------- hyperbolic force prfile Book I p.    --- }
!Qx=0   Qy=0           {driving velocity in X & Y directions}
 { ------ exponential force profile Book I p. 92 - 93; 132 - 5 ---- }
  Cv=15.38 {normalized valu - see Book II}
 d=10  d2=2*d  Vo=0.1 k1=Cv*(Vo^2) k2=pi/d 
  Qx=0 Qy=0        {driving velocity in X & Y directions}
  pLint=BINTEGRAL(p, 'Lint')/(2*0.01)
 { --------- steady state p(x,y) analytical solution --------- }
 { Relaxation Initial Condition hyperbolic Book I p.    --- }
 !Kratio=0.02         	{DEP/thermal energy ratio=(k/zeta)/(KbT/zeta)  }
 !al2=1.-Kratio {2*alpha=(k/zeta)/(KbT/zeta)}
 !pa0= (al2*a^(al2-1.) )/ (b^al2-a^al2)
 { ---  steady state: exponential Book I p.    ; Book II p.    }
 beta=k1/(k2*ktz) 
 pss=exp(beta*exp(-k2*y))
 Cexp=EXPINT(beta*exp(-k2*b)) - EXPINT(beta*exp(-k2*a))
 p0= -k2*pss/Cexp 
 pa0=-k2*exp(beta*exp(-k2*a))/Cexp
 Cint=integral(pss)/(2*c)
 !p0=-k2*pss/Cint {report dif btw Cint & -Cexp/k2 } 
 { --- steady state Constant force e.g. gravity --- }
 ! Cexp=exp(-b*Vbgz/ktz)-exp(-a*Vbgz/ktz)
 ! ps=(-Vbgz/ktz)*exp(-y*Vbgz/ktz)/Cexp
 ! pas=(-Vbgz/ktz)*exp(-a*Vbgz/ktz)/Cexp
 { -- parameters common to all steady state -- } 
  pTau=(pas-pa0)*(1.-exp(-1.))+pa0
  IQx=integral(Qx) IQy=integral(Qy)
{ ------------ flux monitoring ------------- }
Fluy1= ktz*dy(p)-p*Qy   {for monitoring Flux in Y direction}
Flux1= ktz*dx(p)-p*Qx   {for monitoring Flux in X direction} 
 J= vector(p*Qx-ktz*dx(p),p*Qy-ktz*dy(p)) 
INITIAL VALUES
! p = (al2*y^(al2-1.) )/ (b^al2-a^al2) {frm hyperbolic steady state}
 p=p0		     {frm exp collection steady state}
EQUATIONS
ktz*div[grad(p)]-dx[p*Qx]-dy[Qy*p]=dt(p)      { 2D FPE equation }
BOUNDARIES
Region 1
natural(p) = 0
start (-c,a) line to (c,a) to (c,b) to (-c,b) to finish {walk outer boundary }
Region 2
start  'Lint' (-0.01,a) line to (0.01,a)
FEATURE
start (-c,a) line to (-c+0.1,a) to (-c+0.2,a) to (-c+0.3,a) to (-c+0.4,a) to (-c+0.5,a)
start (0,a) line to (0.1,a) to (0.2,a) to (0.3,a) to (0.4,a) to (0.5,a)
start(-0.01,a+0.02) line to (0.01,a+0.02)
TIME 0 to 500 by 0.00000001     { establish time range and initial timestep }
MONITORS
for cycle=10     contour (p) report('IQx')  report('IQy')  report('IE2')
elevation (p) from  (0,a) to (0,b)  as "X-sectional"
elevation (LOG10(p)) from(0,a) to (0,b)  as "X-sectional"
elevation (LOG10(ps)) from (0,a) to (0,b)  as "X-sectional"
contour(Qx) report 'IQx' contour(-LOG10(-Qy)) 
contour (ps) contour((p-ps)/ps)
HISTORIES
! history (p) at (0,a) report(Kratio) report(a) report(b) report(pas) report(pa0)
! report(pTau) 
 history (p) at (0,a) report(a) report(b) report(k1) report(k2) report(pa0) report(pas)
 report(pTau) {report(-Cexp/k2) report(Cint) report(beta)}
PLOTS     { plot prints a tag-delimited data list to the file "filename.TXT":}
 grid(x,y)
 history(pLint) at (0,a) export format "#t#b#" file="c:\windows\desktop\Re__0.1538_pi-d_80.txt"
END  26300