sandbox/Antoonvh/isotropicLES.c

    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.

    Setup

    First we include the Sub-Grid-Scale formulation. Furthermore, we include a function to help visualize the \lambda_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 128^3 to 32^3.

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

    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);
  Csmag=0.2;
  molvis=MU;
  fp1=fopen("LES.dat","w");
  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;
  }
  foreach_face()
    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.*Delta);
      vde+=(Evis[]-MU)*vd;
      vdm+=MU*vd;
      vdt+=Evis[]*vd;
    }
  }
  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 %g %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];
  lin++;
  sprintf(name2,"structure%d.dat",lin);
  FILE * fp2 = fopen(name2,"w");
  long2structure(fp2,u,1000,L0/2.,100);
  fflush(fp2);
  fclose(fp2);
}

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

    event viewer (t = 20; t <= 100; t += 0.1) {
  view(phi = 0.5, theta = 0.25);
  clear();
  scatter(loc);
  save("LESiso.mp4");
}

    Results

    Movie before plots:

    Particles!

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

    set xlabel 'Time'
set ylabel 'Kinetic energy'
set logscale  y
plot 'LES.dat' u 1:3 w l lw 6 t 'Present LES' ,\
     'isotropicLES.DNSbasilisk.dat' u 1:3 w l lw 4 t 'DNS Basilisk' ,\
     'isotropicLES.hit3d.dat' u 1:($3*3./2.) w l lw 4 t 'Spectral DNS hit3d'
    Evolution of kinetic energy (script)

    Evolution of kinetic energy (script)

    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.

    set ylabel 'Dissipation'
plot  'LES.dat' u 1:2 w l lw 3 t 'Present LES' ,\
 'isotropicLES.DNSbasilisk.dat' u 1:2 w l lw 3 t 'DNS Baslisk',\
 'isotropicLES.hit3d.dat' u 1:2 w l lw 3 t 'Spectral DNS hit3d'
    Evolution of dissipation (script)

    Evolution of dissipation (script)

    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:

    set xr [10:25]
set ylabel 'Dissipation'
unset logscale y
plot  'LES.dat' u 1:2 w l lw 3 t 'LES total' ,\
'isotropicLES.DNSbasilisk.dat' u 1:2 w l lw 3 t 'DNS Baslisk' ,\
  'LES.dat' u 1:6 w l lw 3 t 'LES molecular' ,\
  'LES.dat' u 1:4 w l lw 3 t 'LES SGS'
    Evolution of the dissipation during the initial instability and its partitioning (script)

    Evolution of the dissipation during the initial instability and its partitioning (script)

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

    reset
set  xr [0.1:5]
set  yr [0.005:0.5]
set ylabel 'dv?'
set xlabel 'Distance'
set logscale y
set logscale x  
set key bottom right
plot  'structure5.dat' u 1:2 t'data t = 200',\
'structure6.dat' u 1:2 t'data t = 250',\
'structure7.dat' u 1:2 t'data t = 300',\
  (x**(1.33))/5 w l lw 3 t 'x^{4/3} scaling'
    We do observe a range with a scaling characteristic. But is is not the 2/3 law I was expexting. (script)

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

    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 128^3 grid to a 32^3 grid results is approximately 16^2=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 \times 10^5 points/sec during the first 20 time units, and with the LES formulation this number was reduced to 3.9\times 10^5 points/sec. Thus a reduction of 15%. These numbers were obtained with no movie output and less frequent calls to the diagnostic event.