sandbox/b-flood/Testcases/tor2fluv-d.c
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
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'