# src/examples/breaking.c

# 3D breaking Stokes wave (multilayer solver)

A steep, 3D, third-order Stokes wave is unstable and breaks. This is the 3D equivalent of this test case. The bathymetry is given by \displaystyle z_b(y) = - 0.5 + \sin(\pi y)/4 i.e. it is shallower toward the back of the domain which causes the wave to break earlier there.

The solution is obtained using the layered model and demonstrates its robustness and a degree of realism even for this complex case.

```
#include "layered/hydro.h"
#include "layered/nh.h"
#include "layered/remap.h"
#include "layered/perfs.h"
#include "view.h"
```

The initial conditions are given by the wave steepness ak and the Reynolds number Re=c\lambda/\nu with c the phase speed of the gravity wave and \lambda its wavelength.

```
double ak = 0.33;
double RE = 40000.;
#define k_ (2.*pi)
#define h_ 0.5
#define g_ 1.
#define T0 (k_/sqrt(g_*k_))
```

The domain is periodic in x and resolved using 256^2 points and 30 layers.

Some implicit damping is necessary to damp fast modes. This may be related to the slow/fast decoupling of the \theta-scheme mentioned by Durran & Blossey, 2012.

```
theta_H = 0.51;
run();
}
```

The initial conditions for the free-surface and velocity are given by the third-order Stokes solution.

We can use a larger CFL, in particular because we are not dealing with shallow-water/wetting/drying.

` CFL = 0.8;`

The layer thicknesses follow a geometric progression, starting from a top layer with a thickness of 1/3 that of the regular distribution.

```
geometric_beta (1./3., true);
foreach() {
zb[] = - 0.5 + sin(pi*y)/4.;
double H = wave(x, 0) - zb[];
double z = zb[];
foreach_layer() {
h[] = H*beta[point.l];
z += h[]/2.;
u.x[] = u_x(x, z);
w[] = u_y(x, z);
z += h[]/2.;
}
}
}
```

We add (an approximation of) horizontal viscosity.

```
event viscous_term (i++)
horizontal_diffusion ((scalar *){u}, nu, dt);
```

We log the evolution of the kinetic and potential energies.

```
set xlabel 't/T0'
plot [0:6]'log' u 1:2 w l t 'kinetic', '' u 1:3 w l t 'potential', \
'' u 1:(($2+$3)/2.) w l t 'total/2'
```

```
event logfile (i++; t <= 8.*T0)
{
double ke = 0., gpe = 0.;
foreach (reduction(+:ke) reduction(+:gpe)) {
foreach_layer() {
double norm2 = sq(w[]);
foreach_dimension()
norm2 += sq(u.x[]);
ke += norm2*h[]*dv();
}
gpe += sq(eta[])*dv();
}
fprintf (stderr, "%g %g %g\n", t/T0, ke/2., g_*gpe/2.);
}
```

And generate the movie of the free surface (this is quite expensive). The movie is 45 seconds at 25 frames/second.

```
event movie (t += 8.*T0/(45*25))
{
view (fov = 17.3106, quat = {0.475152,0.161235,0.235565,0.832313},
tx = -0.0221727, ty = -0.0140227, width = 1200, height = 768);
char s[80];
sprintf (s, "t = %.2f T0", t/T0);
draw_string (s, size = 80);
for (double x = -1; x <= 1; x++)
translate (x) {
squares ("u29.x", linear = true, z = "eta", min = -0.15, max = 0.6);
}
save ("movie.mp4");
}
```

## Parallel run

The simulation was run in parallel on the Occigen machine on 64 cores, using this script

```
local% qcc -source -D_MPI=1 breaking.c
local% scp _breaking.c user@occigen.cines.fr:
```

```
#!/bin/bash
#SBATCH -J breaking
#SBATCH --constraint=HSW24
#SBATCH --ntasks=64
#SBATCH --threads-per-core=1
#SBATCH --time=1:00:00
#SBATCH --exclusive
#SBATCH --mail-type=BEGIN,END,FAIL
#SBATCH --mail-user=popinet@basilisk.fr
module purge
module load openmpi
module load intel gcc
NAME=breaking
mpicc -Wall -std=c99 -O2 _$NAME.c -o $NAME \
-I/home/popinet/local -L/home/popinet/local/gl -L/home/popinet/local/lib \
-lglutils -lfb_tiny -lppr -lgfortran -lm
export LD_LIBRARY_PATH=/home/popinet/local/lib:$LD_LIBRARY_PATH
export PATH=$PATH:/home/popinet/local/bin
rm -f *.ppm
srun --mpi=pmi2 -K1 --resv-ports -n $SLURM_NTASKS $NAME 2> log > out
```

The number of timesteps was 4159 and the runtime was 44 minutes with movie generation.

## References

[durran2012] |
Dale R Durran and Peter N Blossey. Implicit–explicit multistep methods for fast-wave–slow-wave problems. |