/**
# Forced isotropic turbulence in a triply-periodic box

We compute the evolution of forced isotropic turbulence (see [Rosales
& Meneveau, 2005](/src/references.bib#rosales2005)) and compare the
solution to that of the [hit3d](http://code.google.com/p/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.](isotropic/movie.mp4)
*/

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");
}

/**
We can optionally try adaptivity. */

#if TREE
event adapt (i++) {
  double uemax = 0.2*normf(u.x).avg;
  adapt_wavelet ((scalar *){u}, (double[]){uemax,uemax,uemax}, maxlevel);
}
#endif

/**
## Running with MPI on occigen

On the local machine

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

On occigen (using 512 cores)

~~~bash
module purge
module load openmpi
module load intel
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](/src/gl/INSTALL#standalone-installation) in `$HOME/gl`.

With the `run.sh` script

~~~bash
#!/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

~~~bash
%local qcc -source -D_MPI=1 isotropic.c
%local scp _isotropic.c mesu.dsi.upmc.fr:
~~~

On mesu (using 512 cores)

~~~bash
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

~~~bash
#!/bin/bash 
#PBS -l select=22:ncpus=24:mpiprocs=24
#PBS -l walltime=WALLTIME
#PBS -N isotropic
#PBS -j oe  
# load modules 
module load mpt
# 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.

~~~gnuplot Evolution of kinetic energy
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'
~~~

~~~gnuplot Evolution of microscale Reynolds number
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'
~~~

~~~gnuplot Evolution of dissipation
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 128^3^
problem on 512 cores). Note that these were obtained when switching
off movie outputs.

~~~gnuplot Computational speed in points.timesteps/sec/core
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](http://www-hpc.cea.fr/en/complexe/tgcc-Irene.htm) is the
supercomputer at CEA.

On the local machine

~~~bash
%local qcc -source -D_MPI=1 isotropic.c
%local scp _isotropic.c irene.ccc.cea.fr:
~~~

On irene

~~~bash
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

~~~bash
#!/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}
~~~

~~~gnuplot Weak scaling on Irene (16^3^ per core)
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
~~~

~~~gnuplot Weak scaling on Irene (32^3^ per core)
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
~~~

## See also

* [Same example with Gerris](http://gerris.dalembert.upmc.fr/gerris/examples/examples/forcedturbulence.html)
*/