Non-Convergence When Using 'PointSource' Method

General
1

Hi everyone,

I’m using FEniCS 2019.1.0 and I’m experimenting with FEniCS’ ability to represent concentrated point sources (Dirac Delta ‘functions’). The results I’m seeing are rather surprising -

In the program below, I’m using two methods to approximate the delta on the right hand side of

\tag{1} -\nabla^2 u = \delta(\mathbf{x}),

where the analytical solution is given by

u(x,y) = -\frac{1}{2\pi}\ln(||\mathbf{x}||).

I first use a mollified delta defined by the continuous function

\tag{2} \delta_h(x,y) = \begin{cases} \frac{1}{h^2}\biggl( 1 - \frac{|x|}{h} \biggr)\biggl( 1 - \frac{|y|}{h} \biggr), \quad \max{(x,y}) < 1 \\ 0, \quad \mathrm{else} \end{cases}

Conversely, I also compute the numerical solution to (1) using the PointSource method with magnitude 1.

Two things are rather surprising about this -
Firstly, the approximation with the mollified delta is better than the approximation using the PointSource - regardless of the mesh resolution. Another surprising bit is that the error does not decrease to 0 with the mesh size. In both cases, the error decreases asymptotically to about 0.19.

Does anyone have any ideas about what’s going on here? Why the mollified delta performs better than the built-in FEniCS function and why I don’t get convergence to the analytical solution?

import fenics as fe
import mshr as msh
import numpy as np
import matplotlib.pyplot as plt

N = 10
# (a) using FEniCS' built-in 'PointSource' method

domain = msh.Circle(fe.Point(0,0), 1)
mesh = msh.generate_mesh(domain, N)
V = fe.FunctionSpace(mesh, "P", 1)

u_D = fe.Expression("0", degree = 2)

def boundary(x, on_boundary):
    return on_boundary
bc = fe.DirichletBC(V, u_D, boundary)

u = fe.TrialFunction(V)
v = fe.TestFunction(V)
f = fe.Constant(0)

a = fe.dot( fe.grad(u), fe.grad(v) )*fe.dx
L = f*v*fe.dx

u = fe.Function(V)

assembler = fe.SystemAssembler(a, L, bc)
solver = fe.LUSolver("mumps")
A = fe.Matrix()
solver.set_operator(A)
assembler.assemble(A)

b = fe.assemble(L)
delta = fe.PointSource(V, fe.Point(0,0), 1)
delta.apply(b)
bc.apply(b)
fe.solve(A, u.vector(), b)

u_e = fe.Expression('-1/(2*3.14159)*std::log( pow( pow(x[0],2) + pow(x[1],2), 1/2) )', degree = 6)
V = u.function_space()
mesh = V.mesh()
degree = V.ufl_element().degree()
W = fe.FunctionSpace(mesh, "P", degree + 3)
u_e_w = fe.interpolate(u_e, W)
u_w = fe.interpolate(u, W)
e_w = fe.Function(W)
e_w.vector()[:] = u_e_w.vector().get_local() - u_w.vector().get_local()
error = e_w**2*fe.dx

E = fe.errornorm(u_e, u, 'L2')

print("error in exact delta = ", E, "\n" )


#%% Mollified Delta

domain = msh.Circle(fe.Point(0,0), 1)
mesh = msh.generate_mesh(domain, N)
V = fe.FunctionSpace(mesh, "P", 1)

u_D = fe.Expression("0", degree = 2)

def boundary(x, on_boundary):
    return on_boundary
bc = fe.DirichletBC(V, u_D, boundary)

u = fe.TrialFunction(V)
v = fe.TestFunction(V)

h = mesh.hmax()
delta_approx = fe.Expression('(abs(x[0] -r_x)< h ) && abs(x[1]-r_y) < h'\
                 +' ? 1/pow(h,2)*(1 - abs(x[0]-r_x)/h)*(1-abs(x[1]-r_y)/h) : 0',degree=2,r_x =0,r_y=0,h=h )   


a = fe.dot( fe.grad(u), fe.grad(v) )*fe.dx
L = delta_approx*v*fe.dx

u = fe.Function(V)

assembler = fe.SystemAssembler(a, L, bc)
solver = fe.LUSolver("mumps")
A = fe.Matrix()
solver.set_operator(A)
assembler.assemble(A)

b = fe.assemble(L)
bc.apply(b)
fe.solve(A, u.vector(), b)

# Compute error


u_e = fe.Expression('-1/(2*3.14159)*std::log( pow( pow(x[0], 2) + pow(x[1],2), 1/2) )', degree = 6)
V = u.function_space()
mesh = V.mesh()
degree = V.ufl_element().degree()
W = fe.FunctionSpace(mesh, "P", degree + 3)
u_e_w = fe.interpolate(u_e, W)
u_w = fe.interpolate(u, W)
e_w = fe.Function(W)
e_w.vector()[:] = u_e_w.vector().get_local() - u_w.vector().get_local()
error = e_w**2*fe.dx

E = fe.errornorm(u_e, u, 'L2')
print("error in moll. delta = ", E )
2

You have to consider the regularity requirements of FE formulations. See e.g., Finite element convergence for singular data | Numerische Mathematik.

1 Like