FiPy Простая конвекция

Question

Я пытаюсь понять, как работает FiPy, работая на примере, в частности, я хотел бы решить следующее простое уравнение конвекции с периодической границей:

$$ \partal_t u + \ частичный_x u = 0$$

Если исходные данные задаются как $u(x, 0) = F(x)$, то аналитическим решением является $u(x, t) = F(x - t)$. Я получаю решение, но оно не правильное.

Что мне не хватает? Есть ли лучший ресурс для понимания FiPy, чем документация? Это очень редко...

Вот моя попытка

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

Answer 1

Как указано в списке рассылки FiPy, FiPy не очень хорошо справляется с конвекционными PDE (отсутствует диффузия, чисто гиперболический), поскольку в ней отсутствуют схемы конвекции более высокого порядка. Для этого класса проблем лучше использовать CLAWPACK.

FiPy имеет одну схему второго порядка, которая может помочь с этой проблемой, VanLeerConvectionTerm, см. Пример.

Если VanLeerConvectionTerm используется в вышеупомянутой проблеме, это делает лучшую работу по сохранению шока.

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