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

232

views

1

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:

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.

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)
a = dot(grad(phi), grad(v))*dx
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

### 1 Answer

4

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)
a = dot(grad(phi), grad(v))*dx
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

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:

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
qelement = FiniteElement(family="Quadrature", cell=mesh.ufl_cell(), degree=qdeg, quad_scheme="default")
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)
a = dot(grad(phi), grad(v))*dx(degree=qdeg)
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

Please login to add an answer/comment or follow this question.