src/test/rt.c

    #include "navier-stokes/centered.h"
    #include "two-phase.h"
    #include "navier-stokes/conserving.h"
    #if REDUCED
    # include "reduced.h"
    #endif
    
    #define LEVEL 8
    
    uf.n[left]   = 0;
    p[left] = neumann(0.);
    uf.n[right]  = 0;
    p[right] = neumann(0.);
    
    int main() {
      size (4);
      origin (-2, -2);
      init_grid (1 << LEVEL);
      mu1 = mu2 = 0.00313;
      rho1 = 1.225, rho2 = 0.1694;
    #if REDUCED
      G.y = -9.81;
    #else
      const face vector g[] = {0,-9.81};
      a = g;
    #endif
      DT = 5e-3;
      TOLERANCE = 1e-6;
      run();
    }
    
    event init (t = 0) {
    #if TREE  
      mask (x < -0.5 ? left : x > 0.5 ? right : none);
    #endif
      fraction (f, 0.05*cos (2.*pi*x) + y);
    }
    
    event logfile (i++) {
      stats s = statsf (f);
      printf ("%g %d %g %g %g %g %d %d %d\n", 
    	  t, i, dt, s.sum - 2., s.min, s.max - 1., mgp.i, mgpf.i, mgu.i);
      assert (s.min >= -1e-10 && s.max <= 1. + 1e-10);
      assert (fabs (s.sum - 2.) < 9e-7);
    }
    
    event interface (t = {0,0.2,0.4,0.8}) {
      output_facets (f, stderr);
      FILE * fp = fopen ("levels", "w");
      scalar l[];
      foreach()
        l[] = level;
      output_field ({l}, fp);
      fclose (fp);
    }
    
    #if 0
    event gfsview (i += 10) {
      static FILE * fp = popen("gfsview2D -s rt.gfv", "w");
      output_gfs (fp);
    }
    #endif
    
    #if TREE
    event adapt (i++) {
      coord momb = {0};
      foreach()
        foreach_dimension()
          momb.x += dv()*rho(f[])*u.x[];
      
      adapt_wavelet ({f}, (double[]){5e-3}, LEVEL);
    
      coord moma = {0};
      foreach()
        foreach_dimension()
          moma.x += dv()*rho(f[])*u.x[];

    We check that each component of the momentum is conserved.

      foreach_dimension()
        assert (fabs(momb.x - moma.x) < 1e-10);
    }
    #endif