Дирак дельта источник термин в фипы - PullRequest
Новая
       97

Дирак дельта источник термин в фипы

1 голос
/ 21 сентября 2019

Я хотел бы знать, как представить дельта-функцию Дирака в качестве исходного термина в Fipy.Я хочу решить следующее уравнение enter image description here

Я пробовал следующий код 

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)

Код работает, но мне нужна точная дельта-функция Дирака.Я посмотрел модуль цифр, но не смог найти такую ​​функцию.Sx1 и Sy1 являются константами.Использую python 2.7

1 Ответ

1 голос
/ 24 сентября 2019

Вероятно, неплохо сгладить дельта-функцию Дирака, как это делается с помощью методов диффузионного интерфейса (см. Уравнения 11, 12 и 13 здесь ).Таким образом, это один из вариантов 

def delta_func(x, epsilon):
    return ((x < epsilon) & (x > -epsilon)) * \
        (1 + numerix.cos(numerix.pi * x / epsilon)) / 2 / epsilon

2 * 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 просто нужно умножить функции.

...