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:
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