diff options
| -rw-r--r-- | implosion.py | 2 | ||||
| -rw-r--r-- | ldc_2d.py | 16 | ||||
| -rw-r--r-- | ldc_2d_gl_interop.py | 1 | ||||
| -rw-r--r-- | ldc_3d.py | 2 | ||||
| -rw-r--r-- | notebook/ldc_interactive.ipynb | 167 | ||||
| -rw-r--r-- | shell.nix | 5 | ||||
| -rw-r--r-- | simulation.py | 11 | 
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 @@ -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 @@ -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 +} @@ -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( | 
