# src/examples/swasi.c

# A Shallow Water Analogue for the Standing Accretion Shock Instability

We use the Saint-Venant solver to reproduce the experimental setup of Foglizzo et al. 2012. We use radial coordinates and Basilisk View to visualise the results.

```
#include "grid/multigrid.h"
#include "radial.h"
#include "saint-venant.h"
#if dimension > 1
# include "view.h"
#endif
```

We use seven levels of refinement (i.e. ${128}^{2}$ grid points) by default. The geometrical parameters and other inputs are those of the experimental setup. Note that we chose centimetres and seconds as length and time units.

```
int LEVEL = 7;
#define Ri 4. // inner radius (cm)
#define Ro 32. // outer radius (cm)
#define Ho 0.074 // injection layer thickness (cm)
#define Q 1e3 // flow rate (cm^3/s)
#define R_45 5.6 // hyperbolic surface scaling (cm)
```

We double the equivalent viscosity compared to that given in Foglizzo et al. 2012 (0.03 cm^2/s), to obtain a more stable spiral mode. The “tube height” is set (by trial and error using 1D runs) so that the shock radius obtained numerically is close to 20 cm.

```
#define NU 0.06 // equivalent viscosity (cm^2/s)
#define STEP 6.98 // height of exit tube (cm)
```

The injection velocity is the flow rate divided by the surface area of the outer injection ring.

`#define Uo (Q/(2.*pi*Ro*Ho))`

The maximum runtime (seconds).

`#define TMAX 150.`

## Boundary conditions

The outer injection ring corresponds to the right boundary of the radial coordinate system.

Note that this code can also be run in one (radial) dimension by changing the first line to “grid/multigrid1D.h”.

```
u.n[right] = - Uo;
#if dimension > 1
u.t[right] = dirichlet(0.);
#endif
h[right] = Ho;
```

The inner outflow is the left boundary.

`u.n[left] = dirichlet(- Ro*Uo*Ho/(STEP*Ri));`

## Main program

The level of refinement can be given as a command-line argument.

```
int main (int argc, char * argv[])
{
if (argc > 1)
LEVEL = atoi (argv[1]);
dtheta = 2.*π;
G = 981.; // acceleration of gravity (cm.s^-2)
size (Ro - Ri);
origin (Ri, 0.);
#if dimension > 1
periodic (top);
#endif
init_grid (1 << LEVEL);
run();
}
```

## Initial conditions

We setup the hyperbolic bottom profile and set an initial uniform fluid layer and radial velocity identical to the injection conditions.

```
event init (i = 0)
{
foreach() {
zb[] = - sq(R_45)/r;
h[] = Ho;
u.x[] = - Uo;
}
zb[left] = dirichlet(- sq(R_45)/r);
zb[right] = dirichlet(- sq(R_45)/r);
}
```

## Perturbation

At $t=40$ seconds the shock is close to being stationary, we perturb the velocity field to trigger the instability.

We add a systematic mode 1 perturbation of the radial component of the velocity, tapering off at ${R}_{i}$ and ${R}_{o}$. This gives a clean growth of the mode and allows a simple estimate of the frequency and growth rate (see below).

Optionally we can also add a small, constant azimuthal perturbation which will trigger the “spiral” instability (see second movie below).

```
event perturb (t = 40) {
foreach() {
u.x[] += 1e-2*Uo*((r - Ri)/(Ro - Ri))*(1. - (r - Ri)/(Ro - Ri))*sin(θ);
#if 0
u.y[] += 1e-2*Uo;
#endif
}
boundary ((scalar *){u});
}
```

## Friction

We use a time-implicit discretisation of the linear friction.

```
event friction (i++) {
if (NU > 0.)
foreach() {
double a = h[] < dry ? HUGE : 1. + NU*dt/sq(h[]);
foreach_dimension()
u.x[] /= a;
}
}
```

## Outputs

We log the evolution of the fluid height and norms of the azimuthal (in 2D) or radial (in 1D) component of the velocity as well as the coordinates of the center-of-mass of the fluid.

```
event logfile (i += 10; t <= TMAX) {
stats s = statsf (h);
#if dimension > 1
norm n = normf (u.y);
#else
norm n = normf (u.x);
#endif
double sumx = 0., sumy = 0.;
foreach (reduction(+:sumx) reduction(+:sumy))
sumx += r*cos(θ)*h[]*dv(),
sumy += r*sin(θ)*h[]*dv();
fprintf (ferr, "%g %d %g %g %g %g %g %g %g\n", t, i, s.min, s.max, s.sum,
n.rms, n.max, sumx/s.sum, sumy/s.sum);
}
```

The evolution with time of the y-coordinate of the center of mass is illustrated below. It clearly shows the exponential growth of the oscillation of this position followed by non-linear saturation (see also the animation below).

This can be used to estimate the period of oscillation (2.74 s) and the initial growth rate (0.07 ${s}^{-1}$) as illustrated below

While we run, we display the evolution of a vertical cross-section of the flow.

```
event profiles (i += 100) {
static FILE * fp = popen ("gnuplot", "w");
if (i == 0)
fprintf (fp,
"set term x11 noraise\n"
"set grid\n"
"set xrange [%g:%g]\n", X0, X0 + L0);
fprintf (fp,
"set title 't = %.2f'\n"
"p '-' u 1:3:2 w filledcu lc 3 t ''\n", t);
foreach_leaf()
#if dimension > 1
if (y < Δ)
#endif
fprintf (fp, "%g %g %g\n", x, η[], zb[]);
fprintf (fp, "e\n\n");
fflush (fp);
}
```

We save the stationary profile at $t=50$.

```
event stationary (t = 50) {
FILE * fp = fopen ("stationary", "w");
foreach_leaf()
#if dimension > 1
if (y < Δ)
#endif
fprintf (fp, "%g %g %g %g\n", x, η[], zb[], u.x[]);
fclose (fp);
}
```

## Animation

We create an animation of the free surface height $\eta $. We first define a coordinate mapping function which will be used by Basilisk View.

```
#if dimension > 1
void radial (coord * p) {
double r = p->x, θ = p->y*dtheta/L0;
p->x = r*cos(θ), p->y = r*sin(θ);
}
event movie (t = 60; t += 0.1) {
view (map = radial, fov = 45, width = 600, height = 600, samples = 1);
clear();
squares ("eta", min = -1.4, max = -0.2, linear = true);
save ("eta.mp4");
}
#endif
```

## References

[foglizzo2012] |
Thierry Foglizzo, Frédéric Masset, Jérôme Guilet, and Gilles Durand. Shallow water analogue of the standing accretion shock instability: Experimental demonstration and a two-dimensional model. |