Fipy Verification/Rxn Diffusion Verification Script

from numpy import *
from fipy import *
import matplotlib.pyplot as plt

"""
Comparison of Computed/Analytical Solutions
Steady-State Reaction-Diffusion Solution

Charles Reid
February 2014

Solution to reaction diffusion equation,

u_xx = phi^2 u

u'(x=0) = 0
u'(x=1) = v ( 1 - u(x=1) )

Solution via Wolfram Alpha:

u(x) = v * cosh(ax) / ( v * cosh(a) + a * sinh(a) )

http://wolfr.am/1kddYKK

"""

# the larger v is, the higher the surface concentration is (max of 1)
# (and the longer convergence takes; residual = v)
v = 1.0

# the smaller a is, the further into the particle diffusion goes,
# and the longer convergence takes
# (i.e., very large a means you're just using the skin of the particle)
a = 10.0

x = linspace(0,1,100)

L = 1.0
nx = 100
dx = L / nx

mesh = Grid1D(nx=nx, dx=dx)

u = CellVariable(mesh=mesh)
uexact = CellVariable(mesh=mesh)

u.name = 'u'
uexact.name = 'u exact'

u.faceGrad.constrain([0], mesh.facesLeft )
u.faceGrad.constrain([v*(1.0 - u.faceValue)], mesh.facesRight )

eq = DiffusionTerm(1.0) == ImplicitSourceTerm(a*a)

x = mesh.cellCenters[0]

uexact.setValue( (v*numerix.cosh(a*x))/( v*cosh(a) + a*sinh(a) ) )

restol = 1e-5
res = 1e+10
while res > restol:
    res = eq.sweep(var=u)
    print "Residual = ",res

if __name__ == '__main__':
    viewer = Viewer(vars=(u,uexact))
    viewer.plot()