# Poisson in a hundred ways¶

## First touch¶

Login by SSH to tyche and type: How to login

Open a terminal window by hitting CTRL+ALT+T. Use ssh to connect to a remote system:

ssh -X -C tyche


Note

• -X enable X11 forwarding (allows processes on the remote machine opening windows of graphical applications on the local machine)

• -C enables compression which is mainly beneficial for access from a remote network

• tyche stands here for machine tyche.mathematik.tu-chemnitz.de; username on the local machine is used by default to login to the remote machine; the machine is not accessible from outside the univerity, so one would login through a jump host

luigivercotti@local_machine:~$ssh -X -C user@login.tu-chemnitz.de user@login:~$ ssh -X -C tyche
user@tyche:~\$


Alternatively one can use VPN.

For experts: Kerberos + public key + jump host

Host login
User <user>
ForwardAgent yes
ForwardX11 yes
ForwardX11Trusted yes
GSSAPIAuthentication yes
GSSAPIDelegateCredentials yes

Host tyche
Hostname tyche.mathematik.tu-chemnitz.de
User <user>
ForwardAgent yes
ForwardX11 yes
ForwardX11Trusted yes


to ~/.ssh/config,

• upload public key to login and tyche by

ssh-copy-id login
ssh-copy-id tyche

• install Kerberos

which allows to get Kerberos ticket by

kinit <user>@TU-CHEMNITZ.DE


ssh tyche

source /LOCAL/opt/fenics-2017.2.0/fenics.conf


to prepare environment for using FEniCS. Now fire up interactive Python 3 interpreter:

python3


You should see something like:

Python 3.6.5 (default, Apr  1 2018, 05:46:30)
[GCC 7.3.0] on linux
>>>


Now type:

>>> from dolfin import *
>>> import matplotlib.pyplot as plt
>>>
>>> mesh = UnitSquareMesh(13, 8)
>>> plot(mesh)
[<matplotlib.lines.Line2D object at 0x7fe0003d65c0>, <matplotlib.lines.Line2D object at 0x7fe0003d6748>]
>>> plt.show()


Hint

Click on >>> in the right top corner of the code snippet to make the code copyable.

A graphical plot of the mesh should appear. If any of the steps above failed, you’re not correctly set up to use FEniCS. If everything went fine, close the plot window and hit ^D to quit the interpreter.

## Run and modify Poisson demo¶

Get the Poisson demo from FEniCS install dir and run it:

mkdir -p work/fenics/poisson
cd work/fenics/poisson
cp /LOCAL/opt/fenics-2017.2.0/share/dolfin/demo/documented/poisson/python/demo_poisson.py .
python3 demo_poisson.py


You should see some console output and a plot of the solution.

Now login to tyche from another terminal window and open the demo file using your favourite editor (if you don’t have any you can use gedit, nano, …):

cd work/fenics/poisson
<editor> demo_poisson.py


Now add keyword argument warp='mode' to the plot function call by applying the following diff:

 # Plot solution
import matplotlib.pyplot as plt
-plot(u)
+plot(u, mode='warp')
plt.show()


and run the demo again by python3 demo_poisson.py.

Open Poisson demo documentation on the FEniCS website. Notice that the doc page is generated from the demo file. Go quickly through the docpage while paying attention to

Modify the code to solve the following problem instead:

\begin{align}\begin{aligned}-\Delta u + c u &= f &&\text{in } \Omega,\\u &= u_\mathrm{D} &&\text{on } \Gamma_\mathrm{D},\\\tfrac{\partial u}{\partial\mathbf{n}} &= g &&\text{on } \Gamma_\mathrm{N}\end{aligned}\end{align}

with

\begin{gather} \Omega = (0,1)^2, \qquad \Gamma_\mathrm{D} = \{(x, y), x=1, 0<y<1\}, \qquad \Gamma_\mathrm{N} = \partial\Omega\setminus\Gamma_\mathrm{D}, \\ c = 6, \qquad f(x, y) = x, \qquad u_\mathrm{D}(x, y) = y, \qquad g(x, y) = \sin(5x) \exp(y). \end{gather}

## Semilinear Poisson equation¶

Derive weak formulation for the following semilinear Poisson problem:

(1)\begin{align}\begin{aligned}-\Delta u + u^3 + u &= f &&\text{in } \Omega,\\\tfrac{\partial u}{\partial\mathbf{n}} &= g &&\text{on } \partial\Omega\end{aligned}\end{align}

with

(2)$\Omega = (0,1)^2, \qquad f(x, y) = x, \qquad g(x, y) = \sin(5x) \exp(y).$

Notice that the weak formulation has the form

Find $$u\in H^1(\Omega)$$ such that
$F(u; v) = 0 \qquad \text{for all } v\in H^1(\Omega)$

with certain $$F$$ depending on $$u$$ in nonlinear fashion but being linear in test functions $$v$$. One can find the solution iteratively by the Newton method:

1. Choose $$u_0\in H^1(\Omega)$$,

2. For $$k=1,2,\ldots$$ do

1. Find $$\delta u\in H^1(\Omega)$$ such that

(3)$\frac{\partial F}{\partial u}(u_k; v, \delta u) = -F(u_k; v) \qquad \text{for all } v\in H^1(\Omega),$
2. Set $$u_{k+1} = u_k + \delta u$$.

3. Check certain convergence criterion and eventually stop iterating.

Here Jacobian $$\frac{\partial F}{\partial u}(u; v, \delta u)$$ is Gâteaux derivative of $$F$$. It is generally nonlinear in $$u$$, but linear in $$v$$ and $$\delta u$$. Hence with fixed $$u_k\in H^1(\Omega)$$ the left-hand side and the right-hand side of (3) are a bilinear and linear form respectively and (3) is just ordinary linear problem.

Modify the previous code to adapt it to problem (1), (2). Define $$F$$ by filing the gaps in the following code:

u = Function(V)
v = TestFunction(V)
f = Expression(...)
g = Expression(...)

F = ...


If in doubts, peek into Nonlinear Poisson demo documentation.

Look into documentation of solve function, read section Solving nonlinear variational problems. Now you should be able to call the solve function to obtain the solution.

## Nonlinear Dirichlet problem¶

Modify the code to solve the following Dirichlet problem:

\begin{align}\begin{aligned}-\operatorname{div}(c\nabla u) + 10 u^3 + u &= f &&\text{in } \Omega,\\u &= u_\mathrm{D} &&\text{on } \partial\Omega\end{aligned}\end{align}

with

$\Omega = (0,1)^2, \qquad f(x, y) = 100 x, \qquad u_\mathrm{D}(x, y) = y, \qquad c(x, y) = \tfrac{1}{10} + \tfrac12(x^2+y^2).$

Hint

Supply instance of SubDomain class to DirichletBC. How do you tell SubDomain to define $$\partial\Omega$$? What do you fill in?

class Boundary(SubDomain):
def inside(self, x, on_boundary):
return ...


on_boundary argument evaluates to True on boundary facets, False otherwise.

## Variational formulation¶

For $$u\in H^1(\Omega)$$ consider functional

$E(u) = \int_\Omega \bigl( \tfrac12|\nabla u|^2 + \tfrac14u^4 + \tfrac12u^2 - fu \bigr) \,\mathrm{d}x - \int_{\partial\Omega} gu \,\mathrm{d}s.$

Convince yourself that minimization of $$F$$ over $$H^1(\Omega)$$ is equivalent to problem (1).

By filling the following code:

u = Function(V)
f = Expression(...)
g = Expression(...)

E = ...


define $$E(u)$$ for data (2). Remember that functionals (zero-forms) do not have any test and trial functions.

Obtain $$F(u;v):=\frac{\partial E}{\partial u}(u; v)$$ using derivative:

F = derivative(E, u)


and run the solver like in Task 5. Check you get the same solution.

## Yet another nonlinearity¶

Consider quasilinear equation in divergence form

(4)\begin{align}\begin{aligned}-\operatorname{div}(\mathcal{A}\nabla u) + u &= f &&\text{in } \Omega,\\\tfrac{\partial u}{\partial\mathcal{A}^\top\mathbf{n}} &= 0 &&\text{on } \partial\Omega,\\\mathcal{A} &= \begin{bmatrix} \tfrac{1}{10} + u^2 & 0 \newline 0 & 1 + u^2 \end{bmatrix} &&\text{in } \Omega\end{aligned}\end{align}

with data

(5)$\Omega = (0,1)^2, \qquad f(x, y) = \tfrac12(x+y).$

Derive weak formulation for the problem (4).

Solve the problem (4), (5) using FEniCS. Employ as_matrix function to define $$\mathcal{A}$$:

u = Function(V)
v = TestFunction(V)

A = as_matrix((
(..., ...),
(..., ...),
))


## Reference solution¶

Show/Hide Code

Download Code

from dolfin import *
import matplotlib.pyplot as plt

"""Return solution of Task 3 on space V"""

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x, on_boundary):
return on_boundary and x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
uD = Expression("x[1]", degree=1)
bc = DirichletBC(V, uD, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("x[0]", degree=1)
g = Expression("sin(5*x[0])*exp(x[1])", degree=3)
L = f*v*dx + g*v*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

return u

"""Return solution of Task 5 on space V"""

# Define variational problem
u = Function(V)
v = TestFunction(V)
f = Expression("x[0]", degree=1)
g = Expression("sin(5*x[0])*exp(x[1])", degree=3)
F = inner(grad(u), grad(v))*dx + (u**3 + u)*v*dx - f*v*dx - g*v*ds

# Compute solution
solve(F == 0, u)

return u

"""Return solution of Task 6 on space V"""

# Define Dirichlet boundary
class Boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary

# Define boundary condition
boundary = Boundary()
uD = Expression("x[1]", degree=1)
bc = DirichletBC(V, uD, boundary)

# Define variational problem
u = Function(V)
v = TestFunction(V)
c = Expression("0.1 + 0.5*(x[0]*x[0] + x[1]*x[1])", degree=2)
f = Expression("100*x[0]", degree=1)

# Compute solution
solve(F == 0, u, bc)

return u

"""Return solution of Task 7 on space V"""

# Define variational problem
u = Function(V)
f = Expression("x[0]", degree=1)
g = Expression("sin(5*x[0])*exp(x[1])", degree=3)
E = ( grad(u)**2/2 + u**4/4 + u**2/2 - f*u )*dx - g*u*ds
F = derivative(E, u)

# Compute solution
solve(F == 0, u)

return u

"""Return solution of Task 8 on space V"""

# Define variational problem
u = Function(V)
v = TestFunction(V)
f = Expression("0.5*(x[0] + x[1])", degree=1)
A = as_matrix((
(0.1 + u**2,        0),
(         0, 1 + u**2),
))

# Compute solution
solve(F == 0, u)

return u

if __name__ == '__main__':

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

# Solve all problems

# Compare solution which should be same
err = ( grad(u7 - u5)**2 + (u7 - u5)**2 )*dx
err = assemble(err)
print("||u7 - u5||_H1 =", err)

# Plot all solutions into separate figures
plt.figure()
plot(u3, title='u3', mode="warp")

plt.figure()
plot(u5, title='u5', mode="warp")

plt.figure()
plot(u6, title='u6', mode="warp")

plt.figure()
plot(u7, title='u7', mode="warp")

plt.figure()
plot(u8, title='u8', mode="warp")

# Display all plots
plt.show()