real_time_fvm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Real Time 2D Finite Volume Code - Explicit Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from OpenGL.GLUT import *\n",
    "from OpenGL.GLU import *\n",
    "from OpenGL.GL import *\n",
    "from matplotlib import cm\n",
    "import numpy as np\n",
    "import numexpr as ne\n",
    "import networkx as nx\n",
    "\n",
    "\n",
    "name = 'Explicit 2D FVM'\n",
    "\n",
    "def get_cmap_from_matplotlib(cmap=cm.Reds):\n",
    "    ncol = cmap.N\n",
    "    cmap_rgba = []\n",
    "    for i in range(ncol):\n",
    "        b, g, r, _ = cmap(i) # Not sure why this is inverted, I was expecting r, g, b.\n",
    "        cmap_rgba.append(int(255.0) << 24 | \n",
    "            (int(float(r) * 255.0) << 16) | \n",
    "            (int(float(g) * 255.0) << 8 ) | \n",
    "            (int(float(b) * 255.0) << 0 ) )\n",
    "\n",
    "    return np.array(cmap_rgba), len(cmap_rgba)\n",
    "\n",
    "\n",
    "###### Flow definition #########################################################\n",
    "Lx, Ly = 1.0, 1.0\n",
    "ni, nj = 200, 200\n",
    "\n",
    "dx, dy = Lx / ni, Ly / nj\n",
    "\n",
    "Ax = dy * 1.0\n",
    "Ay = dx * 1.0\n",
    "A = Ax\n",
    "\n",
    "k = 100.0\n",
    "rho = 1000.0\n",
    "cp = 4200.0\n",
    "dV = dx * dy * 1.0\n",
    "q = 1000.0 / dV\n",
    "\n",
    "# Condition for positive coefficients \n",
    "# dt < rho * c * dx^2 / 2k (for 2D)\n",
    "dt = 0.9 * rho * cp * (dx ** 2) / (4.0 * k)\n",
    "\n",
    "def to_continuous(ij):\n",
    "    i,j = ij\n",
    "    return i + j * ni\n",
    "    \n",
    "zero_I = np.zeros(nj)\n",
    "zero_J = np.zeros(ni)\n",
    "last_I = (ni-1) * np.ones(nj)\n",
    "last_J = (nj-1) * np.ones(ni)\n",
    "I = np.arange(0,ni)\n",
    "J = np.arange(0,nj)\n",
    "\n",
    "left  =  np.array([zero_I,J]).astype(int).T\n",
    "right =  np.array([last_I,J]).astype(int).T\n",
    "north =  np.array([I,last_J]).astype(int).T\n",
    "south =  np.array([I,zero_J]).astype(int).T\n",
    "\n",
    "graph = nx.grid_2d_graph(ni,nj)\n",
    "\n",
    "subgraph_left  = graph.subgraph(map(tuple,  left))\n",
    "subgraph_right = graph.subgraph(map(tuple, right))\n",
    "subgraph_north = graph.subgraph(map(tuple, north))\n",
    "subgraph_south = graph.subgraph(map(tuple, south))\n",
    "\n",
    "T_left  = 10.0\n",
    "T_right = 10.0\n",
    "T_north = 10.0\n",
    "T_south = 10.0\n",
    "\n",
    "initial_temperature = 23.0\n",
    "T     = initial_temperature * np.ones((ni, nj))\n",
    "T_old = np.copy(T)\n",
    "\n",
    "edges = np.array([e for e in graph.edges_iter()])\n",
    "\n",
    "e0 = edges[:,0]\n",
    "e1 = edges[:,1]\n",
    "e0 = [e0[:,0],e0[:,1]]\n",
    "e1 = [e1[:,0],e1[:,1]]\n",
    "\n",
    "nodes = np.array([n for n in graph.nodes_iter()])\n",
    "nodes = [nodes[:,0], nodes[:,1]]\n",
    "\n",
    "nodes_left  = np.array([n for n in subgraph_left.nodes_iter() ])\n",
    "nodes_left = [nodes_left[:,0], nodes_left[:,1]]\n",
    "\n",
    "nodes_right = np.array([n for n in subgraph_right.nodes_iter()])\n",
    "nodes_right = [nodes_right[:,0], nodes_right[:,1]]\n",
    "\n",
    "nodes_north = np.array([n for n in subgraph_north.nodes_iter()])\n",
    "nodes_north = [nodes_north[:,0], nodes_north[:,1]]\n",
    "\n",
    "nodes_south = np.array([n for n in subgraph_south.nodes_iter()])\n",
    "nodes_south = [nodes_south[:,0], nodes_south[:,1]]\n",
    "\n",
    "###### OpenGL Properties #######################################################\n",
    "solid = np.ones((nj, ni))\n",
    "plot_rgba = np.zeros(ni*nj, dtype=np.uint32)\n",
    "cmap_rgba, ncol = get_cmap_from_matplotlib(cmap=cm.jet)\n",
    "\n",
    "iter_per_frame = 1\n",
    "\n",
    "import requests\n",
    "\n",
    "\n",
    "def solve_fvm(sources=[]):\n",
    "    try:\n",
    "        Tpresc = float(requests.get('http://10.0.0.196:5000/').text.split()[0])\n",
    "    except:\n",
    "        Tpresc = 30.0\n",
    "    \n",
    "    for i in range(iter_per_frame):\n",
    "\n",
    "        T_old[:] = T[:]  \n",
    "\n",
    "        T[nodes] = rho * cp * dV * T_old[nodes] / dt\n",
    "\n",
    "        #T_old_n = T_old[nodes]\n",
    "        #T[nodes] = ne.evaluate(\"rho * cp * dV * T_old_n / dt\")\n",
    "\n",
    "        # apply boundary conditions\n",
    "        T[nodes_left ] += -2.0 * k * A / dx * (T_old[nodes_left ] - T_left ) *0\n",
    "        T[nodes_right] += -2.0 * k * A / dx * (T_old[nodes_right] - T_right) *0\n",
    "        T[nodes_north] += -2.0 * k * A / dx * (T_old[nodes_north] - T_north) *0\n",
    "        T[nodes_south] += -2.0 * k * A / dx * (T_old[nodes_south] - T_south) *0       \n",
    "\n",
    "        #T_old_e0 = T_old[e0]\n",
    "        #T_old_e1 = T_old[e1]\n",
    "        #flux = ne.evaluate(\"k * A / dx * (T_old_e0 - T_old_e1)\")\n",
    "\n",
    "        flux = k * A / dx * (T_old[e0] - T_old[e1])  \n",
    "        np.add.at(T, e0, -flux)\n",
    "        np.add.at(T, e1,  flux)\n",
    "        T[sources] += q * dV \n",
    "        T[nodes] *= dt / (rho * cp * dV)\n",
    "        \n",
    "        \n",
    "        T[100,100] = Tpresc\n",
    "        \n",
    "    return (T.T).flatten()        \n",
    "\n",
    "    \n",
    "def display():\n",
    "    '''\n",
    "    Set upper and lower limits for plotting.\n",
    "    Do one FVM step.\n",
    "    \n",
    "    '''   \n",
    "    sources = (1 - solid.T).astype(np.bool)\n",
    "\n",
    "    plotvar = solve_fvm(sources=sources)\n",
    "  \n",
    "    # convert the plotvar array into an array of colors to plot\n",
    "    # if the mesh point is solid, make it black\n",
    "    minvar = 20.0\n",
    "    maxvar = 31.0 #1.1 * np.max(plotvar)\n",
    "\n",
    "    frac = (plotvar[:] - minvar)/(maxvar - minvar)\n",
    "\n",
    "    icol = frac * ncol    \n",
    "    plot_rgba[:] = solid.flatten().astype(np.int) * cmap_rgba[icol.astype(np.int)] \n",
    "\n",
    "    # Fill the pixel buffer with the plot_rgba array\n",
    "    glBufferData(GL_PIXEL_UNPACK_BUFFER, plot_rgba.nbytes, plot_rgba, GL_STREAM_COPY)\n",
    "\n",
    "    # Copy the pixel buffer to the texture, ready to display\n",
    "    glTexSubImage2D(GL_TEXTURE_2D,0,0,0,ni,nj,GL_RGBA,GL_UNSIGNED_BYTE, None)\n",
    "\n",
    "    # Render one quad to the screen and colour it using our texture\n",
    "    # i.e. plot our plotvar data to the screen\n",
    "    glClear(GL_COLOR_BUFFER_BIT)\n",
    "\n",
    "    glBegin(GL_QUADS)\n",
    "\n",
    "    x0, y0 = 0.0, 0.0\n",
    "    x1, y1 = ni, nj\n",
    "    \n",
    "    glTexCoord2f(0.0, 0.0)\n",
    "    glVertex3f(x0, y0, 0.0)\n",
    "    \n",
    "    glTexCoord2f(1.0, 0.0)\n",
    "    glVertex3f(x1, y0, 0.0)\n",
    "    \n",
    "    glTexCoord2f (1.0, 1.0)\n",
    "    glVertex3f(x1, y1, 0.0)\n",
    "    \n",
    "    glTexCoord2f (0.0, 1.0)\n",
    "    glVertex3f(x0, y1, 0.0)\n",
    "\n",
    "    glEnd()\n",
    "    glutSwapBuffers()  \n",
    "\n",
    "\n",
    "\n",
    "def resize(w,h):\n",
    "    '''\n",
    "    GLUT resize callback to allow us to change the window size.\n",
    "    \n",
    "    '''\n",
    "    global width, height\n",
    "    width = w\n",
    "    height = h\n",
    "    glViewport (0, 0, w, h)\n",
    "    glMatrixMode (GL_PROJECTION)\n",
    "    glLoadIdentity()\n",
    "    glOrtho(0., ni, 0., nj, -200. ,200.)\n",
    "    glMatrixMode(GL_MODELVIEW)\n",
    "    glLoadIdentity()\n",
    "\n",
    "    \n",
    "def mouse(button, state, x, y):   \n",
    "    global draw_solid_flag, ipos_old, jpos_old\n",
    "\n",
    "    if button == GLUT_LEFT_BUTTON and state == GLUT_DOWN:       \n",
    "        draw_solid_flag = 0\n",
    "        xx = x\n",
    "        yy = y\n",
    "        ipos_old = int(float(xx) / width * float(ni))\n",
    "        jpos_old = int(float(height - yy) / height * float(nj))\n",
    "        \n",
    "    if button == GLUT_RIGHT_BUTTON and state == GLUT_DOWN:        \n",
    "        draw_solid_flag = 1\n",
    "        xx = x\n",
    "        yy = y\n",
    "        ipos_old = int(float(xx) / width * float(ni))\n",
    "        jpos_old = int(float(height - yy) / height * float(nj))\n",
    "\n",
    "\n",
    "\n",
    "def mouse_motion(x,y):\n",
    "    '''\n",
    "    GLUT call back for when the mouse is moving\n",
    "    This sets the solid array to draw_solid_flag as set in the mouse callback\n",
    "    It will draw a staircase line if we move more than one pixel since the\n",
    "    last callback - that makes the coding a bit cumbersome:\n",
    "    '''\n",
    "    global ipos_old, jpos_old\n",
    "\n",
    "    xx = x\n",
    "    yy = y\n",
    "    ipos = int(float(xx) / width * float(ni))\n",
    "    jpos = int(float(height-yy) / height * float(nj))\n",
    "\n",
    "    if ipos <= ipos_old:\n",
    "        i1 = ipos\n",
    "        i2 = ipos_old\n",
    "        j1 = jpos\n",
    "        j2 = jpos_old\n",
    "    \n",
    "    else:\n",
    "        i1 = ipos_old\n",
    "        i2 = ipos\n",
    "        j1 = jpos_old\n",
    "        j2 = jpos\n",
    "    \n",
    "    jlast = int(j1)\n",
    "\n",
    "    for i in range(i1, i2+1):\n",
    "        if i1 != i2:\n",
    "            frac = (i-i1) / (i2-i1)\n",
    "            jnext = int((frac * (j2-j1)) + j1)\n",
    "        else:\n",
    "            jnext = int(j2)\n",
    "        \n",
    "        if jnext >= jlast:\n",
    "            solid[jlast,i] = draw_solid_flag\n",
    "            \n",
    "            for j in range(jlast, jnext+1):    \n",
    "                solid[j,i] = draw_solid_flag\n",
    "        else:\n",
    "            solid[jlast,i] = draw_solid_flag\n",
    "            for j in range(jnext, jlast+1):\n",
    "                solid[j,i] = draw_solid_flag\n",
    "        \n",
    "        jlast = jnext\n",
    "        \n",
    "    ipos_old=ipos\n",
    "    jpos_old=jpos\n",
    "\n",
    "\n",
    "glutInit(name.encode())\n",
    "glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB)\n",
    "glutInitWindowSize(ni, nj)\n",
    "glutInitWindowPosition(50, 50)\n",
    "glutCreateWindow(name.encode())\n",
    "\n",
    "glClearColor(1.0,1.0,1.0,1.0)\n",
    "glMatrixMode(GL_PROJECTION)\n",
    "glLoadIdentity()\n",
    "glOrtho(0,ni,0.,nj, -200.0, 200.0)\n",
    "\n",
    "glEnable(GL_TEXTURE_2D)\n",
    "\n",
    "gl_Tex = glGenTextures(1)\n",
    "glBindTexture(GL_TEXTURE_2D, gl_Tex)\n",
    "\n",
    "glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP)\n",
    "glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP)\n",
    "glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)\n",
    "glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)\n",
    "glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA8, ni, nj, 0, GL_RGBA, GL_UNSIGNED_BYTE, None)\n",
    "\n",
    "gl_PBO = glGenBuffers(1)\n",
    "glBindBuffer(GL_PIXEL_UNPACK_BUFFER, gl_PBO)\n",
    "\n",
    "#setup callbacks\n",
    "glutDisplayFunc(display)\n",
    "glutReshapeFunc(resize)\n",
    "glutIdleFunc(display)\n",
    "glutMouseFunc(mouse)\n",
    "glutMotionFunc(mouse_motion)\n",
    "\n",
    "glutMainLoop()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:petsc_debug]",
   "language": "python",
   "name": "conda-env-petsc_debug-py"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}