Tormod Landet - FEniCS tag

Articles tagged with FEniCS

A linear 2D sloshing solver in FEniCS

on. 29 januar 2014
In the previous blog post I showed how to solve the Laplace equation in 2D with FEniCS. That example, while illustrative for how to use FEniCS, was not particularly exciting. So, I decided to show the simplest example I could think of that involves an actual example of fluid flow.

Before we start: Writing these posts is a learning exercise for me. If you are learning FEniCS or the finite element method, please do not consider this material authoritative. I still hope this can be useful for someone else, but correctness is not guaranteed. In addition no mesh sensitivity studies or rigorous validation has been performed. This work has been done in evenings when I have felt like it, and then producing something that is much more fun than carefully testing it afterwards. Just so you know before we start on the real content of this post.

Sloshing

Sloshing is the motion of a fluid with a free surface in an enclosed space, such as a tank or a glass of water. The term free surface simply means the top of the fluid is allowed to move freely.

The simplest way to simulate sloshing with FEniCS is with a 2D rectangular tank and incompressible, irrotational and inviscid fluid flow. This is called potential flow and we only need to solve the Laplace equation for the scalar velocity potential, normally called \(\phi\), and not the velocities and the pressure which can be calculated from the derivatives of the velocity potential function. This makes potential flow solvers fast, and hence popular, but of course not applicable to all problems of fluid flow!

\begin{equation*} \newcommand{\dif}{\mathrm{d}} \newcommand{\pdiff}[2]{\frac{\mathrm\partial#1}{\mathrm\partial#2}} \newcommand{\laplace}[1]{\nabla^2#1} \laplace{\phi} = 0 \qquad \text{(the Laplace equation)} \end{equation*}

To further simplify our problem we assume that the motion of the free surface is so small that we do not make a large error by not moving the mesh near the surface when waves travel across it. If the waves are small then the mesh points at the still water surface are sufficiently close to the actual position of the free surface that we can get away by having a constant mesh.

The tank with free surface and still water water surface marked

We need to apply linear free surface boundary conditions to get away with not moving the mesh. A quick overview of the linear free surface boundary conditions can be found for example on the Wikipedia page for Airy wave theory.

The free surface introduces a new unknown function, \(\eta(x,t)\) in addition to the velocity potential \(\phi\). This is the surface elevation, and in our simple problem the surface elevation is only a function of the x-position in the tank and the time, t.

The linear kinematic free surface boundary condition says that the free surface moves up and down with the velocity of the fluid at the free surface. Particles on the free surface stay on the free surface.

\begin{equation*} \newcommand{\dif}{\mathrm{d}} \newcommand{\pdiff}[2]{\frac{\mathrm\partial#1}{\mathrm\partial#2}} \newcommand{\laplace}[1]{\nabla^2#1} \pdiff{\eta}{t} = \pdiff{\phi}{z} \quad \text{at the free surface} \end{equation*}

The z-derivative of the velocity potential is equal to the vertical fluid speed as described in the previous blog post.

The dynamic free surface boundary condition is derived from Bernoulli's equation and makes sure the pressure at the free surface is constant. We assume that the atmospheric pressure is zero:

\begin{equation*} \pdiff{\phi}{t} + g\eta = 0 \quad \text{at the free surface} \end{equation*}

We also need to specify that there is no motion through the floor of the tank:

\begin{equation*} \pdiff{\phi}{z} = 0 \quad \text{at the the bottom} \end{equation*}

And finally we want to move the tank back and forth to produce some free surface motion. We do this by introducing a horizontal velocity \(U_H(t)\) on the walls:

\begin{equation*} \pdiff{\phi}{x} = U_H(t) \quad \text{on the walls} \end{equation*}

One nice property of our equation system is that only our free surface boundary conditions involve time, the Laplace equation in the domain itself is independent of what happened at the last time step.
Solution algorithm

From the above we can decide on the following way to make a simulator for this problem:
1. We need a Laplace solver that can calculate the velocity potential given a position of the free surface and a wall velocity, \(U_H\). We specify the value of the velocity potential at the free surface by using the dynamic free surface boundary condition.
2. We need to update the free surface height function \(\eta(x,t)\) by using the kinematic free surface condition.
3. We need to advance our solution in time and run the above in a loop.
For point 1 we can use the previously shown FEniCS Laplace solver, with some small modifications. For point 2 we need to find a way to get the vertical derivative of the velocity potential at the free surface. For the last point we will use the built in time integration routines in SciPy so that we do not have to worry about selecting an appropriately small time step.

We will not apply any filtering of high frequent noise on the free surface function \(\eta\), something which is common to do in Laplace equation solvers with a free surface. We will hence need a quite fine time step to avoid numerical problems. The SciPy routines have automatic time step control which avoids this problem by reducing the time step when noise is trying to sneak in. This simplification in the implementation makes our code a bit slower that it could have been.

For the purpose of time stepping we define \(2 N\) unknowns, where \(N\) is the number of vertices on the free surface. The first half of the unknowns are the free surface elevations, while the second half is the velocity potential at the free surface:

\begin{align*} \pdiff{\eta}{t} &= \pdiff{\phi}{z} &\qquad \text{At each of the N free surface vertices}\\ \pdiff{\phi}{t} &= - g\eta &\qquad \text{At each of the N free surface vertices} \end{align*}

We take the initial velocity potential to be zero throughout the domain and the surface elevation is also zero at the start of the simulation.
Derivatives

Finding derivatives of the computed unknown function is covered in the FEniCS tutorial. Basically you can find the derivative \(w\) of the unknown function \(u\) by solving:

\begin{equation*} w = \pdiff{u}{t} \end{equation*}

which has the following weak form

\begin{equation*} \int_\Omega w\,v\,\dif \Omega = \int_\Omega \pdiff{u}{t}\,v\,\dif \Omega \end{equation*}

This can be written as follows in FEniCS:
```
# The z-derivative or the velocity potential
u2 = df.TrialFunction(V)
v2 = df.TestFunction(V)
a2 = u2*v2*df.dx
L2 = df.Dx(phi, 1)*v2*df.dx
```
This allows us to solve for the z-derivative of the solution, phi, in the whole domain. We only really need the derivative at the free surface, but for this demo I did not bother to spend time to figure out how do that.
Code

The complete code is available from here. I have tried to comment the code well so I hope it is relatively understandable, at least after reading these two blog posts.

If you look carefully you will see that the mesh has been slightly changed from the last post. The diagonals of the triangular elements are now alternating to left and right. Since the problem is highly symmetric this seemed to give significantly better results, at least for the test case I used while developing the code. Testing with higher order elements and quadrilaterals instead of triangles would be interesting, but has not been done.

At resonance the solution does not converge, the motion amplitudes are growing linearly with time, which is expected. Below I have shown surface elevation amplitudes as function of the excitation frequency. It is simply the elevation amplitude in the last oscillation cycle that is shown, so the results around resonances depend on how many oscillations have been run. Some damping could be added to the free surface boundary condition to get converged results.

Some results

The results for a selected test case can be seen below in terms of surface elevation and the velocity potential:

Animation of the free surface elevation

Elevation at each end of the tank as a function of time

Snapshot of the velocity potential at t=100

A series of simulations were performed and the amplitudes of the free surface motion close to the walls was extracted. The results can be seen below. The resonant periods calculated by analytical formula have also been added for comparison:

Transfer function for free surface elevation at the walls. Elevation amplitude divided by the forced motion amplitude.

Analytical natural frequencies can be found from the following formula [Faltinsen09]:

\begin{equation*} \omega_n = \sqrt{g\,k_n\,\tanh(k_n\,h)} \end{equation*}

where \(k_n = n \pi / L\). The mode number is \(n\) where 1 gives the first natural mode, corresponding to a standing half sine wave. The dimensions involved are the length \(L\) and the height of the tank, \(h\).

The calculations can be seen to be predicting roughly the physics they are intended to, but for high frequent motion the resulting velocity potential is not captured well by the selected discretization and these results are probably not very reliable.

Calculating the above transfer function with 300 excitation periods and 20 oscillations for each period on a 40 by 20 mesh with linear triangle elements took more that 12 hours on my laptop. For this specific problem it would probably take much less time to look up and implement the analytical solution, which is not terribly hard to do. The accuracy would also be much better.

I expect that the results given in the plot above are quite unreliable. This could easily be verified by comparing to the analytical solution, but I have not bothered with this here, as the intention was to test FEniCS for a simple dynamic problem, not promote a hopelessly inefficient method of analyzing linear sloshing in a rectangular tank.

I do, however, believe that the method presented can be made sufficiently accurate with the right mesh and/or higher order elements and perhaps more carefull treatment of the derivative of the velocity potential on the free surface. It is then possible to incorporate nonlinear free surface boundary conditions, find solutions for non-regular tank geometries and in general quite easily find answers to questions that are very hard to find by analytical work only.

I hope you have found this and the previous article usefull. I will probably restrict myself to smaller, bite sized, chunks if I am to continue with these posts. The topic of this article took more than one evening to code and write, and I am not sure I will devode that much time to this blog in the future. There is simply too much else to do and so many things to learn! Explaining something is still a great way to learn, so we will see if there is a follow up. Maybe the same problem, just with a discontinous Galerkin implementation instead of a continous? Time will tell.

[Faltinsen09] Faltinsen, OM and AN Timokha, Sloshing, Cambridge University Press, 2009
Solving the Laplace equation with FEniCS

ma. 13 januar 2014
I have found that explaining something is a good way to learn. As I am currently reading some books on the finite element method to refresh my rusting memories from math courses at university I thought it would be a good idea to write a summary of some of the basics. This post will explain how to use the finite element method and FEniCS for solving partial differential equations. FEniCS is a collection of tools for automated, efficient solution of partial differential equations using the finite element method.

Introduction

I will show how to apply FEniCS to a simple partial differential equation. I will also explain the solution method employed in the finite element method to solve the selected equation, the Laplace equation:

\begin{equation*} \newcommand{\dif}{\mathrm{d}} \newcommand{\pdiff}[2]{\frac{\mathrm\partial#1}{\mathrm\partial#2}} \newcommand{\laplace}[1]{\nabla^2#1} \laplace{\phi} = 0 \end{equation*}

When a fluid is incompressible, inviscid and the flow can be assumed to be irrotational then the Laplace equation shown above describes the fluid flow. The fluid velocities, \({\bf U} = (U_x, U_y, U_z)\), can be calculated from the velocity potential \(\phi\) as follows:

\begin{equation*} U_x = \pdiff{\phi}{x} \qquad U_y = \pdiff{\phi}{x} \qquad U_z = \pdiff{\phi}{x} \end{equation*}

Our domain, \(\Omega\), is a 2D square with \(x \in [0, 2]\) and \(y \in [0, 1]\). For the boundary conditions, let's say that we know the value of \(\phi\) at that the short boundaries and the normal derivative at the long edges is zero. The complete equation system can then be for example:

\begin{align*} &\laplace{\phi} = 0 \\ \ \\ &\phi = -1.0 \quad &\text{on} \quad x = 0.0 \\ &\phi = +1.0 \quad &\text{on} \quad x = 2.0 \\ &\pdiff{\phi}{\bf n} = 42 \quad &\text{on} \quad \Gamma_L \\ \end{align*}

where we define \(\Gamma_L\) to be the long edge boundaries \(y=0\) and \(y=1\). The normal vector on the boundary is \(\bf n\) and \(\pdiff{}{\bf n}\) is the normal derivative at the boundary. It is common to use \(\Omega\) to denote the domain and \(\Gamma\) to denote the boundary.
Discretization

Now, to solve the equation with boundary conditions we need a discretized domain. This can be created by using the FEniCS Python package dolfin like this:
```
import dolfin as df

start_x, end_x, n_elem_x = 0.0, 2.0, 4
start_y, end_y, n_elem_y = 0.0, 1.0, 2

mesh = df.RectangleMesh(start_x, start_y, end_x, end_y, n_elem_x, n_elem_y)
```
Running df.plot(mesh) will show the following:

An example mesh with triangular elements

We discretize our domain \(\Omega\) into discrete elements and obtain \(\hat \Omega\). I will use the "hat" to denote discretized values. We use linear elements which means we will approximate the unknown function with piecewise linear functions, \(\phi \approx \hat \phi\). They look like this for a 1D domain with two elements and three vertices (nodes):

1D domain with three shape functions, \(N_1(x)\), \(N_2(x)\) and \(N_3(x)\)

We have one shape function \(N_i\) for each vertex in the mesh. The shape functions are equal to 1.0 at the corresponding vertex and 0.0 at all other vertices. In 2D a standard bi-linear shape function looks as follows:

2D domain with a bi-linear shape function shown for one of the vertices on the boundary
Marking the boundaries

We know how to make a mesh of the domain. We also need to tell FEniCS about the different parts of the boundary. We do this by giving the boundaries different unique integer numbers and marking the boundaries with a mesh function:
```
# Each part of the mesh gets its own ID
ID_LEFT_SIDE, ID_RIGHT_SIDE, ID_BOTTOM, ID_TOP, ID_OTHER = 1, 2, 3, 4, 5
boundary_parts = df.MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundary_parts.set_all(ID_OTHER)

# Define the boundaries
GEOM_TOL = 1e-6
define_boundary(lambda x, on_boundary: on_boundary and x[0] < start_x + GEOM_TOL, boundary_parts, ID_LEFT_SIDE)
define_boundary(lambda x, on_boundary: on_boundary and x[0] > end_x - GEOM_TOL, boundary_parts, ID_RIGHT_SIDE)
define_boundary(lambda x, on_boundary: on_boundary and x[1] < start_y + GEOM_TOL, boundary_parts, ID_BOTTOM)
define_boundary(lambda x, on_boundary: on_boundary and x[1] > end_y - GEOM_TOL, boundary_parts, ID_TOP)
```
The define_boundary function is just a small utility I made to make it easier to mark the boundaries. Normally you need to create a class for each boundary and implement a method, inside, on this class. My function accepts a lambda expression which is used instead:
```
def define_boundary(defining_func, boundary_parts, boundary_id):
    """
    Define a boundary

    - defining_func is used instead of the normal SubDomain inside method: defining_func(x, on_boundary)
    - boundary_parts is a MeshFunction used to distinguish parts of the boundary from each other
    - boundary_id is the id of the new boundary in the boundary_parts dictionary
    """
    class Boundary(df.SubDomain):
        def inside(self, x, on_boundary):
            return defining_func(x, on_boundary)
    boundary = Boundary()
    boundary.mark(boundary_parts, boundary_id)
    return boundary
```
The solution method

The simplest variant of the finite element method to describe is probably the Galerkin weighted residuals method. The core of the method is described below.

The unknown function \(\phi\) is written as a sum of \(M\) shape functions and a set of unknown coefficients, \(a_i\). There is one unknown coefficient for each shape function:

\begin{equation*} \hat \phi(x,y) = \sum_{i=0}^{M} a_i N_i(x,y) \end{equation*}

Suppose we guess some coefficients \(a_i\). We can then integrate our trial function, \(\hat \phi\) over the domain, but since the discretization is not perfect we get a residual and not zero as the equation was supposed to give:

\begin{equation*} \int_{\Omega} \laplace{\hat \phi} \ \dif \Omega = R \end{equation*}

To be able to fix this and find some coefficients \(a_i\) that gives zero residual we need a set of equations. This is done by multiplying the above integral with a set of weighing functions, \(w_i(x,y)\). The number of weighting functions is equal to the number of unknowns, so we get a fully determined set of linear equations. In the Galerkin method the weighting functions are taken equal to the shape functions, \(N_i\) except that we do not include shape functions corresponding to vertices where we know the solution \(\phi(x,y)\) from the Dirichlet boundary conditions.

We then get a set of \(M\) equations:

\begin{equation*} \int_{\Omega} \laplace{\hat \phi} w_i \ \dif \Omega = 0 \end{equation*}

which is the same as:

\begin{equation*} \int_{\Omega} \laplace{\left(\sum_{j=0}^M a_j N_j\right)} N_i \dif \Omega = 0 \end{equation*}

Since the integral will only be non-zero for vertices that are close in the mesh, i.e. \(N_i\) and \(N_j\) have some overlap, we end up with a quite sparse equation system. Very efficient solvers exist for the types of equation systems that come out of the Galerkin method.
Solving the equation system with FEniCS

In the standard finite element terminology our unknown function \(\phi\) is called the trial function and denoted \(u\) while the weighting function is called the test function and denoted \(v\). We switch to this terminology to be compatible with the FEniCS tutorial. We call this equation the weak form:

\begin{equation*} \int_{\Omega} \laplace{(u)}\,v \ \dif \Omega = 0 \end{equation*}

Before we are done we need to modify the equation a bit. As written it makes little sense since the Laplacian of the bi-linear function \(N_j(x,y)\) will always be zero. We modify our weak form by performing integration by parts. The formula for this in higher dimensions is:

\begin{equation*} \int_\Omega \pdiff{f}{x_i}\,g \ \dif \Omega = \int_\Gamma f\,g\,n_i\ \dif \Gamma - \int_\Omega f \pdiff{g}{x_i} \ \dif \Omega \end{equation*}

where \(n_i\) is the \(i\)-th component of the outwards normal vector on \(\Gamma\). I do not think I ever learned that at school? Maybe it's just my rusty memory, but thanks to Wikipedia we will not let lacking education and/or failing memory stop us, now will we? We get this when applying integration by parts to our integral:

\begin{equation*} \int_\Omega \laplace{(u)}\,v \ \dif \Omega = \int_\Gamma \nabla(u)\,v\,n_i\ \dif \Gamma - \int_\Omega \nabla (u)\cdot\nabla(v) \ \dif \Omega \end{equation*}

I guess now is the time for me to remember to mention that the test function \(v\) is always zero at the Dirichlet boundaries. Dirichlet boundaries are boundaries where we have Dirichlet boundary conditions, which means we know the value of \(\phi\), i.e. the short edges of our mesh in this case. No reason to have a weight function there, right? We now have the final expression for the weak form:

\begin{align*} \int_{\Omega} \nabla u \cdot \nabla v \ \dif \Omega &= \int_{\Gamma_N} \laplace{(u)}\,v\,n_i\ \dif \Gamma_N \\ &= \int_{\Gamma_N} \pdiff{u}{\bf n}\,v\ \dif \Gamma_N \end{align*}

The long edges were we prescribe the normal velocity (normal derivative of the velocity potential \(\phi\)) are now denoted \(\Gamma_N\) since this type of boundary condition is called a Neumann boundary condition.

FEniCS understands equations such as the above and calls them forms. To describe the equations we use UFL, Uniform Form Language which is a Python API for defining forms. When interacting with FEniCS through Python we can use UFL through the dolfin package just like the mesh facilities. UFL can also be directly imported, import ufl.
```
# Define a function space of bi-linear elements on our mesh
V = df.FunctionSpace(mesh, 'Lagrange', 1)

# The trial and test functions are defined on the bi-linear function space
u = df.TrialFunction(V)
v = df.TestFunction(V)

# Our equation is defined as a=L. We only define what is inside the integrals, and the difference
# between integrals over the domain and the boundary is whether we use dx or ds:
#   - the term a contains both the trial (unknown) function u and the test function v
#   - the term L contains only test function v
a = df.inner(df.nabla_grad(u), df.nabla_grad(v))*df.dx
L = df.Constant(42)*v*df.ds(ID_BOTTOM) + df.Constant(42)*v*df.ds(ID_TOP)
```
This is basically all there is to solving the Laplace equation in FEniCS. The complete code can be seen below. The parts I have skipped explaining above are well explained in the FEniCS tutorial.
Complete code

The complete code can be found on Bitbucket. The resulting plot from running the code with FEniCS 1.3 is shown below:

Plot of the resulting \(\phi\) function

If you have any questions to the above and how to use FEniCS you can ask here, but I would really recommend the FEniCS QA Forum where you will probably get better and faster answers.

Creating a mesh in FEniCS

I may end up doing a large project using FEniCS which is a collection of tools for automated, efficient solution of differential equations using the finite element method. While playing with implementing a simple solver for fluid flow in a tank by use of FEniCS to solve the Laplace equation I needed to create my own mesh in code.

FEniCS is quite well documented, but I had to look at the source code for some of the mesh conversion routines to find out how to build a mesh from scratch. So, for posterity, here is a reimplementation of dolfin.RectangleMesh in a much more cumbersome (and flexible) way:

import numpy
import dolfin as df

def create_mesh(length, height, nx, ny, show=False):
    """
    Make a square mesh manually

    Should give exactly the same results as using the built in dolfin.RectangleMesh() class
    """
    # The number of mesh entities
    nvert = nx*ny
    ncell = 2*(nx-1)*(ny-1)

    # Positions of the vertices
    xpos = numpy.linspace(0, length, nx)
    ypos = numpy.linspace(0, height, ny)

    # Create the mesh and open for editing
    mesh = df.Mesh()
    editor = df.MeshEditor()
    editor.open(mesh, 2, 2)

    # Add the vertices (nodes)
    editor.init_vertices(nvert)
    i_vert = 0
    for x in xpos:
        for y in ypos:
            editor.add_vertex(i_vert, x, y)
            i_vert += 1

    # Add the cells (triangular elements)
    # Loop over the vertices and build two cells for each square
    # where the selected vertex is in the lower left corner
    editor.init_cells(ncell)
    i_cell = 0
    for ix in xrange(nx-1):
        for iy in xrange(ny-1):
            # Upper left triangle in this square
            i_vert0 = ix*ny + iy
            i_vert1 = ix*ny + iy+1
            i_vert2 = (ix+1)*ny + iy + 1
            editor.add_cell(i_cell, i_vert0, i_vert1, i_vert2)
            i_cell += 1

            # Lower right triangle in this square
            i_vert0 = ix*ny + iy
            i_vert1 = (ix+1)*ny + iy+1
            i_vert2 = (ix+1)*ny + iy
            editor.add_cell(i_cell, i_vert0, i_vert1, i_vert2)
            i_cell += 1

    # Close the mesh for editing
    editor.close()
    print 'Created mesh with %d vertices and %d cells' % (nvert, ncell)

    if show:
        df.plot(mesh)
        df.interactive()

    return mesh

If anyone finds this by searching the web for how to build a mesh programmatically in FEniCS/dolfin, then I hope the above was understandable. The documentation is quite good when you know that you need to search for the MeshEditor class ...

Articles tagged with FEniCS

Sloshing

Solution algorithm

Derivatives

Code

Some results

Introduction

Discretization

Marking the boundaries

The solution method

Solving the equation system with FEniCS

Complete code