{ "cells": [ { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from pystencils.session import *\n", "import sympy as sp\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "sp.init_printing()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Code Generation for the Heat Equation\n", "\n", "The heat equation which is a simple partial differential equation describing the flow of heat through a homogenous medium. We can write it as\n", "$$\n", " \\frac{\\partial u}{\\partial t} = \n", " \\kappa \\left( \n", " \\frac{\\partial^2 u}{\\partial x^2} + \n", " \\frac{\\partial^2 u}{\\partial y^2}\n", " \\right)\n", "$$\n", "where $\\kappa$ is the medium's diffusion coefficient and $u(x, y, t)$ is the unknown temperature distribution at the coordinate $(x,y)$ at time $t$.\n", "\n", "To discretize this equation using pystencils, we first need to define all the fields and other symbols involved." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "u, u_tmp = ps.fields(\"u, u_tmp: [2D]\", layout='fzyx')\n", "kappa = sp.Symbol(\"kappa\")\n", "dx = sp.Symbol(\"dx\")\n", "dt = sp.Symbol(\"dt\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We define the PDE using the pystencils building blocks for transient and spatial derivatives. The definition is implicitly equalled to zero. We use ps.fd.transient for a first derivative by time and ps.fd.diff to express the second derivatives. ps.fd.diff takes a field and a list of spatial dimensions in which the field should be differentiated." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "-κ⋅(D(D(u[0,0])) + D(D(u[0,0]))) + Transient(u_C)", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUAAAAAZCAYAAABAZvchAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJNUlEQVR4Ae2c/3XUOBDHl7wUkLur4EIH/KjgQgdwVwGkA3j8Bf/lhQ6ACiB0EKiAHx1wHcDRQe770Xr8ZK/slbxSVuvsvKeVJY1GoxnNaCzbe+vq6mpRI7x8+fKJ+PpX+cca+dvztJdATRKQnbxWOq2Jp13g5aBGJqXIh+Lr7t751aidPU+VSuBc9nJZKW/VsnWrtghQSjyStD4pvzsktQbnRO3HSh9V/jaEG1Ofm17MmDcdp4TMS9DcJT1p/tw1HSl/VRPfNeulxgjwQsp7PaRACRPH91wJp/dG6Z7q6DMJctObxMQN61RC5iVo7ppaJAPs4VQ5QUQVUL1eiABrSS9evDhW+m+IH7UdKZ3321X3VOmiX7+unJveuvH27VeLEjIvQXNXdSVZTLKFEvPdBb3UFgE+07bFLjYERH5PAjscfR6q/s5Qx4H63PQGhtlXexIoIfMSND2Wd+rSbKGGKLB6vdTmADnDGDvI/a7293J0v/wl2ZSp4/Y4BXLTSxl7J3El6xMl9DQVSsi8BM2p89tqv8YWOB76e1NGKtX1ptPq9D/slLZYkLAtevsyxIZwxqLDn+p3f6hvqD43vdAYM6wjspgcXZSQeQmaO643bOiR0pi9xEyxOl3HMJ2Ck80BahHy9IkobCoQvfHe3ygNtRN98ITY8D6r7oPKx/RXSoLc9JIGv6HIJWReguYOq4eIeOMIMMf8a9dLNgcoYfHqyjsl9wheObuHOwNQfmb1I0IlejOnFkQTDW6PcZLtC5+69m/JkhxgTnqixXxxzreVvqrc2X1V5kn1Y+WjcxTOrEHzz6pDhJWLpujMRYfcAhOQbBqUbLQWc+iltE6ynAE2THILi+AXKnPNDvRaiXA85msOFh+3sUEQTXc2qLx1fiCqDG2r+xrsHKjMTU9DPBdNnD98nvtDqh7HyEOaoPNTPW1EsEVAtJ8WIZxIVHxk1SHDZ6Y5Fx2aHRVbU+tUn1Evk3Wyjkfac0WALtzWpHkpma84fim3COhBDCPC+V0pGMGJFg6EW+Shl6NxnsB7foRvDuiHin+ozNPlFgrQw+F/bgZgvrYAbUzq3OZgFZY3vBDVtnPX9Sj/1tfPx/qo7ZXSVj+V0vipOkSnrCveawvqPSdN0ZqTDm2jxaauHVL1AoPqQ/DySLlvB5N1YpMWPTZ/7i6xSZPLma65Oz07FAIL7ZMSeSzAqG/QGDhOjx3eor5YWobH+H3HYW04MOj7Y1obOTvdB7WDg6H9UG634ie6vlTyHXFuejgw4w2jRcA+4Lz7dQv1Yc4PlHNg7UDXMfwbekqfC9Hmc6nOZtAh5BWEhx7huw/OqNR+2m9Q+Zvq27n02lNkzsK3scfWZU6ac9LhkB31VLIsblPXcKDx0fEd5a3zW3LmjruS7aqhyfphDT8TXecLGpqMh5/iA4pnOMBfKgR3WOsQkTMY0RfnXDgjzgOZTMqZF3ys7Fiig3Ccg1O+Amo3Q2GyANHTX+5KP2onKsUBHiuxyLPSYxzRhHdyeIG+Rb/UwTt1oWMAeDW+delglH9D6uVr+4gP5IADPFJy/PZodIrCCTm4heqJ8JFlZ1F1OvcKwmX+0ToUPoseZ8pYQchNU/ScTJTPQYdmR1GOUHPemq4b5SJzc3StvqfqRP1YN2+V/jQaLdHlBTbnnO1BryG5qAHMwIkwMLI3SjhUJkU040B1GMEUMGXaLWafBlEAxsLYxkt/J2FxE1UAuektqS5/iX5wss6YmgZnxKpbUbDaT+C7wVvoOoZ/Q3d5Yh/GQi/XDdEyT2CsBE2Gn4MOzdb8dZgg2o1Qo/XC2lVi8yYtuFa6Exg9Wifqjw2te+CIf3Dn0QeBwVKrnEFp4NaQ+wTUhiM0wfSbrQxTpjirW6gv9UFFqg3nwvh227XSvyHETujGz02voW8Zwu87Om69nWw09hMlx6NycPvzWsu/DeTlKX3YRPyjAI9MuctEmUcxUoJmM/AcdGi2xrq/VkjRC7hK7hhDTPKmCLerffuB/2idCBdnCl1ejQsCY1j7QRAjrbI18EA3i8Rua0C7DqC5Kozz3kAjQuqE6aLHTkGYy99mraPNgvAdRW56Iu+gw0fDI3MypSIHc3ooNXaB9vlfjjb+G+rD2Iy7DdhU5iGeS9Ccgw7RMefhttZCsitZl6oXbHkwgFJbik4IiMZodeZ92ClNK+BY+udYUHqs9I+UcH+gHRwfcBKcT62cUamO2+qfSoS2vtPo3+MPKRweWyHmpudNAsW/FX12IZ5AMybHAYT2PI16p2QAT31++2Uft+XfKps8pY8vux6ZssUEmUczUoKmBp+DDllz0U4gWuCRiCl6ES4Oa52zTtEJXPIi+CBoTF47cxHixg5QhIK3VM0Ag2Fonzvhc4aHMRMxrSgvhp5w7PyNHRCH6kOnnJseA4km/NvtuD92qA6HduQjpfBv/RL7IJchR2oki+XilfUQvSZiGMlNU/TmoENsKBSUxIg0C06CXtroT32wB84F+7aaohN8B0FXEESbQKRdgwdBrO1V8iQ56FATWDoTLruKA02Ya16RmWr4uekZXygZh9SH0fE0DyJkd2jsdRzt4+Ex3uju6OEOXbIYSTcepIfqdMj6kGJwKtjSpnAduubLKY6/gOdKU+3UEdAPd548XGx9gDWoDufH+V87RlX/CC3GUByv0PxmTE/J1d8cBLehCJjD1clGm5uezUl0eRLFi76tQmgbG09tTkZC4xWjdicb6wNNQDjB8Zat9fyKTxz1qRKLmPnyyk37bqeuk6EETZgYkqnqB9eg2orpULR54Nh5tzRZWNfYoZEFuuZF6C8qd6K/KayIBpuAyR+7xw8AHKV1/UCJP0LchKb+RPFC6eEmNHalr+Z5orTyB68x/KfKSPj8mexlDO09TvyfBEumVelQ/HxXOt7rME6HB0vHWNUvISyh8OxBuxHnFdzSEvFEg/CJjFbOSdcQYEdkp91DRgnUpEPxQvTHJ4+dO4qM050dqeocoJRHiMo7QRbCzk7o/oQ0T5wSt8GE7bHAQXE3lB/pKVwMgxfV94YxIqepTZLr1nUoHthE+UQ1+gudqfOdU7/qHCDClRI52/qufOUgc07Ct7lonjzmjwbht5/aRXbiX7RTI8ZI0ns0JFCBDllDobcN9goakcD/XTRMxcq6e8AAAAAASUVORK5CYII=\n", "text/latex": "$\\displaystyle - \\kappa \\left({\\partial_{0} {\\partial_{0} {{u}_{(0,0)}}}} + {\\partial_{1} {\\partial_{1} {{u}_{(0,0)}}}}\\right) + \\partial_t u_{C}$" }, "metadata": {}, "execution_count": 5 } ], "source": [ "heat_pde = ps.fd.transient(u) - kappa * ( ps.fd.diff( u, 0, 0 ) + ps.fd.diff( u, 1, 1 ) )\n", "heat_pde" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, the PDE will be discretized. We use the Discretization2ndOrder class to apply finite differences discretization to the spatial components, and explicit euler discretization to the transient components. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": " ⎛-2⋅u_C + u_E + u_W -2⋅u_C + u_N + u_S⎞\nu_C + dt⋅κ⋅⎜────────────────── + ──────────────────⎟\n ⎜ 2 2 ⎟\n ⎝ dx dx ⎠", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnkAAAAyCAYAAAAgCc0nAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAUBElEQVR4Ae2d67XUNhuFB9YpgEAHoQOSVBDoIJcKAh0ki3/5x/rSAUkFEDpIUgEJHYQOODkd8O3HxxKyxx5fxvZY9tZaHsm6v/uV5T26+c7Hjx8PNttE4Oeff/5Dkv0k+902JbRURsAIGAEjYAT2hYDe6fck8e+6vpX75pT0d0zyTsGTb5gUTwP4Q/avc0ihfD9Xvj+VeX8p+5p7+ZtQlqDYMgJ7RMB9wx61bpmXRkDP2SOVyXv+C7lbiZ5J3tKaWaA8KRzFv5cdSNikpSpfCN5L2U9CxnL/T+4fdT2R+8/gb9sIGIH9IOC+YT+6tqSXR0DP2zeqxf9kP2yrzd22APvniYCU/VQ1/1z2LASvRAVC9yxFqCyPfxMQzN0Zyf+4xH53svcR2Pj0QWkTcdw3JGp0u0/AaHEaoxZgengLuzeK9kZ263vXJK8HkLlEkaIZYaOT/XbmOj9W/v+qPNYFpIYRvHtlPVL/PbjBoo7HHuTuK6Px6YtU3vHcN1T153ZfxaPpzhg1odLTT+9bBnQeyWaA58iY5B1BkrUHGy1eSNnvZ5YCMsd08E1LOTy0NkbACOwPAfcN+9O5Jb48AhA9pm2P3r1Xl6+bazAFAiWLvy/7lynyO5WHymgbKWQh6EHhxeYL2TQ4/l2wXuAf3Vc2geieIeYfZLeRRQVv2xijbv0ao26M1hJDunLfMJEy3O67gTRGtxgJB6Zsn+vuaLmER/K621EuMVBuhUQtWXE1MAge08XpWsDn8od0MsJI/aKRP+TvG9mNBE/+hJHfbEb5s1Hk0ma1GK0EH/SzWozGNJ4V4Tqm+oPTSN6s+oYV6We17d4YDX4MeiU4E1fevU+VR+W96ZG8XtCvO1LZMBg1e3HBmjIqx7+JYiRRNh3727I+7MLliJXU4Nd43IrSQgCZDo7TznIHkvhBYQ90n5LJNN/o7kqj8F90vdRV2UQSM5jZoXKnxog28J2uZ8r7iz7VV7xWXBV2UXyov+owKUZ9MKnHUR02h2tdxpnvs+obttjuJVPrc96m+1Pt3hi1odbtPxeuyvdPXbxT0XUcUTfJ69ZJDjEgKRCsxlGxugBlI/tL/ry8+hoOXWwjZS+VCaQsNqzyPsSHeNQJKAu0636Hsm4cwxLzkhvS90F2IJDsZOUMwHiEi8IrZkCa3xWXtQydpJECFA9ZqXvd3MdD4U2E8Z38ozxJQjCbCiPIUKhXL72q7D64XhIfoJoMowT33k5htApce1f4zIiSl7bjvkGnBAiLS/ULaHGydt/zOa+0nJ7tftcYVQDreTMHrrWieaeiFzZAFnzA5+TVEMrttmw0/6jeHIgYCMNiYqhMiEKFlKWFKxziwXTtZ3IXjU42w8n/6jqqs8IgUTTSeNae3P/J72vZUT65+VTLQ9lxtE/30QxJo7jgR/5F/WImAxxKy3lFHF0zeE2k0pyNUahqWY+T5yYlcXvhqjwvig/1nRKjIP8Qe624DpFh6bjCLOu+QfXfRLuXHL2e86b20dXujVETat1+c+GqfPmThr75MEHxLrrbXR3HWDkCLLa8kUIjAVqqvmVDhWjFUSq5ITrpmgDC+FeaEigI0UF+TXV+LP+U4JEXDbdO5siPUZYjo/RD01BeGAU7ym8Bj7MwGlO/gRhdGh9EXByjneA6RszONGpfPOO59w3Zt/uBz3mnXhsiGKMGUCbwGoWr9M17kbRxRskkbwJtXDgLyAlKXdSoMUGwvpJdn+akc79OKgPhqpM5plmLOiv9U12QuINs4tJIU1OEpR6lmzLuN/jjNTQNawdbp35bypjS+1yMxtRlCEaXxgf5LoHRHnAdI+PJNHqOt9I3bKHdD3nOT+q1JdAYtQBzpvc5uP6usuNgy9WZFXHyCyJQdqY8xEyHLmZKMkZDYqEn06upYSQunbJkBC6SsbLOX8ov7ATm3/5NmQEv8pQglt6NFnm2dWCNCeTZlobyKftSZi6MxsjThNGl8UGONWG0JVzHyNKaRs8yz9FW+oYtt/um57xVrycCjNEJcM4IOgfX1yo3rB3/1STvDC2sIGmYYlx6JA9SSWfOmpu6qY/aMdL3mzp/dvx80MXLmp2frBvjCJNXuoK5JweNOzX1+xBGXPJqMkPT9CWWTWVN4XcuRmPqMASjS+ODfGdjpPb2SPmEXYanMGN3clvbOpWOsNxw7ZJnaPiW+oYttPsh7XGorom/W4xm7k9G46p6sXwLvTM7ZZJHK83YoEQUOvaFNEp0ldf6MeR6hmVji2v2kvAmP+S4l8Q5IFuZB6TyXRrWcF8Ej0hD3otimMpRyteER5PfEUZpXn3dAzG6KD7INAVGyoP2M+u0fG649m0vfeNJ/i31Ddm3+4Htsa+a03i7xWjm/uRcXP+WkvhTe7ibasvu4QhI0Xwzjq85/KerPnU5PMNhKZj2RJmbMMKPlzCNu27YFh5GLQ+Kh5sjYwpiJpvt4vURmpNpagVQJrt9zzH8c2r713xOvpW0krMNo0q8+s2ZGGWDD3KPxaiOWZ/7FeDap5rZxzmh05PP+Qr0s0i/gIInxqhvm9lK3zCmHfXFaEy8c3GN7wkfoTIG/oY0esA40oOz5N6EYLkZlWIbPkeMTDpSVObNVulf5X4WyszdlixM9xxNl8k/kDimfBktYIt4Qapk84+Fs734RFqKf2saxY1GaRrLjBFW5miqr/zoFGgHEGDwYF1kerbgaIyaylPeqzZT1dm4rkfNbTqVf+tzrjC3e6lwKEaKf7I/Ca1C8bLqO6l3W53lP7gdBRz62ipjEVxVDkupGHR6ctW3cluIJ8F5AbLrJCz6n0SsMl/yqq+NK8qT//UkBVUzobFgzh2Bus1lPb88aJCVyq5dYVy5T6urMP61cA4fO3ujOZUmRFIciPhB9qQkPOQ/k32EUVn/yTHKFB9gP8JojC5ywlV1naV/G4PbTGkadSq53e4/AT4JRl3tnuIUJ8e+k6pPghEZDTUL4hpm+B7dHVrJzOPTKIuGObEcrPMJa8fSrNv80zhj3YHk5UROOmXVQwBRZvo1yNeZhgiKzwuuTrL7pA0PfJ+4q4izMEbZ4YOSxmI0RsEj294cuM7Vv42BZfI0Y3W6Iv1Mjkk9Q2NUR+T4fmGMjivQ7TNF3xB4wcO9kbxueMfFaCMY+MfpwzRrNbRzyeb9Mr+gzDT7rN3ChpE8pmyHYMQI7c0QwRWfIW2+rpEdhqrz7BjljA/tYCRGQ5pQiDuo7eWOaxD6EvZIne5KP8aou2UugVF3LY5jTNU3KJ/wLrx/RTHy4GXKC4+1TmwiqExn6p5zj1jvFBLqdp9GGDDCBNOGGLA+jNEj1n280HVQOMSOFzCYEpeNGeD3VnZ6ftxfun8V/GQT/7ku9PAi+MvdZoiP2aROJD9r7oKMt5Ke+FXcSps9ETUNeq102eK3AEZZ44Oih2KUNo6+bpUxtO1lj2tfbOaIN1Sne9SPMepueQtg1F2J4xhT9w33io0XErb4ELNs1jRxptlnoWy5iwV8su8Ev1ztUj7+1aVkq7c4SgeBg7DxzdVi9Ec2C0/xj99mlfsgf7BklOgIN/lBXtg0wYYMDhSGJLJTFsLIYkkIDuvMWo3CIZqcM1cptzWBA4yAEdg0AuoT6HNG92+bBsfCGYGdIaD+AI7x/koOCMbbUn7WkNU3CeDXSjjKjuWd7FmmvJTvj7pGkbJSpkks1QFiBsGDgKWy4kb++ojQKdy+o1JKA8GjY+asuzAaQLo+hvrYGAEjYASMgBEwAkagjgBc7h7TtWwYCCQO8lFMOyaxGaWq+xXBSscoH+kj6ZGbESYMU5kPdN+686mIpZ9TaRT2i66XupgC7WWIr4jUu26KdWwteUHUmg6fDXkgFxsCAhkL/pTTtO6uzZ90EDmIHaOA1DXshJGztwmy1MnlyQxUJke92BgBI7BSBPSMHo3+p1VV+Bz928F9Q4qy3UZgfQh09Q0NNb7PSF5BEmRDShgdiiRGfqwpw+9o16LC8Ge6MRIjuSF96dlcj+X3h67W0ameaZj2LKaUlX+nUdxGQij/c6YzIMAVHJQfGIARZC2aNv8Y4ZaAvtY9I4OkZ30eRHn2dY8q5+QLRHWwMQJGYMUI6Bmeo387uG9YsdJdNSMwEoG7STrIWv0YEEgRD38Y6UuiF5sP+EeZGka7IhFSOtwQPYhMm+lMk+QDqVrcqHzK5aqQOd3HaddapYpRxLLeRVCZx0E2WJAXxJXpWg4z5luupIEkF0Z+XbIyFHvoEe82Q/8aASNgBIyAETACe0LgOiV5kI86mWMEriBtIhNPa4QC8pYSukBe4tRtieSNbNb9HRmlH5KmIIxHmSzrUZcNfArMJMs3pTzUKPpzQ5gsZMUcEcBb70+/ig/ZK6ZjP/keucAV0xXvNpZ/jYARMAJGwAgYgb0gADe4uUqkhcBEwiCi8Uj3X+oK07cP5VcQC9kQlkAy5CzMveCo2Yw4xXxrYUPSsDkE8tS0/q2W7bS3yK0LkhmI2kH3EDcwCkSXqetQN+QtCKH8kJEdbyEsEmf5100gkWAdcK/HCfesecS0YXgb6t9dI1C2v+clCKH9zr4sYNegW3gjkAEC7hsyUNL5VXyfkjw2SHB8CtOnEAgIB9OIrIX7UfYrXcHwsiimC4PHCRvCM5SINKW5UT7hJXWiuNmCmM4GH7B4oAs8IGwBn3TqGizZLFJMv8pOdweDRRpXt4X5Qb/fK+5XspvCb2N9+gUPGyPQhQDtM67hkpu29Y8uzsS0MQJGYL8IuG/Ytu7hGteR5KnzhzTETRSJ7E1+JK6TjPp9yIK4YYQq+AV7SJq+pDLkPak9BB/FRV4I4JFRWJs/I31htO8oXYNHwANCbLNBBNRWHkms33Tx54ZDMiNZGyAuyyyK9Z9lGv7E4cch3e8G5OOoRsAIrAQBnl9VxX3DSvSxtmqofcC7MDd3b+3Bv5CYkEmRWJniB2lrGm1rfJkMTEO+bWSxqEOPH+rHtQUTsGjCewvy7V4GPR8c68NoOs9afdNPX3wghn/3jex4WSOwpf4ta0XMXXn3DXMjnH3+gRd8iCN5Q0QqG1jIJE36QjePdRWkTvFwv5FdEBLZvKyey07PzjuZRvGDobx/w80YW+WG9XNjkq8qjWSBAFAnSIDNRhGQjnmGMKPartLX13ZC+tKzMYvM/ZM/AtL1qDaSv+T7lMB9wz713lPqwM/e3+2ZoCkaL4qQURGue9aecQByWKfGVC9rzYIhPlNFxdEsePZIE9IyzTlkOjOk27INeWZzjM12EaDd1482GiWtnjWmeHj2/MdgFIJOZARWhYD7hlWpY1WVCdzs3dUZ1WJtD6MC6agcpK1yn+avMEb4+NZqJHmEn0pThjMCSLwwRcmtze2IaRjpMR7bRAD9nj1Co2eHh55nlu8uM61nYwSMQN4IuG/IW39z1p4NnAVnGk3y9KLgIN9vdXE8SG/ypbhjGmYglHOCkmPerNPifD4+t+YXd44aTOosHQYixvPEDnfIHaNvLGmIBn3rJhyczcPMaDlpv9eFeas4cdS7zJdvLvPP/1DeY/d+bklnYwSMwGUQKJ9Z3oPuGy6jgtxKjTzr7jk1V8NjJO+Z7GKkrWdekMLehERxeZmxO9AvpGOAw4J6FGqTMQJq3+iQo00gY1wsfaBTx9RH8lgOwTedicP5keyye0w62TyLId1BfpA/jk3hSB921EIaiXety8YIGIGVI6Bn1n3DynW0puqVfT7vgWKz3lkkD8GUYev0bJPgil9fCN4ULfXj6Ij6Sy4N361buDD9DWEOIzi7xSJnwaVHHki+Ywy5S//M4GaDDToujNz86YkkTm7CWP4Qnqv7tXCIY3hJ4OZ6muapexsjYARWiICeU/cNK9TLyqsUBn2K2ZyrKSo75wtjzrynkH0FebxWHYpv6K6gLq7COAQgbUy5B6IWcuFhjdOupeffipcSQTZRRCKosMq5lrr/LGRm2wgYgewQcN+QncouXuGwIad4T5w9kndxcVwBRoAgCIG9G5H8EICkV0arpU/+wTPVWjkfT/7F8USJiKR9ldzbaQSMwHYQcN+wHV0uJQkzO/GrWSZ5S8E+Uzl66UMOYOyVEZyZinO2EyNQkjkIXYXM6b4YnS3121iqwiD2pI2jfeRX5tmYxp5GwAjkgUD5HLtvyENdq6il2kzYkBdnhSaZrl2FdPuuBEP6LKxnTdfNvqHIVvp0ChYhGHIvRu2kU/6Z4b7WxcgtGy8g92x8upE7TVs/bFxRbIyAEcgYgfT5Rgz3DRkrc+aqs0eCD1BEHuCRvJkRXyJ7KRTWjlKfL1Gey5gOgfJhhLAxNVsY+UHq2AUbdk8/kR8dPSN3XNe65x8+pC8a+RHGblsbI2AEMkdAzzN9uvuGzPW4VPXL/p/3SGUz7J2PHz8uVQeXMyMCUjDDtMzDc9h0ZPEzFumsJ0JA+oKwcQwKBO2BLtbYoUNGaPHjTEo2VxAPv+LzfrrnGBUIIf/s2TUL+YtTt7q3MQJGIGMEymfefUPGOlyq6morLPlhZqeydMskbykNLFCOlMuLHkJQYfILFO0ijIARMAJGwAgYgQsgoHc+sziQvKNBHk/XXkAhMxb5g/J+KoXHqb8Zy3LWRsAIGAEjYASMwOURYBaPL5AdzeKZ5F1eOZPVQApmcT6fwIrbpyfL3BkZASNgBIyAETACq0JA732W8LCcp3GpjkneqtR1fmWkaD51xbw8ircxAkbACBgBI2AENoiA3vPFZjzZlXV4qagmeSkaG3GXCudbpizKtzECRsAIGAEjYAQ2hIDe75zAwGDO16fEMsk7hU7GYWoAfO7qWdkQMpbEVTcCRsAIGAEjYAQCAnqvhxMZvpb7aB1eiIf9fyRZszVOHmoKAAAAAElFTkSuQmCC\n", "text/latex": "$\\displaystyle {{u}_{(0,0)}} + dt \\kappa \\left(\\frac{- 2 {{u}_{(0,0)}} + {{u}_{(1,0)}} + {{u}_{(-1,0)}}}{dx^{2}} + \\frac{- 2 {{u}_{(0,0)}} + {{u}_{(0,1)}} + {{u}_{(0,-1)}}}{dx^{2}}\\right)$" }, "metadata": {}, "execution_count": 6 } ], "source": [ "discretize = ps.fd.Discretization2ndOrder(dx=dx, dt=dt)\n", "heat_pde_discretized = discretize(heat_pde)\n", "heat_pde_discretized" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It occurs to us that the right-hand summand can be simplified." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": " 2 \nu_C⋅dx + dt⋅κ⋅(-4⋅u_C + u_E + u_N + u_S + u_W)\n───────────────────────────────────────────────\n 2 \n dx ", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiUAAAAvCAYAAAA8abqkAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAQJElEQVR4Ae2d6ZHVOBeGL1QH0AwZMBnAEMHXZMAMEQAZDMUv5h8FGQARsGQAEwFLBpABTGfQ3/uoLbft60WWr+/mV1VqW/vRo3OtY1l2X7u4uFjZdRP4559/TpX6tMhxqzg+VPx5d6nlpYjHI/X6h46fltd799gETMAETGCIgOaHV/KP+/Jd70t0WiDwQhCfFP5PxfyS/2o2VwTE5r5Cd3S0QXKFxWcmYAImYAJ1AsynH+tR9ZCNkjqPttAjQTyrJLzQ+S3F3a7ELfZUHMJKko691u9iAbnjJmACJmACgYDmiR86ea/j311IbJR0kbmKZ7L9chX0WYPAe4VfNeIcNAETMAETMIE1AjJIXivysY7c0K65k7WYI4pQp1nNeCPPXpB3Co++my8AVqlQB3snvlUjd3W+iT7myq624fqHjvdy60gtV/QTRR49hqltOJ8JmIAJmMBWCHAjy9zMloiaO+qVEk1g3+TvqMdYZL3PsWpUOgLFxBj2T3Rk2Xp0Xx+Vdir/XR7jYQ73RJVi9W7DsSLz2zYachsmYAImYAKzEmDeuM8c1WzlqI0SOqtOx/0gkzZhFhM7+0nY0HneBJkSRhZ53lLZqKPeosJmH4nHIGFz7hyOvkw29oYEU/86nz8OlT2W9Ll051j40A8z6h9N8+nnk5tqruPJiRlzKE8b/mqWPurHN0VnebTA45YsQ4I6VJaJnTdwwmOKIkw8m3bGOKxC/KZdVx+74ie3r77Hjb6z7rcp2mHsssdvcmf3o4K5dGc/ercZKcyon6P59PPJTTXXPHLMHTy+qa22H/1KiTrMakFzBUFRaU6TIgYJz794v/p2MUny2GKu1Yc0weq5uvpI/Id61suQ+sEPaYqj7knGXmLjDyRrTWkTyzmbCZiACZjA/hL4LtH+aIp3QkQxQbEU/7v81+YkoDDP8/f+g2GSEwOCRyysYPyUxxjhjv65fOkq/SXurvxDeco+kMd9Vp44mfNNEibw2rdJlL6TDZdDfVQ6xgKyITN9wpBi/OjTSx2j+1fhtzFOR/LzkTj04HmMj5lbjnA7b4nfWJRk4LFN55s9hcwHr7cbA9ZRkTl1gCmizaefD6lmNMwoJ8fCufL4hn2P+HIuuV6AfKpIJiz2BzCpl07xXPTZkFIWKhN1onjSmPxmcao7aT+B8jEZYzjED53Rn9iX5koJH3B5iVeez/LsAj5TmBUQJudYbqW4G/LXml55tu4kw2AfleeTPEtiYTLX+T3C8qVBonP6iLGGUqwU5pxne5RhSa3JS1FrjjpmWy2STOjUuY59j8j2Vm+hJdmTdJe8M7u95bQnjMxnWAH3ltGw6Os59kTvEGzJXOP8UbMfTjQ4TEhMzDj2IMSMIaKIC5NXjIhHlcVgYQm/nDh0Hid0VipuKsxE3+v6yigN4+GVfOfKhNKYIFkNwCApZVGYc97AOdcxOJ0jc5SRONJ4o4bVEtxv8tX0ELnrP2P6WMjKWLaOm+IxQFaqEwOGvjP5x0cklEtxcKqyTikzJg+v/3bqjtJ2rrd0RnJEXVnTd6UN6u4YIDl5Z+DEbw39YXx4s23Q7TOjGfh06kMXKMnQyVRpx6hDoxl1seuLXxrXPhZdaX2Mesp0jt9IfY3zMnNJ6Xh8U/3mBheb2qMOhbk7b8atis6Eu3ClB6c4JvyfCEaEjmfyH+U7JzqlpZThC3Dhc++hofU/QGIJKE6sMQeyx8cwMe6L8lUnUy6speGiNFYZJjnVwYoDbTddgK/0NgMLGfraHtNH2m3re5SH8cAQYWUsro7EtNQjF9KmARvKql7S/pXnmOpYzQlGlI7oRFjp6Sm8U71FrkLOIX0f0t1aF1XnpnVnk5wwBKNeJ43tATDaJJ+Ua1lzvFOYHpMOjWZUA5YYkN4timsillq2REbNMinjl6qvrfPHiQQ7p1UdudhwoSkndsWxrEJc23I+k2Rz4iDuf/LBqTx34hglfJa9agjELBwHyxT1YJRgeAR5qxXoHGOqJiN5FYf8tVdWFd9cPWgzxFQs36mNNqNjpXhWJWBRPkoZ0cqYPrb2vdIWY/1OntUlGLG/hPEZs2+IcQhGlo41p7pIS7qLrhVUQGWDzhXyNJPLcNEG+Xelt8iyCd0t+8SJ+rNR3VF9jAX1TuakOvjtYDyjx6lurxmpLxvjIyCDfW1CS2GqPFxH+65/tWqVd291KIdRrXOJgQVyTSRzlS2F0VXu8mxQx1Vvqr7G+aNmnFwvm7p8NYe7hvAjLeLDxacQvpI1nJ7ReIzUeZhMFG4aH9SH1brmRpahLS6sNac6ThWBrxkfCjOJr5Reyki46pQWL9Tlagr14av5dn1eyDOmj4GTypV9j33SMY4T1izK81oeA4IyWMHBKW5XDJDvrtpHvtIrLtz5FHH8MKJjdWmrekvDkiNyTNH3Vt2NHdjScRKnHBkPjNEkPiP7moPz4HVoC4wWyTWn0zllRo5fir7GOea8Ks9JJcBFtrmKwDJ/mNgkEBPWOx1Z9idvrSKFYwM6rTmsoGgR1RIUGFOGfS/IUxoQjcqakwN5Q38kL8YV58jC6gB3HfSLOwr6Uy3LxqPOvQzKv0tXlRM51vpY9KWMJ5Pi6D9lYbBmsCiu5pSfsYZPc4yr+aiva/yq+Uadq23aDTpXLaj4/4jXsfmIaxd6i2hdfUfHmvo+pLvUN7ebyilHvkNiNJXPmL7msDwGHZqb0VK55vQ7p8yY8UvR13id5JpZuuvlWWOFQxd/7kz/kGciw/2uuPPL0xU/4FpFRXzbgYa7OtOWn7i2MrRNuzVXyMQkVqYpjkkY+XmTBMfeFyZRJmT8L4WRqdYHxZEGzL1ykou+p/YR2eEXDJiinzwyiuOIwbI26VNILpTRkbGO5yGh5Q+c0I9tOcYL33Q1OSU3474rvUW2ZN1tdmTm8FyccsTeR0Zz8Wnraw4zrgHlNS6ngg2U2XdGOV08Zq45PHLKtOl4Clf0mUUB8pbupDzTmys6f6MMLI3zJgEKyLI+qwp/6/hWPrpTndQqaglX89aUOSaMLFMzICp1cMrdM7Ij5015ZGXyjbLHvS9MxuyZwfhYKT9vELDnhXReJ8ZY6VqJocguXWofkZGx5I2l8DhGx+oeFsYu8iBvdA91wofK7urYlh7zxSNGTnjUpTJNXYh5Jh9VN7LEi/GZwu8V/qxj7NMu9JZ+dfUZvk1979Nd6tqGm8opR8ZDYjSVz5i+5rA8Bh3KYqTf+m0Bqz6y7eLH9bz52+vKG+OPgetqZkaR1ZjxS+GKfbF2g1waJeoUDTaXxhGmLY6BPyUxOpShqIMJ5FuML47NcIgeWYZ6WxUuVfYiX20TmOIwQvbVECnwBQMqeXzgqoIYZWtOaV3xozioHh6lIBOrEmuKtdZwZoTaqI1Xs5pChjYdbYvbiN4ig9odo++dutvsz1zhqZxy5DokRlP5jOxrDs6D16FcRirH/NF63coB2Shz8Fzpz8yMArKR45fClblj7Qb4emht5J8CAI023XNFhFUIEpSP8w90pghzZ920eHvLUK5wtMdnaQ/ZMYnjj8XxBs9cF4uNM5LucXEbrbcIorK71t2t6U4Pp94xWQqjHj6917IOPr1MK4mbuP7tgw7NyaiCK/n0WLgmdzg1Y4e+9o5fpe5ertStvKyAMYfU3LWLi4taRGpAlfK2y9pymeKj0cEjID5bzwfNwkSsI0Lw/QpePS1XJ/rKKG9wytPaXkz3cfsE4njqeGP7ree12KVHiu/UW1pS+qJ0t42T4rjQPJbnZgMePEKrfqdlMYza+IjFSvGdeqS0NT6K62VKnTjlO7jrX5fMih/F6JLAuL9qY3FcxxEKOtXLSAzX9JU2+sYvyqA8vfqqdLYW1L5zFsuuMEpy/LNnz87kX2SWvT+mnNo5lf84pozz5o3rWG4al/fyo8ZzbBubzC9Zs/UWOcb2VfkPUnencFoCI/MZvr5sk1HONULyLe63mclp1PU9havyfJe/1SbP9dI6GXkiS4d9BCxpY20lO+XnLmvsHgQsa+7Q7PaPwEOJ9HT/xGqXKFdvqW1JupvLaSmMzKf991WN3TKjatOp5wc5r+RyTYVSzTfH71l1skrCixite0SzjRIEV6UYCjzC4flQquP11PA4J6WA8tIBPqLV2oGUOpxnPgLFWD7XMS7JztfYhmqWrDl6S+uL0t1MTothZD7DP8htMBqWYj2H5DroeSWT6zqI4ZiN/p4lN4sYtX8QuyZC2/LJ2DiWa8aWSc0/Z92pMjhf0lLtI43V2SGxmlu35q5/W6zn7MecdZvP8O/2GBjl9OEY9I5+71s/huRR+quhPNkbXdesG0eYgAmYgAmYgAmYwAQCkx7fTGjXRU3ABEzABEzABEygRuCallLy3gmuVeOACZiACZiACZiACUwj4Mc30/i5tAmYgAmYgAmYwIYI+PHNhkC6GhMwARMwARMwgWkEbJRM4+fSJmACJmACJmACGyJwsqF6XI0JmIAJ9BLQNwr4nlH80B7fK8DxLyeSv1t0WcR/TcAEjpWAjZJjHVn3ywT2j8CL4qNPQTKd8x9Cv8rzP7LsTMAETGDlxzdWAhMwgUECMiBuy3+V/68wJgbLtGR4pLL8m4no+AowX4zkH3/ZmYAJmMDKKyVWAhMwgUECMhy+KdMdHfmEAP8BNMfxef8vOQVdxgRMYBkEbJQsY5zdSxOYTKCyyjH2H2qGtlX+dUMIjJQfhcHTSHLQBExgiQT8+GaJo+4+m0AegXsqhhExeWOq6uCRzX35O3miuJQJmMAxErBRcoyj6j6ZwDwE2A+StUpSFUcGCW/esJ+Ex0GTDZxq3T43ARM4bAJ+fHPY42fpTWAWAhXD4Yca+CmPMcLqxnP50ikfr/nyb+Bxd+UfymN0PJDHfVaeD5enq1VR7xMdWXWJYY60Y2cCJrBwAjZKFq4A7r4JNAnIQGBF5L08KxnBWNAxbm5trpSUr/kqz98q80YeQwTDg1d+WREJRonCGCvEkRbfuGFfyRN5OxMwARPw2zfWARMwgSsCMhZY+cAgwXCorl5w/k1x5eMWnbNCgtERHWnsE2G1BPebfDWdb5JQP8fSqR4MEzsTMAETsFFiHTABE6gRwIg4laHQfFOG1ZPyMUxR4ovyVQ0XNq2WhovS/izyhYPCN6phn5uACZhAk4A3ujaJOGwCyybwl7pfe0QjY4LVDR69xEc4gZDi+XZJ1VH2bTXC5yZgAiYwhoCNkjG0nNcEjphAYXxggNSMD4UxNlZKrxkrxEWnNFZSKFuuplAfPubx0QRMwASGCHij6xAhp5vA8ghUH8nQe96UCasiMjLYM8L5L3n2nrDRFWOFfSHnOq+WfaqwN7EKjJ0JmEAaAa+UpHFyLhM4egIyIM7VSQwMHtUEpziMEN6U+XIZs7pXGB6sjOB/KcxqCEZK6RRH2ucywicmYAImkEDg2sUF/8rCzgRMwATCIxoMjPBar4435dkjgrHCBliMjE8yONjMSj7ivsuvFH4pjwHDqgpv12CslI9yFLYzARMwgUEC/weXMZOqdFgFVAAAAABJRU5ErkJggg==\n", "text/latex": "$\\displaystyle \\frac{{{u}_{(0,0)}} dx^{2} + dt \\kappa \\left(- 4 {{u}_{(0,0)}} + {{u}_{(1,0)}} + {{u}_{(0,1)}} + {{u}_{(0,-1)}} + {{u}_{(-1,0)}}\\right)}{dx^{2}}$" }, "metadata": {}, "execution_count": 7 } ], "source": [ "heat_pde_discretized.simplify()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While combining the two fractions on the right as desired, it also put everything above a common denominator. If we generated the kernel from this, we'd be redundantly multiplying $u_{(0,0)}$ by $dx^2$. Let's try something else. Instead of applying simplify to the entire equation, we could apply it only to the second summand.\n", "\n", "The outermost operation of heat_pde_discretized is a $+$, so heat_pde_discretized is an instance of sp.Sum. We take it apart by accessing its arguments, simplify the right hand summand, and put it back together again." ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgcAAAAvCAYAAAB3wWlJAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAQS0lEQVR4Ae2d6ZUctRqG23McwGBnABnYOII7ZOALEdjOwBz/gn8+kIEhAmxnAERgQwaQAb6Twdz3qSkJlWpTqbt6enl1jloq7Xpq+b7SUn3v5uZmY3M6BL7//vvn6s3fcn87nV65JyZgAiZgArsiIPnwRvbFVHkXU5GOOy4COtlP1eLHVgyO67y5tSZgAiawZwI/SE78OlXnPY8cTOE5njid6Eu19ne5j4+n1W6pCZiACZjAXRCQrGCU+VLuj0P1e+RgiMpxhr1Ts98cZ9PdahMwARMwgX0SkFLwk+p7IZcXy5653wtxQBUBAX6kjD/Lfi77VseT8zlVlYxkUl3U+aXcr0aS7Cy47ScX1N76t7PGuyATMAETMIGUAC+UyK3/poH4PXKQE6k8lrD8U5YhfbSwzlyOwhm6+UsWIb6G+VaFogXuwzBC8WAfFbkOEzABEzCBVQkgN54io/JarBzkRLY4FuCrNnu+U4BwFINPWxQ/lZW5o45CMpW4Nk79e1mb91TycY5l4W0zQsCMRsC0weYzzac21lyXkxOza+X6U/brPLeVg5zIdscM67ONEOCpGQtP01T5VRfTGZiPt846v2099Cvv2zoVHm6pl2oa1macgBmNsyHGfKb51Maaax05ZIenFerYFedihCAfNSAz4e/x5EZCd1tBQ9lDCkle1bbH36it+5q62Latzm8CJmACJlBG4C8l+zJPej8P8HEZAQlKpgl+kP1b9h9ZlALe4l/LbhSP0GbRHsKftI8Uxnz9B7np1hG2H/4SwuSS/pUsQ9evQ7j8Y+aJIq7HIncRrjYwnTC6E6JtM+39QvYPHXeUCB3T72dyV22n6jhoY07Tp8d8pvkQa0bzjGpSnDlXphVYF4eNz2hPK1RcSQKI4P9D9lv5sQh7FAVMM3KgsN9kGapphKr8X3EsGxUD+VEEUCg4ORsd42fuhzwM9QyNQii4YyhjrbUMtAnF5louStCYeaV4+sW6h8ChSatwlAYWvMSLrolofxROHHWsZlT+oayVOFhOB8LIfObvgoNlNN/0fooDue5o2DlzDfKj8xz2yEH/ep0M0cWMMOZNGKUgFZj42bGQC0HWGzTCX25umkUgyoMiwdcNEcLhrbt0W+ID5Uvbkdex7THbFtkNMWgUh0LzoY2kzeFCC+lH+6+8KA5MicT2yx+UC0ZjHup4tO5QAe5UPsX9KPtG9s62X6ruXXPiOuT64fwUffhK6UbZKu5OGa3AZ7SvXC9DRm0YZXrXfGjvITAa4jYXdm5c53gMxU8xGkpPmPKMXuOKW3I/B5mFLInGykFEUezhhDD8EoR4yMhowtC6grFw8iE4UQh44w6jBYQvMTzQcoHc5Fe5xP0ui1tqGN1olBm5CO/R6YS2QIR7UH4QVs20SlIZ/c/DNspDm5rRlJC2re8fuc3oitwr2V9lJxUlxdPOuXzvlI5PhpYqG/SbtuemuYFUzpCigXLYW9jTFrBLTigaoW1F51btOnRGu+RT0tfOeRWfEqandA0tZtQBVnhwblwLsXSSFTLK85Scv9LrdVB+3O/U6IMSAgjAznC/Ti4PaIZkOtsJx8KTSnjAv5VlJIL8rD/gLXrJ/Py10nc0Ph03RmURV/RW2WaJjvLSHpSg+FYfIxNPW8dGLn2BQ1SaFNaUobAOrzY7SlaueBD2nzaeMhlRQTn4XHaqHbP52rJQDugTXCaN0gwJ/43CGeGhPXF6aLKgNlLpmzrlbs1JZaCMoYjQllJz0Ix2yUdAZvuaQythqjRcjydxDdUwypmVHJ8h1xIsnTQljDoZbg9mr/EF12uQHx0lwcrBAPWxIMFG+GE7SoCO4/RAlhdBsOEkhXDKkGW0IAhOtLsQ/5P8/1NatMLw9lwkzEL5O3Rp3xO1B8UlNbxhIRwJ520vvInzxsxxKngb4aWwMLKQlsOoQBTA8gceuRJAedSZhzdlLcwHZ87J0AhPU94efrbiVNO+I2O0FZ+Ffa3BefTX0B4YnSXXmk7X5Fl4/kquV2QaJn12b6wc3EJZ+psLqjivrhOHQOStjjQxnAraOMIRlj3FgTSpUXqUBE5u56SlaeSnvHBys6j6Q9VNvdiOUTjKS1hsmcYh3HMlgP43ZbR94bPSQTHK+zTWB7TZoNmm9QX/knysjaBNd6kcbMsp9HuJe0yMtuWzpK9LGIa0p3ANrc0osFringLXJf3dJu2S81fCNTxfedZGcxF99swSQLApEcKOB1hjFIYywJstuwswzKMH5QHojV9hnFDeuIMAjYKTTJkJ+b9IysqSxENO/pfxaH0P/Ri6OEObmxao3TChXaG/9AV+GPh1LsQmdPgHhkP1Daf+N3QoH/XHc/dv0r361uJU04lDZLQWn6G+1jA75WtoV4zMtYbA9nmGzl/J9cozkZc20kbjkYOIotjDsOfPAvlS7kPZX2QR9MxFEpbOozPkzip5RgA2ctN56sssLUkwz2T54NATuWlZxA0ZhC9TD6tOP6h82hIEK1MC73T8QW7oE32FC3Nh7DTgIc96h8AFTsHQ986FOHCcpu0IjBDRunk5IZo68nylCkkoYw13W041bTomRtvyWdLXGpancA1VMdK9/UjAuL/nDDto8ntvLs8pcN2szCgwXHL+SrjynO6NEls5CLgLXZ18TgwKQm56Ye0NguLQM4obC2fIu3jYW+UwxE+beEvvneBexZUBqiOuDxgqom1Dj4HSDoXx4LhMy1H+sF4BBeTPNG7gOEYvzEfZSx9asa5deLblVNOGY2K0LZ+Ffa3BefTXUC0j5eO+HHxu1YDM8hw9V/qzMqMG2cLzV8IV2dF7Eb1oavPPsRN4qw6sddPunE17A3HR5ua1Apq1GEQoHf733AwcY+RnhCR/e5nNd5u7GfngU6HbGBQx7OpG/eRhPMRpsu5zYTTBZ/J6GOEzyTSJ5HycwjW0JqMEV7H3VLgWd7g04cj1Onn+krInuVK20jIihAzpmHs3NzedAB8cHwGdYE4u2yA/O5bWq63s+OgNPyo8CH6mJvgcMx+bisJY/qavCme7Zxxhmcun9BulGayTuEM1Q21WGDf8C1mUJ3gwtZN+5+FsGA3xEYuNwkevI8X1+ChskillYpTuJK6hti+LGDUAFv6cI9eFiLimJq89xfeu17nzF9qgvJPXq+KZ8u58bybk3aAcnIv97rvvrmSfn2J/1a93sk+PpW/tufihtr1L+6r0l7K/1tZ3V/m24XQOjMxn/vm9T0Y194nad3b3ZiWnRc/3Eq5K85fs50PtuYhawnl4GELBnqJ5pk69OpaOSWNlfQRTBGjNi4zy8Ma8dH0Fb0m8bR+VqeV0LozMZ/5y3jOj+Qb1U5zVvdnv/nzIGvezymTUgAXzcdo2bcm5KQdp30/KrxPM0PtruWGo8OD7p7YirJlaWKqwsSWU/hYZpeUm4GNTgzdBUSF3mEjtruF0NozMZ/7i3Aej+Vb0U6hd53hv9kHMh+z0fhZ3Xso6fwTYawLDCe3ww0u5b2R7w+4KY8j6cmjo4ZjC1Iensi+Pqc1L26r+PZe9WprvLtOrvateW2uXvy92a/ZjzbLNZ37o/xQY1fThFK47+n1o/Zhrj+KR9ZPP3WZBorQI9qKz8OuptAf2qseFbfKHoYd7Pc3iyALa/qGBhb35R9YDN9cETMAETMAE1idwIUHJSki+sodhO1z+0QTC8n3npG0MAld28bxxyD/nqmw+LGRjAiZgAiZgAiawJwL3VU/VX6XSPgluRhXIH+dy5Q9z3mxFe6hjvng2aabyKG7J/1I39SgPH3Rg0VpuHhCg+KGFaVN/t9uUo3ze95kT9bEJmIAJmMDJEYjfOZDgQ5iyJ/Iz+ZvFXnIZEeCDH4/l74we6JhFZExBxC/gyY+ywAr0ZtheLmUyXTH6gZ6SPG057MWcVTRU36hRfqZNPK0wSsgRJmACJmACJrDZpLsVEPLhE7aBDcJ0I6HaUQzaSEYI8k8uEha3mCkf/iu5U9MOs3mScpauam+bascETMAETMAETKCUQKocIMBzJYA3/kbYS0A/l02FM0I/VQTIT3ycYpAfcy3LuoaeUf4leRpFo1eIA0zABEzABEzABHZKgDUHwSDUmzl5AiS4Eej8IcNPHMvEv9xthTpCPzWp4pCGf9JBLDeNkH9JnpL/pc6K96EJmMCpE9DziOdI+AAYLxwYPq+dP6NuY/xrAiYwSyBVDpb8VSo3IEK/xKAYjCkBY/mH8nCjhxt/LJ/DTcAEzo8AW7HjImP5me78Q5b/5rAxAROoIBCVg1bLjosLk7KGwhD2uVaeH4ciSJtPNYS4JXlKlZFQtl0TMIEDJ6DnDiOUP8ui+L9NhfyCpjPlyRcwwzQn65gIeySbT5UuKNZJTeB8CVxUdh1hj9CPRjchYQj7obf7wRt0YR7KHVMyYjtmPLRvTCGZyepoEzCBXRNAeMs+Vrk8T9gtVWMYNfhYk9F5TMAEhgnEkYPh6OHQ9oYeUgJeK8eVbKMMKB3+93IboS6XB8AruemWxMk8Sh8M9W31P+qqN7xZhDLtmoAJ3DEB3Zc8JzBV96fyh3VRt6Xc/sFW+v2WEG7XBEygkMBFYbqhZNx8HQVBx3zfgA8fMQfIlw2ZkniWZCY9w33NFknC5Z/LE7Kzc+J9OLBrAiZwMgS4t/Nt1FWd0/OEaQqeL4xG2JiACVQSqBo5aOtiXo/hvHQUAGHfOU7bpThGFPjIUlQOiJ/K08Yz4kC6bacVKMbGBEzgsAgwclA1apB2Q88HXj54LvHRNk8fpnDsN4GFBKqVA918v8nyl498cbBYaCttzYMgKCILu+fkJmACh0SA54Xaw/3MM4NPrKMU8LbP9GI0SscLAV9cxTyRZQSSvN/IYj4oTRxJbMuNX2Ntjzdyi59Nt8X61wRMAAIX22DQjcfIwQu5zZt9YVkoE8VavdLygGAlsm/yQsBOZgKHSED3MC8GbDFEiGOZUkRRwOQjB0xN8r8qpOEbJ+xouCKfXJ43Id9GYSgNbF98Iz87FFA2SPdJ1sYETKCCwFbKAfXpRhydRhhqj9Lni4eGkqVhbG/KHxxpvP0mYAIHTkD3MAL9nSxKQaro42fHQnxhkJ8Xgij85SeOqcjw7HiQxaNwBMUDP5a1TbFMHduYgAksIFA9rZDWseZNuGbZaR/sNwETWJUAwp4/ZQsCPlSGUI/TA23gR6VLFQgWF0YFQnGdb6/o+LNQmF0TMIHdENh65GA3zXApJmACJ07ga/WvMwIooc5oAlMCne8bKDz/Lgp5f5G1MQET2BMBKwd7Au1qTOBcCbRKAIpARwnQMUJ/o/iO0kBYMIpjZIG8cXSB8rAhjV0TMIHdE9jJtMLum+USTcAETpBAOlVA9/i+QfhgGmsK8LOIkLUJLEhEaWDR87X8ad78Q2pKYmMCJrBLAh452CVNl2UCJtAjIMF+rUAEPVMIjVEYygC7Cj7ehmy+ahUARgqwn3TM6ADKQjQKI47dCzYmYAIrErh3c3OzYvEu2gRMwASaqQMEPdsREewPZVlDgNLAQkXC+G4Kiw5JR1jzqXQds50RRYJRBnYhoDTEKQYd25iACaxA4P8C3RkI5peT+QAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle {{u}_{(0,0)}} + \\frac{dt \\kappa \\left(- 4 {{u}_{(0,0)}} + {{u}_{(1,0)}} + {{u}_{(0,1)}} + {{u}_{(0,-1)}} + {{u}_{(-1,0)}}\\right)}{dx^{2}}$" ], "text/plain": [ " dt⋅κ⋅(-4⋅u_C + u_E + u_N + u_S + u_W)\n", "u_C + ─────────────────────────────────────\n", " 2 \n", " dx " ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "heat_pde_discretized = heat_pde_discretized.args[1] + heat_pde_discretized.args[0].simplify()\n", "heat_pde_discretized" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That looks a lot better! There is nothing left to simplify. The right-hand summand still contains a division by $dx^2$, though. Due to their inefficiency, floating-point divisions should be replaced by multiplication with their reciprocals. Before we can eliminate the division, we need to wrap the equation inside an AssignmentCollection. On this collection we can apply add_subexpressions_for_divisions to replace the division by a factor $\\xi_1 = \\frac{1}{dx^2}$ which in the kernel will be computed ahead of the loop." ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Subexpressions:
 $$\\xi_{0} \\leftarrow \\frac{1}{dx^{2}}$$
Main Assignments:
 $${{u_tmp}_{(0,0)}} \\leftarrow {{u}_{(0,0)}} + dt \\kappa \\xi_{0} \\left(- 4 {{u}_{(0,0)}} + {{u}_{(1,0)}} + {{u}_{(0,1)}} + {{u}_{(0,-1)}} + {{u}_{(-1,0)}}\\right)$$
" ], "text/plain": [ "AssignmentCollection: u_tmp_C, <- f(u_S, u_N, u_E, u_W, dt, u_C, kappa, dx)" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@ps.kernel\n", "def update():\n", " u_tmp.center @= heat_pde_discretized\n", "\n", "ac = ps.AssignmentCollection(update)\n", "ac = ps.simp.simplifications.add_subexpressions_for_divisions(ac)\n", "ac" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our numeric solver's symbolic representation is now complete! Next, we use pystencils to generate and compile a C implementation of our kernel. The code is generated as shown below, compiled into a shared library and then bound to kernel_func. All unbound sympy symbols (dx, dt and kappa) as well as the fields u and u_tmp are arguments to the generated kernel function. " ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
FUNC_PREFIX void kernel(double * RESTRICT const _data_u, double * RESTRICT _data_u_tmp, int64_t const _size_u_0, int64_t const _size_u_1, int64_t const _stride_u_0, int64_t const _stride_u_1, int64_t const _stride_u_tmp_0, int64_t const _stride_u_tmp_1, double dt, double dx, double kappa)\n",
"{\n",
"   {\n",
"      #pragma omp for schedule(static)\n",
"      for (int ctr_1 = 1; ctr_1 < _size_u_1 - 1; ctr_1 += 1)\n",
"      {\n",
"         double * RESTRICT _data_u_tmp_10 = _data_u_tmp + _stride_u_tmp_1*ctr_1;\n",
"         double * RESTRICT _data_u_10 = _data_u + _stride_u_1*ctr_1;\n",
"         double * RESTRICT _data_u_11 = _data_u + _stride_u_1*ctr_1 + _stride_u_1;\n",
"         double * RESTRICT _data_u_1m1 = _data_u + _stride_u_1*ctr_1 - _stride_u_1;\n",
"         for (int ctr_0 = 1; ctr_0 < _size_u_0 - 1; ctr_0 += 1)\n",
"         {\n",
"            _data_u_tmp_10[_stride_u_tmp_0*ctr_0] = dt*kappa*(-4.0*_data_u_10[_stride_u_0*ctr_0] + _data_u_10[_stride_u_0*ctr_0 + _stride_u_0] + _data_u_10[_stride_u_0*ctr_0 - _stride_u_0] + _data_u_11[_stride_u_0*ctr_0] + _data_u_1m1[_stride_u_0*ctr_0])/(dx*dx) + _data_u_10[_stride_u_0*ctr_0];\n",
"         }\n",
"      }\n",
"   }\n",
"}\n",
"