Add interactive 2D LDC notebook, fix material initialization

7 files changed, 191 insertions, 13 deletions

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(