src/test/implicit.c

Implicit Saint-Venant solutions for waves

We solve the Saint-Venant equations either explicitly or implicitly, in one dimension.

#include "grid/multigrid1D.h"
#if EXPLICIT
# if ML
#   include "layered/hydro.h"
# else
#   include "saint-venant.h"
# endif
#else
# if ML
#   include "layered/hydro.h"
#   include "layered/implicit.h"
# else
#   include "saint-venant-implicit.h"
# endif
#endif

We vary the wave amplitude to test both linear and nonlinear regimes.

double amp;

int main()
{
  init_grid (512);
  periodic (right);

We set the acceleration of gravity and the maximum timestep (for the implicit solver).

  G = 10.;
  DT = 5e-3;
  theta = 1.3;

We do several runs with increasing wave amplitude.

  amp = 0.1;
  run();
  amp = 0.4;
  run();
  amp = 1.;
  run();
  system ("cat out-* > log");
}

The initial condition is a simple sinusoidal wave. The water depth is one.

event init (i = 0) {
  foreach()
    h[] = 1. + amp*sin(2.[-1]*pi*x);
}

We log the evolution of the timestep.

event logfile (i++) {
  char name[80];
  sprintf (name, "log-%g", amp);
  static FILE * fp = fopen (name, "w");
  fprintf (fp, "%g %g %g\n", t, dt, statsf(h).sum);
}

We output wave profiles at a few time intervals.

event output (t <= 0.2; t += 0.1) {
  char name[80];
  sprintf (name, "out-%g", amp);
  static FILE * fp = fopen (name, "w");
  foreach()
    fprintf (fp, "%g %g\n", x, h[]);
  fprintf (fp, "\n");
  fflush (fp);
}

For small wave amplitudes, non-linear effects are small and we get the following wave evolution for the explicit and implicit solvers.

set xlabel 'x'
set ylabel 'z'
plot 'out-0.1' w l t 'implicit', '../explicit/out-0.1' w l t 'explicit', \
     '../implicit-ml/out-0.1' w l t 'implicit (ML)', \
     '../explicit-ml/out-0.1' w l t 'explicit (ML)'

Explicit and implicit solutions for a wave amplitude of 0.1 (script)

The solutions are not too far given that the explicit scheme requires 698 time steps to verify the CFL condition, against only 51 steps for the implicit scheme.

For an amplitude of 0.4, the explicit solution is close to a weak shock and requires 840 timesteps against 186 timesteps for the implicit scheme.

plot 'out-0.4' w l t 'implicit', '../explicit/out-0.4' w l t 'explicit', \
     '../implicit-ml/out-0.4' w l t 'implicit (ML)', \
     '../explicit-ml/out-0.4' w l t 'explicit (ML)'

Explicit and implicit solutions for a wave amplitude of 0.4 (script)

For an amplitude of 1., strong shocks develop for both solutions. The implicit scheme is a bit more oscillatory near the shocks but this remains acceptable. The explicit solution requires 1092 steps against 471 steps for the implicit scheme.

plot 'out-1' w l t 'implicit', '../explicit/out-1' w l t 'explicit', \
     '../implicit-ml/out-1' w l t 'implicit (ML)', \
     '../explicit-ml/out-1' w l t 'explicit (ML)'     

Explicit and implicit solutions for a wave amplitude of 1 (script)