{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "nbsphinx": "hidden" }, "outputs": [], "source": [ "import psutil\n", "\n", "from pystencils.session import *\n", "\n", "import shutil" ] }, { "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": { "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", "text/latex": [ "$\\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}}$" ], "text/plain": [ "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 " ] }, "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": { "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", "text/latex": [ "$\\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}}$" ], "text/plain": [ " 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 " ] }, "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": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
FUNC_PREFIX void kernel(double * RESTRICT const _data_u0, double * RESTRICT const _data_u1, double * RESTRICT  _data_u2)\n",
       "{\n",
       "   for (int64_t ctr_0 = 1; ctr_0 < 59; ctr_0 += 1)\n",
       "   {\n",
       "      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",
       "      {\n",
       "         _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",
       "      }\n",
       "   }\n",
       "}\n",
       "
\n" ], "text/plain": [ "FUNC_PREFIX void kernel(double * RESTRICT const _data_u0, double * RESTRICT const _data_u1, double * RESTRICT _data_u2)\n", "{\n", " for (int64_t ctr_0 = 1; ctr_0 < 59; ctr_0 += 1)\n", " {\n", " 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", " {\n", " _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", " }\n", " }\n", "}" ] }, "metadata": {}, "output_type": "display_data" } ], "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": [ "{u0: double[60,70], u1: double[60,70], u2: double[60,70]}" ] }, "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": [ { "name": "stdout", "output_type": "stream", "text": [ "No ffmpeg installed\n" ] } ], "source": [ "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\")" ] }, { "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", "We can also run the same kernel on the GPU, by using the ``cupy`` package." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
FUNC_PREFIX __launch_bounds__(256) void kernel(double * RESTRICT const _data_u0, double * RESTRICT const _data_u1, double * RESTRICT  _data_u2)\n",
       "{\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",
       "      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",
       "   } \n",
       "}\n",
       "
\n" ], "text/plain": [ "FUNC_PREFIX __launch_bounds__(256) void kernel(double * RESTRICT const _data_u0, double * RESTRICT const _data_u1, double * RESTRICT _data_u2)\n", "{\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", " 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", " } \n", "}" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "try:\n", " import cupy\n", "except ImportError:\n", " cupy=None\n", " print('No cupy installed')\n", "\n", "\n", "res = None\n", "if cupy:\n", " gpu_ast = ps.create_kernel(update_rule, target=ps.Target.GPU)\n", " 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", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "if cupy:\n", " def run_on_gpu(timesteps=1):\n", " # Transfer arrays to GPU\n", " gpuArrs = [cupy.asarray(cpu_array) for cpu_array in u_arrays]\n", "\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", " cpuArr[:] = gpuArr.get()\n", "assert np.isfinite(np.max(u_arrays[2])) " ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "if cupy:\n", " run_on_gpu(400)" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.5" } }, "nbformat": 4, "nbformat_minor": 2 }