Skip to content
Snippets Groups Projects
Normal view History Permalink
demo_wave_equation.ipynb 52.8 KiB
Newer Older
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
1 2 3 4 5 6 7 8 9 10 
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
Jan Hönig's avatar
`create_kernel` API Update  
Jan Hönig committed
11 12 
    "import psutil\n",
    "\n",
Michael Kuron's avatar
skip more tests based on dependencies  
Michael Kuron committed
13 14 15 
    "from pystencils.session import *\n",
    "\n",
    "import shutil"
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demo: Finite differences - 2D wave equation\n",
    "\n",
    "In this tutorial we show how to use the finite difference module of *pystencils* to solve a 2D wave equations. The time derivative is discretized by a simple forward Euler method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\frac{\\partial^2 u}{\\partial t^2} = \\mbox{div} \\left( q(x,y) \\nabla u \\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin by creating three *numpy* arrays for the current, the previous and the next timestep. Actually we will see later that two fields are enough, but let's keep it simple. From these *numpy* arrays we create *pystencils* fields to formulate our update rule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "size = (60, 70) #  domain size\n",
    "u_arrays = [np.zeros(size), np.zeros(size), np.zeros(size)]\n",
    "\n",
    "u_fields = [ps.Field.create_from_numpy_array(\"u%s\" % (name,), arr)\n",
    "            for name, arr in zip([\"0\", \"1\", \"2\"], u_arrays)]\n",
    "\n",
    "# Nicer display for fields\n",
    "for i, field in enumerate(u_fields):\n",
    "    field.latex_name = \"u^{(%d)}\" % (i,)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*pystencils* contains already simple rules to discretize the a diffusion term. The time discretization is done manually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
72 
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAAmCAYAAAAoXfRlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAABJ0AAASdAHeZh94AAAK8ElEQVR4nO2df6wcVRXHP6U1BcFAwYD8MPyB/FTTfYJAocgDoYIUQwxEIRYRlRiCSRMwEdAcvkDwJ9CGQIxVqBSMoSAVIojyQxAVAfUVihgkRkILotBiQBR4+Pzj3n1v3r7d2Zl9O7M7u+eTNN29c++8M98557z7Zu6ZmTMxMYHjOI7jOM4wMi9PZ0nbAeuAY81sQtJS4HJgK+AbZvY9SQuAG8zshK5bWzGSegE/BkaBe8zs5ESfgderQYc9gDXAzsA4cImZrR0GHXpJk9i9lQZ/HJZz4FoEsugQ+7kWDIcOnTAI2m2Vs//ngJviwc4DrgCOBkaAL0naycy2AJskLeqyrVVkUi9gJXB6Y4ch0Supwziw3MwOAJYAKyRtOyQ69JLkOYAm/jhE58C1CLTVAVyLOkOiQydUXru8E6HTgJ/EzwcDT5jZJjN7FbiT8IsN4Dbg1O6YWGkm9TKzXwKvtOg36HoldXjezMbi578DLwI7xn6DrkMvScZumj8OwzlwLQJZdQDXos6g69AJldcu860xSfOBXczshdi0G7Ap0WUTsHv8/Afgom4YWDaSzgc+DuwLvA48BJxvZhty7qdRrzQqq1c70nSQdCAw18yejU0Dq0M3iL55GXC1mZ2TY5z7YsS1COTUAVyLOgOrQycMinZ5rgjtBGzJ2PefwK75zekLRoFrgMMIt/3Ggbsl7Zg2qAmV1UvSakkXdWl3TXWIel4PnJVo7isd+glJhxK0eqyD4ZX1RSjHH1swyFrk0QH6TAv3ic5x7WYyeUVI0jHAHcB2ZvZGbNsDeBbYH/gHsHVi7HNMXQEifn44ft4a+E9xZheHmX0k+V3SMuBfwOHA7bEtVSsz+zPwX6brlUZl9erAb+p/RawDvm5mv0lsqqwORSJpe+BG4EzAmmzPfQ5SqOw58Licoss+Aa5FncrqkJdhiqfkFaER4E/1A060vQY8ZWabgW3iImkIk573Sdo9rho/HrgrbnsP8GSxppfGOwg6JWe9qVoBNNErjSrrlctvJM0BVgP3mtmahn1VWYci+S5ws5nd12J73thNo8rnwONyim76BLgWdaqsQ16GJp6SxteAsYbtI8AGM/tf/H4/cAjwazMbl3QucB9hovBNM3sp9juSsHgayFZ238rAVn1LLMdbSdDlt4m2Gu21goReku4GFgLbStoInGJm9X0WqlfBWtXI4TeEK2ufAB6TdFLcvszMHqdBB5hRen8CGXWIYyv/eAdJnyckkE+ldKuR7xyQ4o9VPgc1csZltDGTFlnLzZvRgzxWo3s+Ad3NUWWXV9fo0/jIS8na1SgxntoZU2S8JSdCI8Cqhv3VmC7ENcAZxAM2s9sIK8EbORFIGtqs7P4owi2n30u6NTGJSh5My75mtkXSJkmLEsGaHHspcGET25IcFVe4N0XSFcBiYLGZvZXYlEUrSOhlZsek2FG4Xm20ugC4INE0H5iQdF6i7Xgz+1UT23P5jZk9SOu1aY06QNQCmNvq2JrtqFMtiiaPX0ral7A4erGZvZnSv5PYbeWPPT8Hs/DH3HEJubRoViZ8LfDpFuOBWeexIrXIqgN0KUdFZuhWoA7Qp/HRIWVqV3Y8taOweNsqDtwG2IeZB/hBYH39i5k9BDwcb2+0MmIBcJWF5wbUyVp230i7vmnleCsI93/T/j3cYiySroz7PtrM/ppoz6QV9J1eaVp9h+Dg9X+3NWl7tIntRfsNTGmRRwfa9O9lGecKsvvlIuCdwBOSxiWNE/6qOjt+nz+g5yC3P/YgLtuVmyeZTR4rTIssOsT9dTNHdVpePYg5KjdladeLeGpHkfFWvyK0F2E2+5eE4R8iPAV4rMGYa9sYu4Vwuau+nzxl942069uyHM/MXiQ8oyY3klYSbt8cZWFBWJLMWkU7+kWvNK02A5sTNrwCbDazp9Nsp0C/ifua1ELSEWTXAfr08Q45/XIdM5P7dQS9LwPeAN7LgJ2DDv2x7LjMw2zyWKFatNMh9ulmjkpj2HJUNxmmeGpH7nirT4ReBCaAg4CNkg4iXBKboLNy3SR5yzTz0PVyPElXA8uAk4Atkt4VN70aZ5dFagXF6VVE6aJrUSBm9jLwcrJN0r8JSWxD/O7nIFBVHcC1qDOUOUrSGM2f6bfEzJ7LuJvK+1CXdMjCDK3mQXjCb7yP+H1JVwEPAD8CTou//GdDY3ldWtl9I+36FlGOd3b8/56GdgEXFawVFKdX17UqWYs8OrTr37dlnHnxcxDoQVzmodQ81uc5Ko1hy1F1G2tdsKPqPtQtHaCDeJtTxtvnJT0D7GWh0mweoYRulLiQCTjMwgLKe4DTzWxTHNeyb9w+AlxsZicWfhAlUoReVdWqrkX8muYLA69Fr/BzEEjGZaJtFDjHpr+kdODzWKc5KjF+lJm6VU4H6Dw+ZvHzRhkA7ZrFU4YxoxQQb3nfNdYp9fI64kHXy+7HgMtjwMwhlAlP3sts1Tex3xmljANCEXpVVav7gUPSjm2ItOgVfg4Ck3EJk2XCa4GPStooadEQ5bGOchQ01y1uqqIO0GF8dMKAaTctntpRZLxlfZjWbMlSdr8fcIuZTbtk1aJvnWaljINAEXpVVatk6X2rYxsWLXqFn4NA2zJhSfszHHlsNjkqT2l6Feg4PvIyYNpN86F2FBlvpdwaA5B0JnCdZX9mQLv9LQCONLN13dhfv9FNvaquVS+1kLQc2CHHjxirqs5puD8GPI9N4T4xRbf9IufPrqx2ZevWSqvSJkKOU0Uk/Q3YM8eQH5jZGcVY4ziO43Qbnwg5juM4jjO0lLVY2nEcx3Ecp+8oa7G041QSXyPkOI4z2MwDkOT3x5yhxsxavSdnOTnXCNHwGP46HmeO4zj9h68R6kMkvRtYA+wMjAOXmNna3lrlOI7TGs9bTlXxNUL9yTiw3MwOILw1d4WkbXtsk+M4Thqet5xK4leEeoCkLwDnmtneGfuvB5aa2bPFWuY4jjNF3lzVMNbzllMJ/IpQb6gRHv2NpG9JuqtVR0kHAnM9mTiO0wNqxFyVB89bTpXwiVBvqAF/jJ8PpsUbiiXtCFwPnFWOWY7jONOoMZWrMuF5y6kafmusYCTVgCuBQ4Gngc8SXgZ3KnAz8LZE9yfj/XUkzQd+AawyszVl2uw4zvCRkqtOMbM7JJ0M3AjsY2bPxDErgaWEt3u/4HnLqSJ+RahAJO1NeMPu74D3A18GbgLeDjwK1N8cfAiwK3B4HDcHWA3c68nEcZyiaZOrxmK3W4DHga/EMecR/qA7Lk6CPG85lcSvCBVIXPvzkpmdlmhbBXzMzHaRtBT4IbB98qVzkhYDDwCPJXa3zMweL8l0x3GGiHa5KtG2BPgp8FXgAuDDZvZI3OZ5y6kk/mTpgojP1FhCWAOU5E2m/sIaAdY3vnnXzB7Er9Y5jlMCGXMVAGb2c0mPAJcCJ9YnQXGb5y2nkrjTFscI8BawvqH9A0wllxo5FyI6juN0mSy5CgBJRwMLgTnAC2UY5zhF4xOh4pgA5gLz6w2SjiCsB6pPfhYy/TKy4zhO2WTJVUhaCNwKfJHwGpmvlWql4xSET4SK41HgdeDbkvaSdAJwQ9w2Fv+fB+wnaTdJO5RvouM4TvtcJWlP4E7gcjO7FjDgWEmjpVvrOF3GJ0IFYWbPA58BjgM2ABcSKipeA56K3S4EPglsxP+6chynB7TLVfG5QD8Dbjezi+OYDcBaPG85A4BXjTmO4ziOM7T4FSHHcRzHcYaW/wMNNksUDH5eSgAAAABJRU5ErkJggg==\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
73 
      "text/latex": [
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
74 
       "$\\displaystyle \\frac{{u^{(0)}}_{(0,0)} - 2 {u^{(1)}}_{(0,0)} + {u^{(2)}}_{(0,0)}}{dt^{2}} = \\frac{- 4 {u^{(1)}}_{(0,0)} + {u^{(1)}}_{(1,0)} + {u^{(1)}}_{(0,1)} + {u^{(1)}}_{(0,-1)} + {u^{(1)}}_{(-1,0)}}{dx^{2}}$"
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
75 76 
      ],
      "text/plain": [
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
77 78 79 80 
       "u_0_C - 2⋅u_1_C + u_2_C   -4⋅u_1_C + u_1_E + u_1_N + u_1_S + u_1_W\n",
       "─────────────────────── = ────────────────────────────────────────\n",
       "            2                                 2                   \n",
       "          dt                                dx                    "
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "discretize = ps.fd.Discretization2ndOrder()\n",
    "\n",
    "def central2nd_time_derivative(fields):\n",
    "    f_next, f_current, f_last = fields\n",
    "    return (f_next[0, 0] - 2 * f_current[0, 0] + f_last[0, 0]) / discretize.dt**2\n",
    "\n",
    "rhs = ps.fd.diffusion(u_fields[1], 1)\n",
    "\n",
    "wave_eq = sp.Eq(central2nd_time_derivative(u_fields), discretize(rhs))\n",
    "\n",
    "wave_eq = sp.simplify(wave_eq)\n",
    "wave_eq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The explicit Euler scheme is now obtained by solving above equation with respect to $u_C^{next}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
117 
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAAAmCAYAAAD+8btsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAABJ0AAASdAHeZh94AAAVq0lEQVR4nO2debRcVZWHv0AgBFAkkSGICg1ExCgvogSJkUCUAIIzKrZokBYRcdYWAq4fO9DQogKxmUQEFHDZRE0Ew6AGCK0MQSCEIEpARAiDzAookJj+Y59K7qt3q+reml69x/7WeqteneHeffc991Sd2uf8zohVq1YRBEEQBEEQBEEQwMjBNiAIhhNmtiEwD3iHpFVmNheYCiyQ9IFUZmPgAknvHDRDS2JmrwTOBzYFVgDHSpozCHYMS/9C+LjThH97j7gng0uv+L8RvWRnTlvZF/g2sBbwDUlnD8e20gxD/blaa7ANCIJhxn8AF0mqhGZnAx/LFpD0BLDczN7SbeNaYAXwBUk7AHsCp5jZBoNgx3D1L4SPO034t/eIezK49Ir/G1HYTjN7a4dtWd1WzGwkcBKwBzAR+KqZjR2mbaUZevq5atRWYoAUBO3lI8DPK28kXQ38PafcxcABXbKpIWZ2qJktq5Uv6UFJi9P/DwGPAmO6ZF6WYelfCB+3SrTh3iLa/ODTy89EkfZRoaidZrY38EI77cwh21Z2Bm6XtFzS08Bl+AAOhlhb6RC9/lxtbGY71cqMKXZBkMHMjgSOB06TdHjJuqOAzSQ9XKD4zcAx5S3sGH3A4sobM/sm8AZJ06sLpg5lbUn3dc06Xjz+Tfnh4/L0EW24l+gj2vxg00fvPhP9bCtKLTtTNGc/SYe1xbr8c1e3lS2A5Zkiy4FXpP+HWltZTfoe9D7gNcBzwPXAkZKWljhGzz9Xki4xs7OAQ/LyI4IUBAkz2wV/UJY0eYixwBMFyz4CjGvyPJ2gD7gl835nYFF1ITMbA/yQGh1KTvnzzOyYNtgHLwL/Qvi4BfqINtxL9BFtfrDpowPPRKrT6j2otq3IOevZ+RFqtK82MmTbSsn7NRU4HdgVnz64Avh18n9Rhoqv7jWz3fMyIoIUBICZbQRcCHwCUE7+24FLgQ0lPZ/StgTuA14r6Q/AP4H1Cp5yPeAfbTC9NGbWB5wM7ALcBRwMvB6YZWbrAk8D6wBvM7OjgTsk7ZB+EZoH/LekaztgV10fA39lGPs31e2Yj6MNRxvuFNHmB59efSYa2ZbyP4B//o6XdG9Kmw3sC+wq6eECdr4f+HILNjbTVh5gTcSI9H9lkNazbaUR1RFGMzsQeAqYDFyS0oZLX/cb4IPAVdUZLUeQzGxjM3vYzLYpUWeOmTXdkIOgA5wF/ETSgIckMRH4faUjyKQ9C9wJIOlxYHQK9TdiW+COFuxtCjPbDlgI3IB/QB0BXASsj093WAFUFktOwn/VmWxmI4DzgCslnd8h8+r6eDj7N9XttI+jDUcbbjvR5gefXn4mCtgG8FPgNuDoVOcr+JqUvdLgqK6d6d5NAu5uwdRm2soiYIKZvcJcsW1v4IqU15NtpUlego8XshGh4dLXLcbv2wDaEUGaCVwqaXXDLDB/cRaw0MzOlvRUG2wIgqYxs0/iD+hH6xTrY+B86YnAUkn/yqQtxDvq36Zj/xrYEdjAzO4H9pd0HbAbvqCzYkND6dA69peRGT0VmC/piPT+LjN7D/AuSQ+k443DF1LeWFGfMVd7+RCwJJUHOFDSbbXsaoI+Gvu4Kf+msg0lR2tRwsdN+Tcxmc76uI/ebcO596KdPh6mbbiwj8uSd0+izddv842MKdPnpPId7Xe69Ew0tC215ZnAfDO7G/9eOU1SRcShUTvZHHi2yD2oQx8l24qkFenH/qvwe3SipMdSucH6jO8Es3HfXJdJ66NLfR3wTkr0c2X8JekpM9vMzEZKWpHNa2mAZGbr4zJ++1VlTcXnL94IjMAHRL82sx0kPS7pNjP7E/6F9LRWbAiCPMzsOOCoBsV2Bx7ERRneKqme+s1E4HtVaX0M7CBOB2aQOgRJb69xvP2A7AdknnTo7nhY+yYzm5vpeFdTr6yZLTezt6QOCPO9JPbE56FneaHqOiYCt2Y/bCT9hgIR5/QhNzOTNApYlX4RrLC3pP/LqV7Ex836F/IlR88BPl7jGEBxH7fi33QtnfZxT7bhRO69kPREu3w83NpwEz4uy4B70s77kY43bNp8/StYTaE+B7rT7xT1f7Kn9D0oYRuSfmlmNwLH4WILN5awc1Pgbzk2F/oeIFdXa7atXIyrsFXT9c/4qnqt9GPZ45wEvBX/jrQyk9WVvg5YmxL9XJP+ehrYBP8+uJpWp9jtA6wiXXwFSdMlnStpaRrhH5hOPjlTLCQQg05yCj4Ptt7fInzawcuB281shZmtwH/NOCy9H2Vmo4HxDOwo3wzcmk2QdD2wyHxKQC7pl4z/kev/VygqHVpNGZnRicDKapuBN1ZdWx8lF89mODPVr/xdnJP2u+pKRX3cgn+huORoNUV93A3/QhM+7vE23OhedNvHQ6UNl/JxWerck2jza+jXpzSiRJ8Dw6DfKWEbZrYHHl0YARRRPsuyLj6NsJpTKPA9oMf7x2alxJvqx6qu4+R0/D0k/SmT3s2+rmw/14y/VuBtqB+tTrGbAtxU4NeTvPmLi4CjzWy0pCG5kK0WZrYVcA/wA0kz2nC8zwGHAlvji9m+KOmUVo87nJH0KL5XQl3MbB4DO4lzgWV4ZOl54HX4rxjLMvXeBmxJjkyppHMa2PYEHjquHKuMdGg1ZWRGV6XrGJWuCzObgofAT8mU25GqcHdR5POOH6+8N7O/A49LuqtB1W0o6OOy/k3HKiM5Wk1RH3fcv9C0jwv7N52jm224EV318RBqw+30cRmizdNyn1KE4dDvFLLNzHYE5gKfxadTnQDkysHX4BHgpTk2F/0eMIHe7R+bkhJvoR+rlJ+NT2vcXS5QkaVrfV1qL2X6uWb89RJcVKIfrQ6QXo2reDQib/7iA7iiyha0trCu65jZHcBzkvq6cK4P4/67Be9QngOub/cg7MWKpCeBJ7NpZvYM3pEsTe8fxTv6NwH3m9mb8NDyKpqXBM9SRg6zDNXSmb/D28+3zOxEYHs8BA79O7WRwPZmtgU+r/vJDthWzXDwcfi3M4SPnbb42MwWk//Zv6fSmpAGxP1wBtyPNvi2KEPlHjS0zcxejQ/cvi3pHDNbhK81mpoibkV4kJwBUgmGxLObQ0fksc3sNHzm13uAJ8xs85T1dIrKDBt/mSs8PpcXqGl1it1oXPawJpn5i++vmr9YMWZ0izYMBnOBHc1s6y6ca9/Kq6QjJR2TwpZBl5Dv3D0T+L6Z3Qd8EfgxcHfqLFqliHRorQ/XemX7SWdKehA4CNgLWIrPzT6PjEpP4ijgw8D9+C95HWcQfFyGQj4O/zbdhhsRPnba4mNJfZIm5PwVvT9xP5wBfUobfJtlOPQ7dW0z31fncuASSbNSnaXAnDI2SnoWuNvMNm3Szl5+dgt/xreRw/CoygJ88Fn5+wp03V9l+7my/nodcGXegUasWtW86IeZXQisI+mDNfJPxh/IASE6M5uEq9tt3sEQdUcwszfjUwS/LOmknPytaFN0x8yuxP03oiq9becIBh8zuxfYRq6KMxKXvJxKWmSI7wXxmJktAD4maXmqV6/sRGCWpGoRlRclWR9n0qYCh6u/clr4uAmabcOZ+lOpuhcpPXycaNXHTZxvKgOfj7gfibw+pUCdqeS38+h3WsBcPXmJpPmDbUse8Rlfjoq/0tvc60/lWvKXmR2CR1QvqLahXyg4NbDj8WjPz6ryXgXcC8yT9N6UfAuuUJF3cfXmLwJMAJYPtcFR4nf4rzDvw9UySpEGh1/FI2tj8AWJlwKmNZKcx5DZsNTMsiNZy+R93MyyijgHSTov1ZmBq4NMxMOKL+B7DZxR3RiyAy68DRyLq4C8HN9J+c/18rOhcDP7IHA4Pud5XXxjuB8BJ0l6LpXZEJ8fe6OkyZm6o/HQ6ii80Z+fyfs0Hp4/uNG81iHIQhpIh5ovdNyWzLziWmVT9gDpzBc5q30M+ZKj+I824ePmaKoNQ135VwgfZ2nax2Wpc0/ifqyhX5/SiFo+jb69LXwX/27UkwMk4jO+LAuBSZJ+W+v62+SvXfA1/gOoniu7U3rNU7Z4c3q9OZN2BfANMxubMaDI/EVwgYcrGIJozR4qnzGzUos0zewT+Kakz+GKGvcB25Hk0s1sF0l/Aa5OVWbga70sc5irgZcBn8cVQ+Zl8hZn/j8DuB24Bg+PjsWVB883s9dI+nqOidvgG7rdie9sPZr+8pl1883seOBIfI7qj3D5xL3xQdV0M9tT0vOSnk5zjSeZ2UskVZR9JuODI4BpQHZTuGnpdUGO3UOdItKh2wM/VdVc2RplIV8688VMQ8lRM3st4eNmaaUN15J/hfBxlqZ9XJY69yTuxxr63Y9G1PFp9O0tIulxM/ujmY1LU/t6jfiML8dqf9W5/pb8ZWZvAH6l/pvdrqbfFDszuwfYUNIm1QXN7AR8B+T9JP0ik34dvvnSaZm0WvP2TNIxZrYeHjWZriG6nsbMdsfnLX5K0llVeVuRM/3NzMbj83D/AuyWnf5gZtOAXwIXZyJ0mNnVqWzpKXZmto0yG/imtHXxEfTbgK0yYcnK8QBOkDSzql7d/FTmLcC1+KBvZ/k81UrIcy6+nuooScen9FnA1/H1VfNT2gl4dG0hMF7SK1P6WvgCuyclbcMwJA2ez1Vrm91VjrUx3m7mtWzYMCJ83Fna6d90vPBxFe32cclzx/2ootv3I+5BbdL3hP0l/e9g25JHfP6Uo9P+MrMPS/pxrTojMwXHAFtRO6rzpvR6c1W6AbPN7EwlEYbqL/M5HAzcUGtwlC7kD/i8wa4o3JnZHOB6Sd8uWOUa4DHgvXhEqAifxpX7Pq+queGSFpjZxXgUKRtRaZo830l6PkX49sAjMj+sKvIw/aNV1dTL/0R6Pa4yOErnrIQ898EjZcenrAX4AGkaa8Li0/A5oz8DTjWz8ZLuxDX7xwA/rWPbkKad0waVI50ZhI87TbunvoaPB9JuH5c8d9yPKrp9P+Ie1EbSv8zsosG2oxbx+VOOLvir7kA6O8Xujem11sZROwF/VZUai6TL0xfuLfE1SkV4Ade7r8VM4NLsF3wzOwyPLIzDp419QQ12AM5SoP4sYKGZnS3pqUbHk7TSzC4BPmJmGxWpg29KCrCbudBDNZvi2vLj8UFCS6R1Y1/DBx2vYqBiYJ6W/K1K64RqUC+/0oYGKIJIujPNud4646/rcEWRacnejdIxTswcYxo+nW+PWscOgiAIgiAYjMhqMDRp1FayMt+V9UcDvpib2b8BGzMwelQ5yXckFR0cIeksSX/MyzOz9fEow/czaR/C9wI6HhccuBa4LA0AGlKkvqTbgD8BHy16HXiUY118Y7MijE2vX8VFFqr/dk35G5awIZd0z27GF589BJwNHIdHf36Qio3KqfpQTlrR/I3Sa635v5X0l4FHs4DfAK83s01w1ZG1gQWS7kjlK+uOpuEa+zFACoIgCIIgCDpGNoI0Mb3mDYL2Sa+3dNac1edaRf9Fj18CzpP0vfT+s2a2Fz5l7cgCxyxa/2LgAOA0ivEr4Bl8mt2PCpSvRJk2kvS3uiVb50v4gGy1ql0FMzsA+HheJdz39aiXX7m+zcnf/HdcVTnwAc878AHQrrj+/W8zeXub76o8Bbhd0oDdjoMgCIIgCIKgXWQHSNsDL1RHgtKX00+lt7kRpDYzBbipEvpKogI7Ad+qKvdL1kRcalKy/iLgaDMbXUQBSNI/zewyYC8zW09S3U1zcQnhnfBrbFWKsrLp7to18rdNr3lrdnZr8dy1uAWfIjeVqgGSmW2LT8O8R/138K4o0k3DpyBem/HjAuDf8YHsBhRUr6sjEhIEQRAEQRAENZE0IjtAeh5Yx8y2k7QMwMw2wAUIJqQy3YggvZr+u+S+HB8EVEtpPwzUk4Ztpv4DuIjCFuRHQPKYi8sG7km+rGCWU4FDgJPNbFkSH1hNGsxNKri26gk8mlNrmuGf0+tU4JLMOabjUxg7wTm4AMfRZnaxpEfSOdfGB6hrkZk6mbgZjyi9G9iE/pG4ynS6I6ve16WASEgQBEGQMLNX4lsqbAqsAI6VNGdwrQqCIBg8sgOkK/C9jq5Je/xsiP+qvwRfC7I+vkan04xm4GCmW1SiRtViBvX4BT64fC8NBkiS/pBkC88Bbjezy3EBgnXwgc4UXMp6+0Ynle8jdAMwxcwuTMdZicuEL8E15A8C5pjZT/DB3wRgL+AifBPftiLpWjM7EfhPYGk67zP4PkgT8PVG36yqszJJmb87JS3I5N1rZnfjey+txKW/gyAIgvayAhcuWpz2LbzJzC6V9MxgGxYEQTAYZEUa/gv4Dh6VmAHsgCu7HYCvKVncJXWQR3FBiOz7lcBmVeU2o7GgQNn6Y9LrI4UsBdJaoitxee5a092y5S/Ap9ldCLwBOBwXhtgW+AlwWNFz45vxzscHPQKOJSnJpUHS7rggxTvxaWovBd4HnFniHKWQ9DW8zSwDPgZ8Dm9nRwPvUP6GXJVB0d8YqKJYybupoFJgEARBAJjZoWa2rFE5SQ9KWpz+fwj/3BxTt1IQBMEwpt9Gsb2AmX0FmCFpQibtBlxe+pBM2p34DroNRRqK1jezg/HNbLcsafMhwHeBPSRdVaZuEARBEHQCMzsTGCtp/xJ1dsI3IJ/QsHAQBMEwZWTjIl3nCuAbZjZW0mMp7STgfDNbhCucHYqvEyoaCSlafwq1N8qtx8+BM/DoTAyQgiAIgl6gj8ZrY1eTNoz/IfDJThkUBEEwFOi5CBKAmV0HXCDptEzaYfjalnHAUuCLkq7J5M8AzgW2lvTnnGM2qr8evvZpuqTrO3BZQRAEQdARzKwPOBnYBbgLF8y5CtgfX0N8ITC+olRrZrOBfYFdJT2cFGt/BXxP0vndv4IgCILeoRcjSOCbmc42szMlrQSQdDouPFCLrYHfA/fnZRaofzBwQwyOgiAIgqGEmW2Hi9icgUd/XoOL8awPLMaFlo7A14J+Mk1lPwCYnAZHI4DzgCtjcBQEQdCjAyRJl5vZafi+Ofc2Kp/YB/iMpBVNnvYF4LNN1g2CIAiCweJUYL6kI9L7u8zsPcC7JD0AYGYzgflJGXQmMK2ypQcwGVc2XZLqARwo6bZuXUAQBEEv0ZNT7IIgCIIgaEzaw+gvwM6Sbsyknw5sI2l6Ju1aYGdgP0mXdd3YIAiCIcJajYsEQRAEQdCjTMS3sri1Kv2N+PQ6AMxsD2BHYASDt9dgEATBkCAGSEEQBEEwdFkFrA2MqiSY2RRgEnBLer8jMBefRj4POKHrVgZBEAwhYopdEARBEAxRzGwccA/wA+BEYHtckOhVwGuBfwDXAWdKmmVmE4Al+L59Vw+K0UEQBD1ORJCCIAiCYIgi6UHgIGAvfAuLo3BFumeBvwKXA5dImpXKLwXmEFGkIAiCmkQEKQiCIAiCIAiCIBERpCAIgiAIgiAIgsT/A1EcQf0BL8loAAAAAElFTkSuQmCC\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
118 
      "text/latex": [
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
119 
       "$\\displaystyle {u^{(2)}}_{(0,0)} \\leftarrow \\frac{- 4 {u^{(1)}}_{(0,0)} dt^{2} + {u^{(1)}}_{(1,0)} dt^{2} + {u^{(1)}}_{(0,1)} dt^{2} + {u^{(1)}}_{(0,-1)} dt^{2} + {u^{(1)}}_{(-1,0)} dt^{2} + dx^{2} \\left(- {u^{(0)}}_{(0,0)} + 2 {u^{(1)}}_{(0,0)}\\right)}{dx^{2}}$"
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
120 121 
      ],
      "text/plain": [
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
122 123 124 125 126 
       "                     2           2           2           2           2     2                   \n",
       "         - 4⋅u_1_C⋅dt  + u_1_E⋅dt  + u_1_N⋅dt  + u_1_S⋅dt  + u_1_W⋅dt  + dx ⋅(-u_0_C + 2⋅u_1_C)\n",
       "u_2_C := ──────────────────────────────────────────────────────────────────────────────────────\n",
       "                                                    2                                          \n",
       "                                                  dx                                           "
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u_next_C = u_fields[-1][0, 0]\n",
    "update_rule = ps.Assignment(u_next_C, sp.solve(wave_eq, u_next_C)[0])\n",
    "update_rule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before creating the kernel, we substitute numeric values for $dx$ and $dt$. \n",
    "Then a kernel is created just like in the last tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
156 157 158 159 160 161 162 163 
       "<style>pre { line-height: 125%; }\n",
       "td.linenos .normal { color: inherit; background-color: transparent; padding-left: 5px; padding-right: 5px; }\n",
       "span.linenos { color: inherit; background-color: transparent; padding-left: 5px; padding-right: 5px; }\n",
       "td.linenos .special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; }\n",
       "span.linenos.special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; }\n",
       ".highlight .hll { background-color: #ffffcc }\n",
       ".highlight { background: #f8f8f8; }\n",
       ".highlight .c { color: #3D7B7B; font-style: italic } /* Comment */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
164 165 166 
       ".highlight .err { border: 1px solid #FF0000 } /* Error */\n",
       ".highlight .k { color: #008000; font-weight: bold } /* Keyword */\n",
       ".highlight .o { color: #666666 } /* Operator */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
167 168 169 170 171 172 
       ".highlight .ch { color: #3D7B7B; font-style: italic } /* Comment.Hashbang */\n",
       ".highlight .cm { color: #3D7B7B; font-style: italic } /* Comment.Multiline */\n",
       ".highlight .cp { color: #9C6500 } /* Comment.Preproc */\n",
       ".highlight .cpf { color: #3D7B7B; font-style: italic } /* Comment.PreprocFile */\n",
       ".highlight .c1 { color: #3D7B7B; font-style: italic } /* Comment.Single */\n",
       ".highlight .cs { color: #3D7B7B; font-style: italic } /* Comment.Special */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
173 174 
       ".highlight .gd { color: #A00000 } /* Generic.Deleted */\n",
       ".highlight .ge { font-style: italic } /* Generic.Emph */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
175 
       ".highlight .gr { color: #E40000 } /* Generic.Error */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
176 
       ".highlight .gh { color: #000080; font-weight: bold } /* Generic.Heading */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
177 178 
       ".highlight .gi { color: #008400 } /* Generic.Inserted */\n",
       ".highlight .go { color: #717171 } /* Generic.Output */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
179 180 181 182 183 184 185 186 187 188 189 190 
       ".highlight .gp { color: #000080; font-weight: bold } /* Generic.Prompt */\n",
       ".highlight .gs { font-weight: bold } /* Generic.Strong */\n",
       ".highlight .gu { color: #800080; font-weight: bold } /* Generic.Subheading */\n",
       ".highlight .gt { color: #0044DD } /* Generic.Traceback */\n",
       ".highlight .kc { color: #008000; font-weight: bold } /* Keyword.Constant */\n",
       ".highlight .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */\n",
       ".highlight .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */\n",
       ".highlight .kp { color: #008000 } /* Keyword.Pseudo */\n",
       ".highlight .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */\n",
       ".highlight .kt { color: #B00040 } /* Keyword.Type */\n",
       ".highlight .m { color: #666666 } /* Literal.Number */\n",
       ".highlight .s { color: #BA2121 } /* Literal.String */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
191 
       ".highlight .na { color: #687822 } /* Name.Attribute */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
192 193 194 195 
       ".highlight .nb { color: #008000 } /* Name.Builtin */\n",
       ".highlight .nc { color: #0000FF; font-weight: bold } /* Name.Class */\n",
       ".highlight .no { color: #880000 } /* Name.Constant */\n",
       ".highlight .nd { color: #AA22FF } /* Name.Decorator */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
196 197 
       ".highlight .ni { color: #717171; font-weight: bold } /* Name.Entity */\n",
       ".highlight .ne { color: #CB3F38; font-weight: bold } /* Name.Exception */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
198 
       ".highlight .nf { color: #0000FF } /* Name.Function */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
199 
       ".highlight .nl { color: #767600 } /* Name.Label */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 
       ".highlight .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */\n",
       ".highlight .nt { color: #008000; font-weight: bold } /* Name.Tag */\n",
       ".highlight .nv { color: #19177C } /* Name.Variable */\n",
       ".highlight .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */\n",
       ".highlight .w { color: #bbbbbb } /* Text.Whitespace */\n",
       ".highlight .mb { color: #666666 } /* Literal.Number.Bin */\n",
       ".highlight .mf { color: #666666 } /* Literal.Number.Float */\n",
       ".highlight .mh { color: #666666 } /* Literal.Number.Hex */\n",
       ".highlight .mi { color: #666666 } /* Literal.Number.Integer */\n",
       ".highlight .mo { color: #666666 } /* Literal.Number.Oct */\n",
       ".highlight .sa { color: #BA2121 } /* Literal.String.Affix */\n",
       ".highlight .sb { color: #BA2121 } /* Literal.String.Backtick */\n",
       ".highlight .sc { color: #BA2121 } /* Literal.String.Char */\n",
       ".highlight .dl { color: #BA2121 } /* Literal.String.Delimiter */\n",
       ".highlight .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */\n",
       ".highlight .s2 { color: #BA2121 } /* Literal.String.Double */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
216 
       ".highlight .se { color: #AA5D1F; font-weight: bold } /* Literal.String.Escape */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
217 
       ".highlight .sh { color: #BA2121 } /* Literal.String.Heredoc */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
218 
       ".highlight .si { color: #A45A77; font-weight: bold } /* Literal.String.Interpol */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
219 
       ".highlight .sx { color: #008000 } /* Literal.String.Other */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
220 
       ".highlight .sr { color: #A45A77 } /* Literal.String.Regex */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 
       ".highlight .s1 { color: #BA2121 } /* Literal.String.Single */\n",
       ".highlight .ss { color: #19177C } /* Literal.String.Symbol */\n",
       ".highlight .bp { color: #008000 } /* Name.Builtin.Pseudo */\n",
       ".highlight .fm { color: #0000FF } /* Name.Function.Magic */\n",
       ".highlight .vc { color: #19177C } /* Name.Variable.Class */\n",
       ".highlight .vg { color: #19177C } /* Name.Variable.Global */\n",
       ".highlight .vi { color: #19177C } /* Name.Variable.Instance */\n",
       ".highlight .vm { color: #19177C } /* Name.Variable.Magic */\n",
       ".highlight .il { color: #666666 } /* Literal.Number.Integer.Long */</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 
       "<div class=\"highlight\"><pre><span></span><span class=\"n\">FUNC_PREFIX</span><span class=\"w\"> </span><span class=\"kt\">void</span><span class=\"w\"> </span><span class=\"n\">kernel</span><span class=\"p\">(</span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\"> </span><span class=\"k\">const</span><span class=\"w\"> </span><span class=\"n\">_data_u0</span><span class=\"p\">,</span><span class=\"w\"> </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\"> </span><span class=\"k\">const</span><span class=\"w\"> </span><span class=\"n\">_data_u1</span><span class=\"p\">,</span><span class=\"w\"> </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\">  </span><span class=\"n\">_data_u2</span><span class=\"p\">)</span><span class=\"w\"></span>\n",
       "<span class=\"p\">{</span><span class=\"w\"></span>\n",
       "<span class=\"w\">   </span><span class=\"k\">for</span><span class=\"w\"> </span><span class=\"p\">(</span><span class=\"kt\">int64_t</span><span class=\"w\"> </span><span class=\"n\">ctr_0</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"mi\">1</span><span class=\"p\">;</span><span class=\"w\"> </span><span class=\"n\">ctr_0</span><span class=\"w\"> </span><span class=\"o\">&lt;</span><span class=\"w\"> </span><span class=\"mi\">59</span><span class=\"p\">;</span><span class=\"w\"> </span><span class=\"n\">ctr_0</span><span class=\"w\"> </span><span class=\"o\">+=</span><span class=\"w\"> </span><span class=\"mi\">1</span><span class=\"p\">)</span><span class=\"w\"></span>\n",
       "<span class=\"w\">   </span><span class=\"p\">{</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\">  </span><span class=\"n\">_data_u2_00</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"n\">_data_u2</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mi\">70</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span><span class=\"p\">;</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\"> </span><span class=\"n\">_data_u1_01</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"n\">_data_u1</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mi\">70</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mi\">70</span><span class=\"p\">;</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\"> </span><span class=\"n\">_data_u1_00</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"n\">_data_u1</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mi\">70</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span><span class=\"p\">;</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\"> </span><span class=\"n\">_data_u1_0m1</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"n\">_data_u1</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mi\">70</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span><span class=\"w\"> </span><span class=\"o\">-</span><span class=\"w\"> </span><span class=\"mi\">70</span><span class=\"p\">;</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\"> </span><span class=\"n\">_data_u0_00</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"n\">_data_u0</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mi\">70</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span><span class=\"p\">;</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"k\">for</span><span class=\"w\"> </span><span class=\"p\">(</span><span class=\"kt\">int64_t</span><span class=\"w\"> </span><span class=\"n\">ctr_1</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"mi\">1</span><span class=\"p\">;</span><span class=\"w\"> </span><span class=\"n\">ctr_1</span><span class=\"w\"> </span><span class=\"o\">&lt;</span><span class=\"w\"> </span><span class=\"mi\">69</span><span class=\"p\">;</span><span class=\"w\"> </span><span class=\"n\">ctr_1</span><span class=\"w\"> </span><span class=\"o\">+=</span><span class=\"w\"> </span><span class=\"mi\">1</span><span class=\"p\">)</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"p\">{</span><span class=\"w\"></span>\n",
       "<span class=\"w\">         </span><span class=\"n\">_data_u2_00</span><span class=\"p\">[</span><span class=\"n\">ctr_1</span><span class=\"p\">]</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"mf\">-1.0</span><span class=\"o\">*</span><span class=\"n\">_data_u0_00</span><span class=\"p\">[</span><span class=\"n\">ctr_1</span><span class=\"p\">]</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mf\">0.25</span><span class=\"o\">*</span><span class=\"n\">_data_u1_00</span><span class=\"p\">[</span><span class=\"n\">ctr_1</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mi\">1</span><span class=\"p\">]</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mf\">0.25</span><span class=\"o\">*</span><span class=\"n\">_data_u1_00</span><span class=\"p\">[</span><span class=\"n\">ctr_1</span><span class=\"w\"> </span><span class=\"o\">-</span><span class=\"w\"> </span><span class=\"mi\">1</span><span class=\"p\">]</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mf\">0.25</span><span class=\"o\">*</span><span class=\"n\">_data_u1_01</span><span class=\"p\">[</span><span class=\"n\">ctr_1</span><span class=\"p\">]</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mf\">0.25</span><span class=\"o\">*</span><span class=\"n\">_data_u1_0m1</span><span class=\"p\">[</span><span class=\"n\">ctr_1</span><span class=\"p\">]</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mf\">1.0</span><span class=\"o\">*</span><span class=\"n\">_data_u1_00</span><span class=\"p\">[</span><span class=\"n\">ctr_1</span><span class=\"p\">];</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"p\">}</span><span class=\"w\"></span>\n",
       "<span class=\"w\">   </span><span class=\"p\">}</span><span class=\"w\"></span>\n",
       "<span class=\"p\">}</span><span class=\"w\"></span>\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
256 257 258 
       "</pre></div>\n"
      ],
      "text/plain": [
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
259 
       "FUNC_PREFIX void kernel(double * RESTRICT const _data_u0, double * RESTRICT const _data_u1, double * RESTRICT  _data_u2)\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
260 
       "{\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
261 
       "   for (int64_t ctr_0 = 1; ctr_0 < 59; ctr_0 += 1)\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
262 
       "   {\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
263 264 265 266 267 268 
       "      double * RESTRICT  _data_u2_00 = _data_u2 + 70*ctr_0;\n",
       "      double * RESTRICT _data_u1_01 = _data_u1 + 70*ctr_0 + 70;\n",
       "      double * RESTRICT _data_u1_00 = _data_u1 + 70*ctr_0;\n",
       "      double * RESTRICT _data_u1_0m1 = _data_u1 + 70*ctr_0 - 70;\n",
       "      double * RESTRICT _data_u0_00 = _data_u0 + 70*ctr_0;\n",
       "      for (int64_t ctr_1 = 1; ctr_1 < 69; ctr_1 += 1)\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
269 
       "      {\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
270 
       "         _data_u2_00[ctr_1] = -1.0*_data_u0_00[ctr_1] + 0.25*_data_u1_00[ctr_1 + 1] + 0.25*_data_u1_00[ctr_1 - 1] + 0.25*_data_u1_01[ctr_1] + 0.25*_data_u1_0m1[ctr_1] + 1.0*_data_u1_00[ctr_1];\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
271 272 273 274 275 276 
       "      }\n",
       "   }\n",
       "}"
      ]
     },
     "metadata": {},
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
277 
     "output_type": "display_data"
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 
    }
   ],
   "source": [
    "update_rule = update_rule.subs({discretize.dx: 0.1, discretize.dt: 0.05})\n",
    "ast = ps.create_kernel(update_rule)\n",
    "kernel = ast.compile()\n",
    "\n",
    "ps.show_code(ast)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[_data_u0, _data_u1, _data_u2]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ast.get_parameters()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
316 
       "{u0: double[60,70], u1: double[60,70], u2: double[60,70]}"
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ast.fields_accessed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To run simulation a suitable initial condition and boundary treatment is required. We chose an initial condition which is zero at the borders of the domain. The outermost layer is not changed by the update kernel, so we have an implicit homogenous Dirichlet boundary condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "X,Y = np.meshgrid( np.linspace(0, 1, size[1]), np.linspace(0,1, size[0]))\n",
    "Z = np.sin(2*X*np.pi) * np.sin(2*Y*np.pi)\n",
    "\n",
    "# Initialize the previous and current values with the initial function\n",
    "np.copyto(u_arrays[0], Z)\n",
    "np.copyto(u_arrays[1], Z)\n",
    "# The values for the next timesteps do not matter, since they are overwritten\n",
    "u_arrays[2][:, :] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One timestep now consists of applying the kernel once, then shifting the arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run(timesteps=1):\n",
    "    for t in range(timesteps):\n",
    "        kernel(u0=u_arrays[0], u1=u_arrays[1], u2=u_arrays[2])\n",
    "        u_arrays[0], u_arrays[1], u_arrays[2] = u_arrays[1], u_arrays[2], u_arrays[0]\n",
    "    return u_arrays[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets create an animation of the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
384 385 386 387 388 
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No ffmpeg installed\n"
     ]
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
389 390 391 
    }
   ],
   "source": [
Michael Kuron's avatar
skip more tests based on dependencies  
Michael Kuron committed
392 393 394 395 396 
    "if shutil.which(\"ffmpeg\") is not None:\n",
    "    ani = plt.surface_plot_animation(run, zlim=(-1, 1))\n",
    "    ps.jupyter.display_as_html_video(ani)\n",
    "else:\n",
    "    print(\"No ffmpeg installed\")"
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.isfinite(np.max(u_arrays[2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Runing on GPU__\n",
    "\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
414 
    "We can also run the same kernel on the GPU, by using the ``cupy`` package."
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
415 416 417 418 
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
419 
   "execution_count": 12,
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
420 421 422 423 424 
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
425 426 427 428 429 430 431 432 
       "<style>pre { line-height: 125%; }\n",
       "td.linenos .normal { color: inherit; background-color: transparent; padding-left: 5px; padding-right: 5px; }\n",
       "span.linenos { color: inherit; background-color: transparent; padding-left: 5px; padding-right: 5px; }\n",
       "td.linenos .special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; }\n",
       "span.linenos.special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; }\n",
       ".highlight .hll { background-color: #ffffcc }\n",
       ".highlight { background: #f8f8f8; }\n",
       ".highlight .c { color: #3D7B7B; font-style: italic } /* Comment */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
433 434 435 
       ".highlight .err { border: 1px solid #FF0000 } /* Error */\n",
       ".highlight .k { color: #008000; font-weight: bold } /* Keyword */\n",
       ".highlight .o { color: #666666 } /* Operator */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
436 437 438 439 440 441 
       ".highlight .ch { color: #3D7B7B; font-style: italic } /* Comment.Hashbang */\n",
       ".highlight .cm { color: #3D7B7B; font-style: italic } /* Comment.Multiline */\n",
       ".highlight .cp { color: #9C6500 } /* Comment.Preproc */\n",
       ".highlight .cpf { color: #3D7B7B; font-style: italic } /* Comment.PreprocFile */\n",
       ".highlight .c1 { color: #3D7B7B; font-style: italic } /* Comment.Single */\n",
       ".highlight .cs { color: #3D7B7B; font-style: italic } /* Comment.Special */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
442 443 
       ".highlight .gd { color: #A00000 } /* Generic.Deleted */\n",
       ".highlight .ge { font-style: italic } /* Generic.Emph */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
444 
       ".highlight .gr { color: #E40000 } /* Generic.Error */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
445 
       ".highlight .gh { color: #000080; font-weight: bold } /* Generic.Heading */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
446 447 
       ".highlight .gi { color: #008400 } /* Generic.Inserted */\n",
       ".highlight .go { color: #717171 } /* Generic.Output */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
448 449 450 451 452 453 454 455 456 457 458 459 
       ".highlight .gp { color: #000080; font-weight: bold } /* Generic.Prompt */\n",
       ".highlight .gs { font-weight: bold } /* Generic.Strong */\n",
       ".highlight .gu { color: #800080; font-weight: bold } /* Generic.Subheading */\n",
       ".highlight .gt { color: #0044DD } /* Generic.Traceback */\n",
       ".highlight .kc { color: #008000; font-weight: bold } /* Keyword.Constant */\n",
       ".highlight .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */\n",
       ".highlight .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */\n",
       ".highlight .kp { color: #008000 } /* Keyword.Pseudo */\n",
       ".highlight .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */\n",
       ".highlight .kt { color: #B00040 } /* Keyword.Type */\n",
       ".highlight .m { color: #666666 } /* Literal.Number */\n",
       ".highlight .s { color: #BA2121 } /* Literal.String */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
460 
       ".highlight .na { color: #687822 } /* Name.Attribute */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
461 462 463 464 
       ".highlight .nb { color: #008000 } /* Name.Builtin */\n",
       ".highlight .nc { color: #0000FF; font-weight: bold } /* Name.Class */\n",
       ".highlight .no { color: #880000 } /* Name.Constant */\n",
       ".highlight .nd { color: #AA22FF } /* Name.Decorator */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
465 466 
       ".highlight .ni { color: #717171; font-weight: bold } /* Name.Entity */\n",
       ".highlight .ne { color: #CB3F38; font-weight: bold } /* Name.Exception */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
467 
       ".highlight .nf { color: #0000FF } /* Name.Function */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
468 
       ".highlight .nl { color: #767600 } /* Name.Label */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 
       ".highlight .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */\n",
       ".highlight .nt { color: #008000; font-weight: bold } /* Name.Tag */\n",
       ".highlight .nv { color: #19177C } /* Name.Variable */\n",
       ".highlight .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */\n",
       ".highlight .w { color: #bbbbbb } /* Text.Whitespace */\n",
       ".highlight .mb { color: #666666 } /* Literal.Number.Bin */\n",
       ".highlight .mf { color: #666666 } /* Literal.Number.Float */\n",
       ".highlight .mh { color: #666666 } /* Literal.Number.Hex */\n",
       ".highlight .mi { color: #666666 } /* Literal.Number.Integer */\n",
       ".highlight .mo { color: #666666 } /* Literal.Number.Oct */\n",
       ".highlight .sa { color: #BA2121 } /* Literal.String.Affix */\n",
       ".highlight .sb { color: #BA2121 } /* Literal.String.Backtick */\n",
       ".highlight .sc { color: #BA2121 } /* Literal.String.Char */\n",
       ".highlight .dl { color: #BA2121 } /* Literal.String.Delimiter */\n",
       ".highlight .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */\n",
       ".highlight .s2 { color: #BA2121 } /* Literal.String.Double */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
485 
       ".highlight .se { color: #AA5D1F; font-weight: bold } /* Literal.String.Escape */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
486 
       ".highlight .sh { color: #BA2121 } /* Literal.String.Heredoc */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
487 
       ".highlight .si { color: #A45A77; font-weight: bold } /* Literal.String.Interpol */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
488 
       ".highlight .sx { color: #008000 } /* Literal.String.Other */\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
489 
       ".highlight .sr { color: #A45A77 } /* Literal.String.Regex */\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 
       ".highlight .s1 { color: #BA2121 } /* Literal.String.Single */\n",
       ".highlight .ss { color: #19177C } /* Literal.String.Symbol */\n",
       ".highlight .bp { color: #008000 } /* Name.Builtin.Pseudo */\n",
       ".highlight .fm { color: #0000FF } /* Name.Function.Magic */\n",
       ".highlight .vc { color: #19177C } /* Name.Variable.Class */\n",
       ".highlight .vg { color: #19177C } /* Name.Variable.Global */\n",
       ".highlight .vi { color: #19177C } /* Name.Variable.Instance */\n",
       ".highlight .vm { color: #19177C } /* Name.Variable.Magic */\n",
       ".highlight .il { color: #666666 } /* Literal.Number.Integer.Long */</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
510 511 512 513 514 515 516 517 518 519 520 521 522 523 
       "<div class=\"highlight\"><pre><span></span><span class=\"n\">FUNC_PREFIX</span><span class=\"w\"> </span><span class=\"nf\">__launch_bounds__</span><span class=\"p\">(</span><span class=\"mi\">256</span><span class=\"p\">)</span><span class=\"w\"> </span><span class=\"kt\">void</span><span class=\"w\"> </span><span class=\"n\">kernel</span><span class=\"p\">(</span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\"> </span><span class=\"k\">const</span><span class=\"w\"> </span><span class=\"n\">_data_u0</span><span class=\"p\">,</span><span class=\"w\"> </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\"> </span><span class=\"k\">const</span><span class=\"w\"> </span><span class=\"n\">_data_u1</span><span class=\"p\">,</span><span class=\"w\"> </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\">  </span><span class=\"n\">_data_u2</span><span class=\"p\">)</span><span class=\"w\"></span>\n",
       "<span class=\"p\">{</span><span class=\"w\"></span>\n",
       "<span class=\"w\">   </span><span class=\"k\">if</span><span class=\"w\"> </span><span class=\"p\">(</span><span class=\"n\">blockDim</span><span class=\"p\">.</span><span class=\"n\">x</span><span class=\"o\">*</span><span class=\"n\">blockIdx</span><span class=\"p\">.</span><span class=\"n\">x</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"n\">threadIdx</span><span class=\"p\">.</span><span class=\"n\">x</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mi\">1</span><span class=\"w\"> </span><span class=\"o\">&lt;</span><span class=\"w\"> </span><span class=\"mi\">59</span><span class=\"w\"> </span><span class=\"o\">&amp;&amp;</span><span class=\"w\"> </span><span class=\"n\">blockDim</span><span class=\"p\">.</span><span class=\"n\">y</span><span class=\"o\">*</span><span class=\"n\">blockIdx</span><span class=\"p\">.</span><span class=\"n\">y</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"n\">threadIdx</span><span class=\"p\">.</span><span class=\"n\">y</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mi\">1</span><span class=\"w\"> </span><span class=\"o\">&lt;</span><span class=\"w\"> </span><span class=\"mi\">69</span><span class=\"p\">)</span><span class=\"w\"></span>\n",
       "<span class=\"w\">   </span><span class=\"p\">{</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"k\">const</span><span class=\"w\"> </span><span class=\"kt\">int64_t</span><span class=\"w\"> </span><span class=\"n\">ctr_0</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"n\">blockDim</span><span class=\"p\">.</span><span class=\"n\">x</span><span class=\"o\">*</span><span class=\"n\">blockIdx</span><span class=\"p\">.</span><span class=\"n\">x</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"n\">threadIdx</span><span class=\"p\">.</span><span class=\"n\">x</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mi\">1</span><span class=\"p\">;</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"k\">const</span><span class=\"w\"> </span><span class=\"kt\">int64_t</span><span class=\"w\"> </span><span class=\"n\">ctr_1</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"n\">blockDim</span><span class=\"p\">.</span><span class=\"n\">y</span><span class=\"o\">*</span><span class=\"n\">blockIdx</span><span class=\"p\">.</span><span class=\"n\">y</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"n\">threadIdx</span><span class=\"p\">.</span><span class=\"n\">y</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mi\">1</span><span class=\"p\">;</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\">  </span><span class=\"n\">_data_u2_10</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"n\">_data_u2</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"n\">ctr_1</span><span class=\"p\">;</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\"> </span><span class=\"n\">_data_u1_10</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"n\">_data_u1</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"n\">ctr_1</span><span class=\"p\">;</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\"> </span><span class=\"n\">_data_u1_11</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"n\">_data_u1</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"n\">ctr_1</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mi\">1</span><span class=\"p\">;</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\"> </span><span class=\"n\">_data_u1_1m1</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"n\">_data_u1</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"n\">ctr_1</span><span class=\"w\"> </span><span class=\"o\">-</span><span class=\"w\"> </span><span class=\"mi\">1</span><span class=\"p\">;</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"kt\">double</span><span class=\"w\"> </span><span class=\"o\">*</span><span class=\"w\"> </span><span class=\"n\">RESTRICT</span><span class=\"w\"> </span><span class=\"n\">_data_u0_10</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"n\">_data_u0</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"n\">ctr_1</span><span class=\"p\">;</span><span class=\"w\"></span>\n",
       "<span class=\"w\">      </span><span class=\"n\">_data_u2_10</span><span class=\"p\">[</span><span class=\"mi\">70</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span><span class=\"p\">]</span><span class=\"w\"> </span><span class=\"o\">=</span><span class=\"w\"> </span><span class=\"mf\">-1.0</span><span class=\"o\">*</span><span class=\"n\">_data_u0_10</span><span class=\"p\">[</span><span class=\"mi\">70</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span><span class=\"p\">]</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mf\">0.25</span><span class=\"o\">*</span><span class=\"n\">_data_u1_10</span><span class=\"p\">[</span><span class=\"mi\">70</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mi\">70</span><span class=\"p\">]</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mf\">0.25</span><span class=\"o\">*</span><span class=\"n\">_data_u1_10</span><span class=\"p\">[</span><span class=\"mi\">70</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span><span class=\"w\"> </span><span class=\"o\">-</span><span class=\"w\"> </span><span class=\"mi\">70</span><span class=\"p\">]</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mf\">0.25</span><span class=\"o\">*</span><span class=\"n\">_data_u1_11</span><span class=\"p\">[</span><span class=\"mi\">70</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span><span class=\"p\">]</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mf\">0.25</span><span class=\"o\">*</span><span class=\"n\">_data_u1_1m1</span><span class=\"p\">[</span><span class=\"mi\">70</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span><span class=\"p\">]</span><span class=\"w\"> </span><span class=\"o\">+</span><span class=\"w\"> </span><span class=\"mf\">1.0</span><span class=\"o\">*</span><span class=\"n\">_data_u1_10</span><span class=\"p\">[</span><span class=\"mi\">70</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span><span class=\"p\">];</span><span class=\"w\"></span>\n",
       "<span class=\"w\">   </span><span class=\"p\">}</span><span class=\"w\"> </span>\n",
       "<span class=\"p\">}</span><span class=\"w\"></span>\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
524 525 526 
       "</pre></div>\n"
      ],
      "text/plain": [
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
527 
       "FUNC_PREFIX __launch_bounds__(256) void kernel(double * RESTRICT const _data_u0, double * RESTRICT const _data_u1, double * RESTRICT  _data_u2)\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
528 529 530 531 532 
       "{\n",
       "   if (blockDim.x*blockIdx.x + threadIdx.x + 1 < 59 && blockDim.y*blockIdx.y + threadIdx.y + 1 < 69)\n",
       "   {\n",
       "      const int64_t ctr_0 = blockDim.x*blockIdx.x + threadIdx.x + 1;\n",
       "      const int64_t ctr_1 = blockDim.y*blockIdx.y + threadIdx.y + 1;\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
533 534 535 536 537 538 
       "      double * RESTRICT  _data_u2_10 = _data_u2 + ctr_1;\n",
       "      double * RESTRICT _data_u1_10 = _data_u1 + ctr_1;\n",
       "      double * RESTRICT _data_u1_11 = _data_u1 + ctr_1 + 1;\n",
       "      double * RESTRICT _data_u1_1m1 = _data_u1 + ctr_1 - 1;\n",
       "      double * RESTRICT _data_u0_10 = _data_u0 + ctr_1;\n",
       "      _data_u2_10[70*ctr_0] = -1.0*_data_u0_10[70*ctr_0] + 0.25*_data_u1_10[70*ctr_0 + 70] + 0.25*_data_u1_10[70*ctr_0 - 70] + 0.25*_data_u1_11[70*ctr_0] + 0.25*_data_u1_1m1[70*ctr_0] + 1.0*_data_u1_10[70*ctr_0];\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
539 540 541 542 543 
       "   } \n",
       "}"
      ]
     },
     "metadata": {},
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
544 
     "output_type": "display_data"
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
545 546 547 548 
    }
   ],
   "source": [
    "try:\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
549 
    "    import cupy\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
550 
    "except ImportError:\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
551 552 
    "    cupy=None\n",
    "    print('No cupy installed')\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
553 554 555 
    "\n",
    "\n",
    "res = None\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
556 
    "if cupy:\n",
Jan Hönig's avatar
`create_kernel` API Update  
Jan Hönig committed
557 
    "    gpu_ast = ps.create_kernel(update_rule, target=ps.Target.GPU)\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
558 559 560 561 562 563 564 565 566 567 568 569 570 571 
    "    gpu_kernel = gpu_ast.compile()\n",
    "    res = ps.show_code(gpu_ast)\n",
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The run function has to be changed now slightly, since the data has to be transfered to the GPU first, then the kernel can be executed, and in the end the data has to be transfered back"
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
572 
   "execution_count": 13,
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
573 574 575 
   "metadata": {},
   "outputs": [],
   "source": [
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
576 
    "if cupy:\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
577 578 
    "    def run_on_gpu(timesteps=1):\n",
    "        # Transfer arrays to GPU\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
579 
    "        gpuArrs = [cupy.asarray(cpu_array) for cpu_array in u_arrays]\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
580 581 582 583 584 585 586 
    "\n",
    "        for t in range(timesteps):\n",
    "            gpu_kernel(u0=gpuArrs[0], u1=gpuArrs[1], u2=gpuArrs[2])\n",
    "            gpuArrs[0], gpuArrs[1], gpuArrs[2] = gpuArrs[1], gpuArrs[2], gpuArrs[0]\n",
    "\n",
    "        # Transfer arrays to CPU\n",
    "        for gpuArr, cpuArr in zip(gpuArrs, u_arrays):\n",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
587 
    "            cpuArr[:] = gpuArr.get()\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
588 589 590 591 592 
    "assert np.isfinite(np.max(u_arrays[2]))            "
   ]
  },
  {
   "cell_type": "code",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
593 
   "execution_count": 14,
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
594 595 596 
   "metadata": {},
   "outputs": [],
   "source": [
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
597 
    "if cupy:\n",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
598 599 600 601 602 603 604 
    "    run_on_gpu(400)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
Markus Holzer's avatar
Remove LLVM  
Markus Holzer committed
605 
   "display_name": "Python 3 (ipykernel)",
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
606 607 608 609 610 611 612 613 614 615 616 617 618 
   "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",
Markus Holzer's avatar
Replace PyCuda with CuPy  
Markus Holzer committed
619 
   "version": "3.9.5"
Martin Bauer's avatar
Separated modules into subfolders with own setup.py
Martin Bauer committed
620 621 622 623 
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
Markus Holzer's avatar
Remove LLVM  
Markus Holzer committed
624 
}