Large Eddy Simulation of isotropic Turbulence.

On this page we aim to validate the Large Eddy Simulation (LES) formulation as defined in the SGS.h header file. Fortunately, benchmark results from a Direct Numerical Simulation (DNS) are provided under /src/examples here. We only need to slightly alter the formulation presented there to turn the DNS into an LES. This page will highlight these changes and present the results.


First we include the Sub-Grid-Scale formulation. Furthermore, we include a function to help visualize the λ2 contours.

#include "grid/multigrid3D.h"
#include "navier-stokes/centered.h"
#include "SGS.h"
#include "lambda2.h"
#include "structurefunction.h"
#include "view.h"
#define BVIEW 1
#include "particles.h"
#define MU 0.01

face vector av[];
FILE * fp1;

Second, Facilitated by our filtered approach, we reduce the resolution from 1283 to 323.

int maxlevel = 5;
int lin = 0;
int main (int argc, char * argv[]) {
  if (argc > 1)
    maxlevel = atoi(argv[1]);
  L0 = 2.*π;
  X0 = Y0 = Z0 = -L0/2.;
    periodic (right);
  a = av;

Third, a relevant value for the SGS-tuning constant is set and we tell the simulation that the molecular viscosity is the MU value. Additionally, we open a file to we can display the results on this page.

event init (i = 0) {
  init_grid(1 << (maxlevel - 1));
  init_particles_in_cells(); //16x16x16 particles
  init_grid(1 << maxlevel);
  N = (1 << maxlevel);
  foreach() {
    u.x[] = cos(y) + sin(z);
    u.y[] = sin(x) + cos(z);
    u.z[] = cos(x) + sin(y);
  boundary ((scalar *){u});

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

Fourth, The dissipation should be calculated with the eddy viscosity.

event logfile (i+=10; t <= 300) {
  coord ubar;
  foreach_dimension() {
    stats s = statsf(u.x);
    ubar.x = s.sum/s.volume;
  double vd;
  double ke = 0., vdt = 0., vdm = 0.,vde = 0., vol = 0.;
  foreach(reduction(+:ke) reduction(+:vdt) reduction(+:vde) reduction(+:vdm) 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;
  vdm /= vol;
  vde /= vol;
  vdt /= vol;

  if (i == 0)
    fprintf (fp1,
	     "t dissipationt energy disse dissm perf.t perf.speed\n");
  fprintf (fp1, "%g %g %g %g %ld %g %g\n",
	   t, vdt, ke, vde, vdm, perf.t, perf.speed);

As a concsistency check, we evaluate the second-order structure function at fixed time intervals. Or is it really the structure funtion we aim for?

event structurefun(t+=50){
  char name2[100];
  FILE * fp2 = fopen(name2,"w");

For our visual statisfaction, we also generate an mp4 movie of particle advection.

event viewer (t = 20; t <= 100; t += 0.1) {
  view(φ = 0.5, θ = 0.25);


Movie before plots:


Evolution of the \lambda_2 contours.

Evolution of the λ2 contours.

We also look at some statistics. The resolved kinetic energy is shown below:

Evolution of kinetic energy

Evolution of kinetic energy

Compared to the spectral results, the LES ‘suffers’ from the early onset of turbulence, similar as the DNS run with Basilisk. The fluctuating kinetic energy after this initial transition (i.e. t>50) appears to be quite well captured. So atleast the implementation of the SGS formulation did not totally ruin the results. However, one may argue that energy is a rather simple quantity to get correct as it may very well be a result of the construct of the case definition. Therefore, we also look at the resolved dissipation of the simulations.

Evolution of dissipation

Evolution of dissipation

Seems OK. However, the initial transition shows a bit of a deveation, even between both Basilsik-based runs. We look at it in a bit more detail:

Evolution of the dissipation during the initial instability and its partitioning

Evolution of the dissipation during the initial instability and its partitioning

The second-order longitudional structure function was also computed and it is plotted below.

We do observe a range with a scaling characteristic. But is is not the 2/3 law I was expexting.

We do observe a range with a scaling characteristic. But is is not the 2/3 law I was expexting.

We conclude that the SGS-model is not able to properly capture the dissipation dynamics within the transition. This is not a show stopper perse, more that this serves as a cautionary note; for an accurate representation of certain flow statistics an LES is not always a suitable approach. The fact that the SGS model was able to preserve large-eddy quantities, like the kinetic energy, is motivating for its future usage. Furthermore, after the transition the turbulent strucures become larger and the spectrum more resolved. The SGS controbutions also vanish in that case. That is a nice feature.

Finally, note that the reduction from a 1283 grid to a 323 grid results is approximately 162=256 times less effort. Meaning that instead of using 512 cores for the original example, this simulation is reasonably fast on a single processor. So in that respect, well done LES-techniques! Also I checked the overhead of the LES formulation on my own system, using a single core. In DNS mode, I obtained an averaged performance figure of 4.6×105 points/sec during the first 20 time units, and with the LES formulation this number was reduced to 3.9×105 points/sec. Thus a reduction of 15%. These numbers were obtained with no movie output and less frequent calls to the diagnostic event.