FiPy 简单对流
FiPy Simple Convection
我试图通过示例来理解 FiPy 的工作原理,特别是我想求解以下具有周期边界的简单对流方程:
$$\partial_t 你 + \partial_x 你 = 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))
如 FiPy mailing list 所述,FiPy 不擅长处理仅对流的 PDE(没有扩散,纯双曲线),因为它缺少高阶对流方案。对于这个 class 的问题,最好使用 CLAWPACK。
FiPy 确实有一个可能有助于解决此问题的二阶方案,即 VanLeerConvectionTerm,see an example。
如果在上述问题中使用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))
我试图通过示例来理解 FiPy 的工作原理,特别是我想求解以下具有周期边界的简单对流方程:
$$\partial_t 你 + \partial_x 你 = 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))
如 FiPy mailing list 所述,FiPy 不擅长处理仅对流的 PDE(没有扩散,纯双曲线),因为它缺少高阶对流方案。对于这个 class 的问题,最好使用 CLAWPACK。
FiPy 确实有一个可能有助于解决此问题的二阶方案,即 VanLeerConvectionTerm,see an example。
如果在上述问题中使用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))