Fipy中的Dirac delta源项
Fipy中的Dirac delta源项
提问于 2019-09-21 22:28:08
我想知道如何在Fipy中将狄拉克增量函数表示为源项。我想解下面的方程式
我已经尝试了以下代码
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)代码是有效的，但我想要精确的狄拉克增量函数。我查找了numerix模块，但没有找到这样的函数。Sx1和Sy1是常量。我使用的是python 2.7
回答 1
Stack Overflow用户
回答已采纳
发布于 2019-09-25 00:21:15
这可能是一个好主意，使狄拉克增量函数平滑，就像用扩散界面方法所做的那样(参见公式11，12和13 here)。所以，这是一种选择
def delta_func(x, epsilon):
return ((x < epsilon) & (x > -epsilon)) * \
(1 + numerix.cos(numerix.pi * x / epsilon)) / 2 / epsilon2 * epsilon是狄拉克增量函数的宽度，并被选择为几个网格间距。您也可以只使用1 / dx并选择离狄拉克增量函数位置最近的网格点。然而，我认为这变得更加依赖于网格。这是一个1D的工作代码。
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和增量函数以L / 2为中心。2D只需要将函数相乘即可。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/58041222
复制相关文章
相似问题
腾讯云开发者
