{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pystencils as ps\n", "import sympy as sp\n", "from lbmpy.session import *\n", "from lbmpy.moments import is_shear_moment, get_order" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Demo: Thermalized (Randomized) LBM\n", "\n", "This demo notebook shows how to modify the LB collision operator to account for microscopic fluctuations. This technique is also called thermalized or randomized LBM. In these methods a random fluctuation is added to the equilibrium moments. In this simple example we implement a thermalized model by writing our own simple linear congruent random number generator, the draws indepedent random numbers on each cell. A seed is stored for each cell since all cells are processed in parallel." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1) Creating the method definition\n", "\n", "In the thermalized LBM, the equilibrium values of moments are altered by small, random fluctuations. We achieve this by creating a new class to represent this randomized equilibrium, inheriting from ContinuousHydrodynamicMaxwellian, but overriding the computation of raw moments to insert our random numbers. The overridden equilibrium looks like this:\n", "python\n", "class ThermalizedEquilibrium(ContinuousHydrodynamicMaxwellian):\n", " def __init__(self, random_number_symbols, **kwargs):\n", " super().__init__(**kwargs)\n", " self.random_number_symbols = random_number_symbols\n", "\n", " def moment(self, exponent_tuple_or_polynomial):\n", " value = super().moment(exponent_tuple_or_polynomial)\n", " if is_shear_moment(exponent_tuple_or_polynomial, dim=self.dim):\n", " value += self.random_number_symbols[0] * 0.001\n", " elif get_order(exponent_tuple_or_polynomial) > 2:\n", " value += self.random_number_symbols[1] * 0.001\n", " return value\n", "\n", "We use the low-level function create_from_equilibrium to set up a method using this altered equilibrium. This requires also constructing a DensityVelocityComputation instance, and the full collision info dictionary from scratch." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
 \n", " Moment-Based Method\n", " Stencil: D2Q9 Zero-Centered Storage: ✗ Force Model: None
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
 \n", " Continuous Hydrodynamic Maxwellian Equilibrium\n", " \n", " $f (\\rho, \\left( u_{0}, \\ u_{1}\\right), \\left( v_{0}, \\ v_{1}\\right)) \n", " = \\frac{3 \\delta_{\\rho} e^{- \\frac{3 v_{0}^{2}}{2} - \\frac{3 v_{1}^{2}}{2}}}{2 \\pi} + \\frac{3 e^{- \\frac{3 \\left(- u_{0} + v_{0}\\right)^{2}}{2} - \\frac{3 \\left(- u_{1} + v_{1}\\right)^{2}}{2}}}{2 \\pi}$\n", " Compressible: ✗ Deviation Only: ✗ Order: ∞
\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", " \n", "
Relaxation Info
MomentEq. Value Relaxation Rate
$1$$\\rho$$1.8$
$x$$u_{0}$$1.8$
$y$$u_{1}$$1.8$
$x^{2}$$0.001 rand_{0} + \\frac{\\rho}{3} + u_{0}^{2}$$1.8$
$y^{2}$$0.001 rand_{0} + \\frac{\\rho}{3} + u_{1}^{2}$$1.8$
$x y$$0.001 rand_{0} + u_{0} u_{1}$$1.8$
$x^{2} y$$0.001 rand_{1} + u_{0}^{2} u_{1} + \\frac{u_{1}}{3}$$1.8$
$x y^{2}$$0.001 rand_{1} + u_{0} u_{1}^{2} + \\frac{u_{0}}{3}$$1.8$
$x^{2} y^{2}$$0.001 rand_{1} + \\frac{\\rho}{9} + u_{0}^{2} u_{1}^{2} + \\frac{u_{0}^{2}}{3} + \\frac{u_{1}^{2}}{3}$$1.8$
" ], "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from lbmpy.fluctuatinglb import ThermalizedEquilibrium\n", "from lbmpy.methods import DensityVelocityComputation, CollisionSpaceInfo, create_from_equilibrium\n", "from lbmpy.moments import get_default_moment_set_for_stencil\n", "\n", "random_number_symbols = sp.symbols(\"rand_:3\")\n", "\n", "stencil = LBStencil(Stencil.D2Q9)\n", "equilibrium = ThermalizedEquilibrium(random_number_symbols, dim=stencil.D, compressible=False, c_s_sq=sp.Rational(1,3))\n", "cqc = DensityVelocityComputation(stencil, False, False)\n", "\n", "# SRT-Type Collision Info\n", "moments = get_default_moment_set_for_stencil(stencil)\n", "r_rate = 1.8\n", "r_rate_dict = {m : r_rate for m in moments}\n", "c_space = CollisionSpaceInfo(CollisionSpace.RAW_MOMENTS)\n", "\n", "thermalized_method = create_from_equilibrium(stencil, equilibrium, cqc, r_rate_dict, collision_space_info=c_space)\n", "thermalized_method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2) Creating the kernel equations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we have to define rules how to compute the quantities $rand_0$ and $rand_1$. \n", "We do this using a linear congruent RNG on each cell. We pass in a seed array, and in each time step this seed array is updated with the new random numbers." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAACAgAAAA/CAYAAABns2jbAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2dSdLcxLqGC4cXAGYFB+/ANivAnt/BAVZg2AEEI4i4AwfsADy9E5odGCLO3M0O4KzA5t8B931VmfpTKjWZVaoqNU9GqDKVyvbJ71MqG6ne+/bbbx/vdrsXOtrmr+++++5+25NzCEAAAhCAAAQgAAEIQAACEIAABCAAAQhAAAJrI6C50O9Dnd7K/lDnX6+tjtQHApcgIN15X/l8puNLuR9eIk/ygMCaCAQd+ibU6aNgP5X/zZrqSV0gMFcC0rVVPBOqHq/F+EEH50/vJp5PEred71rnnEIAAhCAAAQgAAEIQAACEIAABCAAAQhAAAIQWB0BTaB+oUq9lf2DKyf7sY4XOtpzpqurOxWCwJQEpDNeiPBLiTbeKICBAATKCXwvXfoyRpP7R7m90MdLvREKNgTORED6tqZnwqfCdK+FqvpowHvxCwKq8HutAJxCAAIQgAAEIAABCEAAAhCAAAQgAAEIQAACEFg9Ac2N/q1KfiL7Tays3P/IfV/2X9EPGwIQyCMgvfm3QnqRkwXNPGSEgkBNIPQ/T2T/bk/Z/orAnzoeyl33U76GgQAEpiUgHVv1M2Go39M702IjNQhAAAIQgAAEIAABCEAAAhCAAAQgAAEIQAACyyEQFl78pnN7I4A/5dz1WdblVI6SQgACEIDAEgn46wGvllhwygyBJRPY0jNh+hcDS24zyg4BCEAAAhCAAAQgAAEIQAACEIAABCAAAQhA4BgCfZ9B91+wtj/Lekz6xIEABCAAAQhkE9Ai5U+twN4w8Jf8+XpACwynEJiYwGaeCdkgMLHkkBwEIAABCEAAAhCAAAQgAAEIQAACEIAABCAwbwJaZPGXAZ7r8Geb/6Ojy3hzQN9EcVd4/CCwOQItXfpF5/X/pm8OBhWGwAkE+nQp+PsvOx6ekDxRIQCBHgIt3ftPT7DVPRPyFwM9LY03BCAAAQhAAAIQgAAEIAABCEAAAhCAAAQgsE4Cmgx+o8OLLd4A8EdPLX2t/bcDPUHxhsA2CbR06cU2KVBrCJxOoEuX5OdNbN/reCj3zem5kAIEINAm0NK9zTwT8gWBtiRwDgEIQAACEIAABCAAAQhAAAIQgAAEIAABCKyegCaEH4dK/p/s/9XhhZj255vb5yEKFgQgEAkkuvR79MOGAATKCaS6JLf7pK9lP3FK4dw2G9fK0RIDAoMEpFebeybkCwKDIsFFCEAAAhCAAAQgAAEIQAACEIAABCAAAQhAYKUEvOji/3T2W5nPdMTJYS/A2P2bbBZiBAIDgRECqS6NBOUyBCAwQKDSJV3358x/9KF+6IEPub/W8U4HBgIQmJ5A2o9t4pnw7vQMSRECEIAABCAAAQhAAAIQgAAEIAABCEAAAhCAwOwJeBNA9cazFl9+0PG9D/m91XFfx1MdGAhAYJxAqkt+6/lLHfb7KOqUdWw8GUJAYPMEoi69Fgn/zY3t2kiPrFsYCEBgegJR93bur3Ss/pnw7vQMSRECEIAABCAAAQhAAAIQgAAEIAABCEAAAhCAQDcBT7rqihfhbbwQ708oN/5bWedeGPlGx1ThvGjpfP1FAKfpjQF+I9NviVVGefrtTAwEFkNgjrqkMlnH0KXFSBEFNYG56ZLKw0YARHMTBOamexG6yrX6fowNArG1sSEAAQhAAAIQgAAEIAABCEAAAhCAAAQgAIGzEtCE6wtl4A0Bb5yRbL+x5YX7ejEk+PnTyp8m4byY7w0D9YRtYbhfFfeh4lR/GSDb5bDhP9P3HPhdGIEgw+jSwtqN4s6PALo0vzahRNsggO5dt53vXDd7cocABCAAAQhAAAIQgAAEIAABCEAAAhCAAAS2QEATwX6L/7HsanNAqLMX/OtPKIcwXrz3p13TcN5IEBf1dwXh/CUCbw7wQmq1OUBuG7vfyK/x5YLqCj8QmDmBIP/o0szbieLNnwC6NP82ooTrJIDuXb9d+YLA9duAEkAAAhCAAAQgAAEIQAACEIAABCAAAQhAYDMENCnshX4v2v8u95NWxf3lAJv3de2rvXP3oewXOk/f9s8N568TOK2fQlrR8oaD3+IJNgSWSEByjS4tseEo8+wIoEuzaxIKtBEC6N71GpoNAtdjT84QgAAEIAABCEAAAhCAAAQgAAEIQAACENgMAU0C/6XjB1XYC/9eoN/5XEf9twHB/7cQzkH6TLXAnxHuMyWQbizYKY6/KuCvGdRfJEgz0fV/69xfF0i/OJAGOcmtdL/SYQ4YCBxFwLIZZAhdQpeOkiEi7QksQZdURvokBHZ1BJage4a+Zv1jg8Dq1IoKQQACEIAABCAAAQhAAAIQgAAEIAABCEBgngQ00erNAP7cvxf47fZi+TMd6af+X2aWfjCc0vRGAB/tjQDeNLDT9cbGgeD3hWwvvtabA+T2Vwhs3ur4UOfphobqQvtnKI6ueVPEjzq+bMfjHAK5BCQ/6BK6lCsuhBsgMGddUtnokwbajkvLJjBn3TPZtevfnWWLD6WHAAQgAAEIQAACEIAABCAAAQhAAAIQgAAE5k5Ak6zf6/gzllPurr8X8GUvzPsvBQ6M4vgtymhywzl8vdgfIvtvDd7Y7TR1+GsCdnszwRPZ9cYBub0481a2NzX4rX//1UF7w4G8b01mnF8VLm48uI2MCwIjBCw3OtClW07o0i0LXAUE5q5LKh99UkF7EnQ5BOaueya5Bf1bxQaBIEx+MPKOYx6sl3MfoKQdBJDnDih4QeCMBKRz/i/KL3S8PmM2JA2BixBAni+CmUwgUBEI+uYxiA9PSvrwBA4GAosjYNnVgTwvruUoMAQWR8BfDPgxLbXuPV7w/0n2TeLvuT2HrU24TzlutagfLoyGC+l6sb/aAOB48nOeD3S88rmMNwTEDQROs1FGndsv3TBg92PFqdPUeduMxlH8mA7PD216nI8RQJcS/UOXxsSF6wME5q5L9EkDjcelRROYu+4Z7ur17+6iRUiF1wNA3MVb/W+Xzv2A7p283gmMgcCiCCDPi2ouCrsCAtI5Twr5gcSGSZk9B34XSgB5XmjDUewlE/Biav1ZYLm9mODNZveXXCnKvlkCyPNmm56KQ+CiBJ4qN8/beS4vjr9u0v7UpdG5NwxUG5d0Gt+S9nnd75aEU9hPdTxXfP9Xu79M8LMOzxv63me/dEOAy1fnI7c3AbiscQOBnJW50a/Hk23/XWGcapOA0vmtSpUfCOQRQJcO9Q9dypMdQjUJzF2X6JOa7cXZegjMXfdMevX6t/gNAmok7+L4xK1lo4dwf57MGwQ+0nHwkL4PxS8EZksAeZ5t01CwNRJQP+G3T97ITj9TucaqUqcNEECeN9DIVHFuBPz1GX81wJORNn6Os9+DoI+VJz8QWAgB5HkhDUUxIbBkAqF/TL8A0Fsdha1eBOoNEC7khFMYL+Z7k0DbNPwUzpsBHDY13hzQZd7J817XBfmVxHmp8N6swAaBHph4HxKQrFZzGYdXDn1ydMSxcsIpDLp0iBifBROQTM9Wl1Q2+qQFyxZFHyYwZ91zybeif3eGm2neV0Mj+aG7vRHADyvexYuBwGIIIM+LaSoKCgEIQAACEIAABEzAbxfGTxNDBAJLJ4A8L70FKT8EIDAFAS/GeOE/x3hzQN9GgL74XXE8h+l8MRBYEwF0aU2tSV2uRQA9uhZ58oXA/tls9c+Edxbe0n0P4m64vl28C6/y8ouvhfCDzRvy86fijh4QnSPNK5BGnq8AnSwhAIFlETjH/f4caS6LKqWFAASOIaB7R/v/kr3A+pf8/RYKZoYEznG/P0ea10CHPF+DOnlCAAIzJOB5GS/Yp6Z9Hq85bPuFpXitJE7u5HNMGxsCSyBg/WjrQfs81sNh0aVIAxsCtwTQo1sWuCBwaQKb0L/F/cVAmIB5LmnwYvJ/eqSia0duT9D1egdWv6qGD+Xuewjb6Zo/h5qan+V3MLGZGy4mpPDeCPCl7Pq/28K1P+RnBYt52G3zcG/d/iqc2/nr4PNItgdOX8s/xg2XdtlpxgjXtlUH84my/Ivc6X/epcXbjDwHJoMyqzBRXt8Kkv/j1/8Z2DmQGEtP1y1fUT4thz53evFTwTq9NfJP28jt8lR+B7olv+wyxtQVp09fdrqWpQeF4bLrHcu4BjtwrvVO55HDGqo3aR0Cq0F9dIYKF+U95k8fEkmc0Q7tgywnjHNkNpFX+pDWs5TYlPQ18d452ncmTbR455DehWv+u5qD59nFV/yICgQei+1DcvUhoFn8OET1jTpdVWmL8pwjswoTn3ku2oco36wxiBuvpIxBfh2HcUiEcSE7MOc57kK8C7Px3EKco6qiqr28+c/jfj8rteei2ufHxHG6nXMaVWL8dBJAjzqxzMkTXZpTawyUBV0agHP9S+hRZht0ybH83J97TFqtKenc6wgYCOQS2IT+LXGDgB++veD9j+w/dPxPR4ta+Tf5cB1ufB5ovtPhm58HGp1GYX3NN0kvuFcLorJ97qO+YeaGU5y2cTpdgyWXzcYTEW4n/8/aM+XjAVdtQr4/yvZ/sVVGbk/KvLafjnQRNyvNkMwsLJXfbKIsv5C7Uf+kkKuWZ3Fw/XJl9rXCWlaq/+YLcaM8VDqfm14IZ9mvJ0Tl9sT+C9mf6qj//09ul9F5Wx6r/0GUbfn9r2y3YX2/kXu0jIrXZTr1RelZT0f1oCCc65JV765CLt1PnNp6t/QqTVp+8SnRR/qQSemXJYYs73kVyuzo/Tk3vRAu614awtKHlIn4LEOrLTv7EPn7fuhn1MFNubOs1ISFCrKe+0w32z4ktOfos1eCbg3jkLo6W5LnQpm9eB8SypfVf7gBFX60jHVDNx2MQ5o8zn6mtursT86eMRmMEnDb6HAf1TbP5PFYh9tupzB2/yY7nYf4RufxBRcHG4zjAME4vz/jCXYeAbFGj/JQXSWU20cHunQV+mWZoktlvC4ZGj3Kp90lx/K7UQpeQ/IzMn8NmI+TkCJgmdKx+n7szhJbWw3jB3Gb/9NhRe9qqOqh3YG2ZMTmRocXN73o+fNI3T0R4Dc904V2Lw7Vi50hfm64Ojul+VV9cuiwcn2g4z0d93V4gt3t2DaeaK0Xb33RYWU5rMuUmtw00zhXd6s+UZZ/l9vcFynProeOL44BqnhZMhvS919RpAv35uXz+q2a3PQUx+X9QuG9KSCaqAvfRI9ge7L7nsJWmwPsJ7fvMX64SPN2mqNlVJiGUVpD+pKrB7nhSurdKOdaTsS71ru11GmqeohNlj6G/HL7htxwdTVGdCL3fp+rE843N826jHNwrEWWXQ8d9CF7BvQhc1CugTK09U7nHof4WdYTD76HfhT8BlJZ56VQ/zWMQ0r6DzfmavqQILuLkmeVec19SNYYxEIoDoxDDGJBxrIbihvHoAsq/eqL6i8GuH+vjc49F/ChbH9x0OP3T3U8rQPs5yUb8wsZcWJ0vxRTz3NET+xxAmKMHo1jumYIdOma9AvyRpcKYF0+KHqUyXxAjv2S34vMZAgGgZTA6vXvblrbBbn98OzG8SRc7o7c2VRPZfbg/ZXsq21iUN5eFPXN8ZMUjPzrt/XtnxuulYbT9cKtj1OMH/T/VBm8mSBNywPof8vPE7DtzQyn5HeNuLUsh8wXJ8+h3O/L9nFO4wF4l868lP9XkgUvqqRyMlYWp+XwdRzH19EVz/rSlbf90ryLy6j8xvQlVw9yw5XUu4vFGvzaereoOklm6EPyWixXJ/JSm2eoRctygtT9h49zmuL780hhSu6l9CEjMBd2udY73Y+9eOCNgl5QdX9u482t6duDledcfuhDslpiC/2HQdSy7JMlyrPLLbPmPiS3/zCH4n4u3Lc8FvLRZXJ1ITdcSd/ZVZ61+TV0cMrKzeFeP2V9rpBW3CjW6M/FtXGelkvXLN+eu7Le1mYojgPpevUMLHvpc1t1nS/sOJsehfa5+tj7wjynzg5dmpro+dI7my7p/oYendZu6FE+vwM5lvz5OdVmcxsy0b19w5/4u3r9u3MioGtFt2JXSi1BH9vFe60yDuXricNrP/y7DF4M7ZsMiOXPDRfD2/5c6f6UehzpdhtXG0F64lcDqZ5rS/GuZdkFXqg8X4q1Wb3ryCzqUuzwO4Iceom1v9rgAXz9gCB3HMynXwWIctaV99uQ8qNgH1PGMX3J1YOscLn1PiS2Kp9a78TDG43c2furEZVb50NfdJgDCPqQvFbI0om8pGYbqpZll3Ch8nwpuMfcn3vLJtb0IYGOWBz0nb3g1nEh1Tt/qtDntuPhtwfHnq+vSYI+ZJz+FvoPU0hl2edLlGeX+xLm4n2I7iMlYxAzOKaMjEMuIT39ebR1sD9k+ZU53OvLSz2TGH7OU1H8AkLjKwJjxVP4Y9o0TjyPJc/1bgLHMO9OqdsXXermkuWLLmVhmkugc+oSenRCK6NHRfC65Phg00BRissOjO6d2H5b0L+7JYwExA+ub0Mc/0e935ZpTIDp3ANZL7RMFc4P5M7Xi4BO0w/qfkPnmY7KKM/eXbwxzFxsldU3qqFF70sV1QuaLodZfq7DbN2mv8rPjKPJDVeFV1wvrNWLqzGRtq1w3j1oWflQh9vY/yvvHde10bnfgugy1RtaHeFH04yJKe4sZdnlU9kWI8+R57ltMYkTZENZ3Ru6OHZNeVg3q8G53PUGF7njVwW60rf82nhhubiMijOqLwqTpQe54fbFvf1VvM5634aYzqW8Zqd3KpP7lsXoXGgv+hD6kL7noUXJ83R3l+GUpDfF9+fhFA+vBt2kDzlEM5mPGM+xD/lgsgpeICH6kD1kcRgcM+h61rNX2mRjabbCzk6WXT7VYVHynDI9p1tcrtKHKN+sMUhou+IyKv3Vj0NUxzjX8UBu/w2M5x383H+jw18qrMd8gWN77smffPeXC9O/mfPYyel+rOOpDo8RPfnqZzNPQL9Iw+u8MiFv637vvFYIerKlvOYy53RyXa6ZgDh+qcN/J+C5KstMjvG8QDqfNhhHYd0feQ7OcjFLo7JtUo/cGKo7ujSBVIojurSXp03qEno0gRLt5Qc9aqGUbLWf2w7WDEMU38v9l3DxZTA/w/kvtw/+2kd+fvnB1/1cZ+Mwj+XfeGasrsz8R2WmD5uojcRy1fp3N5eTQPh/OrwhoFrEDUJWTYTGNIKfF4f935MxnDtAbxioF2AKw/m/ix8qTvXALDv+X0j2Q7fiz8Ko7B64m1njs/5XKlycRHikcqVt87fOn+qIN8nccDvF8Y3ZExljgxun+YvC3bjuIZ4H6B60D7arrluenE9dZqchk52m0kCW98yW9HsvFLaSmZ6CWwaKTZApd5p+APB961VHItUDQYe/5dHGeReVUfnm6kuVQfozoAdpsN1QuHBtrN6N9E45UX7o3SkAFVcMLWf0IQX3+zbyIPf0IW0w6z8vuj+X4AgyNXYvpQ8pgdoRVpzpQzq4lHjRh9S03JcWj0MG+g8nnJ0msly3w5Ic1+xDcvoPsywqo+RwK+MQfyHBc1hetPc8lecqqgk+uf1MXU/2yt99edfcUzUprOs3um7jOYsqTbmf63jpc1+Q7fQ8n/KbjnpORO6+tB1tcP7DAUqM8prTeKGk6LMMK55uazPNMgpby1RWhKQ/ygx/jWCb0yNDDu0+l7H3Ndp90jzFE13af2236j/Egz5pUgnbRmLo0W07i0XWs5XCuQ9/EGJ6vcvPgvbz89p9Henzmu/5fu576PCy/bz8WofXCUr7d0W5nlHZeR6cGL+YrrYfu5vDKiiEFSR9w9uDIA+gKhPCePLOA640nBXW/pUpCGdBdvqGXyurzu32rp8b2YsxKq9vRh5AVjejaxZcZTFbG++kb9/gfpH/c/mnA9XBcAob28Jt3164rzJKfxTGO+tro3O/Dev8/IDkrxgMGcuEB9z1Ln4H1nlWmgrnmzuybGjrM/Ft/qKaSSZ8v6ruWXJ7p+Br2d7k5Am5aPx2yH/lZ9mpdEO2dTrKfnqPinG67LSMWfrSlYj8OvWgI2xvOJU/p94dSZZ7KS/0rhxbI0aQN/oQURGLrPt9A+DtSadO5KapcMjyLcu1udL7c3bdJBM591L6kGyihwHRu0MmpT5iyDgkQBOLY/uQzv7DyeamqXD0IaXCu5zw5+pDpuo/TDIt4+rHIUHf/gwidE+23+yOcw2eD6nnQeTvc+t419xT/eUuhUvntqzPvpbOSzgfG6dXmZG0D+a1FN6T0zZ+c+1DnY/Or1Sh9aOws7nXxzKtwRbXOOafvDrnTHuKwqp8lvMl6pF18DMdvtc9LGWhOOhSKbSM8OeU93OmnVG10SAq3+J0SWW2HvmFTxuX3yZ7TUPx0aM9s0l/zynr50x7SghBNvue29rPVn52s/kk1s92cHsNoHqO07nd/sJA/UU1+fk5z3pQr2vKPXujMqN7Z2olsV3lM2HWBoHIVBCsEFZA/+9qY3JFfl7ctXlf1+InOzwI9W7tdLE5N5wHRk6rHrg5cRkrdrpwV3n6R2GtzL4R5C7W1XFzHEr3Kx3pADAn2k5xXBczsRD582GyjjIesHbW/ajU9pst2lG9M8qfOXukwzukbLp4puEsD44T27aKVPjjPB4rHQ/cu/Lbyd/p++b8aWbavWkqjVnLsuunMs5OnlUmt0HsXNNmuBfK7M8rto11MrfN2nHj+bvo6LCrvOXvSZSTjMrpzSdRTz8I7p1tHf9S4v7EoDta39tehsPtZFnLLqPSOFpfFDdLD3LDqdw7he2st69NaZTPrPVO5Zudzpm/ykUfMi6Ivff7GFUcs3Qnhpfdm6bSmrUsuw4q4+zkObQBfQh9SKJm+c65690cdS7cC+hDxsWs934fGJb2H47Wm+bcZTnUmT7EIPYm+xk/RjjGllwcPI/LL2cM4uyyy6g0tzIO8YTuL6EtPH6rx6liULvD9dy5p1ehTZy2F0va6TgfG+t/NLlp70LbvJUdJ6sfy+35tPbcW0y7thVmrvf6uow4FklgiXpkPYzjDZe/yKBLRbgInE9gcbqkqnn+s+7n5PbzsNcC4ma73tqjR71ouDANgexnK2XnZ6j2pgGXwjqZbp59rvP4LO7rO8lxfK47WIvTtdmNlUKZeR6sWo+fEgJ3cwJL6OPOaC9yVw9a8vtBR7qb2f5WpLEF9Nxwnym9dGPBTmlbeT0Q88R8w+iaB7ouZz0Yk9tKYeOFw6zd10NxdM11/lFH3UFWqY/8KLwX9r3r1vHqHUsj0c52WWXxRIPTvxnIpPrvtMxwbhNv5qjZ96WrMG67e7L7dvG6jQ+Mwrt9He9gcFySpsLOXpZdeZVzlvKscnXKvvzdMVpmxvTf1Ss2SjfKbJd8RL9R+UszVppVRy/7TeovtzfG+D5V3aviNZdB7kb95RfvMZarrDIqXLa+xLyjrbi9ehDD2B4Kp2tF9U7TPdatPGevd4EZfUhGI4tVlHXrRJ+hDxm/H57lecgNMld5Vrka99AoPPKnD7nVq9inRTy2o5/vUfQhAiIOjENSCRlwixXjkMBHLIrHIYoz+OxVkqbCzv55yKhCnWf3TKRybbIPUb0HxyChzeKzWewvgtRXVvTbVB8ibtUYT7aft2ziZoH9WfM3a+4ptIVjVmnqvDFfJf/PdbQnobPSdqIyHlvWf0fp9HV4g0DvSxRVLP0ozEXv9crvn5g39vwIqH3em6JUSmdxehTKbD30+KLYKB66VExtvREkD5vVJbXqF6q/X3KMfZ37KPv5K8PtedSGEOg6etQgsu0TycMkepRQLHm2asztOw2VJz4b+sU/n3ue3M/LP/s8MQ7nZ+zGeoPOZ7lm43KrbOhe0oBbd0oesnTvbi4oC5jCWsisHHb7bfpnOjxgjaZSrHgyYA+GU5pWSh/tjQC+Aex0PXZOPvW5w/q/4Oq3lOW2shbtvs6M487Ru+jSzREuxqBR+J90eOHxDx19i+ODaUx80Qw9ydxn4s0vJ5zT+Vj1+7WVmG+wHsza3xMSZvZIxzsdbXPPHgpz8JAhPz/Y+39h0vatyi4/l7MoTcWZrSwHBrOXZ5fzCqZPFivZUXka94WM8nnn607yUH8pICNOO4hl3BM3N+FCThktr7n6UuenPMb0oAqbEW6KetflynWoXLPVO5Vt9jqnMtKH0Iek6tb5POQAS5DntCIXdOfcn0uKM8W9lD4kk7jkmj5kz4pxSFNmPA442zhEcpfz7MU4pNkm1ZnY5YzFj5Lnjuwu4TW3PqTdf5hBThm3OA7xHILnIuJ4rSEv8vc4wEfW3FOI/Lnsg3kL+fmeUW1kCek6eFbaCu97mcPGeRjHtbnR4fZu+/tawyiNi40XlFfWhGOjgJwsmcAi9GgqwOjSVCRJp4PAknTJ/Vn8snBHVYa90KNhPlw9joDkys9KWc9WIQc/X7XXIeP6kp+dbfx8bNN+tvNLqjFMFSDkzxpkRaOaf+R5MLBYsnVnrPASfC+Gx/+c2snd9fcCTsYDlvTTHHXSiuOBUjS54Ry+PQiqPgviC05Th5XcxrvY/Kmb1NivVmKFtftxEicNG92jcZJ0fDMqMorrG42/QOB8rm3MK/JLy+LNC94dFdmNhnNYHf7P9sbhdHTEa57YtfGgtetzRN54EvOsAvpHYT0Y9mJqjB+vWabehZOsNJXGEmTZVVqEPAf2l7S80SR22mm+ltn2mxrp9T635bPx+aAQMOZRy6Nkx/ebv3XUeh/cccNUzGO0jIoXdSJHX6p0FSdHD3aZ4bLr7cyVZl3nWMkSW/GXoHeL0DmxpA+hD4nq1/c85OuLkOdYkQvao/fnwrJk30ulu/QhhXBjcPqQ5thFPOJ4prhvpg+ppCprzOCQ4pX17KWgWWkqvSU8D7nq9CGmcGiu0odIbnL7D5d4tIxKb1PjkNCMHq8dfBo2XEut3LknxzmYu3BbhcTilwq+SRLPSbvvvu45j3tJWoNOlWNO44XBsnLxcgQkF33ylVuIpehRbn1Gw6FLo4g2GWBLuqS6+hnXY95ovI8wo40AABaNSURBVGHAG+7cz2SZEHYu6yBZZSbQZQhINk7tl3KerWJl6rAhX79w86XcqXzvdN4O576vvYF0EWMl1YXnwdj62KME7oyG2A9+GovvEjIPftodhRXEilMbhXtfh+OmncdoOMWxgnoCrF7ADnl6oibuXvNunai4XvhPF/MczzeaeF3Oyjhdp3FgFL8kTjU5d5BIhofy+UnB4uAxI8ZJQeKGjYMBpcrhQbInCNwelZHbzHyTfLr3qSbHssLF8C3b6flIjR8M2vL0VQgQd3BVpwrnNvFERyVHjhcP+aU38tw0LZ/tvOcmyyribjHy7MJObIZk1rrzTjJQ64/cBzLbKk9vegrnTSeNjj6k7TRT+XKSlsV3diTGsulw9f1N7mPKGJPs0ped0szSg9xwyiy73krTZfLGiOpN2VjQQnsJercYnVNb0IfQh/ge2Pc8ZPVcjDwX3ktygvfe84Pu0If0PEuJT1Zfo0agD9lLIuOQcY0c0ses8YXkMitcT1G6nquyxgwF+uCss9JUuCU8D7k+9CEdi7FX7EN8bx4dg7jhjiyjo9p06YvTzOobcsMpn4v2IaFcrkNjzOcKR6MwN3Jnzz2FNM3r55hGsJ1PtWldYfwFjZ9L026lF089l+P8so3yveR4IbtcBLwOAcmD5efoMb3iW7aXrkdHwUeXjsK22khb1iXV3fMPnod4WNrA6FEpsfWHP0WXFLfouU00PZ58lFD1XP4vQS6jd7W5U37u66J5Hhxx42f0X8xYKdSxXkOJFcCGQJvA3bZHx7kXiy38HuTEgYnfMG/8B6CFTocXcr3gHL844POjwikNLxY/V3wvHnuCyQMwvzHnty/sVy30ym3l9c0hNbGcqZ/dHlwfLJaHQCVxXiqOy+KbzDHGby03bijHJNIXR2n7ZmfjiSgbf77RmyX8/3UeLFZGbk9ummdcNDebT3ReL3g6oM6zwlWJ7sM7vXhTdT1dnpey/X+x3m3YztPt8i/5t9vRA3mnY9lrm7qMBWnOWpZdQdVlifLcbpvic9U7S2aVsB9GLT8fy36rw3aXzI6mpzR8z7J8RvlXUpW8Wd7rDUf21Lll1/ehb2THe4XL0QjnsDJZZdwHrdLu1ZcQJksPFDYrnMpcUm/f633v8L3cR1tHYzWG7Fnrneq0RJ2jD7nVW/db9CFBAxcqz0P3j6xrqvfoPT8klHV/zklPYUrupfQh9CEWQd+rGIcEZZQOZY0vcsOFZHcK3/tcpWu545CsZyrnWZDmrJ+HQl2W+EwUm/5oW2042z5EZSvpP8wgq5+LsJR+r76EMLm6kBVO+ZX0nVOMQyzT1vuuMVvEYDtr7ilEiGnW8xHB3/Ms/vKh5y68ITFez027b5zlsWf7xZeQ5aB11vHCYM5cnBUByeKpuhRlfgl6dA726NI5qC4wza3qkurte4DXeh6awZFNhx4dCW6N0SbQpdxnK+PzGMzri/6atZ+pvEZWr43pfBfK80ROz/O/lO35f4/d/QxZy7zc1oX6XG4bp9llNjH276o4fssj8N63337rRWQvHC/yPyNUbpf/a9lW5MrI/UAOv/Ha+G9x+f8tv2eyf6gCJj8lcRT2IM8kKZwQOJpAl2yVyGbMuCROV54xnVxbaXhH2keyD3QrNw3CzZdAaF9/caT9IDTfQmeWTHU6uJ/Ljz4kkx/B5kVgqfKsctOHzEuUJi1NaF/6EMYhk8oViU1PgD5keqakeDqBNfchbTqqq+erGpvf5feP/O7LPmaTQDsLzk8koHbw4oDHj7/IffGxsfJ8oCNuPqlqo3MvTtyTPSgjuu7n7VU+j1UgWj+hvl7s6fqL01ZoTi9NQO2CLl0a+hH5hXbymkv18mc438kevN8ckRVRjiQQ2oR+6Uh+fdHE1euKr2TXLz7LzfxxH7AN+i9Z92JzqQ4eezy9Gz0WbPthuP1g3j6P1XPYvk6sJI53AWEgcA4CltG2LLbPY74OOxd5dhn7yhnLi71cAn4b5tgvpsy91tajtuy2z2MdHHYuOhfLhA2BlMBS5dk616d3af1wL5MAfci+3ehDlim/Wyo1fciWWns5dV1zH9JuhWfy8ORztQAcJqL91mXf+KMd/+rnKqs3Wvtv+OoJdRdK514M9F9M2DzS4Tk1Lzo1Frt9UX7pF/785Z2n8ht8TtT1znydnk1JmgrrN2VT87P8Yjmdj8vnv5dJw0S339b/IJ7IPXW9/1CavlfH8tht46+HjJkt6dIYi9lfVzt3yvQZZKrBoi9fB9K1It1U+El0Selk61FBOdGlRssfngTubnPfqy2PNr63x3t55THnn1Dui/dJZtKXd7iWrUtKZ0iPnFx2v6S0snVJYXPLuGhdCu30jTjWzxrBz1z99YHULHWslNbhIu7AEN0Lz4TiMTfda8jBncbZMk88WIoPxVUNBN1+NzoMv23ig3TDvzCO013MIK1RUU7mTmCR8iz98U70xid65g6a8uURULv6/uq/c1irWaTOrbUxqNfJBBYpz/QhJ7f7bBOgDzloGsYhB0jwmBEB+pAZNQZFqSbX1z4OaTSz+kx/je9D2X7r2X+r6U/otienG3FmeOK/7LiXlkt18fyZF9SriWLZXsx2f/habm+IqIzc/ks7/13pn7LjpLI3TfxX511ze/uI+9+DfO2teNlpOg8dfmPQX1j1gpgXwZxv/BsSOau/OPSn9j330T58D63by+npfOp6v1Oanuv0woX1w5v4Rz/7rbJsRpfMXYcX1LzgU7l1bn1amjmQaddNlZhaptpcuvLN1iMn5nLqmESXcusc8i0pJ7rUbvnDc7eh79G24/GF2sT3oKWYLnm+hB6ZT1fe2TKaqUfO52Mdo/1Sri4pXHYZnbnM0nXJfUXjWUTnz3V0bWJc5FhJdbmG6ZJ/dG//fDvVc/BRutclDHe7PJfkpxvXm3CTaxd7cPe1b3iK8I3sdOfbYJwkAwu0By4YCExKAHmeFCeJTUPAA4DV/nUEOjeNkJDKPAggz/NoB0rRIEAfIhzSTU861G+B6pxxSENMOJkDAfqQObQCZWgRWHUf0qprdSo9TOenuoLM1k9l71uE9WKt3zqtjeup4wt5eAI5vnHvCXl/Kr8ee8rt+b5X8vebjPXfispdG13vy9dhStJ0Wfy1AC+0ROP+2gsCtdH1g3LIz3OEO9npV/fOUW/z8MaRUrMZXRIft9di9cgNqzr0yfQ5ZKqWpYF8S/TI6U2pS7l1dr4l5USXTGzASB7ivXkg1HwvDchzrkyVyFMDxEDeJWlm6ZEzVn45/dK56r10XfLzhY37yQ9lu9/3M0r6LODr5uy6Vv195XH7M7ieqDibGvurvlfpw0Ib9eWN7t3K64HrzoHPMj3+aiuozsd2X1uhrfz+H67KZMSJQX3jTR/8oz82BKYggDxPQZE0JiEQ7ouTpDXjRNC5GTcORSsmgDwXIyPCuQjQh/S+Bco45FxCR7qnEqAPOZUg8ScjsJE+ZDJe10xIbeU32m/C0S6KN8n5qwCeIE+NJ9/9pmKcbPfcXGMxPgR+I/txR/yd/IbydfSsNJWOwzmtxlcR5f/EhxMK5mV0tGx/9aGxCULXz1bvVt6jpypbveliNDABrkpAbTUk02eTqZF8s/TI4JTO1LqUW2dnn11OBz7GqH7o0jHgLhxnRJ5zZeooeRrJOyvNAj0y2dx+6az1Lm3iueiSyuEvIn/q8ujwxgB/wehgc0BSP8ZKCYy2U+yu0oe5HCN5o3vtxkrO7ybuJTvjLqjGLlEJRuM8raCueZDxgWwLSG2G4jiQrleDGtldA5c6HRwQOIEA8nwCPKJC4AgC6NwR0IgyWwLI82ybhoKtlAA6t9KG3Wi1kOeNNjzVhsCJBD73XJqO9iK5k/VE+wNd8waCLuNNAnHzgD+X2jZvg8cj2e1J+958C9N0uW8GylgVQdcPXhSSn++bfnuwbc5Z73ZenK+HQK9Mq4rnlKnOfAv1yK0wtS6N1tmZHlFOR8Osl0CnPIfqjsrUifLUmXdhmll65Poo3dx+6dz1Xq80NWvGWKnJo33WKf8h0LllsDNvdK/dRIfnq9ggoIaOu338P0fZC/cK691TFs4SE28EJXEIC4FsAshzNioCQmASAujcJBhJZCYEkOeZNATF2AwBdG4zTb2JiiLPm2hmKgmBSQnovuHPucZP9B6kret9n8T3W2Y7XffLO7Zt3fNPy/iTvzbxSwPVicKP5esFf4fNSdObD/xWoMv0uQ5vSriv41f59c4Z6prL5M0PBy8nye9c9fZfM3hDhbk4/2fKq2IoN2bBBNSOYzJ9LpnqzVdlKtEj059Ulwr0qLScO6WNLi1YX/qKrnbtlWfH0fWz6FFIuzdv5Vsio0fpUShDZ7905npvRpfEkTVIC1qHEZte+XfwM8tgb97KF93raK/U6056smS3Gtu7q/wZkLjzOKc63lDQt4v5IL7CVv+RJjt7E8JBInhAIIOAZAx5zuBEEAhMRQCdm4ok6cyBAPI8h1agDFsigM5tqbXXX1fkef1tTA0hMBUB3S+8EOGJ16I5MoX3QrzjpgvrfgvSfm3jsDb1XF9BvrlpxrQfKe2vdfhTw56T8QaBf1e5d//4BSIfWUZpnVRvZeJy/hLKZ3Y+XuvcLz9hFkxAbThnXcrVI7fA2XWpR4+cd2k50SVTW5G5lh4ZYWbeuTJ6rB65KNn9Uo8u5ZbRedlsrl8SN9Zs9m1f/2bKfx0+OqaQwcy8c+V6k7p3JzbIGmwJRDq4GK2Swjf+X2w0wv5BvHf3cEZ8gkAgmwDynI2KgBCYhAA6NwlGEpkJAeR5Jg1BMTZDAJ3bTFNvoqLI8yaamUpCYAoCfkmndF7N+f6q4zfFTf/P+6kvyK9e7JbbC+o39pdJNyHk5juapvKIk8H+EkC7Lr8o3+dJmKog/pGfF3Qfyy6ZIzyp3srriY7IYye3mTj/3i846BpmGQRyZbpdm5NkSonl5DuqRy6U5PFSutRVZxchq5yhrOiSQazP5MhzV627ZCpbnkKCOXmPpnmsHrkMilvaL51cb+W5SV1SvVmDbGpSjvw3Y+zPTpZBJZOTN7rXRT/4rWqDgOskBa0flgfqfdSlc6Z9VIGItHoC55S5c6a9+oahgqslcE69OGfaq20QKnYSgXPK3DnTPqnSRIbAFQmcUy/OmfYVkZH1jAmcU+bOmfaMkVI0CKyKgPTYX9gsXphWPMfx5/wbn3kO94V/6dqncn+lw29BerHjpQ6baoOA/LPzzU1zn3xjA0Lw2r2Ww4ue/txz2/gNwnTTQvt641xlOanejcSaJy6Dv45qVpgFEiiR6bR6p8pUbr4K53n2Ud1MytalF5PoUl+dnfcR5UyKXDnRpTaRBZ2r/bP7hrRafTJVIk+5eZekqTKW6pGrld0vTVHvlGPLvQldCu3Zqvo0p+dMe5oS3qaisqJ7C9e9u7fNiQsCEIAABCAAAQhAAAIQgAAEIAABCEAAAhCAQDcBTQZ7Mfp92V0LGN2R5Bsmke/JftIVSP5eiPQCR23kFz/h700FxflmpBn/m3boZaOuxXf/9UBW/VUGT54fXW/DUBovQhoPfd5hvJEBszACatdimXYVFe8kmSrNV+EHdTOU6ay6NFbnWAbZvfeQEAZdMogVmVJ5jlUfk6lMuS/S4bE0fV2Hi1jaJzlOVr+k9E+6fzgjG6WDLu1RbPZXMlAk/xHUFDJYmrfCD/Zjvq7DRdyc7rFBIEomNgQgAAEIQAACEIAABCAAAQhAAAIQgAAEIDBEwBPCH2si1Z+GTY3/EsBvs9vfC/r1J3jl9sLFfdn1lwPkdjo72UML7U7zd4XxxK3f5C/K1+l3mDrNcM2f6a/K0hHWXo3yqRzvh/BvesLX3gp7cr1DYq77uzrhW8c9O5XPaFluo+CaEYEl61Jbj4z1LLp0gh65TO1yokumsi5zFT0KCIvz7kDfltEiPXJ60pGsfukEXWqX0dmiS6awbVMs/xPKYHHeHU3VlutN6h4bBDokAy8IQAACEIAABCAAAQhAAAIQgAAEIAABCECgSUCTu55A9dEw8v9bHl7MrzcBOIDOPQHrhf16w4D9Zbx4/pMdumb3cx3/krt6e0u2Fzwe66jemtd5ab6jaSptG3/+v73Zwf7O1xsT2nX1oohN14L9/op+FW+SeocEf1J6bX6+ZD7t8oUoWHMnEGTroP3kPxtdUlly9ci4J9cl5T+qR864oJzokoGtyKjtS/uGUZnKlaeSvHPTVNOU6pFbc7RfUv6T1TuID7oUQGzVKpF/M5pSBkvyVtjcfmyTuscGga1qMPWGAAQgAAEIQAACEIAABCAAAQhAAAIQgMA0BLyg76M2mpT1G15efPfGAU+8puax/H4IHg7XXnB3vC8VZuzt+IN8S9JU+r/pcPm+11Etwst2mp/peBrSSq1Yx5vUM3Ur/tT1/lFp+qg/ny73VyHPxoaMtBy4F0vgQKbPIFNdcA7yVaBs3VQZJ9Wlgjq7LrnlRJe6Wn6dfgfyXCBTufLUR+4g71wZPUKPXAbnZ3Ozt5q/Z6o3utTEzNktgQP5P5MM3uZ46zrIW5ey9Hmruvfet99+692mLwTgvVuOuCAAAQhAAAIQgAAEIAABCEAAAhCAAAQgAAEI9BPQfKIX/j356vlFm990vJT/Dzr+lNvXuswbXa++DuCLcn8fAsWFjl/ld/B2dUxI13rzTcJkpxnyj3nfUxrP5HewOUF+rs9rHU/ldl0PjPwnr3fIN35FwOXzhoqv5X9zUAA8FklAbdkr0+eQqQhpKF+H0fVsPUrCn6xLJXVO8rUz5t15D1G61mF0yaRWaNS+V9EjoxzKO1zP1iWl5bBRlnv7pJDuYL+ktCbvk5J80SXDwAzK/7lkMGJX+r167zC6ju5FWMEWE3+t6CkbBFpgOIUABCAAAQhAAAIQgAAEIAABCEAAAhCAAAQgAAEIQAACEIAABCAAAQisiUDcIHBnTZWiLhCAAAQgAAEIQAACEIAABCAAAQhAAAIQgAAEIAABCEAAAhCAAAQgAAEIdBNgg0A3F3whAAEIQAACEIAABCAAAQhAAAIQgAAEIAABCEAAAhCAAAQgAAEIQAACqyLABoFVNSeVgQAEIAABCEAAAhCAAAQgAAEIQAACEIAABCAAAQhAAAIQgAAEIAABCHQTYINANxd8IQABCEAAAhCAAAQgAAEIQAACEIAABCAAAQhAAAIQgAAEIAABCEAAAqsiwAaBVTUnlYEABCAAAQhAAAIQgAAEIAABCEAAAhCAAAQgAAEIQAACEIAABCAAAQh0E2CDQDcXfCEAAQhAAAIQgAAEIAABCEAAAhCAAAQgAAEIQAACEIAABCAAAQhAAAKrIsAGgVU1J5WBAAQgAAEIQAACEIAABCAAAQhAAAIQgAAEIAABCEAAAhCAAAQgAAEIdBNgg0A3F3whAAEIQAACEIAABCAAAQhAAAIQgAAEIAABCEAAAhCAAAQgAAEIQAACqyLABoFVNSeVgQAEIAABCEAAAhCAAAQgAAEIQAACEIAABCAAAQhAAAIQgAAEIAABCHQTuBu9v/vuu3+iO9h/ye9+y49TCEAAAhCAAAQgAAEIQAACEIAABCAAAQhAAAIQgAAEIAABCEAAAhCAAARmSkDr/K9VtAddxfMGgVc6nnRcfNfhhxcEIAABCEAAAhCAAAQgAAEIQAACEIAABCAAAQhAAAIQgAAEIAABCEAAAvMl8FRFu9dRvFf/DyUVdn4TxG9hAAAAAElFTkSuQmCC", "text/latex": [ "$\\displaystyle \\left[ {seed}_{(0,0)}^{0} \\leftarrow 1664525 {seed}_{(0,0)}^{0} + 1013904223, \\ {seed}_{(0,0)}^{1} \\leftarrow 1664525 {seed}_{(0,0)}^{1} + 1013904223, \\ {seed}_{(0,0)}^{2} \\leftarrow 1664525 {seed}_{(0,0)}^{2} + 1013904223, \\ rand_{0} \\leftarrow \\frac{{seed}_{(0,0)}^{0}}{4294967295}, \\ rand_{1} \\leftarrow \\frac{{seed}_{(0,0)}^{1}}{4294967295}, \\ rand_{2} \\leftarrow \\frac{{seed}_{(0,0)}^{2}}{4294967295}, \\ {dbg}_{(0,0)} \\leftarrow \\frac{{seed}_{(0,0)}^{0}}{4294967295}\\right]$" ], "text/plain": [ "⎡ \n", "⎢seed_C__0 := 1664525⋅seed_C__0 + 1013904223, seed_C__1 := 1664525⋅seed_C__1 +\n", "⎣ \n", "\n", " seed_C__0 \n", " 1013904223, seed_C__2 := 1664525⋅seed_C__2 + 1013904223, rand₀ := ──────────,\n", " 4294967295 \n", "\n", " seed_C__1 seed_C__2 seed_C__0 ⎤\n", " rand₁ := ──────────, rand₂ := ──────────, dbg_C := ──────────⎥\n", " 4294967295 4294967295 4294967295⎦" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dh = ps.create_data_handling(domain_size=(80, 80))\n", "\n", "seed_type = np.uint32\n", "max_seed_type = np.iinfo(seed_type).max\n", "\n", "# Initialize the seed array\n", "seedField = dh.add_array('seed', dtype=seed_type, values_per_cell=len(random_number_symbols))\n", "for b in dh.iterate():\n", " random_field = np.random.randint(0, high=max_seed_type, dtype=seed_type, size=b['seed'].shape)\n", " np.copyto(b['seed'], random_field)\n", " \n", "debug_output = dh.add_array('dbg')\n", "\n", "linear_congruent_rng_eqs = [ps.Assignment(seedField(i), seedField(i) * 1664525 + 1013904223) \n", " for i, _ in enumerate(random_number_symbols)]\n", "floatEqs = [ps.Assignment(ps.TypedSymbol(s.name, np.float64), seedField(i) / max_seed_type)\n", " for i, s in enumerate(random_number_symbols)]\n", " \n", "rng_eqs = linear_congruent_rng_eqs + floatEqs + [ps.Assignment(debug_output.center, seedField(0) / max_seed_type)]\n", "rng_eqs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These assignments are then added to the LB collision rule." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "collision_rule = create_lb_collision_rule(lb_method=thermalized_method)\n", "collision_rule.subexpressions = rng_eqs + collision_rule.subexpressions" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Subexpressions:
 $${seed}_{(0,0)}^{0} \\leftarrow 1664525 {seed}_{(0,0)}^{0} + 1013904223$$ $${seed}_{(0,0)}^{1} \\leftarrow 1664525 {seed}_{(0,0)}^{1} + 1013904223$$ $${seed}_{(0,0)}^{2} \\leftarrow 1664525 {seed}_{(0,0)}^{2} + 1013904223$$ $$rand_{0} \\leftarrow \\frac{{seed}_{(0,0)}^{0}}{4294967295}$$ $$rand_{1} \\leftarrow \\frac{{seed}_{(0,0)}^{1}}{4294967295}$$ $$rand_{2} \\leftarrow \\frac{{seed}_{(0,0)}^{2}}{4294967295}$$ $${dbg}_{(0,0)} \\leftarrow \\frac{{seed}_{(0,0)}^{0}}{4294967295}$$ $$partial_{m m1 e 0} \\leftarrow f_{3} + f_{5} + f_{7}$$ $$partial_{m 0 e 0} \\leftarrow f_{0} + f_{1} + f_{2}$$ $$partial_{m 1 e 0} \\leftarrow f_{4} + f_{6} + f_{8}$$ $$partial_{m m1 e 1} \\leftarrow f_{5} - f_{7}$$ $$partial_{m 0 e 1} \\leftarrow f_{1} - f_{2}$$ $$partial_{m 1 e 1} \\leftarrow f_{6} - f_{8}$$ $$partial_{m m1 e 2} \\leftarrow f_{5} + f_{7}$$ $$partial_{m 0 e 2} \\leftarrow f_{1} + f_{2}$$ $$partial_{m 1 e 2} \\leftarrow f_{6} + f_{8}$$ $$m_{00} \\leftarrow partial_{m 0 e 0} + partial_{m 1 e 0} + partial_{m m1 e 0}$$ $$m_{10} \\leftarrow partial_{m 1 e 0} - partial_{m m1 e 0}$$ $$m_{01} \\leftarrow partial_{m 0 e 1} + partial_{m 1 e 1} + partial_{m m1 e 1}$$ $$m_{20} \\leftarrow partial_{m 1 e 0} + partial_{m m1 e 0}$$ $$m_{02} \\leftarrow partial_{m 0 e 2} + partial_{m 1 e 2} + partial_{m m1 e 2}$$ $$m_{11} \\leftarrow partial_{m 1 e 1} - partial_{m m1 e 1}$$ $$m_{21} \\leftarrow partial_{m 1 e 1} + partial_{m m1 e 1}$$ $$m_{12} \\leftarrow partial_{m 1 e 2} - partial_{m m1 e 2}$$ $$m_{22} \\leftarrow partial_{m 1 e 2} + partial_{m m1 e 2}$$ $$M_{post 0} \\leftarrow 1.0 m_{00}$$ $$M_{post 1} \\leftarrow 1.0 m_{10}$$ $$M_{post 2} \\leftarrow 1.0 m_{01}$$ $$M_{post 3} \\leftarrow 0.6 m_{00} + 1.8 m_{10}^{2} - 0.8 m_{20} + 0.0018 rand_{0}$$ $$M_{post 4} \\leftarrow 0.6 m_{00} + 1.8 m_{01}^{2} - 0.8 m_{02} + 0.0018 rand_{0}$$ $$M_{post 5} \\leftarrow 1.8 m_{01} m_{10} - 0.8 m_{11} + 0.0018 rand_{0}$$ $$M_{post 6} \\leftarrow 1.8 m_{01} m_{10}^{2} + 0.6 m_{01} - 0.8 m_{21} + 0.0018 rand_{1}$$ $$M_{post 7} \\leftarrow 1.8 m_{01}^{2} m_{10} + 0.6 m_{10} - 0.8 m_{12} + 0.0018 rand_{1}$$ $$M_{post 8} \\leftarrow 0.2 m_{00} + 1.8 m_{01}^{2} m_{10}^{2} + 0.6 m_{01}^{2} + 0.6 m_{10}^{2} - 0.8 m_{22} + 0.0018 rand_{1}$$ $$sub_{k to f 10} \\leftarrow \\frac{M_{post 4}}{2} - \\frac{M_{post 8}}{2}$$ $$sub_{k to f 11} \\leftarrow \\frac{M_{post 2}}{2} - \\frac{M_{post 6}}{2}$$ $$sub_{k to f 12} \\leftarrow \\frac{M_{post 3}}{2} - \\frac{M_{post 8}}{2}$$ $$sub_{k to f 13} \\leftarrow - \\frac{M_{post 1}}{2} + \\frac{M_{post 7}}{2}$$ $$sub_{k to f 14} \\leftarrow - \\frac{M_{post 5}}{4} + \\frac{M_{post 8}}{4}$$ $$sub_{k to f 15} \\leftarrow \\frac{M_{post 6}}{4} - \\frac{M_{post 7}}{4}$$ $$sub_{k to f 16} \\leftarrow \\frac{M_{post 5}}{4} + \\frac{M_{post 8}}{4}$$ $$sub_{k to f 17} \\leftarrow \\frac{M_{post 6}}{4} + \\frac{M_{post 7}}{4}$$
Main Assignments:
 $$d_{0} \\leftarrow M_{post 0} - M_{post 3} - M_{post 4} + M_{post 8}$$ $$d_{1} \\leftarrow sub_{k to f 10} + sub_{k to f 11}$$ $$d_{2} \\leftarrow sub_{k to f 10} - sub_{k to f 11}$$ $$d_{3} \\leftarrow sub_{k to f 12} + sub_{k to f 13}$$ $$d_{4} \\leftarrow sub_{k to f 12} - sub_{k to f 13}$$ $$d_{5} \\leftarrow sub_{k to f 14} + sub_{k to f 15}$$ $$d_{6} \\leftarrow sub_{k to f 16} + sub_{k to f 17}$$ $$d_{7} \\leftarrow sub_{k to f 16} - sub_{k to f 17}$$ $$d_{8} \\leftarrow sub_{k to f 14} - sub_{k to f 15}$$