{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from pystencils.session import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Tutorial 01: Getting Started\n", "\n", "\n", "## Overview\n", "\n", "*pystencils* is a package that can speed up computations on *numpy* arrays. All computations are carried out fully parallel on CPUs (single node with OpenMP, multiple nodes with MPI) or on GPUs.\n", "It is suited for applications that run the same operation on *numpy* arrays multiple times. It can be used to accelerate computations on images or voxel fields. Its main application, however, are numerical simulations using finite differences, finite volumes, or lattice Boltzmann methods. \n", "There already exist a variety of packages to speed up numeric Python code. One could use pure numpy or solutions that compile your code, like *Cython* and *numba*. See [this page](demo_benchmark.ipynb) for a comparison of these tools.\n", "\n", "![Stencil](../img/pystencils_stencil_four_points_with_arrows.svg)\n", "\n", "As the name suggests, *pystencils* was developed for **stencil codes**, i.e. operations that update array elements using only a local neighborhood. \n", "It generates C code, compiles it behind the scenes, and lets you call the compiled C function as if it was a native Python function. \n", "But lets not dive too deep into the concepts of *pystencils* here, they are covered in detail in the following tutorials. Let's instead look at a simple example, that computes the average neighbor values of a *numpy* array. Therefor we first create two rather large arrays for input and output:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "input_arr = np.random.rand(1024, 1024)\n", "output_arr = np.zeros_like(input_arr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We first implement a version of this algorithm using pure numpy and benchmark it." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def numpy_kernel():\n", " output_arr[1:-1, 1:-1] = input_arr[2:, 1:-1] + input_arr[:-2, 1:-1] + \\\n", " input_arr[1:-1, 2:] + input_arr[1:-1, :-2]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7.11 ms ± 137 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" ] } ], "source": [ "%%timeit \n", "numpy_kernel()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now lets see how to run the same algorithm with *pystencils*." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdkAAAAnCAYAAABDqDLfAAAACXBIWXMAAA7EAAAOxAGVKw4bAAANJ0lEQVR4Ae2d79XctBLGN+9JAQE64HaQWwKhAyAVJHQQDt/yLYfbAVBBgA7gVhBuOoAOCG8H731+WsvIsmVbXsuWd0fnePVvNJKekTweSWs/enh4OJmrD4HXr19/p1b91bTsX/K/Udo9cflP5T3X9VThz3V9qvA3usj/Q/Ef5Hec0kJ+0P2utPcdIos4BCKsLsI+4mW4j4yxCKuLcKeaiJ9hP4L9kqwI34vkFfG6KlndLQHXypRFQAPuV9XwVv5/uBT+WRdK0rvnSkepnuR/L+9r+V83mSEd+Z/q+kN57+R7fv9WvEPXlL15Txitgr3hnjeU1sKdWg37POyXUK8lr1uQlSnZJSOsYBkGndg/kx9amSjU/1Ftk4/SxH2sC1qncBV+oiu2YlEav4jmF/nefaUA6R0nmi8a/p30tSLi/WotXiX4NH1fC3vDfaaQVsadWqvBfiYEHTKbJ6dd7k8dIcyMzJHV45m8jGxjBCQ8bhRYsL8p/HlQPYr0pybOsrG3YE+ia8PkK46Vi9J+Q9w7pX/kw95X2kuF/5T/Z5BGXQx4LGWs30knOm8hs9T9ieL+AeCkMJb097o67ZxkujGB2ncR9vRRTb4E9ySGKShUZ1JWyrsJ3MFmD+xTMhlLN3ldLqsxfFN5Y7iPlEnOR/GbnFt3KcaWvg8CEhpKjiXiZ7q4WbPH6oV8Uvi9rntd5OO8wj3Hur8oSJT0fTe5G1M+N2j2dn/zOQqjwFG85HFNOpWB/i/57B/Th1/lo7BC97PS2v6EGXuH1a61sL8E9zkYdqBSu+fI6hZwB5dNse8IYmbE5NUCtVhWLYeMwEzcOxxVZs58HJ1bpmQ7kNYRkWBRUo/UGixYFN8rxWNF96XSsTzv5fdcQ0+ZWMn1aJWA0kOht07lUeYoytaybTPTAfiEipowy69Ydc4p7NPi/niSXX217yLsVZ5+LcZdZScxjAFSnZOyunbcwWQP7GNZzImbvFaR1RyoOzRzcO8UOEcm56P4jt7TTMkOoLpXkoT1nS6/33pCeLrCpeKwaViy4T5rmEfZ+ybB+3G+t4RJRxG2yrFDODOi8ihSlEuslKkfSyt0blCGCXuH1f5VsL8E90wMl0B2tbgDRuXYm7wCBC6RVcCmaDBzPibnlinZomLKZu6XiNuCEvQXivwQDMpTI3yU2pSVyiEoLN6OU3mezj6Q2PAaVMSdQtMRFOyQo56Po4x3iqceHiLSzaJrYr8U9xwMlwBz7biDSa3Ym7z6CCyVVZ9TmZSc+ZicW4/LtM24LkTghcphVfq9UNiw/xofFELBslQ8an1SThcWGkrVW8gMnDdKu5ePg5dTuC62/g8KNh6s1E29NbnVsC+A+xCGS7C7atwBpGLsTV4RAgVkFdVQLDo0H5Nzy5RsMTnkM9ag4287XKNOdChX/vw96UTbnu5NEKMAGSCXuhQP+MdLyCWV+qJ+rI39QtxzMFzSz6vHHVC2xF51PVWVPMROOR5443kwVebq5bVQVqfCuHu55MzHpKxGlWzTkR9VI1bHT4rHFpVvjPnHRYCJjyK8yGls+ENYjJX3EbM4Dk3uDSdiefhoD/dMDJcAYLifUVsNe8mMsV1q68Pkdb5P9O5PhXF3oyRzPiZldTc2U+mILv4fSSen9v8GWan8E138DYVGmKsMAcmFm0S2bBq5xk/w/B+XvU3nREOYF2HECpX6/PL1mfjGfkdwH8Uwgftc9G4ed4DaCfu5Mgrpbl5eI7IKcVolnJhbo/MxqDgpq7uAaDCoiv1Nc3T/b7DwOZHyNODDCI1l7YsAVigyah1xXSjRb3W5sOKvWoKzTF8qjYNZzinMX354AQX7wNBy6Iq9ztjx5J88GR0TT8VVl9/HniKtLX8I9ykMkVOM+5SsfL8Nd4/E+UxDPOazsf+H3fyQxuvm8jrwHAHY3jyZj/Y/lDNwH5pbU2PCV5CcW489xYhPYb8UOEKWzLq0fJLxlhkSEIeRrvWl+ihTtgLa/Vv1F+uzjSvcccrHAv5IfqtkIVA8WabJd0s/ooutW7KXOnj2lpSWMtuwXA936h7DUHk93BssDfc8wa2CfV6VZ+qd5HXUOQJog7LKxX4K96G5RR1Kv2hu3c1oKJboUisW9pRfzWqB4U4OJbSmYtipG/1qNYiQL8v6nSf7PmU3RfRLxoafMF1mNxgz3PcT+sbYL+mozZMGtaWyWgJ6iXtaR8mqApYx3Cui5POWIU7OcQ3uxyrfLQs2tLyT1lspLN/Bh3LcuPkkG/FwuVHJx3BqN8rkEmu++o6qjzxEcAIyxyJkvNzP7ZxoWQ1gHFzlw8pcHEI6YWG4h4BsGN4C+yXdsXnSR22hrPqMplNWv6e1SladQJHwpRdeK+ffPcvTFK5nyYoGBfpWPi9IZt2al9k7esV5UxH7ce5VfQrzXtwvdUF3KKc2o3To14tDNXxBY9XX0WWRmKXo+TN5juOEem8s5TC4RlrDfT+pboD9ks7ZPBlALVdWAywmk1TH6vc0p2TFGEWCkkS5hlYGYfdC+rB1osE6xVplf8g7btDuc2w+QT77sSFNkHUOike1n1dT27Di/6vrhcKzLbZeJw+UULKfJXkfCOLBppbEpiTvwc4cLLE2fGprT03irA2bOe159PDwcBIhFudL+byUvnWK8zcL/oLRsXAUR8mSh1WCcsZyDZWzktyG8WB5l3nOZ/mQZdjWulHYW8+9T6X5crE/VUb52Z9Xa3iyvI1y/RDXmRHnweUa9qQzulyGFDmKMysusfu4SRiSEw+JvVdLxgwsnkbAcE9jU1uOyao2iZxOXsn+raZxchbL0zmFsW5J73wC7ZzrFCjKMNxjZdm4VcYzysP/R9G1N0CF3esE5btlZfncUFFSbbt8/d6fU6bhQz/a9vnyY37Dmz2zzxQuasmK/8NYW24hTxh0HvLm9lnlOOHMXsqi7QjD3s3nbOwN97kjdF064b65rOiBzZNlcnws4FB2XPHhpq9gqfzWyiTundJRWChApwgV5qBU+E5cZ3GE5RXmBKtXVihpt2frecon7TMfp6wuvknKDbRnKTd0k2UaPhzSCuv31SR90fNi/t9FwJLxrI+WJ5lNZKie7IkzwdKyZyJg2M8EamUyw31lQAuzM3ktA9jtyTZFYyXW7qcK3HbfVGGUFcvAzimc+hxbWx5C0Tlr41zK/T6jrI8rzBI0yj5uB0qZvdGeyyxDXUNLjT2+YYLqYE+Z5WaUuTlDwBAwBAwBQ2A2AndSHigxFBBKzjmloRBRbFhxOJZavfJDUXUs0Ia+8zk20bBP5sooH+WJNeoOQcmnLuoNHTRDjn02v+cW5+eUSX6KKGYax9VeTpx1XroQ01jcEDAEDAFDwBCIEXjcJLAvyv4oe6yf6HqrC0vU/Q9WfqhUXyiOFer2TxXG3SvO3mXoWE7GAoTuJD/cL0PJDh1SgTR2KNiUMo1pfXyoDEq9fZDwhBk+B8A61ndGWSM1BAwBQ8AQuEEEnJKV8kABtQeQAhx6aaLFGnUWaUDXC4oOKxZFPeRQmtQZujju86D1VrRP835OmblK3fPu+OpP1qGpTmGLGAKGgCFgCNwkAnc79RqlifJsXaOUUZpD1uagUs8sA9+Usm7bccsB4cmbucJVi6PAwbhJPXAdog8Hxf7wuI8NjoPKJNUlk1UKmcLpfrm4cDVd9hq8/HdxSJm+ESV7vk6pioYwy7ROOcp/ovi38kOrcrSM6L2jvvbAlk80v4MA/3kefKDpUFUW0XhoD9BV1rSc5hwO+yvBfUxGh5NJqjMmqxQy5dPvyleRrKH3+SINBPZtxz6VhqJ8Kbr2ENKMMr4BLF3bSyE8GpEvHMP/PEe5Fi2JgGFfEt1lvE0my3Dbo1TtstrFkm0EwV9iOCwVWqUnAdaJh0JTHlZWLZ9XC5t26LBwfaoOXPVyUq0CMuzrk4zJpD6ZpFp0BFntZskKHJb47PNqqdGzbfpzySP3xdjbtvB6azPs65OtyaQ+maRaVL2sdlOyIKYbO5asfV4tNXw2SJcMWCY+4mGnDdApW4VhXxbfJdxNJktQ26fMUWS1q5JFNAIquTw8JDrR51pc9tmoISDP2LPHzX+c7dR1AqNSycLcsC8F7kK+JpOFwO1Q7Eiy2l3JIh8Bxn5gEVeSd5EGb8uUVYTch5ZtW3i9tRn29cnWZFKfTFItOoysqlCyKRQtvRwCUq68icuWictBnORs2Ceh2S3DZLIb9NkVH01WpmSzRXz8AhqkLFVy6MyWiTcWp2G/MeAzqjOZzACpEpIjysp9T7YS/KwZGyGggcpLPjh0Fjv+f8zSPSe/+R9z1n55zMzifQQM+z4me6eYTPaWwPz6jygrU7Lz5Xv1lBrAf6uTfLqw987qq+/8zh007HcWwED1JpMBUCpNqllWtlxc6aDZqVm8tpLL3PYIGPbbYz5Vo8lkCqF68quVlVmy9QyS3Vqip0AOQLFPyzIyjtdPvlN6+HlCl2E/6yJg2K+L5xrcTCZroLgNjyPI6v9pjxUdD2V9jwAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle {{dst}_{(0,0)}} \\leftarrow \\frac{{{src}_{(1,0)}}}{4} + \\frac{{{src}_{(0,1)}}}{4} + \\frac{{{src}_{(0,-1)}}}{4} + \\frac{{{src}_{(-1,0)}}}{4}$" ], "text/plain": [ " src_E src_N src_S src_W\n", "dst_C := ───── + ───── + ───── + ─────\n", " 4 4 4 4 " ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "src, dst = ps.fields(src=input_arr, dst=output_arr)\n", "\n", "symbolic_description = ps.Assignment(dst[0,0], \n", " (src[1, 0] + src[-1, 0] + src[0, 1] + src[0, -1]) / 4)\n", "symbolic_description" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMQAAADTCAYAAADedbxIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAJ/UlEQVR4nO3cf2jU9x3H8dcnl8SL8aLV09hpbeta/ymMunWyVufagZWIPwgbxT8srLKNVtuZjZU66B+CcxUGXaDgNoaw1spa21VQ8AyDCjZaGLhNGJu/t4Qq/jhqYrTxx10+++OS2zuXXO5yuXy/t93zAYHc3febe6PfZ+5733zv67z3ApBRE/YAQCUhCMCoDXuAauMSyaik5ZIWSoqMsehtSScknfQtcfZrA+J4DxEcl0jOk/SupC+NY7UOST/2LfH05EwFi12mYL2s8cUgSSslfWsSZsEoCCJY3yhxvSfLOgXyIohgNZa43tSyToG8CCJYbsQ9fddrtPnpBWpd8KjOnawvej1MCoIIW7RxQNvfv6glz/aFPQoIInx19dLMZo4gVQiCAAyCAAyCAAxO3agEW1vnqetUVJcu1Gvlhh6t3ngj7JGqFUFUgp37L4Y9AjLYZQIMggAMdpnCtmr2oryPHbp2JsBJIIII39BG37E3pt3b5mjf2fMhT1TV2GWqBANp6djBmGbNTYU9SrUjiErQsbdJS9f0yXEOX9gIImzplNR5IKYV6zm5rwIQRNgO72nSsrV9qhnr49UICkGErft0vT7eN1Ovrl6kK911am+bE/ZI1YyjTGF78Y2kLv+7Sem0tOMFr7b2q2GPVM14hQjbnf4GDQxk9pdef9uJq6CEiiDCdrPnPnk/eHjJS/03p4U7UHUjiGAN//WfTtfoTv9/LzzgfY1u9txXcD1MGoII1q3ht3qnj1ji3t2oUndz39vdGrEcJgVBBOt49jvvpVu9Zncpe7/Tzd4ZOesdC2A2iCCC9pak7uytmtqUIpGUnMvsEkUimdvDHZTUGdiEVY5ruwbMJZJ1ylyJb6GGDnt/cmCx/nJkk7b86gdm0duSTviW+D+Dn7J6EUQFcM61SvrI5+4+IXDsMgEGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQAGQQBG7VgPukQyIsnledj7lni6/CMBE1fqtuu897k/aL6krZKWSWoo8LxJSQlJv/Qt8TvjmhhZzrlWSR957/P9B6IILpH8sqRXJT0laUqBxa9IOijpTRvHsFeIwareljS/yBnikp6X1CjpZ0WuA5SdSySjkt5RZpssRrOk70uql7Rj6M7c9xBPqPgYrFWDAwFhWabiY7DWuUQy20FuEA+UOExU0uwS1wXKYUGJ602XFBu6kRtEZMTifddrtPnpBWpd8KjOnawf4wdzxAphGrn9Fb/tZt86FN6Io40D2v7+RS15tq+kMYGwlLDtFg6irl6a2czhVfzvKWHbZTcHMAgCMAgCMMY8dSNra+s8dZ2K6tKFeq3c0KPVG29M8lxAeYxz2y0uiJ37L45nBudcg6TvSvqppE7v/ebxrA+MxjlXJ+nvg1/tymxbfsyVxrntlnWXyTn3mHPuN5KuSdol6SuSFpfzOVDVIpIekdQq6ZCkLufcT5xzM8v1BMW9QhRys7dJbSsSkh6SVJfzc+udc81leZ7/X9MliX+ngqKSvDK/yKcNfm2XtEPvvPEPPbelQdGp/RN5gsJBrJq9KO9jh66dkST1XJur1L25eZb6mqTLJcxWjfh3Gr+pkqTrVxcreclp/iNnso8Us+3mKBzE0Iode2PavW2O9p09P2KZWXM/05SGc5KeVKZge6Lfp977pwo+TxXj9O/iOOeikm5q+ClGmb9C3//Qn9X8wMJhKxSz7eYo7j3EQFo6djCmWXNToz7eMO0L/fb495TZZdoh6Wp2UKC8nKR7kvolnZD0Q0mz9dyW91Q35d6IpQttuzmKC6Jjb5OWrumTG/sXmPf+svf+55Lul/QdZT489HFRzwEUlpJ0XNKvJT3uvX/Ce/+e9z7/h9OK3HaHFA4inZI6D8S0Yn3Rv/G99wPe+z9571d5718vdj1gLN77lPf+m977Ld77Ud8DDFPCtls4iMN7mrRsbZ9qRp4ZDlS0ErbdwkF0n67XkQ+a9Nq6+brSXaf2tjkTmREITAnbbuGjTC/tTGa/37T8QbW1X53YlEBASth2x/eX6l1Hu8Y/FVABitx2OdsVMAgCMHKDKOqPF3lMZF1goibyMefsH/Ryg/hXiT/wljJnuAJhuVDiekmZsypyg/irpNMl/NA/+pb43RIHAsrhuKTuEtbb51vi2c9UDDvs6lvi3iWSL0h6RZkroTUp/wVjB5Q5OzMh6XclDAKUjW+J33OJ5POSXpa0VJnLq+bbdtOSLilzbdff2wdGXOwYweNs18rBUSbAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAqA17gGrjEskmSSskLZQUkSTt+PAx/e2oXCK51SzaL+mEpE99Szwd+KBVynnvw56harhE8mFJ70qKy3vpctfDGkjXynuXWcBl/jMapvVqZvPVwdU6Jb3oW+L3wpi52rDLFKzNkuKSJJdpIBuD/T4Ssa8IyyQ9E8x4IIhgfX3YrWkzPs++KliN03ty7lkyiTPBIIhgTR12qzF2Y8QSdVP6VVuX+56hYRJngkEQYaqJeKXTffrFRumVb0uXLgwoNuPzUZZ0o9yHSUAQYZvVfF0/etPrq89Icl7Rxi/CHqmaEUTYGqff0fR4SpIUndqXfbONUBBEJZg2uJs0paEv5EmqHn+YqwSxGb2qrWvIOdyKEPAKARi8QlSCra3z1HUqqksX6rVyQ49Wbxx5OBaBIIhKsHP/xbBHQAa7TIBBEIDBLlPYVs1elPexQ9fOBDgJRBDhG9roO/bGtHvbHO07ez7kiaoau0yVYCAtHTsY06y5qbBHqXYEUQk69jZp6RpO26gABBG2dErqPBDTivWctlEBCCJsh/c0adnaPtVEwp4EIojwdZ+u15EPmvTauvm60l2n9rY5YY9UzTjKFLaXdiaz329a/qDa2q+OsTQmGa8QlWTX0a6wR6h2BAEYBBGsUi+CxcWzAkIQwSr10CqngweEIIL1SYnrdZZ1CuRFEMF6S9Lpca7zB0nHJ2EWjIJruwbMJZJO0uPKXOx4rMPetyWd8C3xzwIZDJIIAhiGXSbAIAjA+A9JIWcDKWIk7QAAAABJRU5ErkJggg==\n", "text/plain": [ "