/** # Transcritical flow over a bump with multiple layers This test case is similar to [layered.c]() but using simpler (and more physical) boundary conditions. Both hydrostatic and non-hydrostatic solutions are obtained. The non-hydrostatic solution is compared with [a DNS solution](/src/examples/gaussian-ns.c) using the Navier-Stokes/VOF solver. More details are given in [Popinet (2020)](/Bibliography#popinet2020). */ #include "grid/multigrid1D.h" #if !ML # include "saint-venant.h" #else // ML # include "layered/hydro.h" # define phi q # if !HYDRO # include "layered/nh.h" # endif # include "layered/remap.h" # include "layered/perfs.h" #endif // ML /** The primary parameters are the flow rate, water level, viscosity, bump amplitude and ramping time. */ const double QL = 1., HR = 0.6, NU = 1e-2, BA = 0.4, T0 = 10; int main() { /** The domain is 30 meters long. */ L0 = 30. [1]; G = 9.81; N = 512; // less damping with 1024 /** The viscosity is set to NU and 20 layers are used to discretise vertically. */ nu = NU; nl = 20; // going to 30 changes very little #if ML #if NOMETRIC max_slope = 0.; #endif #if !HYDRO NITERMIN = 2; #endif #endif run(); } /** ## Initialisation and boundary conditions The inflow is a parabolic profile with a total flow rate *Q*. The function below computes the height *zc* of the middle of the layer and returns the corresponding velocity. */ #if !ML double uleft (Point point, scalar s, double Q) { double zc = zb[]; for (int l = 0; l < s.l; l++) zc += layer[l]*h[]; zc += layer[s.l]*h[]/2.; return 3./2.*Q/h[]*(1. - sq(zc/h[] - 1.)); } #else /** For the multilayer solver, we need some trickery to define the inflow velocity profile as a function of $z$. We need to sum the layer thicknesses to get the total depth H and the height zc of the current layer (of index point.l). */ double uleft (Point point, scalar s, double Q) { double H = 0.; double zc = zb[]; for (int l = - point.l; l < nl - point.l; l++) { H += h[0,0,l]; if (l < 0) zc += h[0,0,l]; } zc += h[]/2.; return 3./2.*Q/H*(1. - sq(zc/H - 1.)); } #endif /** We initialise the topography and the initial thickness of each layer *h*. */ event init (i = 0) { double a = 5.; foreach() { zb[] = BA*exp(- sq(x - 10.)/a); #if !ML h[] = HR - zb[]; #else foreach_layer() h[] = (HR - zb[])/nl; #endif } /** The height of the free-surface $\eta$ is imposed on the right boundary. */ eta[right] = dirichlet(HR); /** Boundary conditions are set for the inflow velocity and the outflow of each layer. */ #if !ML for (vector u in ul) { u.n[left] = dirichlet (uleft (point, _s, QL*(t < T0 ? t/T0 : 1.))); u.n[right] = neumann(0.); } h[right] = dirichlet(HR); #else u.n[left] = dirichlet (uleft (point, _s, QL*(t < T0 ? t/T0 : 1.))); u.n[right] = neumann(0); #endif /** In the non-hydrostatic case, a boundary condition is required for the non-hydrostatic pressure $\phi$ of each layer. */ #if !HYDRO && ML phi[right] = dirichlet(0.); #endif } /** We can optionally add horizontal viscosity. */ #if 0 event viscous_term (i++) { // add horizontal viscosity (small influence) #if HYDRO scalar * list = {u.x}; #else scalar * list = {u.x,w}; #endif scalar d2u[]; foreach_layer() for (scalar s in list) { foreach() d2u[] = (u.x[1] + u.x[-1] - 2.*u.x[])/sq(Delta); foreach() u.x[] += dt*nu*d2u[]; } } #endif /** We check for convergence. */ scalar etac[]; event logfile (t += 0.1; i <= 100000) { FILE * fp = fopen ("dh", "w"); foreach() fprintf (fp, "%g %g\n", x, eta[] - etac[]); fclose (fp); double dh = change (eta, etac); if (i > 1 && dh < 5e-5) return 1; } /** Uncomment this part if you want on-the-fly animation. */ #if 0 #include "plot_layers.h" #endif /** ## Outputs We save profiles at regular intervals. */ event profiles (t += 5) { foreach (serial) { #if !ML double H = h[]; #else double H = 0.; foreach_layer() H += h[]; #endif fprintf (stderr, "%g %g %g\n", x, zb[] + H, zb[]); } fprintf (stderr, "\n\n"); } /** For the hydrostatic case, we compute a diagnostic vertical velocity field w. Note that this needs to be done within this event because it relies on the fluxes hu and face heights hf, which are only defined temporarily in the [multilayer solver](hydro.h#update_eta). */ #if HYDRO scalar w = {-1}; event update_eta (i++) { if (w.i < 0) w = new scalar[nl]; vertical_velocity (w, hu, hf); /** The layer interface values are averaged at the center of each layer. */ foreach() { double wm = 0.; foreach_layer() { double w1 = w[]; w[] = (w1 + wm)/2.; wm = w1; } } } #endif // HYDRO /** We also save the velocity and non-hydrostatic pressure fields. */ event output (t = end) { foreach (serial) { double z = zb[]; #if HYDRO printf ("%g %g %g %g\n", x, z, u.x[], w[]); foreach_layer() { z += h[]; printf ("%g %g %g %g\n", x, z, u.x[], w[]); } #elif !ML vector u = ul[0]; scalar w = wl[0]; printf ("%g %g %g %g\n", x, z, u.x[], w[]); int l = 0; for (u,w in ul,wl) { z += layer[l++]*h[]; printf ("%g %g %g %g\n", x, z, u.x[], w[]); } #else // ML printf ("%g %g %g %g %g\n", x, z, u.x[], w[], phi[]); foreach_layer() { z += h[]; printf ("%g %g %g %g %g\n", x, z, u.x[], w[], phi[]); } #endif // ML printf ("\n"); } #if HYDRO delete ({w}); #endif printf ("# end = %g\n", t); } /** ## Results We compare the hydrostatic results obtained with the [layered solver](/src/layered/hydro.h) and those obtained with the [Saint-Venant multilayer code](/src/multilayer.h). ~~~gnuplot Evolution of the free surface (hydrostatic). set yr [0.5:1.2] set xlabel 'x' set ylabel 'z' plot '../gaussian-hydro/log' u 1:2 w l t 'Multilayer', \ '../gaussian-stvt/log' u 1:2 w l t 'De Vita et al.' ~~~ The results which follow are obtained with the non-hydrostatic solver. ~~~gnuplot Evolution of the free surface (non-hydrostatic) plot 'log' u 1:2 w l t '' ~~~ We compare the results obtained with the layered solver to those obtained with the Navier-Stokes VOF solver (and without metric). ~~~gnuplot Evolution of the free surface for both solvers plot 'log' w l t 'Multilayer', \ '../../examples/gaussian-ns/log' w l t 'Navier-Stokes VOF' ~~~ ~~~gnuplot Final free surface profiles plot 'log' index 11 w l t 'Multilayer', \ '../../examples/gaussian-ns/log' index 'prof70' w l t 'Navier-Stokes VOF', \ '../gaussian-nometric/log' index 12 w l t 'no metric' ~~~ ~~~gnuplot Horizontal velocity field { width=100% } set term PNG enhanced font ",15" size 2048,512 set yr [0:1.2] set output 'vel.png' set pm3d set pm3d map interpolate 10,4 unset key # jet colormap set palette defined ( 0 0 0 0.5647, 0.125 0 0.05882 1, 0.25 0 0.5647 1, \ 0.375 0.05882 1 0.9333, 0.5 0.5647 1 0.4392, 0.625 1 0.9333 0, 0.75 1 0.4392 0, \ 0.875 0.9333 0 0, 1 0.498 0 0 ) set contour base set cntrparam levels discrete 0 set cntrlabel onecolor set cntrparam bspline splot 'out' u 1:2:3 lt 3 lc rgb "#ffffff" lw 2 ~~~ ~~~gnuplot Vertical velocity field { width=100% } set term PNG enhanced font ",15" size 2048,512 set output 'w.png' # set cbrange [-0.1:0.1] unset contour splot 'out' u 1:2:4 ~~~ ~~~gnuplot Non-hydrostatic pressure { width=100% } set term PNG enhanced font ",15" size 2048,512 set output 'phi.png' # set cbrange [-0.2:0.35] splot 'out' u 1:2:5 ~~~ ## See also * [Same case with the Navier-Stokes/VOF solver](/src/examples/gaussian-ns.c) */