1

I would like to know how to represent the Dirac delta function as source term in Fipy. I want to solve the following equation enter image description here

I have tried the following code

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)

The code is working but I want the exact Dirac delta function. I looked up the numerix module but couldnt find such function. Sx1 and Sy1 are constants. Am using python 2.7

Mafeni Alpha
  • 308
  • 2
  • 13
  • Would you be cool with a Python answer that isn't built in? `def dirac_delta(x): return 1 if x==0 else 0` gets you 99% of the way there. We could spice it up with `np.isclose()` if you need it. – user1717828 Sep 21 '19 at 14:50
  • Am not sure if defining like this will produce the same desired result as in fipy module they strongly recommended putting sources as they are. But I may also try your suggestion in case no ones knows the good way – Mafeni Alpha Sep 21 '19 at 15:13

1 Answers1

2

It's probably a good idea to smooth the Dirac delta function as is done with diffusion interface methods (see equations 11, 12 and 13 here). So, this is one choice

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

2 * epsilon is the width of the Dirac delta function and is chosen to be a few grid spacings wide. You could also just use 1 / dx and choose the closest grid point to the Dirac delta function's location. However, I think that becomes more grid dependent. Here is a working code in 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')

Here, epsilon = 2 * dx, an arbitrary choice, and the delta function is centered around L / 2. 2D just required multiplying the functions.

wd15
  • 1,068
  • 5
  • 8