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()