Forced isotropic turbulence in a triply-periodic box

We compute the evolution of forced isotropic turbulence (see Rosales & Meneveau, 2005) and compare the solution to that of the hit3d pseudo-spectral code. The initial condition is an unstable solution to the incompressible Euler equations. Numerical noise in the solution eventually leads to the destabilisation of the base solution into a fully turbulent flow where turbulent dissipation balances the linear input of energy.

For the moment we can only use the “multigrid” regular Cartesian grid implementation, as octrees do not support periodic boundaries yet.

#include "grid/multigrid3D.h"
#include "navier-stokes/centered.h"

#define MU 0.01

We need to store the variable forcing term.

face vector av[];

The code takes the level of refinement as optional command-line argument.

int maxlevel = 7;

int main (int argc, char * argv[]) {
  if (argc > 1)
    maxlevel = atoi(argv[1]);

The domain is (2π)3 and triply-periodic.

  L0 = 2.*π;
    periodic (right);

The viscosity is constant. The acceleration is defined below. The level of refinement is maxlevel.

  const face vector muc[] = {MU,MU,MU};
  μ = muc;
  a = av;
  N = 1 << maxlevel;

Initial conditions

The initial condition is “ABC” flow. This is a laminar base flow that is easy to implement in both Basilisk and a spectral code.

For the moment dump() and restore() do not work with MPI.

event init (i = 0) {
  if (!restore (file = "snapshot"))
    foreach() {
      u.x[] = cos(y) + sin(z);
      u.y[] = sin(x) + cos(z);
      u.z[] = cos(x) + sin(y);
  boundary ((scalar *){u});

Linear forcing term

We compute the average velocity and add the corresponding linear forcing term.

event acceleration (i++) {
  coord ubar;
  foreach_dimension() {
    stats s = statsf(u.x);
    ubar.x = s.sum/s.volume;
    av.x[] += 0.1*((u.x[] + u.x[-1])/2. - ubar.x);


We log the evolution of viscous dissipation, kinetic energy, microscale Reynolds number and some performance statistics.

event logfile (i++; t <= 300) {
  coord ubar;
  foreach_dimension() {
    stats s = statsf(u.x);
    ubar.x = s.sum/s.volume;
  double ke = 0., vd = 0., vol = 0.;
  foreach(reduction(+:ke) reduction(+:vd) reduction(+:vol)) {
    vol += dv();
    foreach_dimension() {
      // mean fluctuating kinetic energy
      ke += dv()*sq(u.x[] - ubar.x);
      // viscous dissipation
      vd += dv()*(sq(u.x[1] - u.x[-1]) +
		  sq(u.x[0,1] - u.x[0,-1]) +
		  sq(u.x[0,0,1] - u.x[0,0,-1]))/sq(2.*Δ);
  ke /= 2.*vol;
  vd *= MU/vol;

  if (i == 0)
    fprintf (ferr,
	     "t dissipation energy Reynolds grid->tn perf.t perf.speed\n");
  fprintf (ferr, "%g %g %g %g %ld %g %g\n",
	   t, vd, ke, 2./3.*ke/MU*sqrt(15.*MU/vd), grid->tn, perf.t, perf.speed);

event gfsview (t += 0.5; t <= 150) {
  static FILE * fp =
    popen ("gfsview-batch3D isotropic.gfv | ppm2mpeg > isotropic.mpg", "w");
  output_gfs (fp, translate = true);

event snapshot (i += 100)
  dump (file = "snapshot");

Running with MPI on occigen

On the local machine

%local qcc -source -D_MPI=1 isotropic.c
%local scp _isotropic.c

On occigen (using 512 cores)


with the script

#SBATCH -J basilisk
#SBATCH --nodes=32
#SBATCH --ntasks-per-node=16
#SBATCH --threads-per-core=1
#SBATCH --time=5:00
#SBATCH --output basilisk.output
#SBATCH --exclusive


module purge
module load bullxmpi
module load intel

mpicc -Wall -std=c99 -O2 -D_GNU_SOURCE=1 _isotropic.c -o isotropic -lm
srun --mpi=pmi2 -K1 --resv-ports -n $SLURM_NTASKS ./isotropic $LEVEL \


The two codes agree at early time, or until the solution transitions to a turbulent state. This happens earlier in Basilisk as the symmetry of the base state is not preserved with the same accuracy as in the pseudo-spectral code (the main reason being the tolerance on non-divergence of the incompressible velocity field). Note however that the statistics produced by the two codes agree well after transition to turbulence.

Evolution of kinetic energy

Evolution of kinetic energy

Evolution of microscale Reynolds number

Evolution of microscale Reynolds number

Evolution of dissipation

Evolution of dissipation

The computational speed is respectable (for a relatively small 1283 problem on 512 cores).

Computational speed in points.timesteps/sec/core

Computational speed in points.timesteps/sec/core

See also