src/examples/tangaroa.c

Two-phase flow around RV Tangaroa

This is an improved version of the famous Gerris example, illustrating the combination of complex solid boundaries, air-water turbulent flows and reduced gravity approach.

We use the centered Navier–Stokes solver, two-phase flow and the momentum-conserving option. Note that the momentum-conserving option is crucial to obtain stable solutions for this air-water density ratio configuration.

#include "grid/octree.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "navier-stokes/conserving.h"

We also need to compute distance functions (to describe the ship geometry), use reduced gravity and visualisation functions.

#include "distance.h"
#include "reduced.h"
#include "view.h"
#include "lambda2.h"

On supercomputers we need to control the maximum runtime and we check performances.

#include "maxruntime.h"
#include "navier-stokes/perfs.h"

Importing the geometry

This function computes the solid fraction given a pointer to an STL file, a tolerance (maximum relative error on distance) and a maximum level.

void fraction_from_stl (scalar f, FILE * fp, double eps, int maxlevel)
{

We read the STL file and compute the bounding box of the model.

  coord * p = input_stl (fp);
  coord min, max;
  bounding_box (p, &min, &max);
  double maxl = -HUGE;
  foreach_dimension()
    if (max.x - min.x > maxl)
      maxl = max.x - min.x;

We initialize the distance field on the coarse initial mesh and refine it adaptively until the threshold error (on distance) is reached.

  scalar d[];
  distance (d, p);
  while (adapt_wavelet ({d}, (double[]){eps*maxl}, maxlevel, 5).nf);

We also compute the volume fraction from the distance field. We first construct a vertex field interpolated from the centered field and then call the appropriate VOF functions.

  vertex scalar phi[];
  foreach_vertex()
    phi[] = (d[] + d[-1] + d[0,-1] + d[-1,-1] +
	     d[0,0,-1] + d[-1,0,-1] + d[0,-1,-1] + d[-1,-1,-1])/8.;
  fractions (phi, f);
}

Main function

We can change both the maximum level of refinement and the Froude number at runtime.

RV Tangaroa is 70 metres long. If we assume that it moves at 20 knots (twice its actual cruise speed), this gives a Froude number of approx 0.4.

int LEVEL = 9;
double FROUDE = 0.4, Uship = 1., Lship = 1. [1];

We need additional (fraction) fields for the ship geometry and for the (inflow) boundary condition.

scalar tangaroa[], f0[];

int main (int argc, char * argv[])
{
  maxruntime (&argc, argv);
  if (argc > 1)
    LEVEL = atoi(argv[1]);
  if (argc > 2)
    FROUDE = atof(argv[2]);
    
  init_grid (32);

  rho1 = 1. [0];    // water
  rho2 = rho1/815.; // air

The domain is five times larger than the ship. We change the origin so that the ship is not too close to the inflow.

  size (5.*Lship);
  origin (-L0/2.,-L0/3.,-L0/2.);

We need to tell the code that both tangaroa and f0 are volume fraction fields.

  for (scalar s in {tangaroa,f0})
    s.refine = s.prolongation = fraction_refine;

The acceleration of gravity is given by the Froude number i.e.

  G.z = - sq(Uship/FROUDE)/Lship;
  run();
}

Boundary conditions

The inflow condition fixes the velocity (unity) and the water level (using f0).

u.n[bottom] = dirichlet(Uship);
p[bottom]   = neumann(0.);
pf[bottom]  = neumann(0.);
f[bottom]   = f0[];

Outflow uses standard Neumann/Dirichlet conditions.

u.n[top]  = neumann(0.);
p[top]    = dirichlet(0.);
pf[top]   = dirichlet(0.);

Boundary conditions for the solid and fraction tracers.

tangaroa[back] = 0;
f[back] = 1;

Not sure whether this is really useful.

uf.n[left] = 0.;
uf.n[right] = 0.;

Initial conditions

We can optionally restart, otherwise we open the STL file and initialize the corresponding fraction. We also initialize the f0 field used for the inflow condition and set the initial water level and velocity field.

event init (t = 0) {
  if (!restore (file = "restart")) {
    FILE * fp = fopen ("tangaroa.stl", "r");
    fraction_from_stl (tangaroa, fp, 5e-4[0], LEVEL);
    fclose (fp);
    
    fraction (f0, - z);

    foreach() {
      f[] = f0[];
      u.y[] = Uship;
    }
  }

We setup the default display.

  display ("view (quat = {0.542,0.150,0.209,0.799},"
	   "      fov = 30, near = 0.01, far = 1000,"
	   "      tx = -0.0074, ty = -0.0138, tz = -0.3211);"
	   "draw_vof (c = 'f', color = 'z', min = -0.1, max = 0.1, "
	   "          linear = true, fc = {0.45,0.78,0.92});"
	   "draw_vof (c = 'tangaroa');", true);
}

Boundary conditions on the ship

We use a simple (but crude) imposition of u=0 in the solid.

event velocity (i++) {
  foreach()
    foreach_dimension()
      u.x[] = (1. - tangaroa[])*u.x[];
}

Animations

We generate animations of the ship surface (as represented by the solid fraction) and of the air-water interface, colored with the height.

Several classical features of ship wakes are recognisable: breaking bow wave, breaking stern divergent wave, turbulent boundary layer, Kelvin waves etc…

We also use the \lambda_2 criterion to display the turbulent vortical structures generated by the airflow. The air recirculation at the top of the steep primary Kelvin waves is particularly noticeable.

The computations above were done on the Irene supercomputer using 12 levels of refinement.

event movie (t += 0.01; t <= 10)
{
  view (fov = 5.86528,
	quat = {0.515965,0.140691,0.245247,0.808605},
	tx = -0.07438, ty = -0.0612925,
	width = 1024, height = 768);

  clear();
  draw_vof ("tangaroa");
  scalar Z[];
  Z[back] = dirichlet (z);
  foreach()
    Z[] = z;
  draw_vof ("f", color = "Z", min = -0.1, max = 0.1, linear = true);
  save ("movie.mp4");

  draw_vof ("tangaroa", fc = {0.5,0.5,0.5});
  draw_vof ("f", color = "Z", min = -0.1, max = 0.1, linear = true);
  scalar l2[];
  lambda2 (u, l2);
  isosurface ("l2", -100);
  save ("l2.mp4");
}

#if DUMP
event snapshot (i += 100)
{
  char name[80];
  sprintf (name, "dump-%d", i);
  scalar l2[];
  lambda2 (u, l2);
  dump (file = name);
}
#endif

Mesh adaptation

This computation is only feasible thanks to mesh adaptation, based both on volume fraction and velocity accuracy.

event adapt (i++) {
  double uemax = 0.1;
  adapt_wavelet ({f,tangaroa,u},
		 (double[]){0.01,0.01,uemax,uemax,uemax}, LEVEL, 5);
}

Running in parallel on Irene

Running with MPI-parallelism is a bit more complicated than usual since the distance() function is not parallelised yet. A reasonably simple workaround is to first generate a restart/dump file on the local machine, without MPI, using something like:

CFLAGS=-DDUMP=1 make tangaroa.tst

(and also adjust the maximum level), then kill the running code (using Ctrl-C) and do:

qcc -source -D_MPI=1 tangaroa.c
scp _tangaroa.c popinets@irene.ccc.cea.fr:/ccc/scratch/cont003/gen7325/popinets
scp tangaroa/dump-0 popinets@irene.ccc.cea.fr:/ccc/scratch/cont003/gen7325/popinets/restart

then on irene (to run on 480 cores for 10 hours, with 12 levels of refinement):

ssh popinets@irene.ccc.cea.fr
ccc_mpp -u popinets
sed 's/WALLTIME/36000/g' run.sh | ccc_msub -n 480

with the following run.sh script

#!/bin/bash
#MSUB -r tangaroa
#MSUB -T WALLTIME
#MSUB -@ popinet@basilisk.fr:begin,end
#MSUB -o basilisk_%I.out
#MSUB -e basilisk_%I.log
#MSUB -q skylake
#MSUB -A gen7760
#MSUB -m scratch
##MSUB -w

set -x
cd ${BRIDGE_MSUB_PWD}

mpicc -Wall -std=c99 -D_XOPEN_SOURCE=700 -O2 _tangaroa.c -o tangaroa \
      -L$HOME/gl -lglutils -lfb_tiny -lm
ccc_mprun -n ${BRIDGE_MSUB_NPROC} ./tangaroa -m WALLTIME 12 0.4 \
    2>> log-${BRIDGE_MSUB_NPROC} >> out-${BRIDGE_MSUB_NPROC}

Generating MP4 movies

Note that when running on Irene the ffmpeg MP4 encoder is not available and the logfile will contain the warning:

open_image(): cannot find 'ppm2mp4' or 'ffmpeg'/'avconv'
  falling back to raw PPM outputs

The MP4 files defined above will be renamed movie.mp4.ppm and l2.mp4.ppm. As the extension indicates, these (large) files are now raw (uncompressed) PPM images. To convert them to compressed MP4, you will need to copy them to a machine where ffmpeg (and Basilisk) are installed (i.e. your local machine) and do:

ppm2mp4 movie.mp4 < movie.mp4.ppm
ppm2mp4 l2.mp4 < l2.mp4.ppm