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

The most comfortable solution for password-less logins outside of the university:

• add

  Host login
     Hostname login.tu-chemnitz.de
     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
     ProxyCommand ssh login -W %h:%p
  

  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

and during its validity password-less login from anywhere

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
Type "help", "copyright", "credits" or "license" for more information.
>>>

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

Task 1

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

Task 2

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.

Hint

Constant, Expression, and similar are clickable links leading to their documentation.

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

• definition of weak formulation through forms a and L,
• usage of Constant and Expression classes.

Task 3

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

Task 4

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.

Task 5

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

Task 6

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).

Task 7

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).\]

Task 8

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((
    (..., ...),
    (..., ...),
))
F = inner(A*grad(u), grad(v))*dx + ...

Reference solution

Show/Hide Code

Download Code

from dolfin import *
import matplotlib.pyplot as plt


def solve_task3(V):
    """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)
    a = inner(grad(u), grad(v))*dx + 6*u*v*dx
    L = f*v*dx + g*v*ds

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

    return u


def solve_task5(V):
    """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


def solve_task6(V):
    """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)
    F = c*inner(grad(u), grad(v))*dx + (10*u**3 + u)*v*dx - f*v*dx

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

    return u


def solve_task7(V):
    """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


def solve_task8(V):
    """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),
    ))
    F = inner(A*grad(u), grad(v))*dx + u*v*dx - f*v*dx

    # 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
    u3 = solve_task3(V)
    u5 = solve_task5(V)
    u6 = solve_task6(V)
    u7 = solve_task7(V)
    u8 = solve_task8(V)

    # 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()