tor2fluv-d.c

Transonic shock with Darcy friction

Transonic shock with Darcy friction

Declarations

We use the Saint-Venant solver on a 1D grid and we add the Darcy friction term.

#include "grid/cartesian1D.h"
#include "b-flood/saint-venant-topo.h"
#include "b-flood/darcy.h"

int LEVEL;

scalar e[];
norm nerror;
int ne = 0;
double tmax = 200, q0 = 2, z0;
double dx ;
 
double a1 = 0.674202, a2 = 21.7112, a3 = 14.492, a4 = 1.4305;

// Analytical solution for h(x) and dh/dx
double hex (double x) {
  if (x <= 200/3.)
    return pow(4/G,1/3.)*((4/3.-x/100.)-9*x/1000.*(x/100.-2/3.));
  else
    return pow(4/G,1/3.)*(a1*pow(x/100.-2/3.,4) +
			  a1*cube(x/100.-2/3.) -
			  a2*sq(x/100.-2/3.) +
			  a3*(x/100.-2/3.) + a4);
}

double dhex (double x) {
  if (x <= 200/3.)
    return pow(4/G,1/3.)*(-9*x - 200)/50000.;
  else
    return pow(4/G,1/3.)*(a1*cube(x)/25000000. - a1*sq(x)/200000. +
			  a1*x/7500. + a1/675. - a2*x/5000. + a2/75. + a3/100.);
}

// Darcy friction term in kinematic formulation
double sfd(double x){
  return  -f/(8*G)*sq(q0)/cube(hex(x));
}

// Z and dz/dx
// We use RK4 to solve the topography
double dzex(double x) {
  return (sq(q0)/(G*cube(hex(x))) - 1.)*dhex(x) + sfd(x);
}

double zex(double x, double z) {
  return z + dx/4.*(dzex(x - dx) + 2.*dzex(x - 0.5*dx) + dzex(x));
}

Parameters

Definition of parameters and calling of the Saint-Venant subroutine run().

int main()
{
  f = 0.093;
  L0 = 100.;
  X0 = 0;
  G = 9.81;
  for (LEVEL = 4; LEVEL <= 9; LEVEL++) {  
    N = 1 << LEVEL;
    dx = L0/N;
    run();
    fprintf (stderr, "%d %g %g\n", N, nerror.avg, nerror.rms);
  }
}

Boundary condition

We set h and q (u) at both boundaries (torrential upstream and fluvial downstream).

h[left] = dirichlet(max(hex(0),0));
eta[left] =  dirichlet(max(hex(0)+zb[],zb[]));
u.n[left] = dirichlet(max(q0/hex(0),0));

h[right] = dirichlet(max(hex(100),0));
eta[right] =  dirichlet(max(hex(100)+zb[],zb[]));
u.n[right] = dirichlet(max(q0/hex(100),0));

Initial conditions

event init (i = 0)
{
  // Because the slope is initially dry, we set a maximum time-step. 
  DT = 1e-2;
  z0 = 0;
  foreach(){
    zb[] = zex(x,z0);
    z0 = zb[];
    u.x[] = 0;
  }
  boundary(all);
}

Error norms

We compute the different error norms

event error (i++; t <= tmax) {
  foreach()
    e[] = h[] - hex(x);  
  nerror = normf (e);
}

Gnuplot output

We print the water profile along the channel at final time.

event printprofile (t = tmax)
{
  char name[100];
  FILE * fp;
  sprintf (name, "profil-%i.dat", N);
  fp = fopen(name, "w");
  foreach()
    fprintf (fp, "%g\t%g\t%g\t%g\t%g\n",
	     x, h[], zb[], hex(x), u.x[]);
  fclose(fp);
}

References

[popinet2011]	S. Popinet. Quadtree-adaptive tsunami modelling. Ocean Dynamics, 61(9):1261–1285, 2011. [ .pdf ]
[macdonald1997]	I MacDonald, MJ Baines, NK Nichols, and PG Samuels. Analytic benchmark solutions for open-channel flows. Journal of Hydraulic Engineering, 123(11):1041–1045, 1997. [ DOI ]

Results

set xlabel 'L (m)' 
set ylabel 'Height (m)' 
set xtics 
set ytics
set y2label 'Error (m)'
set y2tics
set key l b
plot [][] './profil-512.dat' u 1:($3+$4) w l lw 0.5 \
                axes x1y1 t 'exact solution :Zb + he', \
          './profil-512.dat' u 1:($2+$3) w l lt 0 lw 7 \
                axes x1y1 t 'N=512 :Zb + h', \
          './profil-32.dat' u 1:($2+$3) axes x1y1 t 'N=32: Zb + h', \
          './profil-512.dat' u 1:3 w l axes x1y1 t 'topo: Zb', \
          './profil-512.dat' u 1:($2-$4) w l \
                axes x1y2 t'error N=512: h - he (right axis)'

reset
set logscale
set xlabel 'Number of cells N' 
set ylabel 'Error norm (m)' 
set xtics 
set ytics
set cbrange [1:2]
ftitle(a,b) = sprintf("order %4.2f", -b)
f1(x) = a1 + b1*x
fit f1(x) 'log' u (log($1)):(log($2)) via a1,b1
f2(x) = a2 + b2*x
fit f2(x) 'log' u (log($1)):(log($3)) via a2,b2
plot exp (f1(log(x))) t ftitle(a1,b1), './log' u 1:2 t 'average error', \
     exp (f2(log(x))) t ftitle(a2,b2), './log' u 1:3 t 'rms error'

Link to the homepage

Homepage