如何在 FiPy 中编写带有自变量的方程项?
How to write equation terms with independent variable in FiPy?
我正在寻找使用 FiPy 求解扩散方程并阅读了他们的一些文档,但似乎找不到任何与编写扩散项相关的内容,该扩散项包含作为自变量函数的附加项 (即 space)。我发现最接近的是 FAQ,他们建议将附加术语重写为 ConvectionTerm
。但是,我认为这仅适用于附加项是解变量而不是自变量的函数的情况。例如,我正在尝试用以下扩散项求解一维扩散方程(其中导数为 w.r.t。自变量 x,y 是解变量):
D * sin(x) * Div_x {sin(x) * Grad_x {y}}
我觉得这是一个非常简单的表达式,但我找不到用 FiPy 表示法表达它的方法。任何帮助将不胜感激!
确切代码:
from fipy import Variable,FaceVariable,CellVariable,Grid1D,ImplicitSourceTerm,TransientTerm,DiffusionTerm,Viewer,ConvectionTerm
from fipy.tools import numerix
D = 1
c0 = 1
ka = 1
r0 = 1
nx = 100
dx = 2*math.pi/100
mesh = Grid1D(nx=nx, dx=dx)
conc = CellVariable(name="concentration", mesh=mesh, value=0.) # This is the "phi" in the docs
valueLeft = c0
valueRight = 0
conc.constrain(valueRight, mesh.facesRight)
conc.constrain(valueLeft, mesh.facesLeft)
timeStepDuration = 0.9 * dx**2 / (2 * D)
steps = 100
show_per_steps = 50
A = 1 / (r0**2 * numerix.sin(mesh.x)[0])
dA = -(numerix.cos(mesh.x)[0])/(r0**2 * numerix.sin(mesh.x)[0]**2)
dsindA = (numerix.cos(mesh.x)[0])**3/(numerix.sin(mesh.x)[0])**2
eqX = TransientTerm() + ImplicitSourceTerm(ka) == DiffusionTerm(D*A*numerix.sin(mesh.x)[0]) - ConvectionTerm(D*dA*numerix.cos(mesh.x)[0])+ D*conc*dsindA
from builtins import range
for step in range(steps):
eqX.solve(var=conc, dt=timeStepDuration)
if __name__ == '__main__' and step % show_per_steps == 0:
viewer = Viewer(vars=(conc), datamin=0., datamax=c0)
viewer.plot()
FiPy 允许项的系数是 space 的函数。例如,以下在 FiPy 中工作,
from fipy import Grid1D, CellVariable, Viewer
from fipy import TransientTerm, numerix, DiffusionTerm
from fipy import LinearLUSolver
m = Grid1D(nx=100, Lx=numerix.pi / 4.)
v = CellVariable(mesh=m)
v[:] = m.x**2
eqn = TransientTerm() == DiffusionTerm(numerix.sin(m.x))
vi = Viewer(v, colorbar=None)
vi.plot()
solver = LinearLUSolver()
for i in range(10):
eqn.solve(v, dt=0.1, solver=solver)
vi.plot()
print('step', i)
input('stopped')
在上面,扩散系数是space的函数。 m.x
是保存单元格中心位置的 CellVariable
。 numerix
用于启用对 FiPy 变量的操作,其方式与 Numpy 对 Numpy 数组的操作方式相同。
现在,在上面的问题中,导数之外有一个 sin(x)
,这在 FiPy 中是不允许的。一切都需要适合导数内部才能与 FiPy 一起工作。因此,我们需要重写该项,以便所有系数都在导数内。对于任何一般情况,我们可以写成
这允许我们使用扩散、对流和源来表示 FiPy 中的术语。如果 f=sin(x)
和 g=sin(x)
那么 FiPy 代码将是
s = numerix.sin(m.cellCenters)
c = numerix.cos(m.cellCenters)
eqn = ... + DiffusionTerm(D * s[0] * s[0]) - ConvectionTerm(D * s * c) + D * y * (c[0] * c[0] - s[0] * s[0])
...
已包括在内,因为我不知道完整的方程式。
我正在寻找使用 FiPy 求解扩散方程并阅读了他们的一些文档,但似乎找不到任何与编写扩散项相关的内容,该扩散项包含作为自变量函数的附加项 (即 space)。我发现最接近的是 FAQ,他们建议将附加术语重写为 ConvectionTerm
。但是,我认为这仅适用于附加项是解变量而不是自变量的函数的情况。例如,我正在尝试用以下扩散项求解一维扩散方程(其中导数为 w.r.t。自变量 x,y 是解变量):
D * sin(x) * Div_x {sin(x) * Grad_x {y}}
我觉得这是一个非常简单的表达式,但我找不到用 FiPy 表示法表达它的方法。任何帮助将不胜感激!
确切代码:
from fipy import Variable,FaceVariable,CellVariable,Grid1D,ImplicitSourceTerm,TransientTerm,DiffusionTerm,Viewer,ConvectionTerm
from fipy.tools import numerix
D = 1
c0 = 1
ka = 1
r0 = 1
nx = 100
dx = 2*math.pi/100
mesh = Grid1D(nx=nx, dx=dx)
conc = CellVariable(name="concentration", mesh=mesh, value=0.) # This is the "phi" in the docs
valueLeft = c0
valueRight = 0
conc.constrain(valueRight, mesh.facesRight)
conc.constrain(valueLeft, mesh.facesLeft)
timeStepDuration = 0.9 * dx**2 / (2 * D)
steps = 100
show_per_steps = 50
A = 1 / (r0**2 * numerix.sin(mesh.x)[0])
dA = -(numerix.cos(mesh.x)[0])/(r0**2 * numerix.sin(mesh.x)[0]**2)
dsindA = (numerix.cos(mesh.x)[0])**3/(numerix.sin(mesh.x)[0])**2
eqX = TransientTerm() + ImplicitSourceTerm(ka) == DiffusionTerm(D*A*numerix.sin(mesh.x)[0]) - ConvectionTerm(D*dA*numerix.cos(mesh.x)[0])+ D*conc*dsindA
from builtins import range
for step in range(steps):
eqX.solve(var=conc, dt=timeStepDuration)
if __name__ == '__main__' and step % show_per_steps == 0:
viewer = Viewer(vars=(conc), datamin=0., datamax=c0)
viewer.plot()
FiPy 允许项的系数是 space 的函数。例如,以下在 FiPy 中工作,
from fipy import Grid1D, CellVariable, Viewer
from fipy import TransientTerm, numerix, DiffusionTerm
from fipy import LinearLUSolver
m = Grid1D(nx=100, Lx=numerix.pi / 4.)
v = CellVariable(mesh=m)
v[:] = m.x**2
eqn = TransientTerm() == DiffusionTerm(numerix.sin(m.x))
vi = Viewer(v, colorbar=None)
vi.plot()
solver = LinearLUSolver()
for i in range(10):
eqn.solve(v, dt=0.1, solver=solver)
vi.plot()
print('step', i)
input('stopped')
在上面,扩散系数是space的函数。 m.x
是保存单元格中心位置的 CellVariable
。 numerix
用于启用对 FiPy 变量的操作,其方式与 Numpy 对 Numpy 数组的操作方式相同。
现在,在上面的问题中,导数之外有一个 sin(x)
,这在 FiPy 中是不允许的。一切都需要适合导数内部才能与 FiPy 一起工作。因此,我们需要重写该项,以便所有系数都在导数内。对于任何一般情况,我们可以写成
这允许我们使用扩散、对流和源来表示 FiPy 中的术语。如果 f=sin(x)
和 g=sin(x)
那么 FiPy 代码将是
s = numerix.sin(m.cellCenters)
c = numerix.cos(m.cellCenters)
eqn = ... + DiffusionTerm(D * s[0] * s[0]) - ConvectionTerm(D * s * c) + D * y * (c[0] * c[0] - s[0] * s[0])
...
已包括在内,因为我不知道完整的方程式。