Hi guys,
I’m using FEniCS 2019.1.0 to test the convergence of finite element methods for problems with singular data. Specifically, I’m looking at the Poisson equation with a delta source - \nabla^2u = \delta(x).
Using a mollified delta, I’m able to achieve convergence to the analytical solution u = -\frac{1}{2\pi}\ln(r) for any point except the origin. However, when using the built-in FEniCS PointSource method, I can’t get convergence at all. Does anyone have any idea why this is? Find attached the code I’m using.
import mshr as mshr
import fenics as fe
import matplotlib.pyplot as plt
import numpy as np
plt.close('all')
class Delta(fe.UserExpression):
def __init__(self, eps, x0, **kwargs):
self.eps = eps
self.x0 = x0
fe.UserExpression.__init__(self, **kwargs)
def eval(self, values, x):
eps = self.eps
values[0] = eps/fe.pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
def value_shape(self): return ()
resolution = [5, 10, 20, 40]
error_moll_delta = np.zeros(4)
error_pt_source = np.zeros(4)
h = np.zeros(4)
order_moll_delta = np.zeros(4)
order_pt_source = np.zeros(4)
for i in range(len(resolution)):
domain_geometry = mshr.Circle(fe.Point(0., 0.), 1 + 0.01)
mesh = mshr.generate_mesh(domain_geometry, resolution[i])
h[i] = mesh.hmax()
V = fe.FunctionSpace(mesh, "Lagrange", 1)
def boundary(x, on_boundary):
return on_boundary
# Define boundary condition
u0 = fe.Constant(0.0)
bc = fe.DirichletBC(V, u0, boundary)
#%% Using Continuous Delta Approximation
# Define variational problem
u = fe.TrialFunction(V)
v = fe.TestFunction(V)
delta = Delta(eps=1E-4, x0=np.array([0., 0.]), degree=5)
f = delta
magnitude = fe.Constant(1)
a = fe.dot( fe.grad(u), fe.grad(v) )*fe.dx
L = magnitude*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)
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 using mollified delta with resolution ",
resolution[i], " is ", E)
error_moll_delta[i] = E
#%% Using PointSource Method
# Define variational problem
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
# Also plot analytical solution and error function
E = fe.errornorm(u_e, u, 'L2')
error_pt_source[i] = E
print("Error using point source with resolution ", resolution[i], " is ", E)
for j in range(1, len(resolution)):
h_ratio = h[j]/h[j-1]
error_ratio_pt_source = error_pt_source[j]/error_pt_source[j-1]
error_ratio_moll_delta = error_moll_delta[j]/error_moll_delta[j-1]
order_moll_delta[j] = np.log(error_ratio_moll_delta)/np.log(h_ratio)
order_pt_source[j] = np.log(error_ratio_pt_source)/np.log(h_ratio)