# Tormod Landet

1. ## Installing R in WSL

Short note for posteriority: the documentation for R tells you to install the key for the Ubuntu 18.04 R repository via the following command:

sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys E298A3A825C0D65DFD57CBB651716619E084DAB9


Using apt-key in this way does not work on Ubuntu 18.04 on Windows 10 (in WSL, Windows Subsystem for Linux) due to use of an unsupported IPC mechanism for communicating with dirmngr (at least with Windows 10 version 1709). The command gives the following error:

Executing: /tmp/apt-key-gpghome.Sf0LlrpiEu/gpg.1.sh --keyserver keyserver.ubuntu.com --recv-keys E298A3A825C0D65DFD57CBB651716619E084DAB9
gpg: connecting dirmngr at '/tmp/apt-key-gpghome.Sf0LlrpiEu/S.dirmngr' failed: IPC connect call failed
gpg: keyserver receive failed: No dirmngr


curl -sL "http://keyserver.ubuntu.com/pks/lookup?op=get&search=0x51716619E084DAB9" | sudo apt-key add


How did I get the number 0x51716619E084DAB9? Add the repository and run apt update:

sudo add-apt-repository 'deb https://cloud.r-project.org/bin/linux/ubuntu bionic-cran35/'
sudo apt update


It will complain:

Err:2 https://cloud.r-project.org/bin/linux/ubuntu bionic-cran35/ InRelease
The following signatures couldn't be verified because the public key is not available: NO_PUBKEY 51716619E084DAB9


Solution:

The full set of commands to install R from the r-project.org repositories:

curl -sL "http://keyserver.ubuntu.com/pks/lookup?op=get&search=0x51716619E084DAB9" | sudo apt-key add
sudo apt-get update
sudo apt-get install r-base


I'm not really much of an R user, but I needed it this one time and I prefer using WSL when I have to use Windows.

2. ## Ocellaris and my PhD project

Ocellaris, the free surface two phase Navier-Stokes solver I have been building as a part of my PhD work at the Department of Mathematics at the University of Oslo is finally ready for the world! It is built using an exactly divergence free unstructured DG FEM numerical method that uses special velocity slope limiters to stabilise the higher order shape functions near the factor 1000 jump in momentum at the air/water free surface, which enables higher order simulations of free surface flows such as ocean waves. I just released version 2019.0.1, which has greatly improved the user guide and general documentation.

You can see an example simulation below, more examples are available on the Ocellaris Blog.

The PhD thesis and the included papers are all soon ready as pre-prints. Looking forward to delivering this winter!

3. ## Cleaning cruft in a Ubuntu installation

I have been running Linux on my laptops for as long as I have been owning laptops. My current one has been through quite a few versions of Ubuntu Linux. There are normally no big changes from version to version, but it steadily improves. In that process some software packages that previously were installed by default are replaced or become obsolete. Most of this software is still updated, but it is no longer present on new installations and it may have no value to me.

Every time I update to a new version of Ubuntu I have a feeling that the cruft is building up, but there is no obvious way that I have found to remove software that has no longer been found worthy of being part of the default install. So, this time I decided to write a script to identify these packages. The script takes a list of software packages to keep, i.e the software that I use directly when I use my machine. From this list it finds all the packages that must be present for my prefered software to work, and lists all the other that can in princible be removed.

In reality, it is hard to come up with a list of all the software that is required for a Ubuntu installation to work properly, so I also made the script read the newest list of default packages. Any of these that happen to be installed are also kept along with their dependencies.

The script itself, unnecessarry-packages.py can be found on Bitbucket. It does not uninstall any packages by itself, but suggest packages that you can uninstall by running apt purge PCKGNAME or similar. I put it online and write about it here in the hope that it will be useful to someone else.

I started out with the following list of packages to keep,

python3 unnecessarry-packages.py ubuntu-desktop krita thunderbird firefox ubuntu-minimal digikam gimp darktable gedit Downloads/ubuntu-16.10-desktop-amd64.manifest -v


and after some trial and error I ended up with the following by noticing packages suggested for removal that I in reality wanted to keep and adding them to the command line invocation,

python3 unnecessarry-packages.py ubuntu-desktop krita thunderbird firefox ubuntu-minimal digikam gimp darktable gedit clang bcmwl-kernel-source systemsettings enfuse ffmpegthumbs flashplugin-installer gstreamer1.0-nice gstreamer1.0-plugins-bad gstreamer1.0-plugins-bad-faad  gstreamer1.0-plugins-bad-videoparsers gstreamer1.0-plugins-ugly gstreamer1.0-plugins-ugly-amr hugin hugin-tools inkscape gimp-ufraw git  gmsh gphoto2 gstreamer1.0-libav gimp-plugin-registry  ttf-mscorefonts-installer mercurial kipi-plugins gimp-gmic sshfs Downloads/ubuntu-16.10-desktop-amd64.manifest -v


This allowed me to get rid of approximately 4 GB of cruft without any noticeable problems afterwards. Note that I did not blindly remove all packages suggested by the above final invocation, but I removed maybe 90% of them. The script is not smart, it can easily suggest that you remove your network driver, so use with caution!

4. ## 3D visualization in IPython notebook

The IPython notebook is a great tool for interactive computing. Below I show a short Python example on how to work interactively with simple 3D models in the IPython notebook.

To begin with, we need a simple way to represent a 3D model in Python. We will here restrict ourselves to geometries made up of plane 2D polygons in 3D:

class PolygonCollection(object):
def __init__(self, nodes, polygons):
"""
A collection of plane polygons in 3D

Nodes is a list of 3-tuples of coordinates
Polygons is a nested list of node indices for each polygon, one
list of indices for each polygon
"""
self.nodes = nodes
self.polygons = polygons

def extents(self):
"""
Calculate the size of the 3D polygon model
"""
xmin, ymin, zmin = 1e100, 1e100, 1e100
xmax, ymax, zmax = -1e100, -1e100, -1e100
for x, y, z in self.nodes:
xmin = min(x, xmin)
ymin = min(y, ymin)
zmin = min(z, zmin)
xmax = max(x, xmax)
ymax = max(y, ymax)
zmax = max(z, zmax)
return xmin, xmax, ymin, ymax, zmin, zmax


We then populate our 3D model with some polygons. The polygon <-> node/vertex connectivity is simply represented as indices into the list of nodes/vertices:

nodes = [(0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 1.0, 0.0), (0.0, 1.0, 0.0),
(1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (0.0, 0.0, 1.0),
(0.0, 1.0, 0.0), (0.0, 1.0, 1.0), (0.0, 0.0, 1.0)]
elements = [(0, 1, 2, 3), (0, 4, 5, 6), (0, 7, 8, 9)]
polygons = PolygonCollection(nodes, elements)


IPython notebook allows us to embed images, Latex equations, HTML, Javascript and more. In this case we will use Javascript to render our 3D model. This can be done as follows:

from IPython.display import display, Javascript
js = polygons2js(polygons)
display(Javascript(js))


The end result can be seen below. In a reasonable modern browser you should be able to rotate the 3D model. If it does not work it may also be that the Javascript libraries imported to do the rendering have changed since the time when this was written. I do not believe the libraries have stable APIs just yet. As of May 2014 it should work in the latest versions of Firefox, Chrome and Internet Explorer.

Update March 2017: fixed the javascript (and the below Python code) so that it works with the latest (Nov 2016) version of JSModeler. The python code below should produce the same Javascript code that is running on this page, but it has not been tested, so there might be a typo somewhere.

The "only" thing needed to make this work is a function polygons2js(polygons) that converts the 3D model to a set of javascript commands that will run inside the notebook. This can be done quite easily with the help of JSModeler and three.js. A simple implementation will look something like this:

import random

def js_wrap(js, js_libraries):
"""
Wrap javascript commands in code that preloads needed libraries
"""
lines = ['function getMultipleScripts(scripts, callback) {',
'  if(scripts.length == 0)',
'    callback();',
'  jQuery.getScript(scripts[0], function () {',
'    getMultipleScripts(scripts.slice(1), callback);',
'  });',
'}',
'var scripts = [']
lines += ['"%s",' % script for script in js_libraries]
lines += ['];',
'getMultipleScripts(scripts, function() {',
js,
'});']
return '\n'.join(lines)

def polygons2js(polygons):
# Find the size of the 3D model. We use this later to make the default camera and rotation point work
xmin, xmax, ymin, ymax, zmin, zmax = polygons.extents()
length = max([xmax-xmin, ymax-ymin, zmax-zmin])
xmean = (xmin + xmax)/2
ymean = (ymin + ymax)/2
zmean = (zmin + zmax)/2

# Generate the Javascript code
# see http://kovacsv.github.io/JSModeler/documentation/tutorial/tutorial.html for details
canvas_id = 'js3dcanvas_%d' % random.randint(0, 1e10)
canvas_style = 'style="border: 1px solid black;" width="800" height="600"'
js = ['var widget = jQuery(\'<canvas id="%s" %s></canvas>\');' % (canvas_id, canvas_style),
'element.append(widget);',
'container.show();',
'var viewerSettings = {',
'  cameraEyePosition : [-2.0, -1.5, 1.0],',
'  cameraCenterPosition : [0.0, 0.0, 0.0],',
'  cameraUpVector : [0.0, 0.0, 1.0]',
'};',
'var viewer = new JSM.ThreeViewer();',
'viewer.Start(widget, viewerSettings);',
'var body = new JSM.Body();']

coords = []
for coord in polygons.nodes:
coords.append('[%f, %f, %f]' % ((coord[0]-xmean)/length,
(coord[1]-ymean)/length,
(coord[2]-ymean)/length))
#print coord, coords[-1]
js.append('var coords = [%s];' % ', '.join(coords))
js.append('for (var i = 0, len = coords.length; i < len; i++) {')
'new JSM.Coord(coords[i][0], coords[i][1], coords[i][2])));')
js.append('}')

polys = [repr(list(poly)) for poly in polygons.polygons]
js.append('var elems = [%s];' % ', '.join(polys))
js.append('for (var i = 0, len = elems.length; i < len; i++) {')
js.append('}')
js.append('var meshes = JSM.ConvertBodyToThreeMeshes(body);')
js.append('for (var i = 0; i < meshes.length; i++) { viewer.AddMesh (meshes[i]); }')
js.append('viewer.Draw();')

# Wrap final js code in a library loader to make sure JSModeler and three.js are available
js = '\n'.join(js)
libraries =  ['https://rawgit.com/kovacsv/JSModeler/master/build/lib/three.min.js',
'https://rawgit.com/kovacsv/JSModeler/master/build/jsmodeler.js',
'https://rawgit.com/kovacsv/JSModeler/master/build/jsmodeler.ext.three.js']
return js_wrap(js, libraries)


I have found this to be quite handy for inspecting small 3D models and modifying the interactively. The PolygonCollection class can also be updated to automatically display the 3D model in the notebook if it objcts of the class are left on the last line of an IPython input cell, or if you call display() on them. For this to work you will have to extend the PolygonCollection class with a _repr_javascript_ metod that returns polygons2js(self) (as a string). IPython will automatically look for this method and call display(Javascript(polygons._repr_javascript_())) for you.

5. ## A linear 2D sloshing solver in FEniCS

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!


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.


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

Page 1 / 3 »