src/test/stokes-ns.c

    Breaking Stokes wave

    We solve the two-phase Navier–Stokes equations using a momentum-conserving transport of each phase. Gravity is taken into account using the “reduced gravity approach”.

    #include "navier-stokes/centered.h"
    #include "two-phase.h"
    #include "navier-stokes/conserving.h"
    #include "reduced.h"

    The primary parameters are the wave steepness ak and the Reynolds number.

    double ak = 0.35;
    double RE = 40000.;
    int LEVEL = 9;

    The error on the components of the velocity field used for adaptive refinement.

    double uemax = 0.005;

    The density and viscosity ratios are those of air and water.

    #define RATIO (1./850.)
    #define MURATIO (17.4e-6/8.9e-4)

    The wave number, fluid depth and acceleration of gravity are set to these values. T0 is the wave period.

    #define k_  (2.*pi)
    #define h_   0.5
    #define g_   1.
    #define T0  (k_/sqrt(g_*k_))

    The program takes optional arguments which are the level of refinement, steepness and Reynolds numbers.

    int main (int argc, char * argv[])
    {
      if (argc > 1)
        LEVEL = atoi (argv[1]);
      if (argc > 2)
        ak = atof(argv[2]);

    The domain is a cubic box centered on the origin and of length L0=1, periodic in the x-direction.

      origin (-L0/2, -L0/2);
      periodic (right);

    Here we set the densities and viscosities corresponding to the parameters above.

      rho1 = 1.;
      rho2 = RATIO;
      mu1 = 1.0/RE; //using wavelength as length scale
      mu2 = 1.0/RE*MURATIO;
      G.y = -g_;

    When we use adaptive refinement, we start with a coarse mesh which will be refined as required when initialising the wave.

    #if TREE  
      N = 32;
    #else
      N = 1 << LEVEL;
    #endif
      DT = 1e-2;
      run();
    }

    Initial conditions

    We either restart (if a “restart” file exists), or initialise the wave using the third-order Stokes wave solution.

    #include "stokes.h"
    
    event init (i = 0)
    {
      if (!restore ("restart")) {

    We need to make sure that fields are properly initialised before refinement below, otherwise a -catch exception will be triggered when debugging.

        event ("properties");
        
        do {
          fraction (f, wave(x,y));
          foreach()
    	foreach_dimension()
    	  u.x[] = u_x(x,y) * f[];
          boundary ((scalar *){u});
        }

    On trees, we repeat this initialisation until mesh adaptation does not refine the mesh anymore.

    #if TREE  
        while (adapt_wavelet ({f,u},
    			  (double[]){0.01,uemax,uemax,uemax}, LEVEL, 5).nf);
    #else
        while (0);
    #endif
      }
    }
    
    event profiles (t += T0/4.; t <= 2.5*T0) {
      output_facets (f, stderr);
      fprintf (stderr, "\n");
    }
    
    event logfile (i++)
    {
      double ke = 0., gpe = 0.;
      foreach (reduction(+:ke) reduction(+:gpe)) {
        double norm2 = 0.;
        foreach_dimension()
          norm2 += sq(u.x[]);
        ke += norm2*f[]*dv();
        gpe += y*f[]*dv();
      }
      printf ("%g %g %g\n", t/(k_/sqrt(g_*k_)), rho1*ke/2., rho1*g_*gpe + 0.125);
    }

    Mesh adaptation

    On trees, we adapt the mesh according to the error on volume fraction and velocity.

    #if TREE
    event adapt (i++) {
      adapt_wavelet ({f,u}, (double[]){0.01,uemax,uemax,uemax}, LEVEL, 5);
    }
    #endif