Fipy 中的 Dirac delta 源项
Dirac delta Source Term in Fipy
我想知道如何在 Fipy 中将 Dirac delta 函数表示为源项。我想解下面的方程
我试过下面的代码
from fipy import *
nx = 50
ny = 1
dx = dy = 0.025 # grid spacing
L = dx * nx
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
phi = CellVariable(name="solution variable", mesh=mesh, value=0.)
Gamma=1
delta=1 # I want knowing how to make this right.
eqG = TransientTerm() == DiffusionTerm(coeff=Gamma)+delta
valueTopLeft = 0
valueBottomRight = 1
X, Y = mesh.faceCenters
facesTopLeft = ((mesh.facesLeft & (Y > L / 2)) | (mesh.facesTop & (X < L / 2)))
facesBottomRight = ((mesh.facesRight & (Y < L / 2)) |
(mesh.facesBottom & (X > L / 2)))
phi.constrain(valueTopLeft, facesTopLeft)
phi.constrain(valueBottomRight, facesBottomRight)
timeStepDuration = 10 * 0.9 * dx ** 2 / (2 * 0.8)
steps = 100
results=[]
for step in range(steps):
eqG.solve(var=phi, dt=timeStepDuration)
results.append(phi.value)
代码有效,但我想要精确的 Dirac delta 函数。我查找了 numerix 模块,但找不到这样的功能。 Sx1 和 Sy1 是常数。我正在使用 python 2.7
像使用扩散接口方法那样平滑 Dirac delta 函数可能是个好主意(参见等式 11、12 和 13 here)。所以,这是一个选择
def delta_func(x, epsilon):
return ((x < epsilon) & (x > -epsilon)) * \
(1 + numerix.cos(numerix.pi * x / epsilon)) / 2 / epsilon
2 * epsilon
是 Dirac delta 函数的宽度,选择为几个网格间距宽。您也可以只使用 1 / dx
并选择距离 Dirac delta 函数位置最近的网格点。但是,我认为这变得更加依赖于网格。这是一维的工作代码。
from fipy import *
nx = 50
dx = dy = 0.025 # grid spacing
L = dx * nx
mesh = Grid1D(dx=dx, nx=nx)
phi = CellVariable(name="solution variable", mesh=mesh, value=0.)
Gamma=1
def delta_func(x, epsilon):
return ((x < epsilon) & (x > -epsilon)) * \
(1 + numerix.cos(numerix.pi * x / epsilon)) / 2 / epsilon
x0 = L / 2.
eqG = TransientTerm() == DiffusionTerm(coeff=Gamma)+ delta_func(mesh.x - x0, 2 * dx)
valueTopLeft = 0
valueBottomRight = 1
timeStepDuration = 10 * 0.9 * dx ** 2 / (2 * 0.8)
steps = 100
viewer = Viewer(phi)
for step in range(steps):
res = eqG.solve(var=phi, dt=timeStepDuration)
print(step)
viewer.plot()
input('stopped')
这里,epsilon = 2 * dx
,任意选择,delta函数以L / 2
为中心。 2D只需要乘以函数。
我想知道如何在 Fipy 中将 Dirac delta 函数表示为源项。我想解下面的方程
我试过下面的代码
from fipy import *
nx = 50
ny = 1
dx = dy = 0.025 # grid spacing
L = dx * nx
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
phi = CellVariable(name="solution variable", mesh=mesh, value=0.)
Gamma=1
delta=1 # I want knowing how to make this right.
eqG = TransientTerm() == DiffusionTerm(coeff=Gamma)+delta
valueTopLeft = 0
valueBottomRight = 1
X, Y = mesh.faceCenters
facesTopLeft = ((mesh.facesLeft & (Y > L / 2)) | (mesh.facesTop & (X < L / 2)))
facesBottomRight = ((mesh.facesRight & (Y < L / 2)) |
(mesh.facesBottom & (X > L / 2)))
phi.constrain(valueTopLeft, facesTopLeft)
phi.constrain(valueBottomRight, facesBottomRight)
timeStepDuration = 10 * 0.9 * dx ** 2 / (2 * 0.8)
steps = 100
results=[]
for step in range(steps):
eqG.solve(var=phi, dt=timeStepDuration)
results.append(phi.value)
代码有效,但我想要精确的 Dirac delta 函数。我查找了 numerix 模块,但找不到这样的功能。 Sx1 和 Sy1 是常数。我正在使用 python 2.7
像使用扩散接口方法那样平滑 Dirac delta 函数可能是个好主意(参见等式 11、12 和 13 here)。所以,这是一个选择
def delta_func(x, epsilon):
return ((x < epsilon) & (x > -epsilon)) * \
(1 + numerix.cos(numerix.pi * x / epsilon)) / 2 / epsilon
2 * epsilon
是 Dirac delta 函数的宽度,选择为几个网格间距宽。您也可以只使用 1 / dx
并选择距离 Dirac delta 函数位置最近的网格点。但是,我认为这变得更加依赖于网格。这是一维的工作代码。
from fipy import *
nx = 50
dx = dy = 0.025 # grid spacing
L = dx * nx
mesh = Grid1D(dx=dx, nx=nx)
phi = CellVariable(name="solution variable", mesh=mesh, value=0.)
Gamma=1
def delta_func(x, epsilon):
return ((x < epsilon) & (x > -epsilon)) * \
(1 + numerix.cos(numerix.pi * x / epsilon)) / 2 / epsilon
x0 = L / 2.
eqG = TransientTerm() == DiffusionTerm(coeff=Gamma)+ delta_func(mesh.x - x0, 2 * dx)
valueTopLeft = 0
valueBottomRight = 1
timeStepDuration = 10 * 0.9 * dx ** 2 / (2 * 0.8)
steps = 100
viewer = Viewer(phi)
for step in range(steps):
res = eqG.solve(var=phi, dt=timeStepDuration)
print(step)
viewer.plot()
input('stopped')
这里,epsilon = 2 * dx
,任意选择,delta函数以L / 2
为中心。 2D只需要乘以函数。