### How to solve Poisson equation with singular source at the boundary?

232
views
1
5 months ago by
Hi,

I am new to the finite element method and FEniCS and I try to solve a Poisson equation with a source term which is singular at one boundary.

As a simplified example consider the following equation
$\nabla^2 \phi(x, y) = \frac{2-x^2-y^2}{{(1-x^2-y^2)}^{3/2}}$
on the domain (quarter circle)
$\Omega = \left\{ (x, y)\, | \, x\ge 0,\, y\ge 0,\, x^2+y^2\le 1 \right\} \subset \mathbb{R}^2$
with the boundary conditions
$\phi(x, y) = \begin{cases} -\sqrt{1-x^2} &\text{for } y = 0\\ 0 &\text{for } x^2+y^2 = 1 \end{cases}$
$\frac{\partial \phi}{\partial n} = 0 \quad \text{for} \quad x = 0 \ .$
The solution of this problem is of course just (1/8 of) a spherical shell. The problem with the spherical shell however is, that in the given Cartesian coordinates its derivative is singular at the boundary with $x^2+y^2=1$.

My question is, if there is some nice way to solve problems of this kind in FEniCS without having to cripple the problem or to use an extremely dense mesh at the boundary.

Of course one could trivially solve my example problem by using spherical coordinates, but this is not the solution I am looking for.

My current implementation for the given example is as follows:
from fenics import FunctionSpace, Expression, Constant, DirichletBC, \
TrialFunction, TestFunction, dot, grad, dx, Function, solve, plot
from dolfin import Point, DOLFIN_EPS
from mshr import Circle, Rectangle, generate_mesh

# Parameters for the mesh
r = 1
dens = 50

# Generate the mesh
r1 = Circle(Point(0, 0), r)
r2 = Rectangle(Point(-r, 0), Point(r, -r))
r3 = Rectangle(Point(-r, 0), Point(0, r))
domain = r1-r2-r3
mesh = generate_mesh(domain, dens)

# Function Space
deg = 3
V = FunctionSpace(mesh, 'P', deg)

# Boundary conditions
p1 = Expression('-pow(1-pow(x[0], 2), 0.5)', degree=deg)
p2 = Constant('0')

def boundary_1(x, on_boundary):
return x[1] < DOLFIN_EPS and on_boundary

def boundary_2(x, on_boundary):
return x[0] > DOLFIN_EPS and x[1] > DOLFIN_EPS and on_boundary

bc1 = DirichletBC(V, p1, boundary_1)
bc2 = DirichletBC(V, p2, boundary_2)
bcs = [bc1, bc2]

# Define variational problem
phi = TrialFunction(V)
v = TestFunction(V)
f = Expression('(2-pow(x[0], 2)-pow(x[1], 2)) * ' +
'pow(1-pow(x[0], 2)-pow(x[1], 2), -1.5)', degree=deg)
L = -f*v*dx

# Compute solution
phi = Function(V)
solve(a == L, phi, bcs)

# Plot the solution
plot(phi)​

The basic variables in my code are the radius r and the density of the mesh dens. I already discovered, that if I reduce r I can obtain quite okay results if the mesh is dense enough, for example r=0.99 and dens=150. But setting artificially r<1 seems to be a really ugly solution in my point of view and furthermore my actual source term is computationally very demanding, such that I can not afford to use a very dense mesh. Any help would be greatly appreciated.
Community: FEniCS Project

4
5 months ago by
Your question is related with this: How to compute the integral of function with singularity?. Consider this workaround:

from dolfin import SpatialCoordinate

# Define variational problem
phi = TrialFunction(V)
v = TestFunction(V)
x = SpatialCoordinate(mesh)
f = (2-x[0]**2-x[1]**2)/((1-x[0]**2-x[1]**2)**1.5)
L = -f*v*dx​
Thank you very much! This works like a charm!
written 5 months ago by BennyS
For my test problem your solution is perfect! However, in my real problem the source term f is a complicated numerical Python function. How would I generalize your SpatialCoordinate trick in this case?

For example if I define f as
def F(x, y):
return (2-x**2-y**2)/np.sqrt(1-x**2-y**2)**3

x = SpatialCoordinate(mesh)
f = F(x[0], x[1])​
I get an error message
AttributeError: 'Sum' object has no attribute 'sqrt'​

If I dont use the sqrt function of numpy there is no problem, but in my real problem I use numpy, scipy etc.

written 5 months ago by BennyS
1
Specifically which functions are you looking for?

In the ufl library there are different math functions already implement for ufl objects. For the above code try:

from ufl.mathfunctions import *

def F(x, y):
return (2-x**2-y**2)/Sqrt(1-x**2-y**2)**3

x = SpatialCoordinate(mesh)
f = F(x[0], x[1])​​

To see which functions are implemented run

import ufl
help(ufl.mathfunctions)
help(ufl)  # More detailed description
written 5 months ago by Hernán Mella
Thank you very much for your answer! My function f can not be expressed in analytical form. It contains a call to a Fortran program via numpys f2py, various root finding procedures etc. Therefore my function f is strictly numerical and it can not be expressed in terms of the functions which are available by ufl.
written 5 months ago by BennyS
1
Maybe you can use a Quadrature finite element to define f as an Expression:

# Function Space
deg = 3
V = FunctionSpace(mesh, 'P', deg)

qdeg = 3
F = FunctionSpace(mesh, qelement)

.
.
.

# Define variational problem
phi = TrialFunction(V)
v = TestFunction(V)
f = interpolate(Expression('(2-pow(x[0], 2)-pow(x[1], 2)) * ' +
'pow(1-pow(x[0], 2)-pow(x[1], 2), -1.5)', degree=deg), F)
L = -f*v*dx(degree=qdeg)

# Compute solution
phi = Function(V)
solve(a == L, phi, bcs)​

In the above code, qdeg is the order of the polynomial which is exactly integrated with the quadrature scheme (please someone correct me if I am wrong)

written 5 months ago by Hernán Mella
Yes, this works great! Thank you very much!
written 5 months ago by BennyS