Note
Click here to download the full example code
Dam break¶
In this example the classic test of the dambreak for the Saint-Venant system is illustrated. Comparison between 2D numerical solution and 1D analytical solution is carried out for the wet and dry dambreak cases.
Note
1D analytical solutions are given in: https://hal.archives-ouvertes.fr/hal-00628246/document at section 4.1, paragraph “Dam breaks”.
import os, sys
import argparse
import matplotlib.pylab as plt
import freshkiss3d as fk
import freshkiss3d.extra.plots as fk_plt
parser = argparse.ArgumentParser()
parser.add_argument('--nographics', action='store_true')
args = parser.parse_args()
#sphinx_gallery_thumbnail_number = 3
Riemann problems:¶
Two cases (Riemann problems) are considered in this example:
Dam break on a wet domain without friction (case = ‘wet’)
\[\begin{split}h(x) = \begin{cases} h_l \quad \text{for} \quad 0 \leq x \leq x_d \\ h_r \quad \text{for} \quad x_d < x \leq L \end{cases}\end{split}\]Dam break on a dry domain without friction (case = ‘dry’)
\[\begin{split}h(x) = \begin{cases} h_l > 0 \quad \text{for} \quad 0 \leq x \leq x_d \\ h_r = 0 \quad \text{for} \quad x_d < x \leq L \end{cases}\end{split}\]
CASE = 'wet'
if CASE == 'wet':
H_L = 0.005
H_R = 0.001
elif CASE == 'dry':
H_L = 0.005
H_R = 0.
FINAL_TIME = 10.
X_D = 5
Time loop:¶
simutime = fk.SimuTime(final_time=FINAL_TIME, time_iteration_max=100000,
second_order=True)
# Plot figure scheduler:
WHEN = [FINAL_TIME/2., FINAL_TIME-1.]
create_figure_scheduler = fk.schedules(times=WHEN)
Mesh:¶
\[L = 10 m, l = 0.5 m\]
dir_path = os.path.abspath(os.path.dirname(sys.argv[0]))
triangular_mesh = fk.TriangularMesh.from_msh_file(dir_path+'/inputs/dambreak.msh')
if not args.nographics:
fk_plt.plot_mesh(triangular_mesh, plot_labels=False)
1D Analytic solutions:¶
1D analytic solution is defined in a DamBreak class. Computation of the solution depends on initial height values on each side of the dam.
dambreak_analytic = fk.DamBreak(triangular_mesh,
h_l=H_L,
h_r=H_R,
x_d=X_D,
scheduler=fk.schedules(times=WHEN))
dambreak_analytic(0.)
In the ‘wet’ case analytical solution is given by:
\[\begin{split}h(x,t) = \begin{cases}
h_l \quad & \text{if} \quad x \leq x_A(t) \\
\dfrac{4}{9g} \left( \sqrt{gh_l} - \dfrac{x-x_0}{2t} \right)^2 \quad & \text{if} \quad x_A(t) \leq x \leq x_B(t) \\
\dfrac{c_m^2}{9} \quad & \text{if} \quad x_B(t) \leq x \leq x_C(t) \\
h_r \quad & \text{if} \quad x_C(t) \leq x
\end{cases}\end{split}\]
\[\begin{split}u(x,t) = \begin{cases}
0 \quad & \text{if} \quad x \leq x_A(t)\\
\dfrac{2}{3} \left( \dfrac{x-x_0}{t} + \sqrt{gh_l} \right) \quad & \text{if} \quad x_A(t) \leq x \leq x_B(t) \\
2 \left( \sqrt{gh_l} - c_m\right) \quad & \text{if} \quad x_B(t) \leq x \leq x_C(t)\\
0 \quad & \text{if} \quad x_C(t) \leq x
\end{cases}\end{split}\]
\[\begin{split}\begin{cases}
x_A(t) = x_0 - t \sqrt{g h_l} \\
x_B(t) = x_0 + t \left( 2 \sqrt{gh_l} - 3c_m \right) \\
x_C(t) = x_0 + t \left( \dfrac{2c_m^2 \left( \sqrt{gh_l} - c_m \right)}{c_m^2-g h_r} \right)
\end{cases}\end{split}\]
\[c_m \quad \text{being solution of} \quad -8g h_r c_m^2 (\sqrt{gh_l}-c_m)^2+(c_m^2-gh_r)^2(c_m^2+gh-r)\]
If case is set to ‘dry’ analytical solution is given by:
\[\begin{split}h(x,t) = \begin{cases}
h_l \quad & \text{if} \quad x \leq x_A(t) \\
\dfrac{4}{9g} \left( \sqrt{gh_l} - \dfrac{x-x_0}{2t} \right)^2 \quad & \text{if} \quad x_A(t) \leq x \leq x_B(t) \\
0 \quad & \text{if} \quad x_B(t) \leq x
\end{cases}\end{split}\]
\[\begin{split}u(x,t) = \begin{cases}
0 \quad & \text{if} \quad x \leq x_A(t) \\
\dfrac{2}{3} \left( \dfrac{x-x_0}{t} + \sqrt{gh_l} \right) \quad & \text{if} \quad x_A(t) \leq x \leq x_B(t)\\
0 \quad & \text{if} \quad x_B(t) \leq x
\end{cases}\end{split}\]
\[\begin{split}\begin{cases}
x_A(t) = x_0 - t \sqrt{g h_l} \\
x_B(t) = x_0 + 2t \sqrt{gh_l} \\
\end{cases}\end{split}\]
Warning
Analytic solution is available in DamBreak for single layer, 1D St Venant only
Layers:¶
Number of layers is set to one for comparison with analytical solutions.
NL = 1
layer = fk.Layer(NL, triangular_mesh, topography=dambreak_analytic.topography)
Primitives:¶
def h_0(x, y):
h = H_R
if x < X_D: h = H_L
return h
if not args.nographics:
fk_plt.plot_init_1d(triangular_mesh, h_0, 'x', '$H_0$')
primitives = fk.Primitive(triangular_mesh, layer,
height_funct=h_0,
Uinit=0.,
Vinit=0.)
Boundary conditions:¶
fluvial_heights = [fk.FluvialHeight(ref=1, height=H_L),
fk.FluvialHeight(ref=2, height=H_R)]
slides = [fk.Slide(ref=3), fk.Slide(ref=4)]
Problem definition:¶
problem = fk.Problem(simutime, triangular_mesh,
layer, primitives,
fluvial_heights=fluvial_heights,
slides=slides,
numerical_parameters={'space_second_order':False},
analytic_sol=dambreak_analytic,
custom_funct={'plot':fk_plt.plot_freesurface_3d_analytic},
custom_funct_scheduler=create_figure_scheduler)
===================================================================
| INITIALIZATION |
===================================================================
Problem size: Nodes=704, Layers=1, Triangles=1200,
Iter = 0 ; Dt = 0.0000s ; Time = 0.00s ; ETA = 0.00s
Problem solving:¶
problem.solve()
if not args.nographics:
plt.show()
===================================================================
| TIME LOOP |
===================================================================
Iter = 6 ; Dt = 0.0404s ; Time = 0.24s ; ETA = 0.20s
Iter = 11 ; Dt = 0.0404s ; Time = 0.44s ; ETA = 0.19s
Iter = 16 ; Dt = 0.0404s ; Time = 0.65s ; ETA = 0.19s
Iter = 21 ; Dt = 0.0404s ; Time = 0.85s ; ETA = 0.18s
Iter = 26 ; Dt = 0.0404s ; Time = 1.05s ; ETA = 0.18s
Iter = 31 ; Dt = 0.0404s ; Time = 1.25s ; ETA = 0.17s
Iter = 36 ; Dt = 0.0404s ; Time = 1.46s ; ETA = 0.27s
Iter = 41 ; Dt = 0.0404s ; Time = 1.66s ; ETA = 0.17s
Iter = 46 ; Dt = 0.0404s ; Time = 1.86s ; ETA = 0.16s
Iter = 51 ; Dt = 0.0404s ; Time = 2.06s ; ETA = 0.16s
Iter = 56 ; Dt = 0.0404s ; Time = 2.27s ; ETA = 0.15s
Iter = 61 ; Dt = 0.0404s ; Time = 2.47s ; ETA = 0.15s
Iter = 66 ; Dt = 0.0404s ; Time = 2.67s ; ETA = 0.14s
Iter = 71 ; Dt = 0.0404s ; Time = 2.87s ; ETA = 0.14s
Iter = 76 ; Dt = 0.0404s ; Time = 3.07s ; ETA = 0.14s
Iter = 81 ; Dt = 0.0404s ; Time = 3.28s ; ETA = 0.13s
Iter = 86 ; Dt = 0.0404s ; Time = 3.48s ; ETA = 0.13s
Iter = 91 ; Dt = 0.0404s ; Time = 3.68s ; ETA = 0.12s
Iter = 96 ; Dt = 0.0404s ; Time = 3.88s ; ETA = 0.12s
Iter = 101 ; Dt = 0.0404s ; Time = 4.09s ; ETA = 0.12s
Iter = 106 ; Dt = 0.0404s ; Time = 4.29s ; ETA = 0.11s
Iter = 111 ; Dt = 0.0404s ; Time = 4.49s ; ETA = 0.11s
Iter = 117 ; Dt = 0.0404s ; Time = 4.73s ; ETA = 0.11s
Iter = 122 ; Dt = 0.0404s ; Time = 4.93s ; ETA = 0.10s
===================================================================
| L2 errors at time 5 s:
| H : 0.000136
| QX : 0.000024
| QY : 0.000001
===================================================================
Iter = 127 ; Dt = 0.0404s ; Time = 5.14s ; ETA = 0.10s
Iter = 132 ; Dt = 0.0404s ; Time = 5.34s ; ETA = 0.09s
Iter = 137 ; Dt = 0.0404s ; Time = 5.54s ; ETA = 0.09s
Iter = 142 ; Dt = 0.0404s ; Time = 5.74s ; ETA = 0.08s
Iter = 147 ; Dt = 0.0404s ; Time = 5.95s ; ETA = 0.08s
Iter = 152 ; Dt = 0.0404s ; Time = 6.15s ; ETA = 0.08s
Iter = 157 ; Dt = 0.0404s ; Time = 6.35s ; ETA = 0.07s
Iter = 162 ; Dt = 0.0404s ; Time = 6.55s ; ETA = 0.07s
Iter = 167 ; Dt = 0.0404s ; Time = 6.76s ; ETA = 0.06s
Iter = 172 ; Dt = 0.0404s ; Time = 6.96s ; ETA = 0.06s
Iter = 177 ; Dt = 0.0403s ; Time = 7.16s ; ETA = 0.06s
Iter = 182 ; Dt = 0.0403s ; Time = 7.36s ; ETA = 0.05s
Iter = 187 ; Dt = 0.0402s ; Time = 7.56s ; ETA = 0.06s
Iter = 192 ; Dt = 0.0401s ; Time = 7.76s ; ETA = 0.05s
Iter = 197 ; Dt = 0.0400s ; Time = 7.96s ; ETA = 0.04s
Iter = 203 ; Dt = 0.0400s ; Time = 8.20s ; ETA = 0.04s
Iter = 208 ; Dt = 0.0399s ; Time = 8.40s ; ETA = 0.03s
Iter = 213 ; Dt = 0.0398s ; Time = 8.60s ; ETA = 0.03s
Iter = 218 ; Dt = 0.0398s ; Time = 8.80s ; ETA = 0.02s
Iter = 223 ; Dt = 0.0397s ; Time = 9.00s ; ETA = 0.02s
===================================================================
| L2 errors at time 9 s:
| H : 0.000130
| QX : 0.000024
| QY : 0.000000
===================================================================
Iter = 228 ; Dt = 0.0397s ; Time = 9.20s ; ETA = 0.02s
Iter = 233 ; Dt = 0.0396s ; Time = 9.40s ; ETA = 0.01s
Iter = 238 ; Dt = 0.0395s ; Time = 9.59s ; ETA = 0.01s
Iter = 244 ; Dt = 0.0395s ; Time = 9.83s ; ETA = 0.00s
Iter = 249 ; Dt = 0.0115s ; Time = 10.00s ; ETA = 0.00s
===================================================================
| END |
===================================================================
Problem.solve() completed in 0.4971175193786621s (wall time)
Total running time of the script: ( 0 minutes 1.483 seconds)