Fipy Verification

This page contains notes on verification of Fipy.

What is Verification

Two important activities of model development are verification and validation.

Verification is confirmation that the code implemented and solved the equations correctly, and relates to how well the discrete mathematical representation matches the (exact) mathematical formulation.

Validation is confirming that the correct equations are being solved, and relates to how well the mathematics describe the physics.

Desirable to verify Fipy because we want to make sure it is solving the correct equations. Also, useful for benchmarking/other activities.

Reaction Diffusion Equation: Solution Verification

The Equation

The reaction diffusion equation for a first-order reaction consists of the following equation:

$\nabla ^{2}u=\phi ^{2}u$

which is derived from a non-dimensionalized material balance which contains only a second-order diffusion term (the left-hand side) and a reaction rate term (the right-hand side). $\phi$ is the Thiele modulus and is a ratio of the reaction rate and the diffusion rate. A large Thiele modulus indicates a reaction rate that is very fast relative to the diffusion rate.

The first (symmetry) boundary condition is given by the Neumann boundary condition:

$\nabla u=0{\text{ at }}x=0$

The second (surface) boundary condition relates to the rate of diffusive flux, and is given by the Robin boundary condition:

$\nabla u=v(1-u){\text{ at }}x=1$

where $v$ is the Sherwood number.

The Analytical Solution

Wolfram Alpha was used to determine an analytical solution to this equation. The solution page is available here: http://wolfr.am/1kddYKK

The solution to the equation is given by the function:

$u(x)={\dfrac {v\cosh {ax}}{v\cosh {a}+a\sinh {a}}}$

(Note that the link above gives the solution in terms of exponentials; this can be simplified/rearranged into the above trigonometric form; see http://wolfr.am/1cdxLHK.)

The Numerical Solution

Download the whole script here: Fipy Verification/Rxn Diffusion Verification Script

The numerical solution to the above equation was implemented in Fipy using the following script.

First, import the things we need:

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

Next, set the values of the two numerical parameters, the Sherwood number and the Thiele modulus:

# 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

Set up the numerical solution in Fipy:

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'

Constrain the solution to satisfy the boundary conditions:

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

Now make your equation, and the exact solution:

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

x = mesh.cellCenters[0]
uexact.setValue( (v*numerix.cosh(a*x))/( v*cosh(a) + a*sinh(a) ) )

Iterate to find a solution, and plot it:

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