{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "9f256eca", "metadata": {}, "outputs": [], "source": [ "from itertools import chain\n", "\n", "import pystencils as ps\n", "from lbmpy.session import *\n", "from lbmpy.methods.default_moment_sets import cascaded_moment_sets_literature\n", "from lbmpy.moments import moments_up_to_component_order, exponent_tuple_sort_key, MOMENT_SYMBOLS\n", "from lbmpy.methods import DensityVelocityComputation, create_from_equilibrium, CollisionSpaceInfo\n", "from lbmpy.equilibrium import GenericDiscreteEquilibrium, ContinuousHydrodynamicMaxwellian\n", "\n", "from pystencils.sympyextensions import scalar_product as dot" ] }, { "cell_type": "markdown", "id": "91cca442", "metadata": {}, "source": [ "# Demo: Lattice Boltzmann Methods for the Shallow Water Equations\n", "\n", "In this notebook we present how *lbmpy*'s low-level method definition interface may be used to rapidly set up two lattice Boltzmann models for solving the shallow water equations (SWE). The shallow water equations approximate the Navier-Stokes equations in a regime where the vertical length scale of the flow is much smaller than the horizontal length scales. Velocity is thus averaged across the third dimension, and density is replaced by water column height $h$.\n", "\n", "We reconstruct the central moment-based method presented by De Rosis in https://doi.org/10.1016/j.cma.2017.03.001 and the cumulant-based method of Venturi et al. (https://doi.org/10.1016/j.advwatres.2019.103474)." ] }, { "cell_type": "code", "execution_count": 2, "id": "4187668f", "metadata": {}, "outputs": [], "source": [ "d2q9 = LBStencil(Stencil.D2Q9)\n", "\n", "# gravitational acceleration\n", "g = sp.Symbol('g')\n", "\n", "# water column height\n", "h = sp.Symbol('h')\n", "\n", "# velocity\n", "u = sp.symbols(f'u_:{d2q9.D}')\n", "\n", "# relaxation rate\n", "ω_s = sp.Symbol('omega_s')\n", "\n", "# moment variables\n", "x, y, _ = MOMENT_SYMBOLS" ] }, { "cell_type": "markdown", "id": "56c221b3", "metadata": {}, "source": [ "## Discrete Shallow Water Equilibrium\n", "\n", "The central moment-based method makes use of the discrete shallow water LB equilibrium, which was originally defined by Zhou et al. in https://doi.org/10.1016/s0045-7825(02)00291-8. We implement it by listing equations for its components, and creating an instance of `lbmpy.equilibrium.GenericDiscreteEquilibrium`:" ] }, { "cell_type": "code", "execution_count": 3, "id": "0de48457", "metadata": {}, "outputs": [], "source": [ "def f_eq(ξ):\n", " ξ_sum = sum(abs(ξ_i) for ξ_i in ξ)\n", " if ξ_sum == 0:\n", " return h * (1 - (5 * g * h) / 6 - (2 * dot(u,u)) / 3)\n", " else:\n", " λ = 1 if ξ_sum == 1 else sp.Rational(1, 4)\n", " common = (g * h) / 6 + dot(ξ, u) / 3 + dot(ξ, u)**2 / 2 - dot(u,u) / 6\n", " return λ * h * common\n", "\n", "eq_populations = [f_eq(ξ) for ξ in d2q9]\n", "discrete_swe_eq = GenericDiscreteEquilibrium(d2q9, eq_populations, h, u)" ] }, { "cell_type": "code", "execution_count": 4, "id": "7a2dfa7c", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Instance of GenericDiscreteEquilibrium\n", "
Discrete Populations:
$f_0 = h \\left(- \\frac{5 g h}{6} - \\frac{2 u_{0}^{2}}{3} - \\frac{2 u_{1}^{2}}{3} + 1\\right)$
$f_1 = h \\left(\\frac{g h}{6} - \\frac{u_{0}^{2}}{6} + \\frac{u_{1}^{2}}{3} + \\frac{u_{1}}{3}\\right)$
$f_2 = h \\left(\\frac{g h}{6} - \\frac{u_{0}^{2}}{6} + \\frac{u_{1}^{2}}{3} - \\frac{u_{1}}{3}\\right)$
$f_3 = h \\left(\\frac{g h}{6} + \\frac{u_{0}^{2}}{3} - \\frac{u_{0}}{3} - \\frac{u_{1}^{2}}{6}\\right)$
$f_4 = h \\left(\\frac{g h}{6} + \\frac{u_{0}^{2}}{3} + \\frac{u_{0}}{3} - \\frac{u_{1}^{2}}{6}\\right)$
$f_5 = \\frac{h \\left(\\frac{g h}{6} - \\frac{u_{0}^{2}}{6} - \\frac{u_{0}}{3} - \\frac{u_{1}^{2}}{6} + \\frac{u_{1}}{3} + \\frac{\\left(- u_{0} + u_{1}\\right)^{2}}{2}\\right)}{4}$
$f_6 = \\frac{h \\left(\\frac{g h}{6} - \\frac{u_{0}^{2}}{6} + \\frac{u_{0}}{3} - \\frac{u_{1}^{2}}{6} + \\frac{u_{1}}{3} + \\frac{\\left(u_{0} + u_{1}\\right)^{2}}{2}\\right)}{4}$
$f_7 = \\frac{h \\left(\\frac{g h}{6} - \\frac{u_{0}^{2}}{6} - \\frac{u_{0}}{3} - \\frac{u_{1}^{2}}{6} - \\frac{u_{1}}{3} + \\frac{\\left(- u_{0} - u_{1}\\right)^{2}}{2}\\right)}{4}$
$f_8 = \\frac{h \\left(\\frac{g h}{6} - \\frac{u_{0}^{2}}{6} + \\frac{u_{0}}{3} - \\frac{u_{1}^{2}}{6} - \\frac{u_{1}}{3} + \\frac{\\left(u_{0} - u_{1}\\right)^{2}}{2}\\right)}{4}$
" ], "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "discrete_swe_eq" ] }, { "cell_type": "markdown", "id": "5b748055", "metadata": {}, "source": [ "## Continuous Shallow Water Equilibrium\n", "\n", "The cumulant-based method uses a special form of the continuous Maxwellian equilibrium, which is adapted to the shallow water equations by redefining its squared speed of sound parameter to be $c_s^2 = gh / 2$:" ] }, { "cell_type": "code", "execution_count": 5, "id": "c5fe3ee1", "metadata": {}, "outputs": [], "source": [ "c_s_sq = g * h / 2\n", "cont_sqe_eq = ContinuousHydrodynamicMaxwellian(dim=d2q9.D, rho=h, c_s_sq=c_s_sq)" ] }, { "cell_type": "code", "execution_count": 6, "id": "aa347620", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Continuous Hydrodynamic Maxwellian Equilibrium\n", " \n", " $f (h, \\left( u_{0}, \\ u_{1}\\right), \\left( v_{0}, \\ v_{1}\\right)) \n", " = \\frac{e^{\\frac{- \\left(- u_{0} + v_{0}\\right)^{2} - \\left(- u_{1} + v_{1}\\right)^{2}}{g h}}}{\\pi g}$\n", "
Compressible: ✓Deviation Only: ✗Order: ∞
\n", " " ], "text/plain": [ "ContinuousHydrodynamicMaxwellian(2D, compressible=True, deviation_only:Falseorder=None)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cont_sqe_eq" ] }, { "cell_type": "markdown", "id": "906652aa", "metadata": {}, "source": [ "## Method Definition\n", "\n", "The remaining parts of the method definition are nearly identical between both methods.\n", "Merely the collision spaces must be specified separately: `CollisionSpace.CENTRAL_MOMENTS` vs. `CollisionSpace.CUMULANTS`.\n", "\n", "The on-line version of this notebook runs the following simulation using the central moment-based method.\n", "To activate the cumulant-based method instead, flip the switch in the next cell:" ] }, { "cell_type": "code", "execution_count": 7, "id": "93f6b569", "metadata": {}, "outputs": [], "source": [ "# toggle central moment / cumulant method\n", "\n", "cumulant_method = False" ] }, { "cell_type": "code", "execution_count": 8, "id": "07a9b7c0", "metadata": {}, "outputs": [], "source": [ "polys = (sp.Integer(1), x, y, x * y, x**2 - y**2, x**2 + y**2, x * y**2, y * x**2, x**2 * y**2)\n", "r_rates = [0, 0, 0, ω_s, ω_s, 1, 1, 1, 1]\n", "rr_dict = dict(zip(polys, r_rates))\n", "\n", "cqc = DensityVelocityComputation(d2q9, True, False,\n", " density_symbol=h,\n", " density_deviation_symbol=sp.Symbol('delta_h'))\n", "\n", "if cumulant_method:\n", " swe_eq = cont_sqe_eq\n", " cspace = CollisionSpaceInfo(CollisionSpace.CUMULANTS)\n", "else:\n", " swe_eq = discrete_swe_eq\n", " cspace = CollisionSpaceInfo(CollisionSpace.CENTRAL_MOMENTS)\n", "\n", "swe_method = create_from_equilibrium(d2q9, swe_eq, cqc, rr_dict, cspace)" ] }, { "cell_type": "code", "execution_count": 9, "id": "56279a42", "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Central-Moment-Based Method\n", " Stencil: D2Q9Zero-Centered Storage: ✗Force Model: None
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Instance of GenericDiscreteEquilibrium\n", "
Discrete Populations:
$f_0 = h \\left(- \\frac{5 g h}{6} - \\frac{2 u_{0}^{2}}{3} - \\frac{2 u_{1}^{2}}{3} + 1\\right)$
$f_1 = h \\left(\\frac{g h}{6} - \\frac{u_{0}^{2}}{6} + \\frac{u_{1}^{2}}{3} + \\frac{u_{1}}{3}\\right)$
$f_2 = h \\left(\\frac{g h}{6} - \\frac{u_{0}^{2}}{6} + \\frac{u_{1}^{2}}{3} - \\frac{u_{1}}{3}\\right)$
$f_3 = h \\left(\\frac{g h}{6} + \\frac{u_{0}^{2}}{3} - \\frac{u_{0}}{3} - \\frac{u_{1}^{2}}{6}\\right)$
$f_4 = h \\left(\\frac{g h}{6} + \\frac{u_{0}^{2}}{3} + \\frac{u_{0}}{3} - \\frac{u_{1}^{2}}{6}\\right)$
$f_5 = \\frac{h \\left(\\frac{g h}{6} - \\frac{u_{0}^{2}}{6} - \\frac{u_{0}}{3} - \\frac{u_{1}^{2}}{6} + \\frac{u_{1}}{3} + \\frac{\\left(- u_{0} + u_{1}\\right)^{2}}{2}\\right)}{4}$
$f_6 = \\frac{h \\left(\\frac{g h}{6} - \\frac{u_{0}^{2}}{6} + \\frac{u_{0}}{3} - \\frac{u_{1}^{2}}{6} + \\frac{u_{1}}{3} + \\frac{\\left(u_{0} + u_{1}\\right)^{2}}{2}\\right)}{4}$
$f_7 = \\frac{h \\left(\\frac{g h}{6} - \\frac{u_{0}^{2}}{6} - \\frac{u_{0}}{3} - \\frac{u_{1}^{2}}{6} - \\frac{u_{1}}{3} + \\frac{\\left(- u_{0} - u_{1}\\right)^{2}}{2}\\right)}{4}$
$f_8 = \\frac{h \\left(\\frac{g h}{6} - \\frac{u_{0}^{2}}{6} + \\frac{u_{0}}{3} - \\frac{u_{1}^{2}}{6} - \\frac{u_{1}}{3} + \\frac{\\left(u_{0} - u_{1}\\right)^{2}}{2}\\right)}{4}$
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "
Relaxation Info
Central MomentEq. Value Relaxation Rate
$1$$h$$0$
$x$$0$$0$
$y$$0$$0$
$x y$$0$$\\omega_{s}$
$x^{2} - y^{2}$$0$$\\omega_{s}$
$x^{2} + y^{2}$$g h^{2}$$1$
$x y^{2}$$- \\frac{g h^{2} u_{0}}{2} - h u_{0} u_{1}^{2} + \\frac{h u_{0}}{3}$$1$
$x^{2} y$$- \\frac{g h^{2} u_{1}}{2} - h u_{0}^{2} u_{1} + \\frac{h u_{1}}{3}$$1$
$x^{2} y^{2}$$\\frac{g h^{2} u_{0}^{2}}{2} + \\frac{g h^{2} u_{1}^{2}}{2} + \\frac{g h^{2}}{6} + 3 h u_{0}^{2} u_{1}^{2} - \\frac{h u_{0}^{2}}{3} - \\frac{h u_{1}^{2}}{3}$$1$
" ], "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "swe_method" ] }, { "cell_type": "markdown", "id": "bda115b1", "metadata": {}, "source": [ "## Circular Dam Break Scenario\n", "\n", "We use the method derived above to simulate a circular dam break scenario.\n", "We discretize a square domain of $40 \\times 40$ meters using $100 \\times 100$ lattice cells,\n", "with periodic boundary conditions.\n", "A water column of $2.5$ meters heights, and $2.5$ meters radius, is placed in the middle.\n", "Otherwise, water depth is uniformly $0.5$ meters.\n", "We neglect any friction due to the bed.\n", "We assume default Earth gravitational pull of $9.81 \\mathrm{m}/\\mathrm{s}$.\n", "The simulated time step length is $0.05$ seconds per time step." ] }, { "cell_type": "code", "execution_count": 10, "id": "9623527e", "metadata": {}, "outputs": [], "source": [ "L_lattice = 100 # LU\n", "L_phys = 40 # meters\n", "L_factor = L_phys / L_lattice # meters / lu\n", "\n", "dt = 0.05 # seconds / step\n", "g_phys = 9.81 # meters / second**2\n", "\n", "g_lattice = g_phys / L_factor * dt**2 # lu / step**2\n", "shear_rr = 1.1" ] }, { "cell_type": "code", "execution_count": 11, "id": "f42339a4", "metadata": {}, "outputs": [], "source": [ "dh = ps.create_data_handling((L_lattice, L_lattice), periodicity=True,\n", " default_target=ps.Target.CPU)\n", "h_field = dh.add_array('h', values_per_cell=1)\n", "u_field = dh.add_array('u', values_per_cell=d2q9.D)\n", "pdf_field = dh.add_array('f', values_per_cell=d2q9.Q)\n", "pdf_tmp_field = dh.add_array('f_tmp', values_per_cell=d2q9.Q)" ] }, { "cell_type": "code", "execution_count": 12, "id": "5fbf5896", "metadata": {}, "outputs": [], "source": [ "from lbmpy.macroscopic_value_kernels import (\n", " macroscopic_values_setter, macroscopic_values_getter)\n", "\n", "setter = macroscopic_values_setter(swe_method,\n", " h_field.center, u_field.center_vector,\n", " pdf_field,\n", " set_pre_collision_pdfs=True)\n", "\n", "setter_kernel = ps.create_kernel(setter).compile()\n", "\n", "getter = macroscopic_values_getter(swe_method,\n", " h_field.center, u_field.center_vector,\n", " pdf_field,\n", " use_pre_collision_pdfs=True)\n", "getter_kernel = ps.create_kernel(getter).compile()\n", "\n", "periodic_sync = dh.synchronization_function(pdf_field.name)" ] }, { "cell_type": "code", "execution_count": 13, "id": "1432d1b4", "metadata": {}, "outputs": [], "source": [ "lbm_config = LBMConfig(lb_method=swe_method)\n", "opt = LBMOptimisation(pre_simplification=True, simplification=True,\n", " symbolic_field=pdf_field, symbolic_temporary_field=pdf_tmp_field)\n", "lb_kernel = create_lb_function(lbm_config=lbm_config, lbm_optimisation=opt)" ] }, { "cell_type": "markdown", "id": "174db9f2", "metadata": {}, "source": [ "## Initialization of PDF field" ] }, { "cell_type": "code", "execution_count": 14, "id": "6938e75c", "metadata": {}, "outputs": [], "source": [ "def init():\n", " dh.fill(u_field.name, 0.0)\n", "\n", " h_arr = dh.cpu_arrays[h.name][1:-1,1:-1]\n", "\n", " x_c = y_c = L_phys / 2\n", " R = 2.5\n", " H_phys = 2.5\n", " b_phys = 0.5\n", "\n", " # set up column\n", " for i in range(100):\n", " for j in range(100):\n", " x = L_factor * (i + 0.5)\n", " y = L_factor * (j + 0.5)\n", "\n", " if (x - x_c)**2 + (y - y_c)**2 <= R**2:\n", " h_arr[i,j] = H_phys / L_factor\n", " else:\n", " h_arr[i,j] = b_phys / L_factor\n", " \n", " dh.run_kernel(setter_kernel, g=g_lattice)" ] }, { "cell_type": "markdown", "id": "b08abe59", "metadata": {}, "source": [ "## Time Step" ] }, { "cell_type": "code", "execution_count": 15, "id": "056def94", "metadata": {}, "outputs": [], "source": [ "def step(output=False):\n", " periodic_sync()\n", " if output:\n", " dh.run_kernel(getter_kernel)\n", " dh.run_kernel(lb_kernel, g=g_lattice, omega_s=shear_rr)\n", " dh.swap(pdf_field.name, pdf_tmp_field.name)\n", " \n", " h_arr = dh.cpu_arrays[h.name][1:-1,1:-1]\n", " return L_factor * h_arr" ] }, { "cell_type": "markdown", "id": "27666725", "metadata": {}, "source": [ "## Simulation" ] }, { "cell_type": "code", "execution_count": 17, "id": "2ed47ba7", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "frames = 200\n", "timesteps = 100\n", "steps_per_frame = (timesteps//frames) + 1\n", "\n", "def run():\n", " for _ in range(steps_per_frame - 1):\n", " step()\n", " return step(output=True)\n", "\n", "init()\n", "\n", "animation = plt.surface_plot_animation(run, frames=frames, zlim=(0, 3))\n", "res = display_animation(animation)\n", "res" ] }, { "cell_type": "code", "execution_count": null, "id": "b5dca013", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": ".venv (3.12.12)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.12" } }, "nbformat": 4, "nbformat_minor": 5 }