1

我试图通过一个示例来了解 FiPy 的工作原理，特别是我想解决以下具有周期性边界的简单对流方程：

$$\partial_t u + \partial_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))
4

1 回答 1

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))
于 2016-03-23T15:18:08.043 回答