{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "nbsphinx": "hidden" }, "outputs": [], "source": [ "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": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAAyBAMAAAAzRVApAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEHaZ70QiZs27Mondq1QBvk6oAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAI5klEQVR4Ae1bXYgcRRCu2c3u3e3ubcbom4asiRJIQE+DohjxJ+IPSDw0KAGDAz6IEPEkoIJ5OBBkHxTPH5A8SAYioviQNYIICXgG/0AkywkhBCGnYkQxJhHiQ2JyVnV1z85P93TPRc3NkeJueqenvvq+rtrZma27AZBWp7Gh9hzHgPxC2rhZQG4hbZwtIM+QNkUsIOeQNoUtIERIGycrmreAooa0MVsP1q6F9rjZQXOkCfA+wM+aI/qpJjRvvq6AP0YpSsHE80MJbFFoDxpTBfLmEn40GOq1+3CjPomG2YOwaz1ANTAczkwfhI/hEFSDzAHzRFEKjjQ/lMAWhI4GjXXT4J43l/AfQGscTsJyc1ayRxo9gEcBvH72kHYG/b+E3YGzPwYpSsG880MJbFHoBwC1aXDOm1P4w7AihI0wHPBqnLajoagFvOTkDYD+X8DT487+GLYoBSuZH0pgi0IPcy1c8+YUfhrfsPAVVCYc80puLfzF8wJW046DkT8c95390bsoBauYH0pgi0KnuRaueXMJX8da+LAf8JJhMW/djlC6bMWRanFM7tsG8kcOZ3/0LkrBEoqjGmOMLEpYn+Za2PPmrqw2RefFnYCjxZZD85R02Ywj1YJGFyO/Ji7a1R/di1KwjOKolnoPFoRSvuh6Yc+bu7Jan64Xp6AxyRjz9jqAu+TRLThSLd4zerc78UPkvw1/zf5xb/HaTsGQLFG+sDTRVaoWdsIEFeZN1CInbz/Guezh8X5liu+jKHa+3Q9wJGAX/G4hakEEekvopu8itV4tBLN/JoqdgiEZIouwFFH9ebVuO2GCCvPG54XCpwLjbqIW9vB4XzoNlTH80BvqZYMlZ1714UjIU1fg8Dv+vsC7mm1CN6D/e3/86uf4Z0LYKRiSIbIISxE1h1Uu7YQJKswbtPE3J2+JWuSH91Z14DIfPgNY+xF+bxtPqdTtHvCXXQ/eHdDy4fEDMwA36JzEnNId+T8wN5fnHwsUQSwUDMkQOaEiumexFq6EkirKW3392U5e3mQtXMJXa7fDtQDXsK5hGkSLReoM5RgbamdgbBNUTkBllmcnadCCVIo0/gGB2LRdsHyIRIZyBDeiNEqhvQ7WIp9woFFSOedN1sIl/BPVSdgBMMK6fqFBtKZwrOzp6ZpHrX594iYYmQLAH7T2BG21IKlb498U0RFn6IJpIcSDRh/SKXUORIRBw0/sTAOtCcO6Nek1SqqcvCU5uBb69dDlI5YCf0kH72WhJhLqjeEh2Zqi68AxqAY4A15XWEivYR14iEAYPCL2L6WtBtTsdl++t9vtIzzrf1BEpw9afRdMCxF0IBpAcXVORIyNNdBia/oNa6ElpAwkNA6ocvK2a33UpPuk2/28233FkAJyTIR/zG+cxpmjNF3xcSNbU4ANi6V+tnlUwXrVTsLuEM+IAN3hStroQfI9lPXHvgxFxy/Wpi5YhgIh0uiGPqXOTiSxdH+VXpPXw1o4a1RU5rwlOeRnlHY9KCaRgn1QuR1nBiZbU+CdhiWz2ebRT1APmtNw1QBAr/QgpTvjPxpydHMXTANRjKIWSXV2IgWm5acaaO3v1t22oaMh1GpUVMa88Y1+xCFroQnPjvEUvI0f/c8opTTK1hQ0/oLWRKZ5hG/PSjDShx1DdA5Fpgcp3Rn/loxu7oJpIIpM1CKpzk6kwFSL1WonGpf0QUPIHCmNisqYN05xxCFroQkva4END9UIvBU2jUWnP4qry9YUNE5Cq5NpHm3d+9H3UO3XNuBfRQZmACndGf+tMrq5C6aBKDquRUKdnUiBqRbH1E40LpnKrsmgUVEZ88YpjjhkLbTrEZ9RuweNwOUrj86EkSi8FsjWFL1z8QJNDZq4vTo39zd4M9c8tSY+awAp3Rn/zTK6uQumgShCdV7E1NmJFJiWn14TNHeendAQcgZSGhWVMW9ci4hD1kITfnBemBqBqjVF14ul447NIwNI6VaJiMYtMrprFwyB2DnZvAPNF02HlDo7UY2wr/Hy3RpiBo0mqigFKQ5Zi2jp6oVqTeWlAO/exc0N36k4No8MIC9xTVEqcMQb69h9FC7DbnwvTn7ivGC8am3ZiRQBnRcKpeb0o0GjiWqQArqPGnBcoo9OKRCOnG19CjxuTdVfg1X0wWpuNsU5ioKuABH9Sd+xC4ZcCJGGDSBndRGRAuc30JQXjRHUTeMgBW69MFoPislPAbemGhuhsnfMsXkEsp/lDMJWFkV/Z9a9C4YQNmoAzYOIwQX6VEU1Rnlz45DdvPwUyNZUm8VP8mDbFgRVZjkgD6ILZmNQEPZzVadQTCQ5JuWYPyShdo3/RQpGWKL4Ii6bTfmi6WhR0JQI2Q7E8IvY2jYMkV7O6hJEDOYGmo0uarcFwtOu8b9IAbemYEJIEM0mu2zZz3IHcStrSEQWXTA7B0Okn7O6OJHEuq4pDnXQWDRvTuH5PcfCRbPJnifZz5KODiB5RgiA6ILhK2y8kY0Z6OKQ+RG5o4RnnFBpNGgT0wXzVjR8HvXFYxczsHgzMLcw7EQqw0sXhqz/VUUqBRd0V14v+hdUxEXyixn49zJAz1UsIiv1cui5ikVkpV4OPVexiGpR0uW8v4FqIJ6rKH0xGm/K75llXc45wKeNxHMVpa8FLJ1Vazjuq1clGvF/Tdr05WB/iTQbpT4WVaCUy2lOQ3WSn6swLrE0B/YppfSYSPlseAxaPX6uonzi04pfVBPb1IsSjY/v2T370M43OuK5ihLp1kmt37LmHHgrDz64hx8T0fks4Lna3bDTh+PyuYoFLNRF2qZg5CRcDts7X/NjIi6YBeRz2QTcBPQjnqtYQMLmIcW7hy58H8Kn/sOlXM5tPnwD9LMIrIr/iNgDnx6AKKXdB94p+lkMhv/GuWIcF/JnORfjnaHvFrUT+G4qv63o0YUPxAMQJVyNdw5GJ482J9thCcWnJa+YhUPt4Ad89s3hz/Vp8ALYfx2293vV/rIFIOW8JYx0hr5t1s6NnKhPnHesCxFg+cp3Z8brb5VTfCph3syap672Zq5c9VzqwP+2+w+M+JTHYLF/LQAAAABJRU5ErkJggg==\n", "text/latex": [ "$$\\frac{{{u^{(0)}}_{(0,0)}} - 2 {{u^{(1)}}_{(0,0)}} + {{u^{(2)}}_{(0,0)}}}{dt^{2}} = \\frac{{{u^{(1)}}_{(-1,0)}} + {{u^{(1)}}_{(0,-1)}} - 4 {{u^{(1)}}_{(0,0)}} + {{u^{(1)}}_{(0,1)}} + {{u^{(1)}}_{(1,0)}}}{dx^{2}}$$" ], "text/plain": [ "u_0_C - 2⋅u_1_C + u_2_C u_1_W + u_1_S - 4⋅u_1_C + u_1_N + u_1_E\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": "iVBORw0KGgoAAAANSUhEUgAAAzsAAAApCAYAAAAS2TciAAAABHNCSVQICAgIfAhkiAAAEepJREFUeJztnXu0HVV9xz8JaUIwigYESnlcJCuhlhCKBJAGcgmtQHlYcVHx1SICroYobVdqMSoGXJa1iq2CSk1cpNGWankFbCnySNatxpTERCFg0GB5ykNAQCAQeXj7x3emZ985M+fMzDln7txzv5+1zrr37D175nf2a+Y3v9/+bTDGGGOMMcYYY0xfMQ1YDUyIvn8beAa4JjhmF+C/KparCvYGhoAtwJ3AqaMkx3hug5C6tEce+rnN6tIO/VrHrt/eU5c6bkWdZEz2hT8GfgrcCyyM0sZqX+gFYX2ljRtwfRlTK84DPhJ8PwY4meaB+zXg7VUJVRG/DcyJ/t8NeBjYaRTkGM9tEFKX9shDP7dZXdqhX+vY9dt76lLHrSgi4149liXsC5OQkrN3JM8WJCuMzb7QC8L6yho3MDr11eu+YsyYZD2wRyJtkOaBexJwWRUCdYlFwD0Fy2wG9umBLO3o1zZIUrRNRqs98jBW28zjoveMlX4+XuoXqqvjMrLFZMm4H/Du0hLlI+wLRwI3BHmfA86O/q9bXxgtkmNnkHRlZzTq60xkVWpiYsWCGFMXJqMB+3iOY3+IJsGxwhzgjuD754EbWxx/KPBb6O1alfRzGyQp0iZVtMcSYBj4csFyY7nNPC56T936eRrjpX6h2jpOypaXVjKeB1zXiVBtSPaFPRNy/Bz4nej/uvWFInwC+AHwHPAk8B/AgSXOU/ex82/Ax9MyrOyY8cquwLM5j32Shil7LJC86cwFNmYcuwvwDeDD6OG3HSuBpR3IFtLPbZAkb5sUbY8yHIHeVm4uUXYst5nHRe/pVT9fiesX6t2Hyyg7rWQ8Gng0Jb2bJPvChJRj4uvXrS9A/jYbBC5HyscC4FXgNmB6wevVfexsR4rzrGSGlR3Tr/wR8DLq+DF7o4lrf+AlYErOc01Bg6iOzALWIPnuAg5Db2zuQL/9ZXTTuAD99ruCslOAVcDFwLouy9Wu/qF/2iBJ2TbpZXvE7AxciR4unknJ75dx43HRe+raz92He1/HrWQDOA34NbBvUOZStB7mzTll/BBwc4dyFu0Lj0T5MXshhSuWd6zcg5IcB/wzcDdqrw+idviD4Jh+mZtuQ79vBN1Qdt4E/IJGZeThGuCvu3BtY7I4GC0ufCWR9jxwH3rQm4wWJLZjRnSuujET2IAmrwOB84GrganopvMaMC869nD0luXo6PsE9FZoDfAvPZCtXf1Df7RBkrJt0uv2iFmO5t81Gfn9MG48LnpPnfu5+3Bv67idbKA55i7gU9H3xcDpwPHojX87GScAx9J53RftCxuA30MP+VOBd9FwExwr96A8vB49/z8dpPXL3PQjFDSh61wCfD2R1s4/8CBUyTv3QiBj0NvrlYm0C4C1wfcVjHyzcRvqry8iX904ksh5NEJQxiTDVbYjK0RjGnlDX94KfDORtgJ4LPj+J2gcJuWcB/wG3Zziz+wcsq0kn9k8T/3H8pZpgzzhQrtFkfCaZdukbHsU4WxgE7pZgUK/JtfsVD1uirSbx0Vvx0Wd+/lK8rtYVd2Hx9vcnkc2gHegB+fzIzkPLSDjHsADOWRpR5m+cDKwFfgZ8NEgvdv3oG6EuF5JOdfDf0dKwQ5BWpm5KWvcQHfrq2hdvUiXPdd2Qv57RyXSb0YmyANRB16FFjSF/oGbgHO7KYwxAVuAv0qkXQd8Jfh+GLAsx7nWIAtmSDJ0aTtahWgMyRv6MjYvz02UXw7cFHy/EPheATmTLAFeCD6vIFN3mJYc/5Cv/qF8G+QNF9oN8obXrKJNlkbXaPUZTCk3C92UDgjShmhWdqocN0XazeMinW6Oizr187L1C9XP/eNpbs8rW8w6tD7k+IKyHUT6msKlFJv/etkXOhlrZUNcdzIuYi5BiumMRHpd56YydfUUDXfJrnAastC0e7s9DZldQ9PSZ4Dvd1MYYyKmogn2mET6wzTCSMacSev+uwt6g5YkLXRpOwZpf0PMG/ryFPQbJzOSTcgHOuYGOgv/OB1NivHn2uh8YdrURJki9Q/l2iBvuNBuMUj78JpVtMmuSGFp9Unbr+IM9CDwavAZRm9XX0X+1VWPmyLt5nHRTC/GxSD16Odl6hdGb+4fZHzM7XllAy2E34bmmEMKyvZ25B2UpMj81+u+0MlYKxviuuy4iPkHtOzkrYn0Os9NZerq58BAmNCpmWce6uTtomWk+QeuR28HWjXMWGUA1cnKLp3vY0ibfSk671926bz9yv7IPPuTIO1otNjwzsSxK2jdf38JXJ9IKxJ+sSh5Q18Oo98YLhacj24qYUSc2ZSLuhXzNDLnx5/nU9JeSpQpUv9QvA2KhAvtNVW3yVOoXlt9Xkwpd3103YODz0bgW9H/L1P9uCnSbh4XzVQ5Lqqu4zL1C/We+/uhD+eVbQ6yCixEbn6fKyjbU+i5MS097/zXy77Q6VgrG+K67LgAKQUfQApNci1NneemMnX1BuTJ8P90quwM0IhU0YovooFwe5D2KIr6sGeHMowG91AunnwZTkdRTLajerwQ1eMA3VWo+omnUN0cEX0/HJnZf8PIiDVlSYZfvDvjU6Zv5w19uQlFu/k8mqhORCE8YWTfnIhcAvYE3lhCnjJUXf+t6qwV3Wi3sdImz9L8O7ehG+XdqL7q3G4eF+3p1rhIw3UsioTeTdIPfTiPbPuitRTxeu4LUDSwdu5VIY/R+ZruXvaFTsda1SGuL0fW/feiOX+P6DMtyq/z3FS0riYhT7JtYWKnys6OtA8tdwnSEE+LBIiJtc+xaNlZhd5c7FfBtU4K/n4C+azennm0Ab09WAJcgUJJnosWVG4l+61HEZLhFw/M+OR5EZAkb+jLR1H44BPQG/q/QTeWF5F/a8wngVPR25Cib9fKUnX9t6qzVnSj3cZKm+Shzu3mcdGebo2LNFzHokjo3ST90IfbyTYd+A7wn4FMm5GLX9LNrRUvoLf5u3Ygay/7QqdjreoQ13+BLGWrkSIZfxZH+XWem4rW1VvRetSuciVwVYv8LP9AkOY4DOzebaEqYC6SPSt89gDds7qsIV3j7eY1THHuI1/4xZBBmv26VzPSfDsJmaLj0JdbaFgaDkaRDc3I+m9VZ8n6Lcsg6T75bpNilG03j4t8dDouBnE/b0fa3D+I5/Zu82mkVNWVTsbaeOwLVc39H0YB0lqyBD1Avyslb98o77ogbTFygUjjMrIVHYCzkMY2VnmY7EgoA7RWRA5HE+PjyFf+YRTdInSfWUp2tJFWeWcE5zgDLWC7D2nWz6GgEB9oI/NMFJrwCWTGHMyRH/KnwHeBX0XXvQtZpULNflr025NBKqYiTX2Y5o2hFkbpZ6bIXzXJcJXtSAvROAF4kGbrZt7Ql+OZPOFCs+q3KEXCa5rWdNJuHhft6aR+3c/zkSdstftw5+yKXPjrSqf3oPHWF6qa+1eQw/p6DXqY3Ccl791R3qeDtNnINS1parwcPVgvoOEbGPoHgvw8r2gnUI35EvrtaZapAbKVnQ+hqBfbkJnw75Fb3GvINBfX/SBSah6goeDEn0G0fmcY+cmGeQcH13oJ+diuRCbk5WgyHgY+myHz99CeAeuBLwBfRYsP2+XH/F103JPAPyE3xnhNwBAjd+ddi0InhgsR/5CG4rYyIWPcP/dl9MkbfrEVvwv8Y4Hj08Kgjlfy1H/R+i2D26QYvWg3t0GDXo0L13ED9+HqWEh913V7rBWjinEzB61zb8v9JCIYBFyMHjRPSqSvAxYl0lpZJEBa23M0FkONRY5Bv+mclLwB0h/WZyJrxs9oNmsuQArPqkT6EOXd2PZPSZuMzISvJGSIzzeMFJas62Xlg95oDQMPMTI05yRkahxG1sOYi6K0E4O0i5EyuIaRETgmokV0/5tx7dGgXfjFbpIVBnU8U2X9p+E2KUc3281t0Ey3x4XruBn34WqYgNYl1RWPtWL0etyclqfgdPTg+Z2M/Fuj/KSWfRwyL+3QVCKbRcAtLfLfhFzg0h7We8U1ZK/BSWMH9PCdtpHWAOmKyBdofrgPWYUe9ENLxxDdX7NzalT2z1LO9zjpJsB2+aANnrIUwJlImbsvSJsfHR9q8huQ1ejcKG9mlH5I9H15xrWNMcYYY4wZQRiNLXZF2phx7CFojUYyYsLNaBfuvQpc92VG+t0lWYKUrvAt/kJkedqOXLOKhDHMU/4i4FPkD3f4GrJWLChQJvaBns9I17P4sxtSomYmC5ZkH7T7bRx7PrbMXBvlpy1QvROFlsyiVX7ch9ak5G1FLnT70QiT+T/I1e7Y6PvO0TlWB+eI8xa0OLcxxhhjjDEt+VuygxO8JcpLs2J0m51QPO5QGXkPcrs6G/n0XYZCE6atLUojb/lNyKKQl1NQvbwvkT5AutXlXrJd/MLP/KDMEOUsO29BlqfXonNchtbpLI3KhG6F4fm+nnG+dvkg97xh4HUZ+bfTvObmFhTkYDfgnVF+rOA8QiPCzU3BccYYY4wxxhTiW2Qv/l5E67Ua3eQ0tOlR6Ne3HrlIhdxL/rjtect/huboYK3YEe1ge3UifYB0RWRjlP6GAtcYopyy82Wao7PFvJdsZSfrfO3yQcriMNnuhw9G+eEGaOdHaacjhWw7qleAf0U78U5Bymk3NrcyxhhjjDHjhDBW/AHI+vFg4pgpwEei/39YgUzzaDw0gxbUvw3t2BtyC3BkjvMVKb8ebdQ1lXwbKW1H7nYnkG+D1dsjWY4Cbsxx/lbEG7RmrZWaEf29NiVvfkpaN/gRckMbpDmQwAzk6ng/I3fSXR39PRa5+X2fRj2uBt6PNsR6XXBsO8ruEm6MMcYYY/qUDYxcEA56wLyShntVFQEDbmCkq9Se0bWPThx3AfDTHOcrUv4giv/O90VlTg7SBki3ghyA1ittJX1dzmSa1xINkf7wPg25df13hlxfTZELFFDiVXpj2TkyOuZ+4M1B+g7A9VHeJxNlJqJQ1k/QHK1tnyjtF9HfU1pc2xhjTGv2RveULWj9ZZ0jXRljTFcILTs3A3PRw/Mq9DB9LLAZeAytpbkveYIekGUhST7wT0hJa0We8rE1p8gmhDciBeZU2u96+xMUem8F8GNkFdqK9p7ZByk6TyKlqB0vIEvUUUgh3YqsPd9GbXY52tPnamTdeQQ4EDgeuAqtY+o269C+QR9He+tcg/YTOiG69lq0705IrLC9M/oeWm8eQhai/dFvy1LsjDHGtOdVtBHfnWj94yZ0H3pxNIUyxpiq2BHtVvsomvg2Ive1N6IH0qGK5LgSPYzHTEYTdDJ+9lfI9/BbpPzhSAFK2yi0FTehYACxS9kAra0gs6O8B1Fks6eRcrCMRtSxmCGylboZSMH6JWqj5BqdI1H0smfQ2qK1KC75IL2x7MScHl3reaS4/hhZdHbMOP6j0bl/RbNb3rIob32O6xpjzHhkEXBPiXKbyR/oxxhjTJdYjB78Q9bTvL/KVooFKMhT/ixkASnKOeiB/JgSZY0xxphO+BrwzYJlDkUK0mhuzGuMMeOS2chladcg7T3IVewsFDr6i8iNKy1yXBp5y38DuKKEzLtHMn+pRFljjDGmEzag7SPysgtat5MnyI8xxpgesA6Z5UMWAg8gt69NNAccOANZVwYyztmu/FTgOeCIUhIbY4wxvWcWck/ejsLxH4Zcz49D7tq/ZuSLvEvRVgtx0JgpwHeBD1YkrzHGmBSOQ25mWWGV07gQrQ2Z1O7ADBahcNTGGGNMHZmJ1jZeitZsnkhj/7LdkUvaRhr7yi1G0SzjCKMTkLvb0sokNsYYk8nHyO+mBvADOlszcw56Y2aMMcbUkVtpXpuzAkVMjXkH2jPvfOStcGiQNw8Fs7kj+MzulbDGGGOMMcYYk4e9kQVnbiJ9OYoIGrIORSE9vgK5jDGm1kwcbQGMMcYY05bfR4Fw7kykvw1ZaGIWAHPQ/f2JakQzxhhjjDHGmPKcjCw7rw/S5kdp8SbRc4BngT8HrqfZ4mOMMcYYY4wxtWNPFIFtGQo4EAYnmIXWuD6CNnAGOAitzzmqckmNMcYYY4wxpiDvBx4CngGGgIuAbcB0tEHossTxVwFrK5TPGGOMMcYYY4wxxhhjjDHGGGOMMcYYY4wxxhhjjDHGGGOMMcYYY4wxxhhjxh7/BzijwUvtrXoMAAAAAElFTkSuQmCC\n", "text/latex": [ "$${{u^{(2)}}_{(0,0)}} \\leftarrow \\frac{{{u^{(1)}}_{(-1,0)}} dt^{2} + {{u^{(1)}}_{(0,-1)}} dt^{2} - 4 {{u^{(1)}}_{(0,0)}} 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", " u_1_W⋅dt + u_1_S⋅dt - 4⋅u_1_C⋅dt + u_1_N⋅dt + u_1_E⋅dt + dx ⋅(-u\n", "u_2_C := ─────────────────────────────────────────────────────────────────────\n", " 2 \n", " dx \n", "\n", " \n", "_0_C + 2⋅u_1_C)\n", "───────────────\n", " \n", " " ] }, "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 (int ctr_0 = 1; ctr_0 < 59; ctr_0 += 1)\n",
"   {\n",
"      double * RESTRICT _data_u2_00 = _data_u2 + 70*ctr_0;\n",
"      double * RESTRICT const _data_u1_00 = _data_u1 + 70*ctr_0;\n",
"      double * RESTRICT const _data_u0_00 = _data_u0 + 70*ctr_0;\n",
"      double * RESTRICT const _data_u1_01 = _data_u1 + 70*ctr_0 + 70;\n",
"      double * RESTRICT const _data_u1_0m1 = _data_u1 + 70*ctr_0 - 70;\n",
"      for (int ctr_1 = 1; ctr_1 < 69; ctr_1 += 1)\n",
"      {\n",
"         _data_u2_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_u0_00[ctr_1] + 1.0*_data_u1_00[ctr_1];\n",
"      }\n",
"   }\n",
"}\n",
"
\n" ], "text/plain": [ "\n", "FUNC_PREFIX void kernel(double * RESTRICT const _data_u0, double * RESTRICT const _data_u1, double * RESTRICT _data_u2)\n", "{\n", " for (int ctr_0 = 1; ctr_0 < 59; ctr_0 += 1)\n", " {\n", " double * RESTRICT _data_u2_00 = _data_u2 + 70*ctr_0;\n", " double * RESTRICT const _data_u1_00 = _data_u1 + 70*ctr_0;\n", " double * RESTRICT const _data_u0_00 = _data_u0 + 70*ctr_0;\n", " double * RESTRICT const _data_u1_01 = _data_u1 + 70*ctr_0 + 70;\n", " double * RESTRICT const _data_u1_0m1 = _data_u1 + 70*ctr_0 - 70;\n", " for (int ctr_1 = 1; ctr_1 < 69; ctr_1 += 1)\n", " {\n", " _data_u2_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_u0_00[ctr_1] + 1.0*_data_u1_00[ctr_1];\n", " }\n", " }\n", "}" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "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, u1, u2}" ] }, "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": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "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": [ "__LLVM backend__\n", "\n", "The next cell demonstrates how to run the same simulation with the LLVM backend. The only difference is that in the create_kernel function the target is set to llvm." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "try:\n", " import llvmlite\n", "except ImportError:\n", " llvmlite=None\n", " print('No llvmlite installed')\n", "\n", "if llvmlite:\n", " kernel = ps.create_kernel(update_rule, target='llvm').compile()\n", " \n", " 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\n", " \n", " def run_LLVM(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]\n", " \n", " assert np.isfinite(np.max(u_arrays[2]))\n", " if shutil.which(\"ffmpeg\") is not None:\n", " ani = plt.surface_plot_animation(run_LLVM, zlim=(-1, 1))\n", " ps.jupyter.display_as_html_video(ani)\n", " else:\n", " print(\"No ffmpeg installed\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Runing on GPU__\n", "\n", "We can also run the same kernel on the GPU, by using the pycuda package." ] }, { "cell_type": "code", "execution_count": 13, "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 const _data_u1_10 = _data_u1 + ctr_1;\n",
"      double * RESTRICT const _data_u0_10 = _data_u0 + ctr_1;\n",
"      double * RESTRICT const _data_u1_11 = _data_u1 + ctr_1 + 1;\n",
"      double * RESTRICT const _data_u1_1m1 = _data_u1 + ctr_1 - 1;\n",
"      _data_u2_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_u0_10[70*ctr_0] + 1.0*_data_u1_10[70*ctr_0];\n",
"   } \n",
"}\n",
"
\n" ], "text/plain": [ "\n", "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 const _data_u1_10 = _data_u1 + ctr_1;\n", " double * RESTRICT const _data_u0_10 = _data_u0 + ctr_1;\n", " double * RESTRICT const _data_u1_11 = _data_u1 + ctr_1 + 1;\n", " double * RESTRICT const _data_u1_1m1 = _data_u1 + ctr_1 - 1;\n", " _data_u2_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_u0_10[70*ctr_0] + 1.0*_data_u1_10[70*ctr_0];\n", " } \n", "}" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "try:\n", " import pycuda\n", "except ImportError:\n", " pycuda=None\n", " print('No pycuda installed')\n", "\n", "\n", "res = None\n", "if pycuda:\n", " gpu_ast = ps.create_kernel(update_rule, 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": 14, "metadata": {}, "outputs": [], "source": [ "if pycuda:\n", " import pycuda.gpuarray as gpuarray\n", "\n", " def run_on_gpu(timesteps=1):\n", " # Transfer arrays to GPU\n", " gpuArrs = [gpuarray.to_gpu(a) for a 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", " gpuArr.get(cpuArr)\n", "assert np.isfinite(np.max(u_arrays[2])) " ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "if pycuda:\n", " run_on_gpu(400)" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "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.7.2" } }, "nbformat": 4, "nbformat_minor": 2 }