{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "nbsphinx": "hidden" }, "outputs": [], "source": [ "from lbmpy.session import *\n", "import pystencils as ps\n", "from pystencils.stencils import visualize_stencil" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Demo: Create lbmpy Method from Scratch\n", "\n", "\n", "\n", "\n", "### Defining transformation to collision space" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA20AAAAUBAMAAADsGchtAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdt3NMolEIma7mVTvEKvunM/GAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFf0lEQVRoBe2az4scVRDHvzuzvTP7o+Oai5CDDkHZQ1AXET0IZs9eskEjshcH/4Esufjj4oJIIAjuQdBcpEFyMIdl8SKIkEXEk+KoeJGQXURQRIka9RCJa1W915np6qoHb4nx4jtkuqvep771Y6en0zM4/DDwAtzVXzddvQ3TzMZ8AmfNYAkNh0gVYmvcWiK/9IMR/YeGeBrobQLfPPNUq3li+75lZsOLNiF7s4jexfMrOGT+cZCGtcoLr3oEF+ItW8MhElkRIV6tk27WS+++oYFw7jbLIsTGRCFz61LbHsOXreaJ7StT71GbwIcjIIu4guJPlMuWCGlIPOXrbpdLDsGFzH+m9tMp22wNh0hkRYR4tUqyWeUenh1qItksiwg2bm+Y2zm6uO2g0M0LtmKxpUfbK5Mof7w6oqAZBH4GPgVO2hohnvK9BjxpE6BCrmzdUPtrm6XhEX5WTARvUybdrEMjLOw0AaSbZRHBxu0Nc1sCZvcwrwsOtlLrsXxn3SSAe2huWcRHwNUVXFY1RY0QT/neBi6NTAJUCKZ0GbXN0vAIPysmgreZVrpZCwMc+r0J0FmqWRYRbNxemVtvANyxh5m/VOBoe0SZ+fSMTcQ+5xDHRzy37pCCqkUa5tz+Bk5vozskr1pciD+37pC8anmEnxUTwdsMlW5W94Y3NzjNsohoI6Lgz7epCtitMPNrM5PatqbMfHrEJmKfswjgwRGKioKqRRrW3MrrNLdFFBV51eJC/LkV7FXLJdyshBBvM9RuaKBTOu3ttK8D/H5DFiFRiJC5zW4AL69i+o9mJrXtFWXm07dsIvY5i8AM6faX2yKkYc2tR9tPrJoEF+LPzdJwCTcrIcTL6Y1XbKBTOu07vTreHI9kblmERCFC5tbZpilUxtyC7b2WHnDMJmKfswh09uiGb9AWIQ1zbtdlbhbBhfhzyyLcrERDvM2MYwOd0mnvA839fCZzyyIkChEyt4WhfdXbDW/9E21BfGITMZUsAhcp/vQvbRHSsOYWr5MWwYX4c8si3KxEQ7yc3nilmwXMVeO99ZHMzWmWTUgUIm7OjT5W5437ErY5gS0iPTebkA8Mq6fO3ED3JZe2zUkvDKkj7v2kpeETFYXyCfY2VyzOaZbcETQBOkvP7Uhrf4xCGjI3fuvPbqKvPzij7TmDp2uYRcRUsogP0FvHzKAt4lwncR54fWQS6eukpeESblZCiLeZcbpZmK/wfBOgM5mb0yyTCFGIkLnxR21/B3PLKnC0OR+cFhFTySGmK0ytm3cZzn0J6P/d79h3Mu5dhrwHc+5L/KxYI3ib3Uo3C08AdzUBOpO5Oc0yiRCFCJlbUVGMY7i8Ul6jg4nFNuACWbRnjWwWEVLJIc6cWvuWruWrtgbH09qd7fLjQGC3+adWVJQWz0jZw7WTNbTHI8ZZWYR4tUayWb37Tm0NdCHJZlmE2GQgMje5gHx953fA0RGVNl5iA31f0PLwm9sitj5/fANZxPH9/d+A7oqtIfFUVuXJu2m3EJ37CRsvLmTup2s/QNnFlkVMZKVisYZ4lT20wyu9u7+/P9DtTTbLIoKNNWRu8rBHqp9fHzehPupRonSj3vTQc66wlF2M+QQ/f8nTCATeD1nU/y7VB8ouZtH4l4n80g9GFPJ9wLm63KI+mHgNT4mVp1/FHcou1nwiPFdWsZIa8Un06kSidHizEGWXXeG5svLcYiK/9IMRYW6z9bvni2Yb5Cx8K6M9b8ad2s7mfCJ8x6JjpTQCUS7GNOJLXYi2s/t2EAcoPb9ZTIS59TZj3RvxdfIldE976u80tZ3JfCL0W8dKaQSiz3ITqy5E23nL7SAOUHp+s5iguR2+9//fKUyMfuLw7MTxxOF//KMO/mXD1NHhPwrjTLx3gTLSAAAAAElFTkSuQmCC\n", "text/latex": [ "$$\\left [ \\left ( 0, \\quad 0\\right ), \\quad \\left ( 0, \\quad 1\\right ), \\quad \\left ( 0, \\quad 2\\right ), \\quad \\left ( 1, \\quad 0\\right ), \\quad \\left ( 1, \\quad 1\\right ), \\quad \\left ( 1, \\quad 2\\right ), \\quad \\left ( 2, \\quad 0\\right ), \\quad \\left ( 2, \\quad 1\\right ), \\quad \\left ( 2, \\quad 2\\right )\\right ]$$" ], "text/plain": [ "[(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2)]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from lbmpy.moments import moment_matrix, moments_up_to_component_order, exponents_to_polynomial_representations\n", "moment_exponents = list(moments_up_to_component_order(2, 2))\n", "moment_exponents" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcgAAAAcBAMAAAAXVzJUAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMkS7zRCZdiKJ71Rmq90icBAQAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADq0lEQVRYCd2YMWgTYRTHX9JemkvSWMS9cRIHaUDUqZhdCjcImUorjXRRzJSiIA26qdBM1TGLOLbo1g7e5CAIgkUQB1HQOoioWBQU4ve978p39937I3djbvjy+P/e/7333XdpkxBlvjbPBpk92OCvXAYQ98GES1lsI9ADyX5/uoVYDv0B7cku3AcTrmSxjeio3ASp9X7hN2I59M802xdtuA8mXMhijrpaLDTFHlCsN8o/IMwOPtBGQ3ThPphwIYs5OqXFyczvsOqBOFVe8SkaAPfBhIewWEV+qLTVzMPV+QnIbEOGfQRwH0y4lsU6+khUzn4sm2iqXLq/jGy4DyZcy2Id9Yi4xw56YqT+tWVJza09Q07cBxOuZTFHG0SVbfJev8uyyZN0HM2VQy8OijOyDffBhCtZzFE9pFm9v/UMmyxc2D0jT5VLPbb7RG6O+2DCE1hsolKTnmuQZZOV0ehbru3Ipq3RSAa4DyZcyWITTQ2po4GwyfYqlR867b1el+6I971989b1rpNNLIoeTLiGWM0h0nymo0m04xcP6LTWhE0OdmjKPTG/OqR1ruEs3qC2X2+JouQx6RLhEmI1l0jzxYwxXPhOC9qd3mR5u0MTTa5sl7ullrkpVjKRH/gH1b6jGlHyYMIlxGoOEeezxjhWH9A+aXd6kx49oiPuAxhUuvSLuzmLR5PuDSEyouTBhOuK1RwizmeNcVz+g06S1JM8O+PsheYCdfbilbohOotF0YMJVxerJYk4nz2XGFYnCR5XUk/mDS4bXzqk/lSJV/qGqDQWRQ8mXF2sliTifPZcYlh9aXqlvenHlSZe0jyXjS/zWhauarBOU4EDIlHwYMIlxGoukeaLGWNYHSofl9qk9zM5Y6lZ/JsSV2gu5DQne277PN1zs41I2iOmM8HV/kOk+aKO2smYS1BtqN5mRFe3FkL6mjwJ701PfTd2xMX71/rGmQSLa+0roZsdiewR0/9TjRuJRk2k+aKOFnMJ9Ymn0uCIKPUvQD3WabGDsrWeKmGSlQcTTgEYl9QucT4ux4vGfE12qTSIYj96jV72aKmhvqMkxOmWd/ijQBJEWZJoPJiwVcKmJiTSfLFhDWZhKaDCMELtWIoKv9A5tSbFejjRlbONmsw2mvFgwlkSxiWZSPMZSwxz+EKtbyMURq/Ry4nejIqSYnHt0mFSEkSqJBoPJmyVMANnAKPxKs2Xwiy8V+vFGBrHsKbfj9Xlcdya3dNtDh9bYRwj9AP9uO31H0gZOpOwIIh4AAAAAElFTkSuQmCC\n", "text/latex": [ "$$\\left ( 1, \\quad y, \\quad y^{2}, \\quad x, \\quad x y, \\quad x y^{2}, \\quad x^{2}, \\quad x^{2} y, \\quad x^{2} y^{2}\\right )$$" ], "text/plain": [ "⎛ 2 2 2 2 2 2⎞\n", "⎝1, y, y , x, x⋅y, x⋅y , x , x ⋅y, x ⋅y ⎠" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "moments = exponents_to_polynomial_representations(moment_exponents)\n", "moments" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAF5CAYAAABqT9akAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl0XmWB+PHvk3RNW8qSUspSoQUXBAb0yCIiKG4BfwPiyDAyqIgL6igj6oz83BiP4CAyqHUZF9zBUWdk9PAjIgdQUFxREXcQKhRa2hRaupckz++PN2+bpHmTd73r93MOp+TNvW+eP27eb+597hJijEiSNJGutAcgScouIyFJqmlaUj8o9A90Ay8ATgDmAaHGohFYB3w/9vXelNDwpJpC/8BM4FTgWGA2brvKidA/0AO8GDgamDXJohFYC9wc+3pvHfMeSc1JhP6By4DTG1zti7Gv94OdGI9Uj5E/bq4Cjmtw1S/Fvt5LOzAkqS4jf9x8GTiywVWXxb7ej1e/SORwU+gfWEzjgQA4J/QP7NHu8UgNOJrGAwHwj6F/YM92D0ZqwIk0HgiA14b+gdnVL5Kak/ibJtfrBg5v50CkBrWy7R7WzoFIDWp2250FPLn6RVKRmJnSulKr3HaVVzNaWHfHtpvYxPWELnzRAdxz1yy6uytf77H3IJ//xX2pjkmayjeX7c4t35zPintmcNwpG7joc6vSHpLUkO9dM4+vX7kXa1dNZ/5eg1zwkVU87aQtEy2abiQAznvfak573fq0hyHVba9Fg5x5wVruuGUO27fWOtNJyqaffLeHr3xwAe/4z4c47NitrHlo0g6kHwkpb577dxsB+POvZ7F2pb9DypdrLu/lZW9ZyxHHbwVg4QGDky2e/sV0V1/ey5kHL+Utz1vML26aPfUKkqSmDA3Cfb+fxfq13bzqqIM4+6lLuPIte7N1c8094nQjce571/CFO+7lq7+9lxeevY5LXr0/D9w9PdUxSVJRrV3VzdAg/Pj6eVx+3f18/Ja/ct/vZ/GlS/eqtUq6kTj8mVuZs1tkxqzIqec+xhOP3MJP+uekOiZJKqqZsytXT5967qMs2G+IPfYe4vTXP8Kvvl/zczf9w01jhIh3pZWkzpi/1zB77D1IqP98i/Qi8dgjXfz4+h62bQkMPg7f/co8/nRHD0c/f1NqY5LqMfg4bNsSGB6C4SF2bMNSHjznpeu57gu7s3ZVN+vXdvGdz+7B05+7sdbi6Z2ZMfh44CuX9fKhN8ykqyuy6KDtvPOzD3Lgof62Kdu+dMle/M8ndh7D/dF1u/HSN63lvIvXpjgqqT6vfPdaHnu0m9cddxDTZ0SO7dvAKy56pNbi6UViz4VDfPIH96f286VmnXexQVB+TZ8BFy5bzYXLVtezeMbmJCRJWZJUJFqZjXYmW2ly21UZ7dh2k4rE5pTWlVrVyvY34b1wpIS0ZdtNKhI/B4abWG8r8Ks2j0VqxE+bXG8bbrtKV7Pb7mPA76tfJBKJ2Nc7ACxrYtUPxb5e/xpTamJf713AN5tY9UOxr9e9YKXpJ0B/g+sMA5fEvt4d93NK7PGlAKF/4BDg2dT/jOvlCQ1NmlToHziCyjOue3DbVU6E/oEAHMXOZ1zX2naHgUepPON6xZj3SDISkqR88RRYqUkhhBtCCJemPQ6pk9yTkJoQQlgC/IHKyRULYozbUx6S1BHuSUjNecPIvwE4Lc2BSJ3knoTUoBDCdGANMH/kpZ/FGI9JcUhSx+R+TyKEsDCEcGMI4eUhhBlpj0el8LeM/d05IoRwcFqDUXmEEOaHEN4aQrh25I+Vjst9JKgcEz4BuAp4OITw/hDCwpTHpGJ7G5XTuKu62Xn4SWq7EMKTQwifA1YCHwIOByZ9NnW75D4SMcb1wJep3NF2d+AdwPIQwv+EEJ6R6uBUOCMT1keNe3k68Br3ZNVOIYSuEMKpIYQfAr8EXgnMpvKH8SUxobmC3EdixBXsrOqskf9OB74fQvidh6LURm9g4t8bJ7DVFtVDSsAK4L+A46nEYfSjHf4rsfEUZeJ6pLbH1/j2RioRWQZ8Isb4cGIDU2FMMGE9nhPYaloI4cnA24GXU7l6v2eCxbYDH4sxviOpcRVlTwLgUioxmMhcPBSl1o2fsB7PCWw1ZNQhpR9RuSFk9ZDSRIGAyu0zPpbU+KBYexJdVHbPFtWx+DCV43r3AcfFGDd0cmwqhhDC7cBxkyzyOLAsxvi2hIakHAshPBW4kcpJEHPrWCUC340xntLRgY1TmD2JGOMwcDn13UO9C5gJ7EGaj3BVbtSYsB7PCWw1IgC7UV8gADYBH+zccCZWmEiM+Dy173I42hDwMHBMjPHRzg5JBVFrwno8J7BVlxjjb4GTqX2YfLyHgR92bkQTK1QkRk6H/SqTnz88OhArJllOAnZMWL8WqGcPYR6VyUdpSjHGnwLPY+pQbCTB015HK1QkRlxB5djwRIapnJ1iINSI/0PticSJPCOEsLRTg1GxjApFrc8tqMxHJHba62iFOx4fY/xTCOGX7Ho67BCVK2P3oTJpLdXrN0z8dLqXj/x7zbjXt1PZW5XqFanMaU1kO/DpGGMqT+kszNlNo4UQTgG+zs4JoeohppOAP4+8tiDGOJD86FQUIYQI3BljPDLtsSi/QghHU3ke9Ubg+VTOeBo9mb0VeGKM8YEUhlfIw00A3wWqp7WOnoO4G5gz8vqaEEJvGoOTJNglELvFGH/C2DmKCNySViCgoJEYdTrsIOMmqWOMmzEUklI2QSAi7DKZvZ0UTnsdrZCRGPF54AtMMEltKCSlqVYgqkaF4hOkcNrraIWck6hXCKGHygUq4ByFGuSchJoxVSCypsh7ElNyj0JSkvIWCCh5JMBQSEpGHgMBRgIwFJI6K6+BACOxg6GQ1Al5DgQYiTEMhaR2ynsgwEjswlBIaociBAKMxIQMhaRWFCUQYCRqMhSSmlGkQICRmJShkNSIogUCjMSUDIWkehQxEGAk6mIoJE2mqIEAI1E3QyFpIkUOBBiJhhgKSaMVPRBgJBpmKCRBOQIBRqIphkIqt7IEAoxE0wyFVE5lCgQYiZYYCqlcyhYIMBItMxRSOZQxEGAk2sJQSMVW1kCAkWgbQyEVU5kDAUairQyFVCxlDwQYibYzFFIxGIgKI9EBhkLKNwOxk5HoEEMh5ZOBGMtIdJChkPLFQOzKSHSYoZDywUBMzEgkwFBI2WYgajMSCTEUUjYZiMkZiQQZCilbDMTUjETCDIWUDQaiPkYiBYZCSpeBqJ+RSImhkNJhIBpjJFJkKKRkGYjGGYmUGQopGQaiOUYiAwyF1FkGonlGIiMMhdQZBqI1RiJDDIXUXgaidUYiYwyF1B4Goj2MRAYZCqk1BqJ9jERGGQqpOQaivYxEhhkKqTEGov2MRMYZCqk+BqIzjEQOGAppcgaic4xEThgKaWIGorOMRI4YCmksA9F5RiJnDIVUYSCSYSRyyFCo7AxEcoxEThkKlZWBSJaRyDFDobIxEMkzEjlnKFQWBiIdRqIADIWKzkCkx0gUhKFQURmIdBmJAjEUKhoDkT4jUTCGQkVhILLBSBSQoVDeGYjsMBIFZSiUVwYiW4xEgRkK5Y2ByB4jUXCGQnlhILLJSJSAoVDWGYjsMhIlYSiUVQYi24xEiRgKZY2ByD4jUTKGQllhIPLBSJSQoVDaDER+GImSMhRKi4HIFyNRYoZCSTMQ+WMkSs5QKCkGIp+MhAyFOs5A5JeREGAo1DkGIt+MhHYwFGo3A5F/RkJjGAq1i4EoBiOhXRgKtcpAFIeR0IQMhZplIIrFSKgmQ6FGGYjiMRKalKFQvQxEMRkJTclQaCoGoriMhOpiKFSLgSg2I6G6GQqNZyCKz0ioIYZCVQaiHIyEGmYoZCDKw0ioKYaivAxEuRgJNc1QlI+BKB8joZYYivIwEOVkJNQyQ1F8BqK8jITawlAUl4EoNyOhtjEUxWMgZCTUVoaiOAyEwEioAwxF/hkIVRkJdYShyC8DodGMhDrGUOSPgdB4RkIdZSjyw0BoIkZCHWcoss9AqBYjoUQYiuwyEJqMkVBiDEX2GAhNxUgoUYYiOwyE6mEklDhDkT4DoXoZCaXCUKTHQKgRRkKpMRTJMxBqlJFQqgxFcgyEmmEklDpD0XkGQs0yEsoEQ9E5BkKtMBLKDEPRfgZCrTISyhRD0T4GQu1gJJQ5hqJ1BkLtYiSUSYaieQZC7WQklFmGonEGQu1mJJRphqJ+BkKdYCSUeYZiagZCnWIklAuGojYDoU4yEsoNQ7ErA6FOMxLKFUOxk4FQEoyEcsdQGAglx0gol8ocCgOhJBkJ5VYZQ2EglDQjoVwrUygMhNJgJJR7ZQiFgVBajIQKocihMBBKk5FQYRQxFAZCaTMSKpQihcJAKAuMhAqnCKEwEMoKI6FCynMoDISyxEiosPIYCgOhrDESKrQ8hcJAKIuMhAovD6EwEMoqI6FSyHIoDISyzEioNLIYCgOhrDMSKpUshcJAKA+MhEonC6EwEMoLI6FSSjMUBkJ5YiRUWmmEwkAob4yESi3JUBgI5ZGRUOklEQoDobwyEhKdDYWBUJ6FpLbX0D8wCzgDOB6YC4Qai0ZgPXAbcG3s6x1MZIASEELoATaNfLkgxjhQc9u96evnMGvOoxz/4utGlh8G1gI3xb7e60fez0AoNaF/YA7wUuA4oIfan7vDwBrge7Gv98Yx75HENhv6BwLwBSoDbcQNsa/3LR0YklTT+FBw/ZrLgWfusuCKe57ItBnb2GfxXyd4m6s4ZcH3MRBKSegf6AauAY5scNWPxr7eT1a/SOpw09/QeCAAXhj6B57Q7sFIkxlz6GnJYWsYGjyh4TfZuul1zJpjIJSm42g8EACvDv0D06tfJBWJw1pY94h2DSKEUGtXSxpjRygWPwlWLl/K0GB33Stv2zKLgZVL2X/pZgyEGtDmz6hmP3fnAQdWv0gqEjNSWpdQ8ZwQwneoTEguauX9VB4xxs286bLzgfpDsW3LLNY8uJgQIlfe8FwDoXqFEJ4FPBpC+EwI4fA2vGVbPnfTPbvpA6/ah3948lLOeMLBnPv0g/j2Z+a3661DCL0hhH8BVgDfAV4MzAZ2a9fPUAnMmb+dfZfcDYwNxfqBLj75r3D+8TM55/Al3PDVeWMCse+Su+melubIlT/7AdOBVwM/CSHcFUJ41cgcWfv99Y/T+dv9DuEDr9pnssXS3YrPuvARFj/xYWbMitz3uxlcdMYBHHLkVg49elszbzeyq3YS8FbgBVRm7GePWmSo5TGrfLq6Kh/6D917CCuXL2XRgX/hoxfuzbTp8JEbt/How6v5wCv3Z4+9A/strSzrkU01ZwjopnIm0mHAMuATIYSrgWUxxrva9pM+8Y6FLHnq1qkWS3dP4uAjtjNjVmV3PIRICPDgvQ3vItXYa5jJ2EBIzauGAmD5H5by8xvncdrrYFYPHHp05IgTAj/5LgZCbTaXSjDau3fxvWvm0bPbEIc/c/NUi6Z/Md1/vHlvTt//EN544kHsvmCQ40/dWM9q4+YaVgAXA/sy+TUYUvOqoXj4fujqgoWLIQ53sebBxRxwcGTVXzcZCHXI+L2LNU3PXWxc18XXrujl/EvX1LN4+gdNL1y2mguuXM1vfjSbO2+bzfSZU070jew1XEBlfmEO9UehCzgphHBQ8wNWqbz7i4dx1Ilj/2rrnjbA7LmVK7KHBqcTQmR+7zq2bOxhy8ady/73x48Np1yxR6LjVZ4dXedyc0f+fTVwdgjhXuCKGOMX61r7qn/r5blnrmefJ9R1oXL6kQDongZHnbiFm76xG9d+anfOvGDdFGtcRmW+odE9oTnAfzY1RpXTT2+AxU8a+9rQIGzZtPPrGAPr1+7BzNmwdtX+O17/5S0fSWiUKo5Gzoar7l08Bfg88MUp1/jTHTO56/YePnnr8np/SDYiUTU8CCuX1zMncTjwT8A5VGIxd/LFd9gAPCPG+KcmR6iSCf0D5wLvHPPi3N1nMzx0AA8/AAsPqLy28r4NPOEpj7P/wQM7lrvyhlfEvt6fJjhc5VgI4e+Bz1K5TqEeG4DtwCeAz9S1xq9u7WHgoem84oilAGzb0sXwMLzhhJl86raJ7hyQ4pzE2lXdfO+aeWzeEBgahB9f38Pt1+/GkSdMOZESY/xtjPF8YAHwJuAuYDPgfZ7UWdu2zGLjugN42knwnc/C0OPbuOdOuOOWeZx8Zl3zaVILtgFbgZuBs4CFMcb3xRgfrGvt0163js/99F4+fstyPn7Lcp739+s48oRNXPLfK2qtkt6eRAjQ/6Xd+fS7FhKHYa9Fg7zq3as58Yy6f9FGror9MvDlEMJhNLd3IdVn9HUQb/3YX7js/IO54AUzmTt/kLPfMY2eeYsZGryH7mnDaQ9VhTNmr6HuKIw3e05k9pydlwLMmjPM9JnD7Lmw5uUB6UViz4VDXHnDA+16uxjjb4HzQwgXAn8HvB1YSuXKwWwdVlP+jL9QLgR402XsuMHf8HAYuY7iYBYdaCjUDtuozFHcDlwB3BBjbO+1XuddvHaqRdI/BbbNYoybY4xfjjEeARwDXEXlUNRGKtdOSI2ZKBDjjb6OYuXygxkaLNzvlhIxk8pew1oqJ+gcHGM8OcZ4fdsDUaekNuRW7l/T9LoTzF18C1jZwlhUNr++9cApA1E1PhR//EW9E5ASwB1ULgZufK5hYm353E0qEq1M6LU8GThq7+IfYoyPtfp+KocQwtFc/6X31hWIqtGh+Mg/3xhC2KvDw1RBxBjviTG+rI17DZumXmTqdZOKxO1NrjdE5aEtUqJ2PFHujz/f3PCtNrq6IosO+hkr7gYYMBRKSbOfuytiX++O02ETiUTs630Q+FijqwEfjH296zswJKmmUY8c3cTAyrmE8NEG32KQ7u73MDxcfWa2oVDiYl/vH6lcZNeIrcD7Rr+Q2DOuAUL/wAFUHgM51bHa9cBtsa93VedHJe00JhAwr/o8iAm33Y++9bPMnvMAr/vA+0deiVQmHG+Nfb2PjLzf6Eeh9sYYpzybRGqn0D9wIJWn1M2ZZLEIrAZ+EPt6xxySTzQSUpbVCsQky0fgzhjjpI+INBTKM0/Tk2g8EI0Y88xsDz0pZ4yESq+TgagyFMorI6FSSyIQVYZCeWQkVFpJBqLKUChvjIRKKY1AVBkK5YmRUOmkGYgqQ6G8MBIqlSwEospQKA+MhEojS4GoMhTKOiOhUshiIKoMhbLMSKjwshyIKkOhrDISKrQ8BKLKUCiLjIQKK0+BqDIUyhojoULKYyCqDIWyxEiocPIciCpDoawwEiqUIgSiylAoC4yECqNIgagyFEqbkVAhFDEQVYZCaTISyr0iB6LKUCgtRkK5VoZAVBkKpcFIKLfKFIgqQ6GkGQnlUhkDUWUolCQjodwpcyCqDIWSYiSUKwZiJ0OhJBgJ5YaB2JWhUKcZCeWCgajNUKiTjIQyz0BMzVCoU4yEMs1A1M9QqBOMhDLLQDTOUKjdjIQyyUA0z1ConYyEMsdAtM5QqF2MhDLFQLSPoVA7GAllhoFoP0OhVhkJZYKB6BxDoVYYCaXOQHSeoVCzjIRSZSCSYyjUDCOh1BiI5BkKNcpIKBUGIj2GQo0wEkqcgUifoVC9jIQSZSCyw1CoHkZCiTEQ2WMoNBUjoUQYiOwyFJqMkVDHGYjsMxSqxUioowxEfhgKTcRIqGMMRP4YCo1nJNQRBiK/DIVGMxJqOwORf4ZCVUZCbWUgisNQCIyE2shAFI+hkJFQWxiI4jIU5WYk1DIDUXyGoryMhFpiIMrDUJSTkVDTDET5GIryMRJqioEoL0NRLkZCDTMQMhTlYSTUEAOhKkNRDkZCdTMQGs9QFJ+RUF0MhGoxFMVmJDQlA6GpGIriMhKalIFQvQxFMRkJ1WQg1ChDUTxGQhMyEGqWoSgWI6FdGAi1ylAUh5HQGAZC7WIoisFIaAcDoXYzFPlnJAQYCHWOocg3IyEDoY4zFPllJErOQCgphiKfjESJGQglzVDkj5EoKQOhtBiKfDESJWQglDZDkR9GomQMhLLCUOSDkSgRA6GsMRTZZyRKwkAoqwxFthmJEjAQyjpDkV1GouAMhPLCUGSTkSgwA6G8MRTZYyQKykAorwxFthiJAjIQyjtDkR1GomAMhIrCUGSDkSgQA6GiMRTpMxIFYSBUVIYiXUaiAAyEis5QpMdI5JyBUFkYinQYiRwzECobQ5E8I5FTBkJlZSiSZSRyyECo7AxFcoxEzhgIqcJQJMNI5IiBkMYyFJ1nJHLCQEgTMxSdZSRywEBIkzMUnWMkMs5ASPUxFJ1hJDLMQEiNMRTtZyQyykBIzTEU7WUkMshASK0xFO1jJDLGQEjtYSjaw0hkiIGQ2stQtM5IZISBkDrDULTGSGSAgZA6y1A0z0ikzEBIyTAUzTESKTIQUrIMReOMREoMhJQOQ9EYI5ECAyGly1DUz0gkzEBI2WAo6mMkEmQgpGwxFFMzEgkxEFI2GYrJGYkEGAgp2wxFbUaiwwyElA+GYmJGooMMhJQvhmJXRqJDDISUT4ZiLCPRAQZCyjdDsZORaDMDIRWDoagwEm1kIKRiMRRGom0MhFRMZQ+FkWgDAyEVW5lDYSRaZCCkcihrKIxECwyEVC5lDIWRaJKBkMqpbKEwEk0wEFK5lSkURqJBBkISlCcURqIBBkLSaGUIhZGok4GQNJGih8JI1MFASJpMkUNhJKZgICTVo6ihMBKTMBCSGlHEUBiJGgyEpGYULRRGYgIGQlIrihQKIzGOgZDUDkUJhZEYxUBIaqcihMJIjDAQkjoh76EwEhgISZ2V51CUPhIGQlIS8hqKUkfCQEhKUh5DUdhIhBCWhhAeCCG8sMb3DYSkxNUbihDC20IIfwghzExudLsqbCSAtwKLgG+ND4WBkJSmqUIRQngb8H5gMfCyhIc3Riji52MIYS7wMNAz8tJm4IwY4w0GQu0SQojAnTHGI9Mei/IphNBD5bMIoDfGuHZUIKqfX3+IMR6aygAp7p7EK4DRH/49VPYo/hUDISkjJtij+DfGBgJgcQjhmMQHN6JwexIhhAD8FThgksW6DITqFUJ4HvB1IIz71h4j/z467vXtwDExxr92emwqhnF7FOMNA9+JMb4kwSHtMC2NH9phJ7Pzl3cim4EXADckMxwVwGNU/rKbVeP747e37cD6jo5IRfMGKp9NPRN8rwt4UQhhUYxxZbLDKubhpv8LzJ3k+9VDTxOe9SRN4OfAQ3UuOwxcG2Nc18HxqEAmmIOYcDHgTcmMaKxCRSKEsAQ4ro5FDYXqNnJo8gpqHw4YbQvwkc6OSEVRZyAAZgL/FEKY0flRjVWoSFA57bW7zmV7gOtCCPt2cDwqjqupb9taTeXkCGlSIYSTgcuZOhBV3aRwOmxhIjFy2uurgel1LL4RWAf8O7Cmk+NSMcQY1wPfonI4qZbNwJWeFKE6/Qb4PJW9z811LD8XeFdHRzSBwkSCXU97HW+YyuGC3wOvBxbGGN8TY3w8icGpED5K5Re6li7gKwmNRTkXY1wTY3wNlYt+3wOsAjZMsVrip8MW4hTYKU573Try7/XAZTHGnyU2MBXKyHZ2D7Bkgm8PA9+MMZ6V7KhUFCGELqCPyt7CkVSOiow/AzUC307ydNii7ElMdNrrRiqnIX4YOCjG+FIDoVZMMYHthLVaEmMcjjH+vxjjM4GnAV9i10NRgZHTYZMaV1H2JG4GnkPlr7mtwHLgUip/2W1PcWgqmBDCfCqHBcZfM3EfsNT5CLXTyPZ2HvAvVK7MngtsAz4cY3x3ImPI+zYdQjgIuBd4HLiOyiElzy5Rx4QQrgbOYuee+GbgnTHGZemNSkUWQuimcijqIuBYKnMXC5KYUy1CJKYDrwH+N42rEVU+IzeJvJmd99zZCizyAjolIYTwFOAZwFeS2HPNfSSkpI2bwHbCWoVWlIlrKTHjJrCdsFahuSchNWFkQnE18CBOWKvAjIQkqSYPN0mSakrseRKhf2AucDZwPJPfyhsqF8HdBlwT+3q3TrGs1FGhf2B3Ktvusew8o6mWR4HvA1+Lfb2DHR6aNKnQP7AI+EcqF+fNnGLxCbfdRCIR+ge6gM8BRzWw2jOp3Pb7tR0ZlFSH0D8wHfgy8KQGVnsWlW39wo4MSqpD6B/YE/galXtD1WuXbTepw01Po7FAVD079A8c3O7BSA14Fo0FouqUkb/ipLScQmOB2LHe6G03qUg8JaV1pVY1u/0FmouL1C5PbnK9MHrdpCJRzzMeakn8SUzSKK1su62sK7WqLdtuYhPXu3jJ4kPGfL19W+D5Z63jnz+6OqURSfXZvjVw5Vv25rc/nsPGx7pZeMB2znnnAMe/uJ7Hm0rpevDeaSx720LuuXM206ZHjnnhBt58xWqmTdyU9CJx7f137/j/zRsDZx96MM8+faoHbkjpGxyE3n0H+fdv38+iJwzyo+vm8OE37cuBh97Hfks8o0nZtuxtC5m/1xBX/+4vbHi0i4vOOIBvfXJ3zrxgwnuPZeM6iVu+OY95ew5y1ImTPfVLyoaeuZHzLl7LfksG6eqGE07bxIJ9t/OnO8bfPlzKnjUrpvPs0zYwc3akd98hjnz2Ju7/U83TY7MRiZu/uRsnnv4YIRvDkRqydmU3q+6fwUGH+uwSZd+pr36UH1w7jy2bAg8/MI1f3zqHp59c81Bp+p/KK5dP44939PDCcx5LeyhSwx7fDh987SKefdpjHPRUI6HsO/KELTxw90xetvQQzn3aEpYctpWTXrKx1uLpR+KGr+7GE4/cwv5LO/7wDKmthofg0vMWMW165IKPPJz2cKQpDQ/Be8/an2NftIFvLb+ba35/DxvXd/GpixbUWiX9SPzg2vk892Xr0x6G1JA4DB86fx/WD0zj4qsfYrpnaisH1q/t5pGHp3HGG9cxY1Zk9wXDPP+sx/jVD2rebibdSNx52yweXT2N557pWU3KlyvevJAVf5nBB76xglk93kpZ+bDH3kP07vc4//vp3Rl8HB57pIubvrEbi5+0rdYq6Z0CC3Dj1+bzjOdtYM5u/pIpPx66bxo3f2M+02ZEzn7qztvGvP6SVbzoHP/gUba966qH+PS79ubbn9mTrq4wIRNOAAABC0lEQVTIU47ewhsvq3l9WrqRePsnPY6r/Nn3oEGuX/PntIchNeVJT9/Gf3z3gXoXT39OQpKUWUlFopXDSR6KUprcdpVXbdl2k4pEK9dAeP2E0tTKHIPbrtLUlm03qUj8iOaq9jjw0zaPRWrEbU2utxH4dTsHIjXoh02utwn4VfWLRCIR+3pXAR8EhhpYbRB4T+zr9WwRpSb29d4DfAwYbmC17cBFsa/XK7CVpluBrze4znbgnaO33RBjcodNQ//AXsAxVJ5xHWosFqk84/rHsa/X3XVlQugfWEBl2+1h8m13HXB77OuteZsDKUmhf+AAKk8HnewGlDW33UQjIUnKF0+BlSTVZCQkSTX9f3sA4pkyxlNTAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from lbmpy.stencils import get_stencil\n", "d2q9 = get_stencil(\"D2Q9\", ordering='walberla')\n", "visualize_stencil(d2q9)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADgCAMAAAADk+xoAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZondIrt8bFiOv0QAAAAJcEhZcwAADsQAAA7EAZUrDhsAAA0ZSURBVHgB7Z1tW9vYDkXDa2daoO29+f+/dXxsdiIlsvfWE6XJDOLDHGMtbx0vQmCCCruH/fz2uOu3cgM/F7e73cP+6Xl6eynv0IG792H2cT8kP7SOaxp4b8nX1Ltkn0p++fG61ZSUcamIJXE1VeRqsa17cZJfH58e9+uSSRltRCyJq6kiV4uRe3GSJ/Z5QzIvoxlJAYZVxEWM3cOVun7GRptsyXC+tkbW1tjpfIS35A1hcymytnFNhLfkDWEteV1O9GAK6FoMDaLUfiTDztoaWVtjp/MR3pI3hM2lyNrGNRHekjeEteR1OdGDKaBrMTSIUvuRDDtra2RtjZ3OR7iX/PT4Y//2+L6WQcq4TMSSuJoqcrXY9r14yWB7LTXQkkt1xmEtOfZSerYll+qMw1py7KX0bEsu1RmHteTYS+nZllyqMw5rybGX0rMtuVRnHNaSYy+lZ73k16eP948/PU4kNhWxUjtFYV7y2zQP9/rz13q2OBCCAAmnTZc0EZtgqauMpe4FsE93kt9/DOjj7Yi6I3EgBNeIOGmKNBHbiV1FDN0vw53k3/P07MvGfEv0aik2EqwKzpvOwSI2WKWrjs3dL8Od5P0s+ft+fYhW3H9mY7zpnCZigxU3KWKZewF7sgkr+XX/NKjv+2cD+8P6jQlNxxZEbN6tuEkRg4ALcCv51/5jJD4sC8LdekEnl3N8R2g6YBGbc8VNihi2egHuJc+P5D8tmTZdJEvY/UsWPiW3Ppyvbz8Pb2/LN9tb+OcjRGg6SBGrkpy8F4KPR/K3/bfljpcvLi+3+MK31XTeHN/bcg/Tf4UP7WBFDLEX4P+z/5zh9++R+P6Hv4WjTefb5HuDDdXeBdYOrdYPbLp9Tt4t3/A/rf3PyJRoL11vcKgoOG86x4nYYJWuOjZ3vwx3knc/x/9W/9h48ULcf2pjtOmSJmK6jmvcC27cf6S95NeP6R+drTsWB0LQScRJU6SJ2E7sKmLofhHuJSOy11IDLblUZxzWkmMvpWdbcqnOOKwlx15Kz7bkUp1xWEuOvZSeHZJfn9f/PXVps68a9mt6gb5/38WVP/r9dHFlwSO+Jf9xyTcZIBGbitgfcJZu4R/JdIDkGnMjtOlyUyI2wbWbFNNgPsKdZDJAIk54iBh2RZomseLhluS9rOBOMh8gEV+DFbFZIG+awQYrdq/F5k2uNHeS+QDJFTbGm877F7HB1m5STJs3udLcShZ+Iix2FLGxMaFpAhvofUsWBkhEeyI2jAhNE9hA710yHSAR7YnYYo82TWADvW/JwmeuaE/EhhGhaQIb6H1Krh9uSUjeiVMrInaXkq8z3JKRLE6tiNhdSrbfXdQNt2Qki1MrInb/ksuGWzKSedMhbqdiEyl2r8XmPa40d4/kHRkgESc8RAzbIk2TWPFwS/JeYtxLxg31WmqgJZfqjMNacuyl9GxLLtUZh7Xk2Evp2ZZcqjMOa8mxl9KzQ3IPt5QqPQ/r4ZZzJ+Vn+umiXOl5YEs+d1J+xksmAySkjM2JGPD//uolkwESUoYtEQN+XKPBkGNVHVo5XFEbR9LQNcKcZDJnQspoI2LAsa4MhohlYIe1No6koesK5iSTORNSRicRA25W8govKZug5ZDwpHwaJ+IR5iSTARJSxq5EDLhZow3qZUMuh7VxJA3dI8xKJj84JmV0ETHgdo02aOqkbMjlkPCkfBon4hFmJZM5E1LGpkQMuF2jDZo6KRtyOSQ8KZ/GiXiEecmbcya/ll9RtPWLXcbGROz0Hsb70QYNR8qGXA4JT8qncSIeYVYy+UQnZWxKxHbkt50gzqzR/g/l2rh8GjYSbXJIVodbxPkSEcOuzBptUC8bcjmsjSNp6B5hmeEWcb5ExLArs0Yb1MuGXA5r40gaukeYfbpgwy3ifImIYVdmjTaolw25HNbGkTR0jzAnmQ2QiL88RcSwreMabfBYZV8XDbkc1saRNHSPMC+ZzJmQMvqIGHCs8WAIqurQisqTboeYzwMRjzEv+TS63y8x0JJLNG6HtORtPyXVllyicTukJW/7Kam25BKN2yEtedtPSXVIPrx2UZLYIWcG3GsXZ9U+UWKgny5KNG6HtORtPyVVL5mMpZAyNiRiSTyZivR7WL1kMpZCyrgfEUviydQpPZozQVNeBklSgJ2u9jInmYylkDLaiFgST6ay3+CyMoaCTWEVMeBYTy5zkslYCimjg4gl8WTqSI9e2kVXXgZJUoCdrvYyJ5mMpZAy2ohYEk+mjnR7n+hmVlIGKWLAsdrLrGTyc2ZSRrqIJfFk6pxu7xPtzErKIEUMOFZ7mZVMxlJIGekilsSTqXO6vU+0MyspgxQx4FjtZV7yDYZbxFkYEcMtzqu9T1dY3iFlXCFiwLHay4Zk/MMc8ilJykgXsSROU/PjKFYDdhOsCkaau3+Ys3xxWf0LQaSMDYpYEk+mjnSih5SxPxEDjtVeZp8udmQshZQRL2JJPJk60u19optZSRmkiAHHai9zkpdv+Ff/LBEpI17EkngydaTb+0Q3s5IySBEDjtVe5iTfZrhFnIURMdzktNr7NKdxSMpJDDhWm+4lk7EUUka+iCXxZCobhonHULCnwypiB/7zwF/mJZ+y/X6JgZZconE7pCVv+ymptuQSjdshLXnbT0m1JZdo3A5pydt+SqpDcg+3lKhcD+nhlnU3ZZV+uihTuR7UktfdlFW8ZDJAQsrYlIgl8WQq0u9h9ZLJAAkp435ELIknU6d0O1+CZmYlZUPOhyIeYU7y8prtx9qfuSdl7ErEkngytWi4BZs8mVbB6dN1BXOSyQAJKaOjiCXxZOpIty/poptZSdmQ86GIR5iTvPwc7fv+5bTB8j4p4yIRS+LJ1JEe3S+68rIh50OSBjzCrGTyE2FSRhcRS+LJVEVLpAO7ClYRjzArmQyQkDL2JWJJPJk6p0f3i7bTSsqGVNKAR6lecg+3wNX5Gtk7p8IP3ZDcwy1nssi0yhl/PBF9LHq45ehn+yiyF1wRYfbpoodbAmmHU5G9Q/F4EGFO8vINfw+3HJWZo8ieKeMwwpzkHm6BqmCN7ImYl0wGSEgZLUUsiSdTi4ZbsEk/rYKzZ2uMeclnF/WJCgMtucIiyWjJRFBFuSVXWCQZLZkIqii35AqLJKMlE0EV5SG5h1sqTG5k9HDLhpyqUj9dVJncyGnJG3KqSl4yGSAhZexJxID/91cvmQyQkDJsiRhwOoUCMJobQc2uIleLYQNRqpNMBkhIGW1EDPjKQAjKWEWMDbUk49Sun7EruJNMBkhIGfsXMeDTesFLtSblcFgbJ6ahe4Q7yWSAhJTRRsSAT2u0MVPGoYgVx6ldP3cZ4VYyGSAhZagQMeBjjTZm65/HIlYcp3bd2KSVTAZISBlWRAz4WMX7ELHiOLXr5w1FuJd8g+GWryL5psMtX0Hy7YdbvtjTxW2GW76a5JsMt3w1ybcZbvla311M850fz8+PD+NrUfhGyrhGxIDHAyGoHlYRY0MtyBPjRGw71X4LB7LXYgMtuVhoFNeSIyvF51pysdAoriVHVorPteRioVFcS46sFJ8bknu4pVjqaVwPt5waucL7/XRxBamnkS351MgV3veSyVgKKWN7IpbEk6lIv4fVSyZjKaSM+xGxJK6nRvMlaGZWEcMVF+BOMhlLIWXsRsSSuJq6Ml+CblhFrAZ3kslYCiljQyKWxBOp0Q+M0c2sIoYrLsCdZDKWQsrYjYgl8USqqEPEsM8LcCuZjKWQMjYjYkk8kyrqEDFs9ALcSiZjKaSMzYhYEs+kijpEDBu9APeSbzDcIv69IRGbjYg6RKxI8re//p6TyKckKWMzIpbEM6mivS0s+ZtbCP7/v3a7h/3nz06XLy7/+j9LtGUPH9tpFTFccQFuny5uM9wi/r0hERtKRB0iVi15+Yb/X/+bW0R7IlYt+TbDLeLfGxKxu38k32a4RZyFEbHi4RY8ji8acnHPyUjstdZAS671Gaa15FBL7cmWXOszTGvJoZbaky251meY1pJDLbUnW3KtzzCtJYdaak+25FqfYVpLDrXUnvSSyQAJKWNnIpbEk6lIv4fVSyYDJKSM+xGxJK6nimMotRjuJkp1kskACSmjjYglcTVVnFqpxXAvK6lOMhkgIWV0ErEknkgVX42vxXA7UaqTTAZISBltRCyJJ1Kj+0Q3s9ZiCI5SrWTyE2FSRhcRS+KZ1Og+0c6stRiCo1QrmQyQkDK6iFgSz6RG94l2Zq3FEBylesk93AJXxzWydqyeHUW4lUw+JUkZ7UQsiWdSo/tEO7PWYgiOUq3k3fLFpYdbIGxeI2sO8O9EuJNMBkhIGc1ELIknUqP7RDez1mIIjlKd5OUb/h5ugbB5jaw5wL8T4U5yD7d4YdeRTAZISBk7FLEkrqaKYyi1GO4lTvWPZLC9lhpoyaU647CWHHspPduSS3XGYS059lJ6dpG8H2+PpcEdNhv4Oaud/6Ty83j73l7qDbzPap93/wAdke0u6RSmcAAAAABJRU5ErkJggg==\n", "text/latex": [ "$$\\left[\\begin{matrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\\\0 & 1 & -1 & 0 & 0 & 1 & 1 & -1 & -1\\\\0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1\\\\0 & 0 & 0 & -1 & 1 & -1 & 1 & -1 & 1\\\\0 & 0 & 0 & 0 & 0 & -1 & 1 & 1 & -1\\\\0 & 0 & 0 & 0 & 0 & -1 & 1 & -1 & 1\\\\0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1\\\\0 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1\\\\0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1\\end{matrix}\\right]$$" ], "text/plain": [ "⎡1 1 1 1 1 1 1 1 1 ⎤\n", "⎢ ⎥\n", "⎢0 1 -1 0 0 1 1 -1 -1⎥\n", "⎢ ⎥\n", "⎢0 1 1 0 0 1 1 1 1 ⎥\n", "⎢ ⎥\n", "⎢0 0 0 -1 1 -1 1 -1 1 ⎥\n", "⎢ ⎥\n", "⎢0 0 0 0 0 -1 1 1 -1⎥\n", "⎢ ⎥\n", "⎢0 0 0 0 0 -1 1 -1 1 ⎥\n", "⎢ ⎥\n", "⎢0 0 0 1 1 1 1 1 1 ⎥\n", "⎢ ⎥\n", "⎢0 0 0 0 0 1 1 -1 -1⎥\n", "⎢ ⎥\n", "⎣0 0 0 0 0 1 1 1 1 ⎦" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "M = moment_matrix(moments, stencil=d2q9)\n", "M" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAByMAAAA0BAMAAADbH4nKAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhAiZrvNmd12RKuJdf+/AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAUVUlEQVR4Ae1da4wkVRU+0z3VPTPdszOKiugijUaj6ypDfMRE43Qi/sAfTq9RIwHZAo0Q0cxgIsRXGN+IDybRqCOKLdHIQ2HwAYqoDeoPNoQZjDHxlR2JmqBENyLra2E891H3Vefc6p5amW22K9muc8859zv3fqdu163qe2cBRseIgREDAzNQufZrsToF5ljVkW3EwIiBbTDwLLg/VqvAHKsKp2z9PWofGUcMPPYZeN7AXbwFFtJIpQJzWDPZVJrK1tYsnHDGK0P7qDxi4PhioNIduL8/hMWNSKUCc67mr5QmOeMVs/CEnHWkGDFwnDHw/O30d74TrVVgDurWV7RifDQkA2pGxeOQgeoTt9Pp2+OVCsxh5Q9oxWhIhsyMyschA/tiU1COj4LJboE5h7pvTqlGQzJHzUhx/DFwxXa6/O54pQJzrnJ1VqmGY0hemmt/eUX2iqs80jGDUF86Sk05N4Lz/8gFPBoBIwmvPhTpMGeabk23ORvqC8xEzW8p3VAMyeYK0YH+VNPf+APnqF9xceYh1N/VR5sFIed8A2rRh6dpnvASuYBf39NmWlg2ICInB3pwXofBl2o+4Y1urJ5ns/yd+YKX5MMVmCMUAMynMtJQDMl78133aIoU9sC9nLW5yVmGVD8dHWa6U4KQ7unQPBTt5D2stUQuKt3pWQ63XECBPIUTv3UOXur5hF+WRiu6RsvfKVtbrkHJBeYYBaC/GIZhSCbb/5mmvhb58ehzeUaHWrOvXdx8Qcgvl66DiWVo7m2x/pNLjKlELuA9UHuYgYVyAQXyBVMrcDnA+AoXAoBN+Pf4OoHF4S+wyGKBOUoB1NX0eRiG5K5Nqvd96caW4OAc57nY4SzDqf9rH80WhLwIdsNYD94BL2cr1LiHqxK5gBNh+kEuYrmAArmDXdoN4zeucSEAuISPs98UOSyHv5wNFQXmKAWQHJGQwzAkzbqlopW8e64IB9nBDqynFHdCV+by4jB3UD9+2A+eZwPtkpDph2ChDT+EmZAuC8ANb5MLIOEtQD5XGLUZNNG6Q5mACvlgZ/wfANXIkOQSPjVrG2K7ZSVrdfmjZhkOvXCmU02LcQpgVboNw5C8OutbwUreSrprJXPV5xsArgpUtjjO3QqsyzBJk5teawk20C4JmVrDp67kHzA259VwCwsdt2Rlkwsa3jrmc4VRJ9esQyCVCaiQb5APyLEhySV8bNM0xnbLSsaIguWPnGU45lf9xa2n5DgF+u4xBEOyZr5aC1byNtLwRgG3gfjm5I4vcoah1C+mXrMJNtAuCZnYhMvrrz0MjSWvhluY7LklI9tc0PDGEfK5anTFvZk7ygRUyJ/CB+R3Ru+SwCR8wfJgu2Ult8mWvw41y3DM8FK3npLjFMCibMcQDEn8atFHwUrexkb4+iD5F5xt6c5gzPliIz0WBPyGdo88G2hVhEwtTz+lgreMRs+t4Mn1Wa+YFWwuSPjMDc/5XO1fTm51HAKxTECFfC3c223FhySTcOfZxnbLSk5LHf6oWYZjBmpIxikAdbcegiGJj+3mKFjJG/7kWz1891dN3bwQXcmfdz/GNbvD9oVsoF0Rktz3k1f/FGcULrNB5fDLTZvdGgS8BxLmav5tB9qeg1coE1Ahn7v3hQ+040OSSfj73ZbYblnJ2H3+crMMx0wOyTgFMCFvPtmQrEVuJqZBrlBP3ZKRB8YxNQEYyMW29bndipQUfu3je/HYMdYirReQ2j6Vl5B+ZWjJAdI81XIz9JANBHIIwW/5mY0ctlE83kiu4OaCgHddIczVdZ41VygR0EWOPUvCWCsXViie7Gptt6xk7A5/1CzDMZND0m2owbRC5ZCQsyF5oTX0Kf2J9Bscx4GhIS+zHpWulSlpT6B0v9QDkyjm3gZJn9qyPG3zY1dKVSxFSw6Q5Kku8+n6hmygzSWEehay1e+3oiM5uQAC3vGEXK5yd3HXG5gt+n0FdJGjQ5JOePIftyW2W1Yydoc/apbhmMkh6TbUYFqhLn+MyYbkZ9Fwesdai6WzSBeBs+2DhnSekQpW8ta7QejFdqDwi/QSlsnU9xqslKxQ/qVoyQGSPNnHPO2fYwP1LiEHIr9LArMQxskFBe+2NMxV/i7uepcJ6CFni7d9cF2iEz5+xHG23bKSNTv8UbMMx0wNSa+hFtRIqiF6SGL45M9/G2hIVnoGygpUN6y1SCIhwS6tKFrJ+3P4hQrxMh3pN/GI1cOUPbyWKJ+Ibi9hK0dLDpDkqbEc+Bk2HL1LSOWa2J9zmk+daka0uQAK3vg5q66rPaWdXnOshLj9gC5y7YZHdDwiBD5Kk1p3AYPtlpVsJZc/Z5ZBXm/2jWtmdhtqMa2kbtd6SDZSNKwPNCSTWYtlJIljSoMKJCTYaQW50NcGGf/0Cz6sShkH1kZK424ujMctRtqWsI+oVY6WHCDJ00zX97Ns+Po+SwtzlKPNRRG8yVU2JCk4V1c2oIvFyXTC685TuO2WlRg0Z5ZBXm+nm3qk2Vgd4REh6yF5kZAHG5LwJVEnOCROoBugSEHCEQNALvQ1Vhjb2jqkSn1ykPzTVjZSbc2I2xImN/LVStKSA6R4mun5bpYNX99naf8S5WhzUQRvctXvkCwbkGpuqCMT7q0pst2yUoiiy84sg7re3vrjB7KKlDmzeed/iZIekgeETA1JYnPTBLrWTgby8R5xmE0y0e0zMUhJIwMqWi32u4TmAg6WsdZ5YuTYKwwL+mh2sQ93/+m3xASUiJTVAsj6kHutgS6C3nPuefudvbChuroPzO+dWkZ/0fA9up57Upd0LIrrzclZLyAc4bKCuqRlDBrA7k0y9oIhOVhAhj8TjBZMDDLhUFHTWdstPgE0PhRcbwXmZUSVlyOcJPD1kPyokKkhSWxueh+6ihWy7xJ1ggNxKuQmGWf7DLFWPwZZE7cyGlQGd5BNY+Ic1MVo3N/DD+qvZU7g7eE18J7e9w1aJjiR7rwzU5pz1of6ilEZAWlJWvXbGyt0PwJgvXfqTFM9E0zDKerl0o9YlAwles56ATMtwk/mQsUgrKiye5OMvWBIDhaQ5s/EYgQTg0w4TMmFlU63WoWb14JA8eutYMSarAI8VeDqIfltIRNDsk5sbuoK32dQb5QAEOdN+OPM5cLFO+z2GWqtflf4MpBy/k+DyggWWRblR5yiKeFTmcMP+2gkq8mPBhqeCad1XmxVWrKR6q3qZmjuCgX2gXoLj7RUOpWHqindDx+4pvZOEWskTcPtmwPbCvksFotiXSNSV9iwFzAjJVFyDpkLFcPRWlFcLnrrl1EWDMmucOw7IM2ficUIXaEXMciEw65DwmS7VZuTm9eEss8jfr0VDEmTVYCbRTw9JE8WMjEkxWYTce06x/QGTAJcArDgKDMRccQmmSNZ2Zzt9hnimo1Cjou7JA0q0R1kdf99zurqzaurcn30ljgk4aYhKPwOYEWtS6CG5NgGhoPPuxW0bCM1NiBcsm76MJ6LB4C0JDC5zPXDB07U3iliQZZpOEW9XKgZiyJ6IQlxPnTPzMn0ghmSIhcqhqniCuJy0Vu/1D6JyurqVSetrm6i04wM63oLecCA5HWgkG2vwgTYGLEhabslEkDO22n+wuvNNkRucXbMNAUmq/0NSbnZxKexCfAqgBY7JHH/Sn6pd3z7TBRSTpZIUNkuB9nef+PfWj0Y70IFayfUxFUMSfIp04m0vw0nyOj2w/SBGZJ4TfbQu4AcCYyBFtrULMQ0nBqSas0yH8U2NSKZXsD+FuGmcsFdr97WL7MbM36XHDQgxR/RUE9lY5AJ1xNXnSBRUyfAA4kW4tdbwV3SZBXgByJKwcQVfxq+KmgM3md/ArUWwG8CgyiKCfB10JwNTfg79uQa7mxp4jdYfmYXhVQbOylQGYREjlPUheYcnI+19RXmN1ZMXMn9I06khRRu9WuB6QP145OgRW2DoPoRAmN5Hf3zs1PTcIp6+XAciRK0ly6aXuhdCYGX3mQrvjDIQ14uqvnmV7v4kBw0IMUf2RRHaWOQCc9e79h9KvjodYdTv1iMX28FQ9JkFfp6vUNsbppoj/8HJlP29Q7gJpnNsBfx7TNRSPVIToHKICRynKIPwUIXPou1yZ+pJpYAXoRfHT+L9GGhk1vBafrAvN6pdtah2SkiRwIjfZfXO8SQNA2nXu/I9zGRKGF3yLLpBchH05yPmFboGDkbKuzepNeZ3ZjxITloQPY6oJqjdTYGmXBoyt8lnW5hAj4Vwcub4tdbwZA0WdWLbfVdco8Is96B5N9ePLXZZP+Kq6w8/ivz952ImmvwX+Av385fCweXwK8D4fYZ3xyFhKeJ6BJUCEXI0iegKGjlqXtOeORUHHgw/TABWOniX42YOFSL9QHvkruDhpg+TLVIWg7OXQkX6H74vXfJkcBy75S8S3INF9QHIGpnjxtF9C0AEKrg8GFML+DgBjrmqotc6BgCJ7A7e5NwBWhDUIwjuCdP5qNUQHUd+BAG2Qq+g+0UmXBQS0udbmECxEMzHkEHlTL/GVxvoUNgDkDN5ajfPukhKWZCN55y9RL8oOPiqc0mje+6uuTKzzQf10PNZ4TW95eT2XP33pSCXwfC7TO+OQoJHxNxJKgQipClT8BB0Mo3/Lt9x27hp/7+kN8WwIln8sDPDlwai4SPfA8HDTF9mGwjcp6Wc+8+5+t4jRaQI4Hl3ik1cfVxTMMF9UGzYWIFlW4ULOYaInXehw9jegGXyQvBjw8yFzqGRPHtzt4ku08iHJLlAlL8ef2RBS4GmXCoym9mp1vJfQfUT5V98Cfjhddb0KLQ7LNmsqrfbegh2Ug1SjUTZBkn1eLI/x6AytqaMFVT/LCHwrkOFX4dockO+Szpm7WJhIRTldVA+FWNGr3MU6rhoPb0Z88JQ4of+aMyK3U+IO6Hzw7f4ERqyDeuvjmrdJEQqil+2MPQiyoB49fMAauK4lmymuJH7lA8+SCgO6OdM8xqmquNqVOsKEsAk7kr1qtpVlZnnQujrKZGREFfLkJl90nYIXn+0z8oTKUCUvwJUDw0vJSZGHTCvWXnsrrtSTVVCu/T409YzPUm3Xhzcs3zU3Spig/iCJedKxfxItIe+HuGOFryM/hQa599f6h3ccdTIr51/Dq7ncpVMRB8s7aSkCC2wmpQ6edX9ZDXNNDr9RneClUxBQlamVkbm1JqZWV1frcp+gYnUrMrfpf0zVmtvUIIAiIt6igiRwErX7FGMsDRKIqnli7pk7Ph1mGLBNCsqIoBTAb6IymE1UUu3MO368tFONh9EkmqK9S6ML+EckuXg1NfAUn+FFAGL0utAFwX6YTnfhq5H/Zv6Bp+B7XS40/ozPUmHXjzvqXak9CFBEV99UFRX98l4SOigMc56qQ/F9tCSHriMzzOkgrfHxCnsTTRC+u4u1LkWv1BIMVKMQUqI/pVXWR7lzRtXZ+D/2IhbKW2qzVoPiCujEu12Td4kfbswefunvbzTmpzVhgwo7eIHBDA6pBrJEMcZZLUh+Gd9ZsOWySAZkWChTA6ukYLqyvGtA+efLu6XJTVvHE1zs0HYWYtvDKMtb+AFH8aQsPLEtMpoBMO3zetUMKt8IlM43dQa13+MkfnzJt/D/BndCRBUa9W9mVD8kINiZe/c6iX7XVHY0V1lfn+ABfCtPzbGn4dfD7zD9+c2WhI8f5WgUo/v6qHLO+/GZg837shn6XDVmqf+VQIPiBO65aFFg/f4EXKm2UN/FADOgyY0dsXORkUnkMcZZI8+a1Dg31J6LBFAmhWJFgORoXQ25jC6uG7dN+uLhcF4OyTUAqYPlm+kioVkOLPh5clJgbQCYebNER2+u2Bdib6HdRal7/M0TnzZnx5eho6kqCoV8uIsiGJLxgHO+op6T8wjoPCQNZnHZ+YSO6VkxNXplb2wBWYLwjKAxUvIb3L0JIDZHi6OOfIK2KsyFpTK2TlvnMBzj4JizTPX2WlA2KUCLxqA5Nw9SrLNrNYKuCPMeMbaVhv8+hjy8KWDUneb8ctiZxhb7cZZ7f4mifypiG07J/rv9ExViTKGD14yuUi+RbfwqMQMAavAjMJXxiAOQVUwB9nPgnvkhs8B4s9YRuCIQl38Z0otJx/W4f1yS8yYl2HwdBc6buVMVYUyB0MVplc1J77NQYV1eUDRuFlYC7hY5t8u0hLAX+sGWetX14iEaVyPRWnYRiSY23R0m0elZPZiswXM+t/rBue0H8DI6woEOZ+AqVyAafyt4ijETACL3vFJXzqUP/MKc8C/jjzZK962hwf6+PSNAxDsr7C96LY8rgO51PmG5/D3En9uwYIzrMiQdg7brlcjD2Ra+JRCcjDq7BcwtWiHq5ppL6AP8783p+e1ibxhDI5Ik3DMCS5vxvP9s0a3gJwJfe1JH8isq7DL032+uxDjBUFcV7KQX2RMxTqaylMsK8FygeMwsvG8Qm/rbDxnkMBf1Hz39g7BDTVhr+hGJJntz1GBihsdfgheTY3VgfAP6Zck+/02ZwYKxIikVtXSLTt52LmIX5IHoWAMXjVEz7hA75yLeAvar6aZFUqxzblaSiGZO0Wvh9xy5MBfpwyLh9l9MOrvojratClGCvSVeyGYY7t5wJf8TcOM6hHIWAMXkXlEz7WZdpFqwv4480X5/8jKSeC+tV0KF7vALzZafhA4qVQVRP0fK1kI68bck2t3V8HIqwogDdGcLadi/oSzPcY4KMQMAYvo0YS7ixGZBroqQv44813dc7iv+zgCyrIUNwlPT4GKkxff+1jb+QNxADlvEOs/Pr6T1KtOVq6MvB6NPTZlAL+eHPl+m/yIar63ddjfEjyBIwsIwYsA4ttK++UtE/fPEZDcqcyMIp7DDFgtvXtYJuu0LFHQ3IHkzAKfcwwcOOOt6S+rJuAQ/KEM1654+0ZNWDEwI4ysKu3o+Ex+K9UA5IzXjELp2xRfztxp1s4ij9i4NFkIPLe5VFpRvJHFaaytTX7PwyjZRSUjSEVAAAAAElFTkSuQmCC\n", "text/latex": [ "$$\\left [ \\left ( 1, \\quad \\rho, \\quad \\omega\\right ), \\quad \\left ( y, \\quad \\rho u_{1}, \\quad \\omega\\right ), \\quad \\left ( y^{2}, \\quad \\rho u_{1}^{2} + \\frac{\\rho}{3}, \\quad \\omega\\right ), \\quad \\left ( x, \\quad \\rho u_{0}, \\quad \\omega\\right ), \\quad \\left ( x y, \\quad \\rho u_{0} u_{1}, \\quad \\omega\\right ), \\quad \\left ( x y^{2}, \\quad \\frac{\\rho u_{0}}{3}, \\quad \\omega\\right ), \\quad \\left ( x^{2}, \\quad \\rho u_{0}^{2} + \\frac{\\rho}{3}, \\quad \\omega\\right ), \\quad \\left ( x^{2} y, \\quad \\frac{\\rho u_{1}}{3}, \\quad \\omega\\right ), \\quad \\left ( x^{2} y^{2}, \\quad \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}, \\quad \\omega\\right )\\right ]$$" ], "text/plain": [ "⎡ \n", "⎢ ⎛ 2 2 ρ ⎞ \n", "⎢(1, ρ, ω), (y, ρ⋅u₁, ω), ⎜y , ρ⋅u₁ + ─, ω⎟, (x, ρ⋅u₀, ω), (x⋅y, ρ⋅u₀⋅u₁, ω),\n", "⎣ ⎝ 3 ⎠ \n", "\n", " ⎛ 2 2 \n", " ⎛ 2 ρ⋅u₀ ⎞ ⎛ 2 2 ρ ⎞ ⎛ 2 ρ⋅u₁ ⎞ ⎜ 2 2 ρ⋅u₀ ρ⋅u₁ \n", " ⎜x⋅y , ────, ω⎟, ⎜x , ρ⋅u₀ + ─, ω⎟, ⎜x ⋅y, ────, ω⎟, ⎜x ⋅y , ───── + ───── +\n", " ⎝ 3 ⎠ ⎝ 3 ⎠ ⎝ 3 ⎠ ⎝ 3 3 \n", "\n", " ⎞⎤\n", " ρ ⎟⎥\n", " ─, ω⎟⎥\n", " 9 ⎠⎦" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from lbmpy.maxwellian_equilibrium import get_moments_of_continuous_maxwellian_equilibrium\n", "\n", "eq_moments = get_moments_of_continuous_maxwellian_equilibrium(moments, order=2, dim=2, \n", " c_s_sq=sp.Rational(1, 3))\n", "omega = sp.symbols(\"omega\")\n", "relaxation_info = [(moment, eq_value, omega) for moment, eq_value in zip(moments, eq_moments)]\n", "relaxation_info" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\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", "
MomentEq. Value Relaxation Rate
$1$$\\rho$$\\omega$
$y$$\\rho u_{1}$$\\omega$
$y^{2}$$\\rho u_{1}^{2} + \\frac{\\rho}{3}$$\\omega$
$x$$\\rho u_{0}$$\\omega$
$x y$$\\rho u_{0} u_{1}$$\\omega$
$x y^{2}$$\\frac{\\rho u_{0}}{3}$$\\omega$
$x^{2}$$\\rho u_{0}^{2} + \\frac{\\rho}{3}$$\\omega$
$x^{2} y$$\\frac{\\rho u_{1}}{3}$$\\omega$
$x^{2} y^{2}$$\\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}$$\\omega$
\n", " " ], "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from lbmpy.methods.creationfunctions import create_generic_mrt\n", "\n", "force_model = forcemodels.Guo(sp.symbols(\"F_:2\"))\n", "method = create_generic_mrt(d2q9, relaxation_info, compressible=False, force_model=force_model, cumulant=False)\n", "method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example of a update equation without simplifications" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Subexpressions:
$$vel0Term \\leftarrow f_{4} + f_{6} + f_{8}$$
$$vel1Term \\leftarrow f_{1} + f_{5}$$
$$\\rho \\leftarrow f_{0} + f_{2} + f_{3} + f_{7} + vel0Term + vel1Term$$
$$u_{0} \\leftarrow \\frac{F_{0}}{2} - f_{3} - f_{5} - f_{7} + vel0Term$$
$$u_{1} \\leftarrow \\frac{F_{1}}{2} - f_{2} + f_{6} - f_{7} - f_{8} + vel1Term$$
$$forceTerm_{0} \\leftarrow \\left(- \\frac{\\omega}{2} + 1\\right) \\left(- \\frac{4 F_{0} u_{0}}{3} - \\frac{4 F_{1} u_{1}}{3}\\right)$$
$$forceTerm_{1} \\leftarrow \\left(- \\frac{\\omega}{2} + 1\\right) \\left(- \\frac{F_{0} u_{0}}{3} + \\frac{F_{1} \\left(2 u_{1} + 1\\right)}{3}\\right)$$
$$forceTerm_{2} \\leftarrow \\left(- \\frac{\\omega}{2} + 1\\right) \\left(- \\frac{F_{0} u_{0}}{3} + \\frac{F_{1} \\left(2 u_{1} - 1\\right)}{3}\\right)$$
$$forceTerm_{3} \\leftarrow \\left(- \\frac{\\omega}{2} + 1\\right) \\left(\\frac{F_{0} \\left(2 u_{0} - 1\\right)}{3} - \\frac{F_{1} u_{1}}{3}\\right)$$
$$forceTerm_{4} \\leftarrow \\left(- \\frac{\\omega}{2} + 1\\right) \\left(\\frac{F_{0} \\left(2 u_{0} + 1\\right)}{3} - \\frac{F_{1} u_{1}}{3}\\right)$$
$$forceTerm_{5} \\leftarrow \\left(- \\frac{\\omega}{2} + 1\\right) \\left(\\frac{F_{0} \\left(2 u_{0} - 3 u_{1} - 1\\right)}{12} + \\frac{F_{1} \\left(- 3 u_{0} + 2 u_{1} + 1\\right)}{12}\\right)$$
$$forceTerm_{6} \\leftarrow \\left(- \\frac{\\omega}{2} + 1\\right) \\left(\\frac{F_{0} \\left(2 u_{0} + 3 u_{1} + 1\\right)}{12} + \\frac{F_{1} \\left(3 u_{0} + 2 u_{1} + 1\\right)}{12}\\right)$$
$$forceTerm_{7} \\leftarrow \\left(- \\frac{\\omega}{2} + 1\\right) \\left(\\frac{F_{0} \\left(2 u_{0} + 3 u_{1} - 1\\right)}{12} + \\frac{F_{1} \\left(3 u_{0} + 2 u_{1} - 1\\right)}{12}\\right)$$
$$forceTerm_{8} \\leftarrow \\left(- \\frac{\\omega}{2} + 1\\right) \\left(\\frac{F_{0} \\left(2 u_{0} - 3 u_{1} + 1\\right)}{12} + \\frac{F_{1} \\left(- 3 u_{0} + 2 u_{1} - 1\\right)}{12}\\right)$$
Main Assignments:
$$d_{0} \\leftarrow f_{0} + forceTerm_{0} + \\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right) - \\omega \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right) - \\omega \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right) + \\omega \\left(- f_{0} - f_{1} - f_{2} - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho\\right)$$
$$d_{1} \\leftarrow f_{1} + forceTerm_{1} - \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{2} + \\frac{\\omega \\left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{1}\\right)}{2} - \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{2} + \\frac{\\omega \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right)}{2}$$
$$d_{2} \\leftarrow f_{2} + forceTerm_{2} + \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{2} - \\frac{\\omega \\left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \\rho u_{1}\\right)}{2} - \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{2} + \\frac{\\omega \\left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{1}^{2} + \\frac{\\rho}{3}\\right)}{2}$$
$$d_{3} \\leftarrow f_{3} + forceTerm_{3} + \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{2} - \\frac{\\omega \\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0}\\right)}{2} - \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{2} + \\frac{\\omega \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right)}{2}$$
$$d_{4} \\leftarrow f_{4} + forceTerm_{4} - \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{2} + \\frac{\\omega \\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \\rho u_{0}\\right)}{2} - \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{2} + \\frac{\\omega \\left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \\rho u_{0}^{2} + \\frac{\\rho}{3}\\right)}{2}$$
$$d_{5} \\leftarrow f_{5} + forceTerm_{5} + \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{4} - \\frac{\\omega \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right)}{4} - \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{4} + \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4}$$
$$d_{6} \\leftarrow f_{6} + forceTerm_{6} + \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{4} + \\frac{\\omega \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right)}{4} + \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{4} + \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4}$$
$$d_{7} \\leftarrow f_{7} + forceTerm_{7} - \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{4} + \\frac{\\omega \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right)}{4} - \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{4} + \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4}$$
$$d_{8} \\leftarrow f_{8} + forceTerm_{8} - \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{4} - \\frac{\\omega \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right)}{4} + \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{4} + \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4}$$
" ], "text/plain": [ "Equation Collection for d_0,d_1,d_2,d_3,d_4,d_5,d_6,d_7,d_8" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "collision_rule = method.get_collision_rule()\n", "collision_rule" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generic simplification strategy - common subexpresssion elimination" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
NameRuntimeAddsMulsDivsTotal
OriginalTerm- 293 261 0 554
sympy_cse61.27 ms 114 67 0 181
" ], "text/plain": [ ".Report at 0x7fd077bef438>" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "generic_strategy = ps.simp.SimplificationStrategy()\n", "generic_strategy.add(ps.simp.sympy_cse)\n", "generic_strategy.create_simplification_report(collision_rule)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A custom simplification strategy for moment-based methods" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
NameRuntimeAddsMulsDivsTotal
OriginalTerm- 293 261 0 554
expand32.60 ms 116 227 0 343
replace_second_order_velocity_products14.12 ms 126 235 0 361
expand10.95 ms 118 227 0 345
factor_relaxation_rates15.49 ms 118 184 0 302
replace_density_and_velocity5.22 ms 118 184 0 302
replace_common_quadratic_and_constant_term21.72 ms 106 148 0 254
factor_density_after_factoring_relaxation_times15.14 ms 106 136 0 242
subexpression_substitution_in_main_assignments12.94 ms 102 132 0 234
add_subexpressions_for_divisions4.48 ms 102 132 0 234
cse_in_opposing_directions12.80 ms 102 116 0 218
sympy_cse42.44 ms 86 73 0 159
" ], "text/plain": [ ".Report at 0x7fd077c7d080>" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "simplification_strategy = create_simplification_strategy(method, cse_pdfs=True, cse_global=True)\n", "simplification_strategy.create_simplification_report(collision_rule)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Seeing the simplification in action" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Initial Version
$$d_{7} \\leftarrow f_{7} + forceTerm_{7} - \\frac{\\omega \\left(- f_{5} - f_{6} + f_{7} + f_{8} + \\frac{\\rho u_{1}}{3}\\right)}{4} + \\frac{\\omega \\left(f_{5} - f_{6} - f_{7} + f_{8} + \\rho u_{0} u_{1}\\right)}{4} - \\frac{\\omega \\left(f_{5} - f_{6} + f_{7} - f_{8} + \\frac{\\rho u_{0}}{3}\\right)}{4} + \\frac{\\omega \\left(- f_{5} - f_{6} - f_{7} - f_{8} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}\\right)}{4}$$
expand
$$d_{7} \\leftarrow - f_{7} \\omega + f_{7} + forceTerm_{7} + \\frac{\\omega \\rho u_{0}^{2}}{12} + \\frac{\\omega \\rho u_{0} u_{1}}{4} - \\frac{\\omega \\rho u_{0}}{12} + \\frac{\\omega \\rho u_{1}^{2}}{12} - \\frac{\\omega \\rho u_{1}}{12} + \\frac{\\omega \\rho}{36}$$
replace_second_order_velocity_products
$$d_{7} \\leftarrow - f_{7} \\omega + f_{7} + forceTerm_{7} + \\frac{\\omega \\rho u_{0}^{2}}{12} - \\frac{\\omega \\rho u_{0}}{12} + \\frac{\\omega \\rho u_{1}^{2}}{12} - \\frac{\\omega \\rho u_{1}}{12} + \\frac{\\omega \\rho \\left(u0Pu1^{2} - u_{0}^{2} - u_{1}^{2}\\right)}{8} + \\frac{\\omega \\rho}{36}$$
expand
$$d_{7} \\leftarrow - f_{7} \\omega + f_{7} + forceTerm_{7} + \\frac{\\omega \\rho u0Pu1^{2}}{8} - \\frac{\\omega \\rho u_{0}^{2}}{24} - \\frac{\\omega \\rho u_{0}}{12} - \\frac{\\omega \\rho u_{1}^{2}}{24} - \\frac{\\omega \\rho u_{1}}{12} + \\frac{\\omega \\rho}{36}$$
factor_relaxation_rates
$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\frac{\\rho u0Pu1^{2}}{8} - \\frac{\\rho u_{0}^{2}}{24} - \\frac{\\rho u_{0}}{12} - \\frac{\\rho u_{1}^{2}}{24} - \\frac{\\rho u_{1}}{12} + \\frac{\\rho}{36}\\right)$$
replace_density_and_velocity
$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\frac{\\rho u0Pu1^{2}}{8} - \\frac{\\rho u_{0}^{2}}{24} - \\frac{\\rho u_{0}}{12} - \\frac{\\rho u_{1}^{2}}{24} - \\frac{\\rho u_{1}}{12} + \\frac{\\rho}{36}\\right)$$
replace_common_quadratic_and_constant_term
$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\frac{f_{eq common}}{36} + \\frac{\\rho u0Pu1^{2}}{8} - \\frac{\\rho u_{0}}{12} - \\frac{\\rho u_{1}}{12}\\right)$$
factor_density_after_factoring_relaxation_times
$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\frac{f_{eq common}}{36} + \\rho \\left(\\frac{u0Pu1^{2}}{8} - \\frac{u_{0}}{12} - \\frac{u_{1}}{12}\\right)\\right)$$
subexpression_substitution_in_main_assignments
$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\frac{f_{eq common}}{36} + \\rho \\left(\\frac{u0Pu1^{2}}{8} - \\frac{u0Pu1}{12}\\right)\\right)$$
add_subexpressions_for_divisions
$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\frac{f_{eq common}}{36} + \\rho \\left(\\frac{u0Pu1^{2}}{8} - \\frac{u0Pu1}{12}\\right)\\right)$$
cse_in_opposing_directions
$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(- f_{7} + \\rho \\left(- \\xi_{72} + \\xi_{73}\\right) + \\xi_{71}\\right)$$
sympy_cse
$$d_{7} \\leftarrow f_{7} + forceTerm_{7} + \\omega \\left(\\rho \\left(- \\xi_{72} + \\xi_{73}\\right) + \\xi_{71} + \\xi_{76}\\right)$$
" ], "text/plain": [ ".IntermediateResults at 0x7fd077d6cef0>" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "simplification_strategy.show_intermediate_results(collision_rule, symbols=[sp.Symbol(\"d_7\")])" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "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.6.6" } }, "nbformat": 4, "nbformat_minor": 1 }