src/test/ponds-ml.c

    Stress test for wetting and drying

    Bumps uniformly covered with water create ponds and streams (i.e. a lot of wetting/drying going on).

    set term @PNG enhanced size 1024,512 font ",8"
    set output 'eta.png'
    
    ! awk '{ if ($1 == "file:") file = $2; else print $0 > file; }' < out
    
    unset key
    set hidden3d
    unset xtics
    unset ytics
    unset border
    set palette defined ( 0 0 0 0.5647, 0.125 0 0.05882 1, 0.25 0 0.5647 1, \
                              0.375 0.05882 1 0.9333, 0.5 0.5647 1 0.4392, \
                          0.625 1 0.9333 0, 0.75 1 0.4392 0, \
                          0.875 0.9333 0 0, 1 0.498 0 0 )
    
    dry=1e-3
    set view 28,55
    
    set multiplot layout 2,2 scale 1,1.3
    splot './eta-1' u 1:2:4 every 3:3 w l, \
          './eta-1' u 1:2:($3>dry?$3+$4:1e1000):(sqrt($5**2+$6**2)) w pm3d
    splot './eta-2' u 1:2:4 every 3:3 w l, \
          './eta-2' u 1:2:($3>dry?$3+$4:1e1000):(sqrt($5**2+$6**2)) w pm3d
    splot './eta-3' u 1:2:4 every 3:3 w l, \
          './eta-3' u 1:2:($3>dry?$3+$4:1e1000):(sqrt($5**2+$6**2)) w pm3d
    splot './eta-8' u 1:2:4 every 3:3 w l, \
          './eta-8' u 1:2:($3>dry?$3+$4:1e1000):(sqrt($5**2+$6**2)) w pm3d
    unset multiplot
    Free surface coloured with the norm of the velocity. (script)

    Free surface coloured with the norm of the velocity. (script)

    set output 'level.png'
    dry = 0.
    set multiplot layout 2,2 scale 1,1.3
    splot './level-1' u 1:2:4 every 3:3 w l, \
          './level-1' u 1:2:($3>dry?$3+$4:1e1000):5 w pm3d
    splot './level-2' u 1:2:4 every 3:3 w l, \
          './level-2' u 1:2:($3>dry?$3+$4:1e1000):5 w pm3d
    splot './level-3' u 1:2:4 every 3:3 w l, \
          './level-3' u 1:2:($3>dry?$3+$4:1e1000):5 w pm3d
    splot './level-8' u 1:2:4 every 3:3 w l, \
          './level-8' u 1:2:($3>dry?$3+$4:1e1000):5 w pm3d
    unset multiplot
    Evolution of the level of refinement. (script)

    Evolution of the level of refinement. (script)

    reset
    set term @PNG enhanced size 1024,512 font ",8"
    set output 'vectors.png'
    
    set xrange [0:1000]
    set yrange [0:1000]
    set size ratio -1
    unset key
    dry=1e-3
    plot './eta-8' u 1:2:($3>dry?$5*30.:0):($3>dry?$6*30.:0) every 2:2 w vec lc 0
    
    ! rm -f eta-? level-?
    Velocity field at the end of the simulation. (script)

    Velocity field at the end of the simulation. (script)

    // use KDT database rather than analytical function
    #define KDT 1
    
    #if ML
    # include "layered/hydro.h"
    # if IMPLICIT
    #   include "layered/implicit.h"
    # else
    #   include "layered/nh.h"
    # endif
    #elif IMPLICIT
    # include "saint-venant-implicit.h"
    #else
    # include "saint-venant.h"
    #endif
    
    #if KDT
    # include "terrain.h"
    #endif
    
    #define LEVEL 7
    
    int main()
    {
      init_grid (1 << LEVEL);
      size (1000.);
      G = 9.81;
    #if IMPLICIT
      DT = 10.;
    #endif
      run();
    }
    
    #define zb(x,y) ((cos(pi*x/L0)*cos(pi*y/L0) + \
    		  cos(3.*pi*x/L0)*cos(3.*pi*y/L0)) - 2.*x/1000.)
    
    #if !KDT
    void refine_zb (Point point, scalar zb)
    {
      foreach_child()
        zb[] = zb(x,y);
    }
    #endif
    
    event init (i = 0)
    {
    #if !KDT
      zb.refine = refine_zb; // updates terrain
      foreach() {
        zb[] = zb(x,y);
        h[] = 0.1;
      }
    #else
      FILE * fp = popen ("xyz2kdt ponds", "w");
      for (double x = 0.; x <= 1000; x += 1.)
        for (double y = 0.; y <= 1000.; y += 1.)
          fprintf (fp, "%g %g %g\n", x, y, zb(x,y));
      pclose (fp);
      terrain (zb, "ponds", NULL);
      foreach()
        h[] = 0.1;
    #endif
    }
    
    event friction (i++) {
    #if IMPLICIT && !ML
      // quadratic bottom friction, coefficient 1e-4 (dimensionless)
      foreach() {
        double a = h[] < dry ? HUGE : 1. + 1e-4*dt*norm(q)/sq(h[]);
        foreach_dimension()
          q.x[] /= a;
      }
      boundary ((scalar *){q});
    #else
      // quadratic bottom friction, coefficient 1e-4 (dimensionless)
      foreach() {
        double a = h[] < dry ? HUGE : 1. + 1e-4*dt*norm(u)/h[];
        foreach_dimension()
          u.x[] /= a;
      }
      boundary ((scalar *){u});
    #endif
    }
    
    event logfile (i += 10) {
      stats s = statsf (h);
    #if IMPLICIT && !ML
      scalar u[];
      foreach()
        u[] = h[] > dry ? q.x[]/h[] : 0.;
      norm n = normf (u);
    #else
      norm n = normf (u.x);
    #endif
      if (i == 0)
        fprintf (stderr, "t i h.min h.max h.sum u.x.rms u.x.max dt\n");
      fprintf (stderr, "%g %d %.4g %g %g %g %g %g\n", 
    	   t, i, s.min, s.max, s.sum, n.rms, n.max, dt);
    #if !IMPLICIT
      assert (s.min > 0.);
    #endif
    }
    
    event outputfile (t <= 1200.; t += 1200./8) {
    #if ML || !IMPLICIT 
      static int nf = 0;
      printf ("file: eta-%d\n", nf);
      output_field ({h, zb, u}, stdout, N, linear = true);
    
      scalar l[];
      foreach()
        l[] = level;
      printf ("file: level-%d\n", nf);
      output_field ({h, zb, l}, stdout, N);
    
      nf++;
    #endif
    }
    
    // int event (t += 100)
    //  output_matrix (h, stdout, N, true);
    
    event adapt (i++) {
      // we do this so that wavelets use the default bilinear
      // interpolation this is less noisy than the linear + gradient
      // limiters used in Saint-Venant not sure whether this is better
      // though.
      scalar h1[];
      foreach()
        h1[] = h[];
      boundary ({h1});
      adapt_wavelet ({h1}, (double[]){1e-2}, LEVEL, 4); 
    }