.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples_tutorials/example_initialization.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_tutorials_example_initialization.py: ================================================================================ Initialization ================================================================================ Simple 2-layers case with initialization examples for topography, primitives and tracer. .. GENERATED FROM PYTHON SOURCE LINES 9-14 .. code-block:: default import numpy as np import matplotlib.pylab as plt import freshkiss3d as fk .. GENERATED FROM PYTHON SOURCE LINES 15-16 functions used to plot examples: .. GENERATED FROM PYTHON SOURCE LINES 16-69 .. code-block:: default def PLOT_VELOCITY_TOP(triangular_mesh, primitives): x = np.asarray(triangular_mesh.triangulation.x) y = np.asarray(triangular_mesh.triangulation.y) trivtx = np.asarray(triangular_mesh.triangulation.trivtx) plt.rcParams["figure.figsize"] = [10, 4] fig = plt.figure() #Subplot 1 ax1 = fig.add_subplot(121) ax1.triplot(x, y, trivtx, color='k', lw=0.5) im1 = ax1.tricontourf(x, y, trivtx, primitives.U[:, 1], 30, cmap=plt.cm.jet) fig.colorbar(im1, ax=ax1) ax1.set_title("Velocity x component on layer {}".format(1)) #Subplot 2 ax2 = fig.add_subplot(122) ax2.triplot(x, y, trivtx, color='k', lw=0.5) im2 = ax2.tricontourf(x, y, trivtx, primitives.V[:, 1], 30, cmap=plt.cm.jet) fig.colorbar(im2, ax=ax2) ax2.set_title("Velocity y component on layer {}".format(1)) plt.show() def PLOT_U_LAYERS(triangular_mesh, primitives): x = np.asarray(triangular_mesh.triangulation.x) y = np.asarray(triangular_mesh.triangulation.y) trivtx = np.asarray(triangular_mesh.triangulation.trivtx) plt.rcParams["figure.figsize"] = [10, 8] fig = plt.figure() #Subplot 1 ax1 = fig.add_subplot(221) ax1.triplot(x, y, trivtx, color='k', lw=0.5) im1 = ax1.tricontourf(x, y, trivtx, primitives.U[:, 0], 30, cmap=plt.cm.jet) fig.colorbar(im1, ax=ax1) ax1.set_title("Velocity x component on bottom layer") #Subplot 2 ax2 = fig.add_subplot(222) ax2.triplot(x, y, trivtx, color='k', lw=0.5) im2 = ax2.tricontourf(x, y, trivtx, primitives.U[:, 1], 30, cmap=plt.cm.jet) fig.colorbar(im2, ax=ax2) ax2.set_title("Velocity x component on top layer") #Subplot 3 ax3 = fig.add_subplot(223) ax3.triplot(x, y, trivtx, color='k', lw=0.5) im3 = ax3.tricontourf(x, y, trivtx, primitives.V[:, 0], 30, cmap=plt.cm.jet) fig.colorbar(im3, ax=ax3) ax3.set_title("Velocity y component on bottom layer") #Subplot 2 ax4 = fig.add_subplot(224) ax4.triplot(x, y, trivtx, color='k', lw=0.5) im4 = ax4.tricontourf(x, y, trivtx, primitives.V[:, 1], 30, cmap=plt.cm.jet) fig.colorbar(im4, ax=ax4) ax4.set_title("Velocity y component on top layer") plt.show() .. GENERATED FROM PYTHON SOURCE LINES 70-78 Time loop: -------------------- The simulation ends at **final_time** or when time iteration reach **time_iteration_max**. vtk_scheduler schedules saving of the solution in VTK file. Here 10 VTK files are set for the whole simulation time. .. GENERATED FROM PYTHON SOURCE LINES 79-85 .. code-block:: default simutime = fk.SimuTime(final_time=0.1, time_iteration_max=1000, second_order=True) vtk_scheduler = fk.schedules(count=10) .. GENERATED FROM PYTHON SOURCE LINES 86-88 Mesh: -------------------- .. GENERATED FROM PYTHON SOURCE LINES 89-92 .. code-block:: default triangular_mesh = fk.TriangularMesh.from_msh_file('../simulations/inputs/square.msh') .. GENERATED FROM PYTHON SOURCE LINES 93-100 Topography: -------------------- There is several ways to define topography, given the fact that it is a table of size [NC] conaining z_b values at each nodes of the mesh. One can provide a table ``topography`` from file (``np.loadtxt(str(topo_file_path))``) or user defined: .. GENERATED FROM PYTHON SOURCE LINES 101-104 .. code-block:: default topography = np.zeros(triangular_mesh.NV) .. GENERATED FROM PYTHON SOURCE LINES 105-106 .. seealso:: We could also provide a function **topography_functor** that returns z_b(x,y): .. GENERATED FROM PYTHON SOURCE LINES 107-114 .. code-block:: default def topo(x, y): x_0 = 5 y_0 = 5 topo = 0.05*(x-x_0) + 0.02*(y-y_0)**2 return topo .. GENERATED FROM PYTHON SOURCE LINES 115-120 Layers: -------------------- Based on toporaphy initialization, **topography** or **topography_functor** should be used as argument in layer call .. GENERATED FROM PYTHON SOURCE LINES 121-126 .. code-block:: default NL = 2 layer = fk.Layer(NL, triangular_mesh, topography=topography) #layer = fk.Layer(NL, triangular_mesh, topography_functor=topo) .. GENERATED FROM PYTHON SOURCE LINES 127-140 Primitives: -------------------- There is several ways to initialize these variables: - For H we can provide water **height** or **free_surface** - For velocities (default=0.) we can either provide **QXinit**, **QYinit** or **Uinit**, **Vinit**. .. note:: For each primitive you can either provide a constant 'Pinit' which will be set within the entire domain or specify a function of **Pinit = f(x,y)** (To call functions you need to add the suffix **_funct** to **Pinit**). If you want to set different initial values on each layer you can also provide a matrix that contains value on every node of the mesh: **Pinit(NC,NL)**. Since topography is flat in this example defining **height** or **free_surface** is the same. First we can set **U** and **V** to simple values: .. GENERATED FROM PYTHON SOURCE LINES 141-147 .. code-block:: default primitives = fk.Primitive(triangular_mesh, layer, free_surface=1., Uinit=0.3, Vinit=0.2) .. GENERATED FROM PYTHON SOURCE LINES 148-149 If we want to use a function to define **U** we can do: .. GENERATED FROM PYTHON SOURCE LINES 150-174 .. code-block:: default def U_0(x, y): r = np.sqrt((x-5.0)**2+(y-5.0)**2) if x < 5.: u = np.sinc(r) else: u = np.sinc(r + np.pi) return u def V_0(x, y): r = np.sqrt((x-5.0)**2+(y-5.0)**2) if y < 5.: v = np.sinc(r) else: v = np.sinc(r + np.pi) return v primitives = fk.Primitive(triangular_mesh, layer, free_surface=1, Uinit_funct=U_0, Vinit_funct=V_0) #Test plot: PLOT_VELOCITY_TOP(triangular_mesh,primitives) .. image-sg:: /auto_examples_tutorials/images/sphx_glr_example_initialization_001.png :alt: Velocity x component on layer 1, Velocity y component on layer 1 :srcset: /auto_examples_tutorials/images/sphx_glr_example_initialization_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 175-177 lets say we want to define specific value for each layer, we will provide the matrix **Uinit(NC,NL)** and **Vinit(NC,NL)**: .. GENERATED FROM PYTHON SOURCE LINES 178-200 .. code-block:: default NC = triangular_mesh.NV U_0 = np.zeros((NC, NL)) V_0 = np.zeros((NC, NL)) for C in range(NC): x = triangular_mesh.vertices[C].x y = triangular_mesh.vertices[C].y r = np.sqrt((x-5.0)**2+(y-5.0)**2) U_0[C, 0] = np.sinc(r) U_0[C, 1] = -np.sinc(r) V_0[C, 0] = -np.sinc(r) V_0[C, 1] = np.sinc(r) primitives = fk.Primitive(triangular_mesh, layer, free_surface=1, Uinit=U_0, Vinit=V_0) #Test plot: PLOT_U_LAYERS(triangular_mesh,primitives) .. image-sg:: /auto_examples_tutorials/images/sphx_glr_example_initialization_002.png :alt: Velocity x component on bottom layer, Velocity x component on top layer, Velocity y component on bottom layer, Velocity y component on top layer :srcset: /auto_examples_tutorials/images/sphx_glr_example_initialization_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 201-202 .. seealso:: We have shown how to define **U**, **V**. We also could have provided **QX** and **QY** instead. .. GENERATED FROM PYTHON SOURCE LINES 205-212 Tracer: -------------------- Like velocity there is 3 ways of defining **Tinit**: constant value, function (**Tinit_funct**) or table of size (NC,NL). Here is an example with a simple function to define **Tinit** and a 'salinity' compressible law: .. GENERATED FROM PYTHON SOURCE LINES 213-223 .. code-block:: default def T_0(x, y): if x < 5.: value = 30.*(x-5)**2+(y-5)**2 else: value = 0. return value tracer = fk.Tracer(triangular_mesh, layer, primitives, Tinit_funct=T_0) .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 1.123 seconds) .. _sphx_glr_download_auto_examples_tutorials_example_initialization.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_initialization.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_initialization.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_