sandbox/b-flood/Testcases/fluv2tor-m.c
Transcritical test case with Manning friction
Declarations
We call the Saint-Venant solver on a 1D grid and we add the Manning friction term.
#include "grid/cartesian1D.h"
#include "b-flood/saint-venant-topo.h"
#include "b-flood/manning.h"
int LEVEL;
scalar e[];
norm nerror;
double tmax = 1000, q0 = 2, z0;
double dx;
// Analytical solution for h(x) and dh/dx
double hex (double x) {
if (x <= 500)
return pow(4/G,1/3.)*(1 - 1/3.*tanh(3*(x/1000.-1/2.)));
else
return pow(4/G,1/3.)*(1 - 1/6.*tanh(6*(x/1000.-1/2.)));
}
double dhex (double x) {
if (x <= 500)
return -pow(4/G,1/3.)/(1000*sq(cosh(3/2.-3*x/1000.)));
else
return -pow(4/G,1/3.)/(1000*sq(cosh(3-3*x/500.)));
}
// Manning's friction term
double sfm (double x) {
return -sq(n)*sq(q0)/pow(hex(x),10/3.);
}
// Z and dz/dx
// We use Runge-Kuta 4 algo to fix the topography
double dzex (double x) {
return (sq(q0)/(G*cube(hex(x)))-1)*dhex(x) + sfm(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-topo subroutine run().
int main()
{
n = 0.0218;
L0 = 1000.;
X0 = 0;
G = 9.81;
tmax = 1000;
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 fix h and u at the left boundary (fluvial).
h[left] = dirichlet(max(hex(0),0));
eta[left] = dirichlet(max(hex(0)+zb[],zb[]));
u.n[left] = dirichlet(max(q0/hex(0),0));
We fix a free exit condition on the right boundary (torrential).
u.n[right] = neumann(0);
h[right] = neumann(0);
eta[right] = neumann(0);
Initial conditions
event init (i = 0)
{
// Because the slope is initially dry, we set a maximum time-step.
DT = 1e-2;
// the topography start at the altitude z = 0 at the left of the domain
z0 = 0;
foreach() {
zb[] = zex(x,z0);
z0 = zb[];
u.x[] = 0;
h[] = 0;
}
boundary (all);
}
Error norms
We compute the differents error norms.
Gnuplot output
We print the water profile along the channel at the final time for each resolution.
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
[macdonald1997analytic] |
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
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'
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'