sandbox/M1EMN/Exemples/slope_inc.c

    Une vague se propage en eau peu profonde de profondeur constante et arrive sur une plage en pente droite:

    la plage n’est pas parallele à la vague initiale. Résolution de Saint Venant simple en 2D.

    avec au temps t=0, h=1+aexp(-x^2) et un fond z=0 pour x<x_0 qui remonte avec un angle \alpha ensuite.

    On utilise le solveur de Saint-Venant et on utilise l’adaptativite:

    #include "saint-venant.h"
    
    double tmax,Cf;
    #define LEVEL 9

    CL

    u.n[left] = neumann(0);
    u.t[left] = neumann(0);
    h[left] = h[];
    u.n[top] = neumann(0);
    u.t[top] = neumann(0);
    h[top] = h[];

    fond plat puis en pente avec un angle de côté, une vague arrive

    event init (i = 0)
    { double xi;
      foreach(){
         xi=x+y/4.+1;   // fond incline
         zb[] =   (xi>-1)*(xi+1)/7.;    
         xi=x+7;
         h[]=fmax(1+.2*exp(-2*xi*xi)-zb[],0);
         u.x[]=.2*exp(-2*xi*xi);     
             }
         boundary (all); 
    }

    en fait 0

    event friction (i++) {
      // Poiseuille friction (implicit scheme)
      foreach()
      {double ff = h[] < dry ? HUGE : (1. + Cf*dt/h[]/h[]);
       foreach_dimension()
       u.x[] /= ff;
        }
    }

    adapatif

    event adapt (i++) {
      astats s = adapt_wavelet ({h}, (double[]){1e-3}, LEVEL);
      fprintf (stderr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);
    }

    film et plot

    event movies (t <= tmax;t+=.1) {
      static FILE * fp = popen ("ppm2mpeg > eta.mpg", "w");
      scalar m[], etam[];
      
      foreach() {
        etam[] = eta[]*(h[] > dry);
        m[] = etam[] - zb[];
      }
      boundary ({m, etam});
      output_ppm (etam, fp, mask = m, min = 0.9, max = 1.2, n = 512, linear = true);}
    
    
     
    event cut(t <= tmax;t+=.5) {
    	double dx=8*L0/pow(2.,LEVEL),dy=dx;
    	fprintf (stderr,"--- %lf \n",t); 
    	  printf("set hidden3d; sp[][][0:1.25]  '-' u 1:2:($3+$4) not w l linec 3 ,''u 1:2:4 not w l linec 2\n");
        for(double x=-L0/2.;  x<L0/2.;x+=dx)
         {for(double y=-L0/2.;y<L0/2.;y+=dy)
          { 
         printf (" %g %g %g  %g \n", x, y, interpolate (h, x, y),  interpolate (zb, x, y) );}
         printf ("\n");
         } 
         printf ("e\n");
    }

    paramètres

    int main() {
      L0 = 20.;
      X0 = -L0/2.;
      Y0 = -L0/2.;
      G = 1;
      N = 1 << LEVEL;  // ie 2^LEVEL
      tmax = 40;
      // friction nulle
      Cf = 0.000000;
      
      run();
    }

    Pour exécuter, le plus simple:

    qcc -g -O2 -DTRASH=1 -Wall  -o slope_inc slope_inc.c -lm
    ./slope_inc | gnuplot

    ce qui fait apparaître dans la fenêtre la simulation au cours du temps

    si on veut reproduire:

    plage

    plage

    alors il faut faire

    qcc -g -O2 -DTRASH=1 -Wall  -o slope_inc slope_inc.c -lm
    ./slope_inc > dmp
    echo "set term gif animate; set output 'slope_inc.gif'" > slope_inc.out
    cat dmp >> slope_inc.out
    rm dmp
    cat slope_inc.out | gnuplot