sandbox/eddie/refract.c
Waves refraction on a shallow elliptical reef platform
#include "grid/multigrid.h"
#include "green-naghdi.h"
#define HS 1.5
#define TS 12
#define DEPTH 1.5
#define XYR 0.70
#define XW 708
#define RROTATE 0
int main()
The domain is 1536m long (reef flat is 1000m x 700m) This is a low resolution run (N=256 gives 6m cells) for higher resolution change N to 512 (3m) or 1024 (1.5m)
{
L0 = 1536;
X0 = -L0/2.;
Y0 = -L0/2.;
G = 9.81;
N = 256;
breaking = 1;
run();
}
monochromatic wave field is generated at the left boundary
u.x[left] = - radiation ((1*HS)*sin(2.*pi*t/TS));
u.x[right] = + radiation (0);
Define average and maximum stats variables
scalar maxa[];
scalar SH[];
scalar SHc[];
scalar AH[];
scalar Vel[];
scalar SVel[];
scalar AVel[];
scalar vmax[];
scalar vmin[];
scalar Sux[];
scalar Suy[];
scalar Aux[];
scalar Auy[];
event init (i = 0)
{
bathymetry is an elongate reef platform with ~35 degree slope and flat 1.5m deep reef flat
double h0 = 100;
double cosa = cos (RROTATE*pi/180.), sina = sin (RROTATE*pi/180.);
foreach() {
double xr = x*cosa - y*sina, yr = x*sina + y*cosa;
double z0 = 0;
xr >= -5.82 ? (5.82 + xr)/50. : 0.;
double zs = sq(xr/XW) + sq(yr/(XW*XYR)) <= 1. ?
1 + 50*(sqrt(1 - sq(xr/XW) - sq(yr/(XW*XYR)))) : 0.;
zb[] = min(-h0+(max((z0 + (zs*5.5) - h0), 0)),0) - DEPTH;
h[] = max (0., -zb[]);
eta[] = h[] > dry ? h[] + zb[] : 0;
maxa[] = 0.;
}
}
quadratic bottom firction with extra friction added near right boundary to prevent reflection
event friction (i++) {
foreach()
{
double a = x < 650 ? h[] < dry ? HUGE : 1. + 0.01*dt*norm(u)/h[] :
h[] < dry ? HUGE : 1. + 2*(x - 650.)*dt*norm(u)/h[];
foreach_dimension()
u.x[] /= a;
Vel[] = norm(u);
}
}
event logfile ( t+= 1) {
stats s = statsf (h);
norm n = normf (u.x);
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 %g %g %g %g %g %g\n", t, i, s.min, s.max, s.sum,
n.rms, n.max, dt);
}
MOVIES
event movies (t += 0.5) {
static FILE * fp = popen ("ppm2mpeg > Vel.mpg", "w");
scalar m[], etam[];
foreach() {
etam[] = Vel[]*(h[] > dry);
m[] = etam[] - zb[];
}
boundary ({m, etam});
output_ppm (etam, fp, mask = m, min = 0, max = 2, n = 512, linear = true);
}
max and time-averaged stats are calculated after the wave field is developed on the reef platform
Set amax as highest wave amplitude
if (fabs(eta[]) > maxa[])
maxa[] = fabs(eta[]);
Set SH as the sum water level
if (h[] > dry)
SH[] = SH[] + (h[] + zb[]);
Set SHc as the number of sum water level
if (h[] > dry)
SHc[] = SHc[] + 1;
Set AH as the average water level
if (h[] > dry)
AH[] = SH[] / SHc[];
Set SVel as the sum velocity
if (h[] > dry)
SVel[] = SVel[] + (Vel[]);
Set AH as the average velocity
if (h[] > dry)
AVel[] = SVel[] / SHc[];
Set vmax as highest velocity
if (h[] > dry && Vel[] > vmax[])
vmax[] = Vel[];
Set vmax as highest velocity
if (h[] > dry && Vel[] < vmin[])
vmin[] = Vel[];
Set Sux as the sum u.x
if (h[] > dry)
Sux[] = Sux[] + u.x[];
Set Suy as the sum u.y
if (h[] > dry)
Suy[] = Suy[] + u.y[];
Set Aux as the average u.x
if (h[] > dry)
Aux[] = Sux[] / SHc[];
Set Auy as the average u.y
if (h[] > dry)
Auy[] = Suy[] / SHc[];
}
}
grid data is output every 300 seconds and at the end
event out (t += 300) {
char name[20]; sprintf (name, "outzxy-%g.out", t);
FILE * fp = fopen (name, "w");
output_field ({eta,zb,maxa,AH,AVel,Vel,vmax,vmin,u,Aux,Auy}, fp, linear = true);
}
event end (t = 1209) {
FILE * fp = fopen ("end", "w");
output_field ({eta,zb,maxa,AH,AVel,Vel,vmax,vmin,u,Aux,Auy}, fp, linear = true);
}
set term @PNG enhanced size 640,640 font ",8"
set pm3d map
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 )
set size ratio -1
set xlabel 'x (m)'
set ylabel 'y (m)'
splot [-768:768][-768:768]'end' u 1:2:4 w pm3d t ''
a0 = 1
splot [-768:768][-768:768]'end' u 1:2:($3/a0) w pm3d t ''
a0 = 1
splot [-768:768][-768:768]'end' u 1:2:($5/a0) w pm3d t ''
a0 = 1
splot [-768:768][-768:768]'end' u 1:2:($6/a0) w pm3d t ''
a0 = 1
splot [-768:768][-768:768]'end' u 1:2:($7/a0) w pm3d t ''
a0 = 1
splot [-768:768][-768:768]'end' u 1:2:($8/a0) w pm3d t ''