lamb-dipole2nd.c

Lamb-Chaplygin dipole collision with a wall using second order accurate interpolation methods

Lamb-Chaplygin dipole collision with a wall using second order accurate interpolation methods

We can test the effect of applying second order accurate interpolation strategies using a simulation case nearly identical to an earlier version. You can view these methods here.

Numerical Set-up

Only minor things need to be changed compared to the version using the default settings. First we need to include the file that contains the definitions for the new attributes.

#include "navier-stokes/centered.h"
#include "additional_attributes.h"

scalar omg[],r[],s[],psi[];
int maxlevel = 10;
double xo= 1, yo= 0.00312;
double k=3.83170597;
double temp= 12;
double vis = 1./500.;
const face vector muc[] = {vis,vis}; 
u.t[right]=dirichlet(0);

int main()
{
  L0=10;
  X0=Y0=-L0/2;
  init_grid(1<<(8));
  run();
}

Initialization

We can set the restriction, coarsen, prolongation and refine attributes of the fields we solve for to the second order accurate formulations.

event init(t=0)
{
  # if 1
  u.x.restriction=quadratic_coarsening;
  u.x.coarsen=quadratic_coarsening;
  u.x.prolongation=refine_quadratic;
  u.x.refine=refine_quadratic;
  
  u.y.restriction=quadratic_coarsening;
  u.y.coarsen=quadratic_coarsening;
  u.y.prolongation=refine_quadratic;
  u.y.refine=refine_quadratic;
  
  p.restriction=quadratic_coarsening;
  p.coarsen=quadratic_coarsening;
  p.prolongation=refine_quadratic;
  p.refine=refine_quadratic;
  # endif
  mu=muc;
  refine(pow(pow((x-xo),2)+(pow((y-yo),2)),0.5) < 1.5 && level<11);
  foreach() {
    r[] = pow(pow((x-xo),2)+(pow((y-yo),2)),0.5);
    s[] = (y-yo)/r[];
    psi[] = ((r[]>1)*((1/r[]))*s[]) + ((r[]<1)*((-2*j1(k*r[])*s[]/(k*j0(k)))+(s[]*r[])));
  }
  boundary({psi});
  foreach() {
    u.x[] = ((psi[0,1]-psi[0,-1])/(2*Delta));
    u.y[] = -(psi[1,0]-psi[-1,0])/(2*Delta);
  }
  boundary(all);
}

Adaptation

We do not wish to use the second order accurate formulations for the “adapt_wavelet” function as that this would result in a different grid structure. Therefore, we define two dummy scalar fields uu[] and vv[] that use the default attributes such that the wavelet error estimation is not altered. Note that potentially the second order accurate attributes can be usefull for error estimation as well.

event adapt(i++)
{
  scalar uu[],vv[];
  foreach()
    {
      uu[]=u.x[];
      vv[]=u.y[];
    }
  boundary({uu,vv});
  adapt_wavelet((scalar *){uu,vv},(double []){0.01,0.01},maxlevel); 
}
double aa=0.1;
event gfsview(t+=aa;t<=temp)
{
  double pz=0;
  double en=0;
  int n=0;
  foreach() 
    omg[]=(u.x[0,1]-u.x[0,-1] - (u.y[1,0]-u.y[-1,0]))/(2*Delta);
  
  foreach(reduction(+:pz) reduction(+:en) reduction(+:n)) {
    pz+=pow(omg[],2)*Delta*Delta;
    n++;
    foreach_dimension()
      en+=sq(u.x[])*Delta*Delta;
  }
  fprintf(ferr,"%d\t%g\t%g\t%g\t%d\n",i,t,pz,en,n);
  
  #if 0
  static FILE * fp =
    popen ("gfsview-batch2D lamb.gfv| ppm2gif > lamb4000.gif", "w");
  output_gfs (fp);
  fprintf (fp, "Save stdout { format = PPM width = 640 height = 480}\n");

  static FILE * fp2 =
    popen ("gfsview-batch2D lambinzoom.gfv | ppm2gif > lamb4000inz.gif", "w");
  output_gfs (fp2);
  fprintf (fp2, "Save stdout { format = PPM width = 1000 height = 500}\n");
  #endif
    aa=0.1/(pz/100);
}

Results and Performance

It appears that the results become more accurate using the new second order accurate formulated attributes. We can for example compare the total resolved enstrophy (\Omega).

Evolution of enstrophy and number of grid points

The performance appears to be affected as this is what I retrieved from compiling with the “-DTRACE=2” tag

Using the default attributes:

# Quadtree, 3193 steps, 154.649 CPU, 154.7 real, 3.99e+05 points.step/s, 20 var
# 4 procs, MPI: min 38 (25%) avg 41 (27%) max 45 (29%)
   calls    total     self   % total   function
   83053    39.79    28.41     18.4%   boundary():/home/antoon/basilisk/src/grid/cartesian-common.h:274
    6388    48.14    27.48     17.8%   project():/home/antoon/basilisk/src/poisson.h:418
    3194    35.90    22.76     14.7%   viscous_term():/home/antoon/basilisk/src/navier-stokes/centered.h:295
 4062230    16.04    16.04     10.4%   mpi_boundary_level():/home/antoon/basilisk/src/grid/tree-mpi.h:421
 4033794    12.94    12.94      8.4%   mpi_boundary_restriction():/home/antoon/basilisk/src/grid/tree-mpi.h:428
    3202    19.77    12.41      8.0%   balance():/home/antoon/basilisk/src/grid/balance.h:390
    6416     7.87     7.87      5.1%   mpi_boundary_update_buffers():/home/antoon/basilisk/src/grid/tree-mpi.h:929
   31662     8.82     7.56      4.9%   mpi_boundary_coarsen():/home/antoon/basilisk/src/grid/tree-mpi.h:1107
    3194    38.90     6.94      4.5%   advection_term():/home/antoon/basilisk/src/navier-stokes/centered.h:252
   77728     3.95     3.95      2.6%   mpi_all_reduce0():/home/antoon/basilisk/src/common.h:580
    3194    43.54     2.96      1.9%   adapt_wavelet():/home/antoon/basilisk/src/grid/tree-common.h:293
    3194    28.46     1.66      1.1%   projection():/home/antoon/basilisk/src/navier-stokes/centered.h:352
    3202     1.72     1.56      1.0%   z_indexing():/home/antoon/basilisk/src/grid/tree-mpi.h:1389
    .....
    And some more lines with less than 1% total time

And With the second order accurate attributes:

# Quadtree, 3266 steps, 183.536 CPU, 183.6 real, 3.61e+05 points.step/s, 20 var
# 4 procs, MPI: min 41 (22%) avg 44 (24%) max 46 (25%)
   calls    total     self   % total   function
    6534    64.33    39.21     21.4%   project():/home/antoon/basilisk/src/poisson.h:418
   87943    46.46    34.36     18.7%   boundary():/home/antoon/basilisk/src/grid/cartesian-common.h:274
    3267    45.19    29.45     16.0%   viscous_term():/home/antoon/basilisk/src/navier-stokes/centered.h:295
 4681147    17.31    17.31      9.4%   mpi_boundary_level():/home/antoon/basilisk/src/grid/tree-mpi.h:421
 4652085    15.27    15.27      8.3%   mpi_boundary_restriction():/home/antoon/basilisk/src/grid/tree-mpi.h:428
    3276    19.98    12.83      7.0%   balance():/home/antoon/basilisk/src/grid/balance.h:390
   32372     8.95     7.96      4.3%   mpi_boundary_coarsen():/home/antoon/basilisk/src/grid/tree-mpi.h:1107
    6563     7.83     7.83      4.3%   mpi_boundary_update_buffers():/home/antoon/basilisk/src/grid/tree-mpi.h:929
    3267    40.08     7.23      3.9%   advection_term():/home/antoon/basilisk/src/navier-stokes/centered.h:252
   81029     3.50     3.50      1.9%   mpi_all_reduce0():/home/antoon/basilisk/src/common.h:580
    3267    44.86     3.05      1.7%   adapt_wavelet():/home/antoon/basilisk/src/grid/tree-common.h:293
    3267    45.33     1.76      1.0%   projection():/home/antoon/basilisk/src/navier-stokes/centered.h:352
    3276     1.77     1.59      0.9%   z_indexing():/home/antoon/basilisk/src/grid/tree-mpi.h:1389
    .....
    And some more lines with a low impact on performance

These results were obtained using my laptop.