.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples_simulations/example_dambreak.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_simulations_example_dambreak.py: ================================================================================ 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". .. GENERATED FROM PYTHON SOURCE LINES 14-27 .. code-block:: default 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 .. GENERATED FROM PYTHON SOURCE LINES 28-49 Riemann problems: -------------------- Two cases (Riemann problems) are considered in this example: - Dam break on a wet domain without friction (case = 'wet') .. math:: 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} - Dam break on a dry domain without friction (case = 'dry') .. math:: 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} .. GENERATED FROM PYTHON SOURCE LINES 50-63 .. code-block:: default 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 .. GENERATED FROM PYTHON SOURCE LINES 64-66 Time loop: -------------------- .. GENERATED FROM PYTHON SOURCE LINES 67-76 .. code-block:: default 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) .. GENERATED FROM PYTHON SOURCE LINES 77-83 Mesh: -------------------- .. math:: L = 10 m, l = 0.5 m .. GENERATED FROM PYTHON SOURCE LINES 84-91 .. code-block:: default 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) .. image-sg:: /auto_examples_simulations/images/sphx_glr_example_dambreak_001.png :alt: Mesh :srcset: /auto_examples_simulations/images/sphx_glr_example_dambreak_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 92-98 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. .. GENERATED FROM PYTHON SOURCE LINES 99-107 .. code-block:: default 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.) .. GENERATED FROM PYTHON SOURCE LINES 108-135 - In the 'wet' case analytical solution is given by: .. math:: 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} .. math:: 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} .. math:: \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} .. math:: 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) .. GENERATED FROM PYTHON SOURCE LINES 138-161 - If case is set to 'dry' analytical solution is given by: .. math:: 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} .. math:: 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} .. math:: \begin{cases} x_A(t) = x_0 - t \sqrt{g h_l} \\ x_B(t) = x_0 + 2t \sqrt{gh_l} \\ \end{cases} .. GENERATED FROM PYTHON SOURCE LINES 164-167 .. warning:: Analytic solution is available in DamBreak for single layer, 1D St Venant only .. GENERATED FROM PYTHON SOURCE LINES 170-175 Layers: -------------------- Number of layers is set to one for comparison with analytical solutions. .. GENERATED FROM PYTHON SOURCE LINES 176-180 .. code-block:: default NL = 1 layer = fk.Layer(NL, triangular_mesh, topography=dambreak_analytic.topography) .. GENERATED FROM PYTHON SOURCE LINES 181-183 Primitives: -------------------- .. GENERATED FROM PYTHON SOURCE LINES 184-198 .. code-block:: default 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.) .. image-sg:: /auto_examples_simulations/images/sphx_glr_example_dambreak_002.png :alt: Initial condition :srcset: /auto_examples_simulations/images/sphx_glr_example_dambreak_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 199-201 Boundary conditions: --------------------- .. GENERATED FROM PYTHON SOURCE LINES 202-207 .. code-block:: default 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)] .. GENERATED FROM PYTHON SOURCE LINES 208-210 Problem definition: -------------------- .. GENERATED FROM PYTHON SOURCE LINES 211-221 .. code-block:: default 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) .. rst-class:: sphx-glr-script-out .. code-block:: none =================================================================== | INITIALIZATION | =================================================================== Problem size: Nodes=704, Layers=1, Triangles=1200, Iter = 0 ; Dt = 0.0000s ; Time = 0.00s ; ETA = 0.00s .. GENERATED FROM PYTHON SOURCE LINES 222-224 Problem solving: ----------------- .. GENERATED FROM PYTHON SOURCE LINES 225-230 .. code-block:: default problem.solve() if not args.nographics: plt.show() .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples_simulations/images/sphx_glr_example_dambreak_003.png :alt: Free surface at time = 5.0s :srcset: /auto_examples_simulations/images/sphx_glr_example_dambreak_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples_simulations/images/sphx_glr_example_dambreak_004.png :alt: Free surface at time = 9.0s :srcset: /auto_examples_simulations/images/sphx_glr_example_dambreak_004.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none =================================================================== | 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) .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 1.483 seconds) .. _sphx_glr_download_auto_examples_simulations_example_dambreak.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_dambreak.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_dambreak.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_