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

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

We use the \lambda_2 criterion and Basilisk View for visualisation of vortices.

#include "lambda2.h"
#include "view.h"

We monitor performance statistics and control the maximum runtime.

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

const double 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 (as well as an optional maximum runtime).

int maxlevel = 7;

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

The domain is (2\pi)^3 and triply-periodic.

  L0 = 2.*pi [1];
foreach_dimension()
periodic (right);

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

  const face vector muc[] = {MU,MU,MU};
mu = muc;
a = av;
N = 1 << maxlevel;
run();
}

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

event init (i = 0) {
double u0 = 1., k = 1.;
if (!restore (file = "restart"))
foreach() {
u.x[] = u0*(cos(k*y) + sin(k*z));
u.y[] = u0*(sin(k*x) + cos(k*z));
u.z[] = u0*(cos(k*x) + sin(k*y));
}
}

## 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;
}
double a = 0.1;
foreach_face()
av.x[] += a*((u.x[] + u.x[-1])/2. - ubar.x);
}

## Outputs

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

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.*Delta);
}
}
ke /= 2.*vol;
vd *= MU/vol;

if (i == 0)
fprintf (stderr, "t dissipation energy Reynolds\n");
fprintf (stderr, "%g %g %g %g\n",
t, vd, ke, 2./3.*ke/MU*sqrt(15.*MU/vd));
}

We generate a movie of the vortices.

Animation of the \lambda_2 isosurface (a way to characterise vortices) and cross-sections of velocity and vorticity.

event movie (t += 0.25; t <= 150)
{
view (fov = 44, camera = "iso", ty = .2,
width = 600, height = 600, bg = {1,1,1}, samples = 4);
clear();
squares ("u.y", linear = true);
squares ("u.x", linear = true, n = {1,0,0});
scalar omega[];
vorticity (u, omega);
squares ("omega", linear = true, n = {0,1,0});
scalar l2[];
lambda2 (u, l2);
isosurface ("l2", -1);
save ("movie.mp4");
}

#if TREE
double uemax = 0.2*normf(u.x).avg;
}
#endif

## Running with MPI on occigen

On the local machine

%local qcc -source -D_MPI=1 isotropic.c
%local scp _isotropic.c occigen.cines.fr:

On occigen (using 512 cores)

module purge
mpicc -Wall -O2 -std=c99 -D_XOPEN_SOURCE=700 _isotropic.c -o isotropic \
-I$HOME -L$HOME/gl -lglutils -lfb_tiny -lm
sed 's/WALLTIME/10:00/g' run.sh | sbatch

Note that this assumes that the Basilisk gl libraries have been installed in $HOME/gl. With the run.sh script #!/bin/bash #SBATCH -J basilisk #SBATCH --nodes=32 #SBATCH --constraint=HSW24 #SBATCH --ntasks-per-node=16 #SBATCH --threads-per-core=1 #SBATCH --time=WALLTIME #SBATCH --output basilisk.output #SBATCH --exclusive LEVEL=7 module purge module load openmpi module load intel srun --mpi=pmi2 -K1 --resv-ports -n$SLURM_NTASKS \
./isotropic -m WALLTIME $LEVEL \ 2> log-$LEVEL-$SLURM_NTASKS > out-$LEVEL-$SLURM_NTASKS ## Running with MPI on mesu On the local machine %local qcc -source -D_MPI=1 isotropic.c %local scp _isotropic.c mesu.dsi.upmc.fr: On mesu (using 512 cores) module load mpt mpicc -Wall -O2 -std=c99 -D_XOPEN_SOURCE=700 _isotropic.c -o isotropic \ -L$HOME/gl -lglutils -lfb_tiny -lm
sed 's/WALLTIME/10:00/g' run.sh | qsub

with the run.sh script

#!/bin/bash
#PBS -l select=22:ncpus=24:mpiprocs=24
#PBS -l walltime=WALLTIME
#PBS -N isotropic
#PBS -j oe
# change to the directory where program job_script_file is located
cd $PBS_O_WORKDIR # mpirun -np 64 !!!! does not work !!!! NP=512 mpiexec_mpt -n$NP ./isotropic -m WALLTIME 2>> log.$NP >> out.$NP

## Results

The two codes agree at early time, or until the solution transitions to a turbulent state. The statistics produced by the two codes agree well after transition to turbulence.

set xlabel 'Time'
set ylabel 'Kinetic energy'
set logscale  y
plot 'isotropic.occigen' u 1:3 w l t 'Basilisk (occigen)', \
'isotropic.mesu' u 1:3 w l t 'Basilisk (mesu)', \
'isotropic.hit3d' u 1:($3*3./2.) w l t 'Spectral' set ylabel 'Microscale Reynolds number' plot 'isotropic.occigen' u 1:4 w l t 'Basilisk (occigen)', \ 'isotropic.mesu' u 1:4 w l t 'Basilisk (mesu)', \ 'isotropic.hit3d' u 1:4 w l t 'Spectral' set ylabel 'Dissipation function' plot 'isotropic.occigen' u 1:2 w l t 'Basilisk (occigen)', \ 'isotropic.mesu' u 1:2 w l t 'Basilisk (mesu)', \ 'isotropic.hit3d' u 1:2 w l t 'Spectral' The computational speed is respectable (for a relatively small 1283 problem on 512 cores). Note that these were obtained when switching off movie outputs. set ylabel 'Speed' unset logscale plot 'isotropic.occigen' u 1:($7/512) w l t 'occigen', \
'isotropic.mesu' u 1:($7/512) w l t 'mesu' ## Scalability on irene Irene is the supercomputer at CEA. On the local machine %local qcc -source -D_MPI=1 isotropic.c %local scp _isotropic.c irene.ccc.cea.fr: On irene mpicc -Wall -O2 -std=c99 -D_XOPEN_SOURCE=700 -xCORE-AVX512 \ _isotropic.c -o isotropic \ -L$HOME/gl -lglutils -lfb_tiny -lm
sed -e 's/WALLTIME/600/g' -e 's/LEVEL/7/g' run.sh | ccc_msub -n 512

with the run.sh script

#!/bin/bash
#MSUB -r isotropic
#MSUB -T WALLTIME
#MSUB -o basilisk_%I.out
#MSUB -e basilisk_%I.log
#MSUB -q skylake
#MSUB -A gen7760
#MSUB -m scratch

set -x
cd ${BRIDGE_MSUB_PWD} ccc_mprun -n${BRIDGE_MSUB_NPROC} ./isotropic -m WALLTIME LEVEL \
2> log-${BRIDGE_MSUB_NPROC} > out-${BRIDGE_MSUB_NPROC}
set xlabel '# of cores'
set ylabel 'Speed (points x timestep/sec/core)'
set logscale x
plot '-' w lp t ''
8     366497
64    273009
512   196904
4096  167429
e
set xlabel '# of cores'
set ylabel 'Speed (points x timestep/sec/core)'
set logscale x
plot '-' w lp t ''
8     281250
64    189218
512   165312
4096  156683
e