.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples_simulations/example_topography_source.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_topography_source.py: ===================== Topography source ===================== This is a simple topography source case showing the ability of freshkiss3d to apply source terms on toporgaphy to simulate bed movement induced phenomena like tsunamis. .. GENERATED FROM PYTHON SOURCE LINES 10-150 .. code-block:: default import os, sys import argparse import numpy as np import random as rd import matplotlib.pylab as plt import freshkiss3d as fk import freshkiss3d.extra.plots as fk_plt from scipy.interpolate import griddata os.system('rm -r outputs') #sphinx_gallery_thumbnail_number = 5 parser = argparse.ArgumentParser() parser.add_argument('--nographics', action='store_true') args = parser.parse_args() def load_landslice(dir_path_z_variations): files_path_z_variations = os.listdir(dir_path_z_variations) times = [] xt_grid = [] yt_grid = [] zbt_grid = [] for ffile in files_path_z_variations: path_ffile = dir_path_z_variations+'/'+ffile x_grid = [] y_grid = [] zb_grid = [] with open(path_ffile, "r") as f: time = f.readline() time = float(time) times.append(time) lines = f.readlines() for line in lines: line = line.split(";") line = [float(i) for i in line] x_grid.append(line[0]) y_grid.append(line[1]) zb_grid.append(line[2]) x_grid = np.array(x_grid, dtype='d') y_grid = np.array(y_grid, dtype='d') zb_grid = np.array(zb_grid, dtype='d') xt_grid.append(x_grid) yt_grid.append(y_grid) zbt_grid.append(zb_grid) xt_grid = [i for _,i in sorted(zip(times, xt_grid))] yt_grid = [i for _,i in sorted(zip(times, yt_grid))] zbt_grid = [i for _,i in sorted(zip(times, zbt_grid))] times = sorted(times) return xt_grid, yt_grid, zbt_grid, times def msh_points_into_data_domain(x_msh, y_msh, max_x_grid, min_x_grid, \ max_y_grid, min_y_grid): x_msh_in_grid = [] y_msh_in_grid = [] C_msh_in_grid = [] C = 0 for x, y in zip(x_msh, y_msh): if min_x_grid < x < max_x_grid: if min_y_grid < y < max_y_grid: x_msh_in_grid.append(x) y_msh_in_grid.append(y) C_msh_in_grid.append(C) C += 1 return x_msh_in_grid, y_msh_in_grid, C_msh_in_grid def load_zb_init(x_msh, y_msh, topo_xy, x_msh_in_grid, y_msh_in_grid, \ zt_msh_in_grid): NC = len(x_msh) zb = np.empty(NC, dtype = 'd') for C in range(NC): zb[C] = topo_xy(x_msh[C], y_msh[C]) zb_0 = {} c = 0 for x_fk, y_fk in zip(x_msh, y_msh): zb_0[(x_fk, y_fk)] = zb[c] c += 1 c = 0 for x, y in zip(x_msh_in_grid, y_msh_in_grid): if (x, y) in zb_0: zb_0[(x, y)] += zt_msh_in_grid[0][c] c += 1 zb_0 = np.array(list(zb_0.values()),dtype='d') return zb_0 def interpolate_topography_data_on_freshkiss_mesh(x_msh, y_msh, xt_data, yt_data, \ zbt_data, interpolate_method='nearest'): if interpolate_method not in ('nearest', 'linear', 'cubic'): raise ValueError( "Choose between these interpolate methods :\n" "'nearest', 'linear', 'cubic'\n" ) zt_msh = [] if interpolate_method=='nearest': for x_data, y_data, zb_data in zip(xt_data, yt_data, zbt_data): zt_msh.append(griddata((x_data, y_data), \ zb_data, (x_msh, y_msh), \ method = 'nearest')) if interpolate_method=='linear': for x_data, y_data, zb_data in zip(xt_data, yt_data, zbt_data): zt_msh.append(griddata((x_data, y_data), \ zb_data, (x_msh, y_msh), \ method = 'linear')) if interpolate_method=='cubic': for x_data, y_data, zb_data in zip(xt_data, yt_data, zbt_data): zt_msh.append(griddata((x_data, y_data), \ zb_data, (x_msh, y_msh), \ method = 'cubic')) print("Using method {} to set data to Freshkiss3d.\n"\ .format(interpolate_method)) return zt_msh .. GENERATED FROM PYTHON SOURCE LINES 151-153 Time loop: -------------------- .. GENERATED FROM PYTHON SOURCE LINES 154-163 .. code-block:: default FINAL_TIME = 15. simutime = fk.SimuTime(final_time=FINAL_TIME, time_iteration_max=20000, \ second_order=True) create_figure_scheduler = fk.schedules(times=[0., 1., 5., 9., 12., 15.]) .. GENERATED FROM PYTHON SOURCE LINES 164-166 Mesh: -------------------- .. GENERATED FROM PYTHON SOURCE LINES 167-187 .. code-block:: default dir_path = os.path.abspath(os.path.dirname(sys.argv[0])) TG, vertex_labels, boundary_labels = fk.read_msh(dir_path + '/inputs/square2.mesh') x = np.asarray(TG.x) y = np.asarray(TG.y) x *= 10. y *= 10. triangular_mesh = fk.TriangularMesh(TG, vertex_labels, boundary_labels) if not args.nographics: fk_plt.plot_mesh(triangular_mesh) x_msh = np.asarray(TG.x) y_msh = np.asarray(TG.y) .. image-sg:: /auto_examples_simulations/images/sphx_glr_example_topography_source_001.png :alt: Mesh :srcset: /auto_examples_simulations/images/sphx_glr_example_topography_source_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 188-213 Topography source term: ------------------------ ``ExternalEffects`` objects are called in ``problem.solve()`` at each time step. In this example external effects are ``TopographySource``. There can be as much external effects as needed as long as they are stored in a dictionary type. Let us consider four cases: * In the first case we define continuous topography source over time to simulate a slowly moving bed. * In the second case we define sparse topography source over time to simulate sudden modification of the bed like a landslide. * In the third case we read several files which contain the topography coordinates (x, y, zb) * for several times for the whole domain. * In the fourth case we read several files which contain the topography coordinates (x, y, zb) * for several times for the moved part of the domain .. note:: Note that in both cases (continuous or sparse over time), the source can either be homogeneous or heterogeneous in space i.e. user can set either a constant set of sources with the ``sources`` parameter or a function (or set of function) of x, y with the ``sources_funct`` parameter. User can also set a vector of size (NC) containing the source term on each cell of the mesh with the ``sources_cell`` parameter. For given coordinates, we can interpolate the data on the mesh with different methods : nearest (neighbour), linear or cubic. .. GENERATED FROM PYTHON SOURCE LINES 214-313 .. code-block:: default CASE = 4 if CASE == 1: # The bed movement rate in set in m/s with ``source_unit='m/s'``. def bed_source_rate(x, y, x_0=90., y_0=100.): return .2*np.exp(-.005*((x-(x_0+5.))**2+(y-y_0)**2)) \ - .2*np.exp(-.005*((x-(x_0-5.))**2+(y-y_0)**2)) if not args.nographics: fk_plt.plot_init_1d(triangular_mesh, bed_source_rate, absc='x', ordo='Bed source rate $m/s$', title='Source term') continuously_moving_bed = fk.TopographySource( times=[1., 10.], sources_funct=bed_source_rate, source_unit='m/s') external_effects = {"bed": continuously_moving_bed} elif CASE == 2: # The landslide in set in m with ``source_unit='m'`` and occurs at :math:`t=1s` def landslide_source(x, y, x_0=90., y_0=100.): return 2.0*np.exp(-.005*((x-(x_0+5.))**2+(y-y_0)**2)) \ - 2.0*np.exp(-.005*((x-(x_0-5.))**2+(y-y_0)**2)) if not args.nographics: fk_plt.plot_init_1d(triangular_mesh, landslide_source, absc='x', ordo='Landslide source $m$', title='Source term') tsunami = fk.TopographySource( times=[1.], sources_funct=landslide_source, source_unit='m') external_effects = {"tsunami": tsunami} elif CASE == 3: dir_path_z_variations = dir_path + "/inputs/landslide" xt_grid, yt_grid, zbt_grid, times = \ load_landslice(dir_path_z_variations) zt_msh = interpolate_topography_data_on_freshkiss_mesh(\ x_msh, y_msh, xt_grid, yt_grid, zbt_grid) Dzt_msh = [] for i in range(len(zt_msh)-1): Dzt_msh.append(zt_msh[i+1] - zt_msh[i]) tsunami = fk.TopographySource( times = times[1::], sources_cell = Dzt_msh, source_unit = 'm') external_effects = {"tsunami": tsunami} elif CASE == 4: dir_path_z_variations = dir_path + "/inputs/landslide_part_domain" xt_grid, yt_grid, zbt_grid, times = \ load_landslice(dir_path_z_variations) max_x_grid = max(xt_grid[0]) min_x_grid = min(xt_grid[0]) max_y_grid = max(yt_grid[0]) min_y_grid = min(yt_grid[0]) x_msh_in_grid, y_msh_in_grid, C_msh_in_grid = \ msh_points_into_data_domain(x_msh, y_msh, max_x_grid, \ min_x_grid, max_y_grid, min_y_grid) zt_msh_in_grid = interpolate_topography_data_on_freshkiss_mesh(\ x_msh_in_grid, y_msh_in_grid, xt_grid, yt_grid, zbt_grid) Dzt_msh = [] for i in range(len(zt_msh_in_grid)-1): Dzt_msh.append(zt_msh_in_grid[i+1] - zt_msh_in_grid[i]) tsunami = fk.TopographySource( times = times[1::], sources_cell = Dzt_msh, cells = C_msh_in_grid, source_unit = 'm') external_effects = {"tsunami": tsunami} else: external_effects = None .. rst-class:: sphx-glr-script-out .. code-block:: none Using method nearest to set data to Freshkiss3d. .. GENERATED FROM PYTHON SOURCE LINES 314-316 Layers: -------------------- .. GENERATED FROM PYTHON SOURCE LINES 317-352 .. code-block:: default FREE_SURFACE = 8.0 X_1 = 10. ; X_2 = 100. def topo(x, y): if x < X_1: topo = 10.0 elif x > X_2: topo = 0.0 else: topo = -(10.0/(X_2-X_1))*(x - X_2) return topo NL = 2 if CASE == 3: NC = triangular_mesh.NV zb_0 = np.zeros(NC, dtype='d') for C in range(NC): zb_0[C] = topo(x_msh[C], y_msh[C]) + zt_msh[0][C] layer = fk.Layer(NL, triangular_mesh, topography=zb_0) elif CASE == 4: zb_0 = load_zb_init(x_msh, y_msh, topo, x_msh_in_grid, y_msh_in_grid, \ zt_msh_in_grid) layer = fk.Layer(NL, triangular_mesh, topography=zb_0) else: layer = fk.Layer(NL, triangular_mesh, topography_funct=topo) .. GENERATED FROM PYTHON SOURCE LINES 353-355 Primitives: -------------------- .. GENERATED FROM PYTHON SOURCE LINES 356-361 .. code-block:: default primitives = fk.Primitive(triangular_mesh, layer, free_surface=FREE_SURFACE) .. GENERATED FROM PYTHON SOURCE LINES 362-364 Boundary conditions: --------------------- .. GENERATED FROM PYTHON SOURCE LINES 365-370 .. code-block:: default slides = [fk.Slide(ref=r) for r in [1,2,3,4]] .. GENERATED FROM PYTHON SOURCE LINES 371-373 Writter: ---------------------- .. GENERATED FROM PYTHON SOURCE LINES 374-380 .. code-block:: default vtk_writer = fk.VTKWriter(triangular_mesh, scheduler=fk.schedules(count=10), scale_h=3.) .. GENERATED FROM PYTHON SOURCE LINES 381-383 Problem definition: -------------------- .. GENERATED FROM PYTHON SOURCE LINES 384-395 .. code-block:: default problem = fk.Problem(simutime, triangular_mesh, layer, primitives, slides=slides, external_effects=external_effects, vtk_writer=vtk_writer, custom_funct={'plot':fk_plt.plot_freesurface_topo_3d}, custom_funct_scheduler=create_figure_scheduler) .. rst-class:: sphx-glr-script-out .. code-block:: none =================================================================== | INITIALIZATION | =================================================================== Problem size: Nodes=3303, Layers=2, Triangles=6408, Iter = 0 ; Dt = 0.0000s ; Time = 0.00s ; ETA = 0.00s .. GENERATED FROM PYTHON SOURCE LINES 396-398 Problem solving: ----------------- .. GENERATED FROM PYTHON SOURCE LINES 399-403 .. 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_topography_source_002.png :alt: Free surface at time = 0.0s :srcset: /auto_examples_simulations/images/sphx_glr_example_topography_source_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples_simulations/images/sphx_glr_example_topography_source_003.png :alt: Free surface at time = 1.0s :srcset: /auto_examples_simulations/images/sphx_glr_example_topography_source_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples_simulations/images/sphx_glr_example_topography_source_004.png :alt: Free surface at time = 5.0s :srcset: /auto_examples_simulations/images/sphx_glr_example_topography_source_004.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples_simulations/images/sphx_glr_example_topography_source_005.png :alt: Free surface at time = 9.0s :srcset: /auto_examples_simulations/images/sphx_glr_example_topography_source_005.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples_simulations/images/sphx_glr_example_topography_source_006.png :alt: Free surface at time = 12.0s :srcset: /auto_examples_simulations/images/sphx_glr_example_topography_source_006.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none =================================================================== | TIME LOOP | =================================================================== Iter = 8 ; Dt = 0.0389s ; Time = 0.31s ; ETA = 3.14s Iter = 16 ; Dt = 0.0389s ; Time = 0.62s ; ETA = 3.07s Iter = 24 ; Dt = 0.0389s ; Time = 0.93s ; ETA = 3.01s Iter = 32 ; Dt = 0.0389s ; Time = 1.24s ; ETA = 2.93s Iter = 40 ; Dt = 0.0388s ; Time = 1.55s ; ETA = 2.85s Iter = 48 ; Dt = 0.0387s ; Time = 1.86s ; ETA = 2.84s Iter = 56 ; Dt = 0.0386s ; Time = 2.17s ; ETA = 2.81s Iter = 64 ; Dt = 0.0385s ; Time = 2.48s ; ETA = 2.76s Iter = 72 ; Dt = 0.0385s ; Time = 2.79s ; ETA = 2.64s Iter = 80 ; Dt = 0.0385s ; Time = 3.10s ; ETA = 2.58s Iter = 87 ; Dt = 0.0386s ; Time = 3.37s ; ETA = 2.53s Iter = 95 ; Dt = 0.0387s ; Time = 3.68s ; ETA = 2.45s Iter = 103 ; Dt = 0.0388s ; Time = 3.99s ; ETA = 2.35s Iter = 111 ; Dt = 0.0389s ; Time = 4.30s ; ETA = 2.34s Iter = 119 ; Dt = 0.0389s ; Time = 4.61s ; ETA = 2.27s Iter = 127 ; Dt = 0.0389s ; Time = 4.92s ; ETA = 2.15s Iter = 135 ; Dt = 0.0389s ; Time = 5.23s ; ETA = 2.07s Iter = 143 ; Dt = 0.0389s ; Time = 5.54s ; ETA = 2.02s Iter = 150 ; Dt = 0.0389s ; Time = 5.82s ; ETA = 1.98s Iter = 158 ; Dt = 0.0389s ; Time = 6.13s ; ETA = 1.91s Iter = 166 ; Dt = 0.0389s ; Time = 6.44s ; ETA = 1.84s Iter = 174 ; Dt = 0.0388s ; Time = 6.75s ; ETA = 1.80s Iter = 182 ; Dt = 0.0388s ; Time = 7.06s ; ETA = 1.70s Iter = 190 ; Dt = 0.0389s ; Time = 7.37s ; ETA = 1.64s Iter = 198 ; Dt = 0.0389s ; Time = 7.68s ; ETA = 1.58s Iter = 206 ; Dt = 0.0389s ; Time = 7.99s ; ETA = 1.51s Iter = 213 ; Dt = 0.0389s ; Time = 8.27s ; ETA = 1.46s Iter = 221 ; Dt = 0.0389s ; Time = 8.58s ; ETA = 1.41s Iter = 229 ; Dt = 0.0389s ; Time = 8.89s ; ETA = 1.33s Iter = 237 ; Dt = 0.0390s ; Time = 9.20s ; ETA = 1.23s Iter = 245 ; Dt = 0.0390s ; Time = 9.51s ; ETA = 1.16s Iter = 253 ; Dt = 0.0390s ; Time = 9.82s ; ETA = 1.13s Iter = 261 ; Dt = 0.0390s ; Time = 10.14s ; ETA = 1.04s Iter = 269 ; Dt = 0.0390s ; Time = 10.45s ; ETA = 0.96s Iter = 276 ; Dt = 0.0390s ; Time = 10.72s ; ETA = 0.90s Iter = 284 ; Dt = 0.0390s ; Time = 11.03s ; ETA = 0.83s Iter = 292 ; Dt = 0.0389s ; Time = 11.34s ; ETA = 0.79s Iter = 300 ; Dt = 0.0388s ; Time = 11.65s ; ETA = 0.72s Iter = 308 ; Dt = 0.0387s ; Time = 11.96s ; ETA = 0.65s Iter = 316 ; Dt = 0.0387s ; Time = 12.27s ; ETA = 0.58s Iter = 324 ; Dt = 0.0387s ; Time = 12.58s ; ETA = 0.52s Iter = 332 ; Dt = 0.0387s ; Time = 12.89s ; ETA = 0.45s Iter = 339 ; Dt = 0.0388s ; Time = 13.17s ; ETA = 0.39s Iter = 347 ; Dt = 0.0388s ; Time = 13.48s ; ETA = 0.33s Iter = 355 ; Dt = 0.0388s ; Time = 13.79s ; ETA = 0.26s Iter = 363 ; Dt = 0.0388s ; Time = 14.10s ; ETA = 0.19s Iter = 371 ; Dt = 0.0388s ; Time = 14.41s ; ETA = 0.13s Iter = 379 ; Dt = 0.0388s ; Time = 14.72s ; ETA = 0.06s Iter = 387 ; Dt = 0.0110s ; Time = 15.00s ; ETA = 0.00s =================================================================== | END | =================================================================== Problem.solve() completed in 6.336383581161499s (wall time) .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 10.386 seconds) .. _sphx_glr_download_auto_examples_simulations_example_topography_source.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_topography_source.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_topography_source.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_