limiters.py

import fipy.tools.numerix as np
from fipy.terms import VanLeerConvectionTerm
import fipy as fp
import pylab as plt

# monkey patch VanLeerConvectionTerm so we can override
# __getGradient in the derived class
VanLeerConvectionTerm._getGradient = VanLeerConvectionTerm._VanLeerConvectionTerm__getGradient
VanLeerConvectionTerm._VanLeerConvectionTerm__getGradient = \
    lambda self, a, b: self._getGradient(a, b)

class MinModConvectionTerm(VanLeerConvectionTerm):
    def _getGradient(self, normalGradient, gradUpwind):
        gradUpUpwind = -gradUpwind + 2 * normalGradient
        min2 = np.minimum(abs(gradUpwind), abs(gradUpUpwind))
        grad = np.where(gradUpwind * gradUpUpwind < 0.,
                             0.,
                             np.where(gradUpUpwind > 0.,
                                           min2,
                                           -min2))
        return grad

class MonotonizedCentralDifferenceConvectionTerm(VanLeerConvectionTerm):
    def _getGradient(self, normalGradient, gradUpwind):
        gradUpUpwind = -gradUpwind + 2 * normalGradient
        avg = 0.5 * (abs(gradUpwind) + abs(gradUpUpwind))
        min3 = np.minimum(np.minimum(abs(2 * gradUpwind),
                                               abs(2 * gradUpUpwind)), avg)
        grad = np.where(gradUpwind * gradUpUpwind < 0.,
                             0.,
                             np.where(gradUpUpwind > 0.,
                                           min3,
                                           -min3))
        return grad

class SuperBeeConvectionTerm(VanLeerConvectionTerm):
    def _getGradient(self, normalGradient, gradUpwind):
        gradUpUpwind = -gradUpwind + 2 * normalGradient
        min1 = np.minimum(abs(2 * gradUpwind), abs(gradUpUpwind))
        min2 = np.minimum(abs(gradUpwind), abs(2 * gradUpUpwind))
        max2 = np.maximum(min1, min2)
        grad = np.where(gradUpwind * gradUpUpwind < 0.,
                             0.,
                             np.where(gradUpUpwind > 0.,
                                           max2,
                                           -max2))
        return grad


terms = [fp.ExplicitUpwindConvectionTerm, fp.VanLeerConvectionTerm,
         MinModConvectionTerm, MonotonizedCentralDifferenceConvectionTerm,
         SuperBeeConvectionTerm]

nx = 100
l = 1.0
dx = l / nx
mesh = fp.Grid1D(nx=nx, dx=dx)
u = 1.0
uv = fp.FaceVariable(mesh=mesh, value=u,  rank=1)

cell_vars = [fp.CellVariable(mesh=mesh, value=0.0) for var in terms]
for cell_var, term in zip(cell_vars, terms):
    cell_var.constrain(1.0, mesh.facesLeft())
    cell_var.equation = (fp.TransientTerm() + term(coeff=uv) == 0)

cfl =  0.75
dt = cfl * dx / u
t, tend = 0, 0.75 * l

while t < tend:
    for cell_var in cell_vars:
        cell_var.equation.solve(var=cell_var,  dt=dt)
    t += dt

plt.title("Convection Term Comparison\n$q_t+q_x=0$ with $q(x,0)=0$ and $q(0,t)=1$")
plt.axvline(0.75 * nx)
for cell_var in cell_vars:
    plt.plot(cell_var.value, "o")
plt.legend(["exact"]+[str(term) for term in terms], loc=3)
plt.show()