How to define product of dirac delta functions as initial condition

General
1

HI ,
I am solving heat equation u_t - u_xx -u_yy = 0 , with u(x,y, 0) = \delta(x)\delta(y), u = 0 on boundary
I tried this code and when I plot the solution I am getting flat surface but i think there should be peak in the solution.

from dolfin import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
nx = 64
# Define the mesh and function space
mesh = RectangleMesh(Point(-10, -20), Point(20, 30), nx, nx)
V = FunctionSpace(mesh, 'P', 1)

# Define the initial condition u(x,y,0) = Ξ΄(x)Ξ΄(y)
class Delta(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)

    def eval(self, values, x):
        values[0] = 1.0 if abs(x[0]) < 1e-12 and abs(x[1]) < 1e-12 else 0.0

u0 = interpolate(Expression('0.0', degree=1), V)


# Define the boundary condition
def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, Constant(0.0), boundary)

# Define the variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.0)

t = 0.0
T = 1.0
dt = 0.01
a = u*v*dx + dt*inner(grad(u), grad(v))*dx
L = u0*v*dx

# Time-stepping loop
while t < T:
    t += dt
    u = Function(V)
    solve(a == L, u, bc)
    u0.assign(u)
#plot comupted solution
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Define a 3D plot
fig = plt.figure()
ax = fig.gca(projection='3d')

# Extract the coordinates and values of the FEM solution
X = mesh.coordinates()[:, 0]
Y = mesh.coordinates()[:, 1]
Z = u.vector()[:]

# Plot the solution as a surface
ax.plot_trisurf(X, Y, Z, cmap=plt.cm.jet)

# Set the axis labels and title
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u')

# Show the plot
plt.show()
2

That initial condition will not work if you do not have a degree if freedom at (0,0), which you do not have with your mesh

3

Sorry, I didn’t get the idea of what you want to say. Can you please explain with an example or provide me any reference? It will be helpful for me.

4

The degrees of freedom has physical coordinates, which you can get through

print(V.tabulate_dof_coordinates)

these are the coordinates sent into

(which you dont call anywhere in your code btw.)

If none of these coordinates are close to (0,0), then you will only get zero values in your initial function (say if you called)

u0 = interpolate(Delta(), V)

Thus you need to make sure your discretization aligns with your point source if you want to use it as an initial condition.

If you don’t want alignment, you need to project Delta into V (but you face the same issue if none of the quadrature points align with the point source).

Alternatively, you can used a smooth approximation of the point source, see: Implement point source in the variational formulation - #2 by MiroK

1 Like
5

I am trying to code a solution for the groundwater pollution model using the FEniCs as given in the paper titled Alternating approaches for groundwater pollution. https://reader.elsevier.com/reader/sd/pii/S2405896315001500?token=971FB3A154F7774257B3FE02E4BC701DA61B47C17ACC23B6F98B9D989ADD185DD73FF6BB2B4B39792471273BB9246755&originRegion=eu-west-1&originCreation=20230418173647

I got the idea that you have explained earlier. Now I have defined the Dirac distribution as the limit of Gauss function, now I am getting the peak in the solution, but I am not sure if I am doing right or not because in defining the initial condition I have considered only one variable x but there should be one more variable y.

The code is given as:

β€œβ€"

from fenics import*

import numpy as np

Define the mesh and function space

nx = 64
mesh = RectangleMesh(Point(-10, -20), Point(20, 30), nx, nx)
V = FunctionSpace(mesh, β€˜P’, 1)

Define the Gaussian function G(x) FOR APPROXIMATION OF dIRAC DELTA FUNCTION

def G(x, v1=0.02, D=0.02, sigma=0.025):
return 1/(4piD*sigma)*np.exp(-((x - v1)**2) / (4 D sigma))

Define the initial condition u(x,y,0) = G(x)G(y)

class InitialCondition(UserExpression):
def init(self, **kwargs):
super().init(**kwargs)

def eval(self, values, x):
values[0] = G(x[0]) * G(x[1])

u0 = interpolate(InitialCondition(degree=1), V)

Define the boundary condition

def boundary(x, on_boundary):
return on_boundary

bc = DirichletBC(V, Constant(0.0), boundary)

Define the variational problem

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.0)
v1 = Constant(0.02)
v2 = 0
e1 = as_vector([v1,v2])
D = Constant(0.02)
t = 0.0
T = 250
dt = 0.025
a = uvdx + dtDinner(grad(u), grad(v))dx + dot(e1, grad(u))vdx
L = u0v*dx

Time-stepping loop

while t < T:
t += dt
u = Function(V)
solve(a == L, u, bc)
u0.assign(u)

import pylab as plt
Q = plot(u)
plt.xlabel(β€˜x’)
plt.ylabel(β€˜y’)
plt.title(β€˜t={}’.format(t))
plt.set_cmap(β€˜jet’)
plt.colorbar(Q)
plt.show()

Define a 3D plot

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

Extract the coordinates and values of the FEM solution

X = mesh.coordinates()[:, 0]
Y = mesh.coordinates()[:, 1]
Z = u.vector()[:]

Plot the solution as a surface

fig = plt.figure()
ax = fig.gca(projection=β€˜3d’)
ax.plot_trisurf(X, Y, Z, cmap=plt.cm.jet)

Set the axis labels and title

ax.set_xlabel(β€˜x’)
ax.set_ylabel(β€˜y’)
ax.set_zlabel(β€˜u’)

Show the plot

plt.show()

β€œβ€β€œβ€

6

Gentle reminder, I am stuck on this point.

7

First of all, you need to format your code with 3x` encapsulation, i.e.

```python
from fenics import * 
# add code here
```

Secondly, you can verify your initial condition by plotting u0. I.e. either use plot(u0) or File("u0.pvd") << u0 and visualize it in paraview.

8 
from dolfin import *
import numpy as np

# Define the mesh and function space
nx = 64
mesh = RectangleMesh(Point(-10, -20), Point(20, 30), nx, nx)
V = FunctionSpace(mesh, 'P', 1)

# Define the Gaussian function G(x) FOR APPROXIMATION OF dIRAC DELTA FUNCTION 
def G(x, v1=0.02, D=0.02, sigma=0.025):
    return 1/(4*pi*D*sigma)*np.exp(-((x - v1)**2) / (4 *D* sigma))

# Define the initial condition u(x,y,0) = G(x)G(y)
class InitialCondition(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)

    def eval(self, values, x):
        values[0] = G(x[0]) * G(x[1])

u0 = interpolate(InitialCondition(degree=1), V)

# Define the boundary condition
def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, Constant(0.0), boundary)

# Define the variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.0)
v1 = Constant(0.02)
v2 = 0 
e1 = as_vector([v1,v2])
D = Constant(0.02)
t = 0.0
T = 250
dt = 0.025
a = u*v*dx + dt*D*inner(grad(u), grad(v))*dx + dot(e1, grad(u))*v*dx
L = u0*v*dx

# Time-stepping loop
while t < T:
    t += dt
    u = Function(V)
    solve(a == L, u, bc)
    u0.assign(u)
     
import pylab as plt 
Q = plot(u)
plt.xlabel('x')
plt.ylabel('y')
plt.title('t={}'.format(t))
plt.set_cmap('jet')
plt.colorbar(Q)
plt.show()
    

# Define a 3D plot
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Extract the coordinates and values of the FEM solution
X = mesh.coordinates()[:, 0]
Y = mesh.coordinates()[:, 1]
Z = u.vector()[:]

# Plot the solution as a surface
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(X, Y, Z, cmap=plt.cm.jet)

# Set the axis labels and title
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u')

# Show the plot
plt.show()
9

Is there a question in this post?

It seems like your scaling of the delta function is wrong

1 Like
10

Thank you for your help. I got the results.