Skip to content

Instantly share code, notes, and snippets.

@jkfurtney
Created March 20, 2013 21:29
Show Gist options
  • Save jkfurtney/5208651 to your computer and use it in GitHub Desktop.
Save jkfurtney/5208651 to your computer and use it in GitHub Desktop.
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment