FiPy Simple Convection

I am trying to understand how FiPy works by working an example, in particular I would like to solve the following simple convection equation with periodic boundary:

$$\partial_t u + \partial_x u = 0$$

If initial data is given by $u(x, 0) = F(x)$, then the analytical solution is $u(x, t) = F(x - t)$. I do get a solution, but it is not correct.

What am I missing? Is there a better resource for understanding FiPy than the documentation? It is very sparse...

Here is my attempt

from fipy import *
import numpy as np

# Generate mesh
nx = 20
dx = 2*np.pi/nx
mesh = PeriodicGrid1D(nx=nx, dx=dx)

# Generate solution object with initial discontinuity
phi = CellVariable(name="solution variable", mesh=mesh)
phiAnalytical = CellVariable(name="analytical value", mesh=mesh)
phi.setValue(1.)
phi.setValue(0., where=x > 1.)

# Define the pde
D = [[-1.]]
eq = TransientTerm() == ConvectionTerm(coeff=D)

# Set discretization so analytical solution is exactly one cell translation
dt = 0.01*dx
steps = 2*int(dx/dt)

# Set the analytical value at the end of simulation
phiAnalytical.setValue(np.roll(phi.value, 1))

for step in range(steps):
    eq.solve(var=phi, dt=dt)

print(phi.allclose(phiAnalytical, atol=1e-1))

As addressed on the FiPy mailing list, FiPy is not great at handling convection only PDEs (absent diffusion, pure hyperbolic) as it's missing higher order convection schemes. It is better to use CLAWPACK for this class of problem.

FiPy does have one second order scheme that might help with this problem, the VanLeerConvectionTerm, see an example.

If the VanLeerConvectionTerm is used in the above problem, it does do a better job of preserving the shock.

import numpy as np
import fipy

# Generate mesh
nx = 20
dx = 2*np.pi/nx
mesh = fipy.PeriodicGrid1D(nx=nx, dx=dx)

# Generate solution object with initial discontinuity
phi = fipy.CellVariable(name="solution variable", mesh=mesh)
phiAnalytical = fipy.CellVariable(name="analytical value", mesh=mesh)
phi.setValue(1.)
phi.setValue(0., where=mesh.x > 1.)

# Define the pde
D = [[-1.]]
eq = fipy.TransientTerm() == fipy.VanLeerConvectionTerm(coeff=D)

# Set discretization so analytical solution is exactly one cell translation
dt = 0.01*dx
steps = 2*int(dx/dt)

# Set the analytical value at the end of simulation
phiAnalytical.setValue(np.roll(phi.value, 1))

viewer = fipy.Viewer(phi)
for step in range(steps):
    eq.solve(var=phi, dt=dt)
    viewer.plot()
    raw_input('stopped')
print(phi.allclose(phiAnalytical, atol=1e-1))