01_tutorial_getting_started.ipynb 223 KB
 Martin Bauer committed Mar 21, 2019 1 2 3 4 5 6 7 8 { "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [  Markus Holzer committed Feb 10, 2022 9 10 11 12 13  "import pystencils as ps\n", "from pystencils import plot as plt\n", "\n", "import numpy as np\n", "import sympy as sp"  Martin Bauer committed Mar 21, 2019 14 15 16 17 18 19 20 21 22 23 24  ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Tutorial 01: Getting Started\n", "\n", "\n", "## Overview\n", "\n",  Martin Bauer committed Apr 28, 2019 25 26 27  "*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",  Martin Bauer committed Mar 21, 2019 28  "\n",  Martin Bauer committed Apr 28, 2019 29  "![Stencil](../img/pystencils_stencil_four_points_with_arrows.svg)\n",  Martin Bauer committed Mar 21, 2019 30  "\n",  Martin Bauer committed Apr 28, 2019 31  "As the name suggests, *pystencils* was developed for **stencil codes**, i.e. operations that update array elements using only a local neighborhood. \n",  Martin Bauer committed Mar 21, 2019 32  "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",  Martin Bauer committed Apr 28, 2019 33  "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:"  Martin Bauer committed Mar 21, 2019 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 72  ] }, { "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": [  Markus Holzer committed Feb 10, 2022 73  "4.65 ms ± 22.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"  Martin Bauer committed Mar 21, 2019 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95  ] } ], "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": {  Markus Holzer committed Feb 10, 2022 96 97 98 99 100 101  "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": [ "Assignment(dst_C, src_E/4 + src_N/4 + src_S/4 + src_W/4)" ]  Martin Bauer committed Mar 21, 2019 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122  }, "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": {  Markus Holzer committed Feb 10, 2022 123 124 125 126  "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMQAAADTCAYAAADedbxIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAJ/UlEQVR4nO3cf2jU9x3H8dcnl8SL8aLV09hpbeta/ymMunWyVufagZWIPwgbxT8srLKNVtuZjZU66B+CcxUGXaDgNoaw1spa21VQ8AyDCjZaGLhNGJu/t4Qq/jhqYrTxx10+++OS2zuXXO5yuXy/t93zAYHc3febe6PfZ+5733zv67z3ApBRE/YAQCUhCMCoDXuAauMSyaik5ZIWSoqMsehtSScknfQtcfZrA+J4DxEcl0jOk/SupC+NY7UOST/2LfH05EwFi12mYL2s8cUgSSslfWsSZsEoCCJY3yhxvSfLOgXyIohgNZa43tSyToG8CCJYbsQ9fddrtPnpBWpd8KjOnawvej1MCoIIW7RxQNvfv6glz/aFPQoIInx19dLMZo4gVQiCAAyCAAyCAAxO3agEW1vnqetUVJcu1Gvlhh6t3ngj7JGqFUFUgp37L4Y9AjLYZQIMggAMdpnCtmr2oryPHbp2JsBJIIII39BG37E3pt3b5mjf2fMhT1TV2GWqBANp6djBmGbNTYU9SrUjiErQsbdJS9f0yXEOX9gIImzplNR5IKYV6zm5rwIQRNgO72nSsrV9qhnr49UICkGErft0vT7eN1Ovrl6kK911am+bE/ZI1YyjTGF78Y2kLv+7Sem0tOMFr7b2q2GPVM14hQjbnf4GDQxk9pdef9uJq6CEiiDCdrPnPnk/eHjJS/03p4U7UHUjiGAN//WfTtfoTv9/LzzgfY1u9txXcD1MGoII1q3ht3qnj1ji3t2oUndz39vdGrEcJgVBBOt49jvvpVu9Zncpe7/Tzd4ZOesdm/zRIBFE0N6S1J29VVObUiSSknOZXaJIJHN7uIOSOgObsMpxbdeAuUSyTpkr8S3U0GHvTw4s1l+ObNKWX/3ALHpb0gnfEv9n8FNWL4KoAM65Vkkf+dzdJwSOXSbAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAIAjAqB3rQZdIRiS5PA973xJPl38kYOJK3Xad9z73B82XtFXSMkkNBZ43KSkh6Ze+JX5nXBMjyznXKukj732+/0AUwSWSX5b0qqSnJE0psPgVSQclvWnjGPYKMVjV25LmFzlDXNLzkhol/azIdYCyc4lkVNI7ymyTxWiW9H1J9ZJ2DN2Z+x7iCRUfg7VqcCAgLMtUfAzWOpdIZjvIDeKBEoeJSppd4rpAOSwocb3pkmJDN3KDiIxYvO96jTY/vUCtCx7VuZP1Y/xgjlghTCO3v+K33exbh8IbcbRxQNvfv6glz/aVNCYQlhK23cJB1NVLM5s5vIr/PSVsu+zmAAZBAAZBAMaYp25kbW2dp65TUV26UK+VG3q0euONSZ4LKI9xbrvFBbFz/8XxzOCca5D0XUk/ldTpvd88nvWB0Tjn6iT9ffCrXZlty4+50ji33bLuMjnnHnPO/UbSNUm7JH1F0uJyPgeqWkTSI5JaJR2S1OWc+4lzbma5nqC4V4hCbvY2qW1FQtJDkupyfm69c665LM/z/2u6JPHvVFBUklfmF/m0wa/tknbonTf+oee2NCg6tX8iT1A4iFWzF+V97NC1M5Kknmtzlbo3N89SX5N0uYTZqhH/TuM3VZJ0/epiJS85zX/kTPaRYrbdHIWDGFqxY29Mu7fN0b6z50csM2vuZ5rScE7Sk8oUbE/0+9R7/1TB56linP5dHOdcVNJNDT/FKPNX6Psf+rOaH1g4bIVitt0cxb2HGEhLxw7GNGtuatTHG6Z9od8e/54yu0w7JF3NDgqUl5N0T1K/pBOSfihptp7b8p7qptwbsXShbTdHcUF07G3S0jV9cmP/AvPeX/be/1zS/ZK+o8yHhz4u6jmAwlKSjkv6taTHvfdPeO/f897n/3BakdvukMJBpFNS54GYVqwv+je+937Ae/8n7/0q7/3rxa4HjMV7n/Lef9N7v8V7P+p7gGFK2HYLB3F4T5OWre1Tzcgzw4GKVsK2WziI7tP1OvJBk15bN19XuuvU3jZnIjMCgSlh2y18lOmlncns95uWP6i29qsTmxIISAnb7vj+Ur3raNf4pwIqQJHbLme7AgZBAEZuEEX98SKPiawLTNREPuac/YNebhD/KvEH3lLmDFcgLBdKXC8pc1ZFbhB/lXS6hB/6R98Sv1viQEA5HJfUXcJ6+3xLPPuZimGHXX1L3LtE8gVJryhzJbQm5b9g7IAyZ2cmJP2uhEGAsvEt8XsukXxe0suSlipzedV8225a0iVlru36e/vAiIsdI3ic7Vo5OMoEGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGAQBGLVhD1BtXCLZJGmFpIWSIpKkHR8+pr8dlUskt5pF+yWdkPSpb4mnAx+0SjnvfdgzVA2XSD4s6V1JcXkvXe56WAPpWnnvMgu4zH9Gw7RezWy+Orhap6QXfUv8XhgzVxt2mYK1WVJckuQyDWRjsN9HIvYVYZmkZ4IZDwQRrK8PuzVtxufZVwWrcXpPzj1LJm8kWAQRrKnDbjXGboxYom5Kv2rrct8zNEziTDAIIkw1Ea90uk+/2Ci98m3p0oUBxWZ8PsqSbpT7MAkIImyzmq/rR296ffUZSc4r2vhF2CNVM4IIW+P0O5oeT0mSolP7sm+2EQqCqATTBneTpjT0hTxJ1eMPc5UgNqNXtXUNOYdbEQJeIQCDV4hKsLV1nrpORXXpQr1WbujR6o0jD8ciEARRCXbuvxj2CMhglwkwCAIw2GUK26rZi/I+dujamQAngQgifEMbfcfemHZvm6N9Z8+HPFFVY5epEgykpWMHY5o1NxX2KNWOICpBx94mLV3DaRsVgCDClk5JnQdiWrGe0zYqAEGE7fCeJi1b26eaSNiTQAQRvu7T9TryQZNeWzdfV7rr1N42J+yRqhlHmcL20s5k9vtNyx9UW/vVMZbGJOMVopLsOtoV9gjVjiAAgyCCVepFsLh4VkAIIlilHlrldPCAEESwPilxvc6yToG8CCJYb0k6Pc51/iDp+CTMglFwbdeAuUTSSXpcmYsdj3XY+7akE74l/lkQcyGDIACDXSbAIAjA+A9JIWcDPN19qQAAAABJRU5ErkJggg==\n", "text/plain": [ "