# src/test/multiriverinflow.c

# Flow rates for multiple rivers

In this example, we impose different flow rates on different rivers situated on the same boundary of a Saint-Venant simulation.

```
#include "saint-venant.h"
#include "discharge.h"
```

The domain is 10 metres squared, centered on the origin. Time is in seconds.

```
#define LEVEL 7
int main()
{
size (10.);
origin (- L0/2., - L0/2.);
G = 9.81;
N = 1 << LEVEL;
run();
}
```

## Initial conditions

We chose a reasonably complicated river bed with two “valleys” (see Figure below).

We allocate a new field *river* which is set to different values (1 and 2) for each of the rivers.

```
scalar river[];
event init (i = 0)
{
```

We start with a dry riverbed, so that the problem does not have a natural timescale the Saint-Venant solver can use. We set a maximum timestep to set this timescale.

```
DT = 1e-2;
foreach() {
zb[] = 0.05*pow(x,4) - x*x + 2. + 0.2*(y + Y0);
river[] = x < 0 ? 1 : 2;
}
boundary ({zb, river});
}
```

## Boundary conditions

We impose inflow/outflow on both the top and bottom boundary. In addition, the tangential velocity on the top boundary *u.t* is set to zero.

```
u.n[top] = neumann(0);
u.t[top] = dirichlet(0);
u.n[bottom] = neumann(0);
```

To impose a given flow rate for the two rivers on the top boundary, we compute the elevation $\eta $ of the water surface necessary to match this flow rate. This gives two elevation values *eta1* and *eta2*, one for each river. The flow rates are set to 4 and 2 m^{3}/sec for river 1 and 2 respectively.

```
double eta1, eta2;
event inflow (i++) {
eta1 = eta_b (4, top, river, 1);
eta2 = eta_b (2, top, river, 2);
```

Once we have the required values for the water surface elevations at the top boundary of both riverbeds, we impose them on both $h$ and $\eta =h+{z}_{b}$.

```
h[top] = max ((river[] == 1 ? eta1 : eta2) - zb[], 0);
η[top] = max ((river[] == 1 ? eta1 : eta2) - zb[], 0) + zb[];
}
```

## Outputs

We compute the evolution of the water volumes in both riverbeds.

```
event volume (i += 10) {
double volume1 = 0, volume2 = 0;
foreach(reduction(+:volume1) reduction(+:volume2)) {
double dv = h[]*sq(Δ);
if (x < 0) volume1 += dv;
else volume2 += dv;
}
fprintf (stderr, "%g %g %g %g %g\n",
t, volume1, volume2, eta1, eta2);
}
```

We use gnuplot to produce an animation of the water surface.

```
event init_animation (i = 0) {
printf ("set view 80,05\n"
"set xlabel 'x'\n"
"set ylabel 'y'\n"
"set zlabel 'z'\n"
"set hidden3d; unset ytics ; unset xtics\n");
}
event animation (t <= 1; i += 10) {
double dx = 2.*L0/N, dy = dx;
printf ("set title 't = %.3f'\n"
"sp [%g:%g][%g:%g][-5:5] '-'"
" u 1:2:($3+$4-.05) t 'free surface' w l lt 3,"
" '' u 1:2:4 t 'topography' w l lt 2\n",
t, X0, -X0, Y0, -Y0);
for (double x = X0; x <= X0 + L0; x += dx) {
for (double y = Y0; y <= Y0 + L0; y += dy)
printf ("%g %g %g %g\n",
x, y, interpolate (h, x, y), interpolate (zb, x, y));
putchar ('\n');
}
printf ("e\n"
"pause %.5lf \n\n", 0.);
}
```