From e096aca2a38141ff0f3e78f3ebadd7a58760f7a6 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Sat, 22 Jun 2019 19:59:36 +0200 Subject: Add interactive 2D LDC notebook, fix material initialization --- implosion.py | 2 +- ldc_2d.py | 16 ++-- ldc_2d_gl_interop.py | 1 - ldc_3d.py | 2 +- notebook/ldc_interactive.ipynb | 167 +++++++++++++++++++++++++++++++++++++++++ shell.nix | 5 ++ simulation.py | 11 +-- 7 files changed, 191 insertions(+), 13 deletions(-) create mode 100644 notebook/ldc_interactive.ipynb diff --git a/implosion.py b/implosion.py index 63c1950..b8cb046 100644 --- a/implosion.py +++ b/implosion.py @@ -2,8 +2,8 @@ import numpy import time import matplotlib -import matplotlib.pyplot as plt matplotlib.use('AGG') +import matplotlib.pyplot as plt from simulation import Lattice, Geometry from symbolic.generator import LBM diff --git a/ldc_2d.py b/ldc_2d.py index 40edf89..8268e10 100644 --- a/ldc_2d.py +++ b/ldc_2d.py @@ -1,15 +1,19 @@ import numpy import time +from string import Template import matplotlib -import matplotlib.pyplot as plt matplotlib.use('AGG') +import matplotlib.pyplot as plt from simulation import Lattice, Geometry from symbolic.generator import LBM import symbolic.D2Q9 as D2Q9 +lid_speed = 0.1 +relaxation_time = 0.52 + def MLUPS(cells, steps, time): return cells * steps / time * 1e-6 @@ -33,16 +37,18 @@ def cavity(geometry, x, y): else: return 1 -boundary = """ +boundary = Template(""" if ( m == 2 ) { u_0 = 0.0; u_1 = 0.0; } if ( m == 3 ) { - u_0 = 0.1; + u_0 = $lid_speed; u_1 = 0.0; } -""" +""").substitute({ + 'lid_speed': lid_speed +}) nUpdates = 100000 nStat = 5000 @@ -58,7 +64,7 @@ lattice = Lattice( geometry = Geometry(256, 256), moments = lbm.moments(optimize = False), - collide = lbm.bgk(f_eq = lbm.equilibrium(), tau = 0.52), + collide = lbm.bgk(f_eq = lbm.equilibrium(), tau = relaxation_time), boundary_src = boundary) diff --git a/ldc_2d_gl_interop.py b/ldc_2d_gl_interop.py index 26a8dea..1c3439f 100644 --- a/ldc_2d_gl_interop.py +++ b/ldc_2d_gl_interop.py @@ -1,5 +1,4 @@ import numpy -import time from string import Template from simulation import Lattice, Geometry diff --git a/ldc_3d.py b/ldc_3d.py index e059b31..34b5136 100644 --- a/ldc_3d.py +++ b/ldc_3d.py @@ -2,8 +2,8 @@ import numpy import time import matplotlib -import matplotlib.pyplot as plt matplotlib.use('AGG') +import matplotlib.pyplot as plt from simulation import Lattice, Geometry from symbolic.generator import LBM diff --git a/notebook/ldc_interactive.ipynb b/notebook/ldc_interactive.ipynb new file mode 100644 index 0000000..0f59727 --- /dev/null +++ b/notebook/ldc_interactive.ipynb @@ -0,0 +1,167 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "from simulation import Lattice, Geometry\n", + "from symbolic.generator import LBM\n", + "import symbolic.D2Q9 as D2Q9" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "lbm = LBM(D2Q9)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "def cavity(geometry, x, y):\n", + " if x == 1 or y == 1 or x == geometry.size_x-2:\n", + " return 2\n", + " elif y == geometry.size_y-2:\n", + " return 3\n", + " else:\n", + " return 1" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "boundary = \"\"\"\n", + " if ( m == 2 ) {\n", + " u_0 = 0.0;\n", + " u_1 = 0.0;\n", + " }\n", + " if ( m == 3 ) {\n", + " u_0 = 0.1;\n", + " u_1 = 0.0;\n", + " }\n", + "\"\"\"" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib\n", + "import matplotlib.pyplot as plt\n", + "import numpy" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [], + "source": [ + "def generate_moment_plot(lattice, m):\n", + " velocity = numpy.ndarray(shape=tuple(reversed(lattice.geometry.inner_size())))\n", + " for x, y in lattice.geometry.inner_cells():\n", + " velocity[y-1,x-1] = numpy.sqrt(m[1,lattice.gid(x,y)]**2 + m[2,lattice.gid(x,y)]**2)\n", + " plt.figure(figsize=(8,8))\n", + " plt.imshow(velocity, origin='lower', vmin=0.0, cmap=plt.get_cmap('seismic'))" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "def test(nX=64, nY=64, nSteps=1000, tau=0.6):\n", + " lattice = Lattice(\n", + " descriptor = D2Q9,\n", + " geometry = Geometry(nX, nY),\n", + "\n", + " moments = lbm.moments(optimize = False),\n", + " collide = lbm.bgk(f_eq = lbm.equilibrium(), tau = tau),\n", + "\n", + " boundary_src = boundary)\n", + " lattice.setup_geometry(cavity)\n", + " \n", + " for i in range(0,nSteps):\n", + " lattice.evolve()\n", + " generate_moment_plot(lattice, lattice.get_moments())" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [], + "source": [ + "import ipywidgets as widgets\n", + "from IPython.display import display\n", + "from ipywidgets import interactive" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "47a0362082b9460fa8d821fa95513742", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(IntSlider(value=64, description='nX', max=1024, min=32, step=32), IntSlider(value=64, de…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "interactive(test, nX=(32,1024,32), nY=(32,1024,32), nSteps=(0,100000,500), tau=(0.515,1.0,0.01))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "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.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/shell.nix b/shell.nix index 118835b..7c2ba9f 100644 --- a/shell.nix +++ b/shell.nix @@ -47,12 +47,17 @@ pkgs.stdenvNoCC.mkDerivation rec { matplotlib Mako pyevtk + + # jupyter, optional jupyterlab + ipywidgets ]); in [ local-python pkgs.opencl-info + # jupyter dependencies + pkgs.nodejs ]; shellHook = '' diff --git a/simulation.py b/simulation.py index 20b87c8..74c59ac 100644 --- a/simulation.py +++ b/simulation.py @@ -68,8 +68,6 @@ class Lattice: properties=[(cl.context_properties.PLATFORM, self.platform)]) self.queue = cl.CommandQueue(self.context) - self.np_material = numpy.ndarray(shape=(self.geometry.volume, 1), dtype=numpy.int32) - self.pop_size = descriptor.q * self.geometry.volume * self.float_type[0](0).nbytes self.moments_size = (descriptor.d+1) * self.geometry.volume * self.float_type[0](0).nbytes @@ -85,7 +83,7 @@ class Lattice: else: self.cl_moments = cl.Buffer(self.context, mf.WRITE_ONLY, size=self.moments_size) - self.cl_material = cl.Buffer(self.context, mf.READ_ONLY | mf.USE_HOST_PTR, hostbuf=self.np_material) + self.cl_material = cl.Buffer(self.context, mf.READ_ONLY, size=self.geometry.volume * numpy.int32(0).nbytes) self.build_kernel() @@ -105,10 +103,13 @@ class Lattice: return z * (self.geometry.size_x*self.geometry.size_y) + y * self.geometry.size_x + x; def setup_geometry(self, material_at): + material = numpy.ndarray(shape=(self.geometry.volume, 1), dtype=numpy.int32) + material[:,:] = 0 + for idx in self.geometry.inner_cells(): - self.np_material[self.gid(*idx)] = material_at(self.geometry, *idx) + material[self.gid(*idx)] = material_at(self.geometry, *idx) - cl.enqueue_copy(self.queue, self.cl_material, self.np_material).wait(); + cl.enqueue_copy(self.queue, self.cl_material, material).wait(); def build_kernel(self): program_src = Template(filename = str(Path(__file__).parent/'template/kernel.mako')).render( -- cgit v1.2.3