FiPy Simple Convection

Question

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

Sorry, I can't seem to fix the latex math... – pdevar Mar 23 '16 at 11:49 — pdevar, Mar 23 '16 at 11:49

score 1 · Accepted Answer · answered Mar 23 '16 at 15:18

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

Thanks for the respons wd15. I have been looking into Clawpack and think that will work better, but I hope to make use of FiPy in the futue. — pdevar, Mar 24 '16 at 16:37

FiPy Simple Convection

1 Answers1