From 46138613ef12acad6c673770797db8ca673a84a7 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Mon, 26 Feb 2018 18:23:20 -0600
Subject: [PATCH] Add PDE-reduced translation notebook
---
.gitignore | 1 +
.../PDE-reduction and translations.ipynb | 305 ++++++++++++++++++
2 files changed, 306 insertions(+)
create mode 100644 contrib/translations/PDE-reduction and translations.ipynb
diff --git a/.gitignore b/.gitignore
index d22b6772..7fcc2940 100644
--- a/.gitignore
+++ b/.gitignore
@@ -15,3 +15,4 @@ a.out
doc/_build
.cache
.DS_Store
+.ipynb_checkpoints
diff --git a/contrib/translations/PDE-reduction and translations.ipynb b/contrib/translations/PDE-reduction and translations.ipynb
new file mode 100644
index 00000000..043f61c4
--- /dev/null
+++ b/contrib/translations/PDE-reduction and translations.ipynb
@@ -0,0 +1,305 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 53,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import pyopencl as cl\n",
+ "import sumpy.toys as t\n",
+ "import numpy as np\n",
+ "import numpy.linalg as la\n",
+ "import matplotlib.pyplot as plt\n",
+ "from sumpy.visualization import FieldPlotter\n",
+ "from pytools import add_tuples\n",
+ "\n",
+ "from sumpy.expansion.local import VolumeTaylorLocalExpansion\n",
+ "from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansion\n",
+ "from sumpy.kernel import (YukawaKernel, HelmholtzKernel, LaplaceKernel)\n",
+ "\n",
+ "order = 4\n",
+ "\n",
+ "knl = LaplaceKernel(2)\n",
+ "\n",
+ "pde = [(1, (2,0)), (1, (0, 2))]\n",
+ "\n",
+ "mpole_expn = VolumeTaylorMultipoleExpansion(knl, order)\n",
+ "local_expn = VolumeTaylorLocalExpansion(knl, order)\n",
+ "\n",
+ "cl_ctx = cl.create_some_context(answers=[\"port\"])\n",
+ "\n",
+ "tctx = t.ToyContext(\n",
+ " cl_ctx,\n",
+ " knl,\n",
+ " mpole_expn_class=type(mpole_expn),\n",
+ " local_expn_class=type(local_expn),\n",
+ " \n",
+ " #YukawaKernel(2), extra_kernel_kwargs={\"lam\": 5},\n",
+ " #HelmholtzKernel(2), extra_kernel_kwargs={\"k\": 0.3},\n",
+ " )\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 54,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pt_src = t.PointSources(\n",
+ " tctx,\n",
+ " np.random.rand(2, 50) - 0.5,\n",
+ " np.ones(50))\n",
+ "\n",
+ "mexp = t.multipole_expand(pt_src, [0, 0], order)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 55,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([ 5.00000000e+01, 2.49673488e+00, 1.47014861e-01,\n",
+ " 2.60554926e+00, -9.96815120e-01, 2.40178914e+00,\n",
+ " 5.42629997e-02, 3.38317711e-03, 2.93648067e-01,\n",
+ " -5.10921327e-03, 3.78790947e-02, -2.80879544e-02,\n",
+ " 1.38914975e-01, -2.60301325e-02, 3.05226127e-02])"
+ ]
+ },
+ "execution_count": 55,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "mexp.coeffs"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 56,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def build_pde_mat(expn, pde):\n",
+ " coeff_ids = expn.get_coefficient_identifiers()\n",
+ " id_to_index = expn._storage_loc_dict\n",
+ " \n",
+ " # FIXME: specific to scalar PDEs\n",
+ " pde_mat = np.zeros((len(coeff_ids), len(coeff_ids)))\n",
+ " \n",
+ " row = 0\n",
+ " for base_coeff_id in coeff_ids:\n",
+ " valid = True\n",
+ " \n",
+ " for pde_coeff, coeff_id_offset in pde:\n",
+ " other_coeff = add_tuples(base_coeff_id, coeff_id_offset)\n",
+ " if not other_coeff in id_to_index:\n",
+ " valid = False\n",
+ " break\n",
+ " \n",
+ " pde_mat[row, id_to_index[other_coeff]] = pde_coeff\n",
+ " \n",
+ " if valid:\n",
+ " row += 1\n",
+ " else:\n",
+ " pde_mat[row] = 0\n",
+ " \n",
+ " return pde_mat[:row]\n",
+ "\n",
+ "pde_mat = build_pde_mat(mpole_expn, pde)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 57,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def find_nullspace(mat, tol=1e-10):\n",
+ " u, sig, vt = la.svd(pde_mat, full_matrices=True)\n",
+ " zerosig = np.where(np.abs(sig) < tol)[0]\n",
+ " if zerosig:\n",
+ " nullsp_start = zerosig[0]\n",
+ " assert np.array_equal(zerosig, np.arange(nullsp_start, pde_mat.shape[1]))\n",
+ " else:\n",
+ " nullsp_start = pde_mat.shape[0]\n",
+ " \n",
+ " return vt[nullsp_start:].T\n",
+ " \n",
+ "nullsp = find_nullspace(pde_mat)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 64,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "7.3852783309125538e-16"
+ ]
+ },
+ "execution_count": 64,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "la.norm(pde_mat @ nullsp)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 61,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def build_translation_mat(mexp, to_center):\n",
+ " n = len(mexp.coeffs)\n",
+ " result = np.zeros((n, n))\n",
+ " \n",
+ " for j in range(n):\n",
+ " unit_coeffs = np.zeros(n)\n",
+ " unit_coeffs[j] = 1\n",
+ " unit_mexp = type(mexp)(\n",
+ " mexp.toy_ctx, mexp.center, mexp.rscale, mexp.order,\n",
+ " unit_coeffs, derived_from=mexp)\n",
+ " \n",
+ " result[:, j] = t.multipole_expand(unit_mexp, to_center).coeffs\n",
+ " \n",
+ " return result\n",
+ "\n",
+ "tmat = build_translation_mat(mexp, np.array([0, 0.5]))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 62,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[ 1. , 0. , 0. , 0. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ],\n",
+ " [ 0. , 1. , 0. , 0. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ],\n",
+ " [ 0.5 , 0. , 1. , 0. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ],\n",
+ " [ 0. , 0. , 0. , 1. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ],\n",
+ " [ 0. , 0.5 , 0. , 0. , 1. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ],\n",
+ " [ 0.125 , 0. , 0.5 , 0. , 0. ,\n",
+ " 1. , 0. , 0. , 0. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ],\n",
+ " [ 0. , 0. , 0. , 0. , 0. ,\n",
+ " 0. , 1. , 0. , 0. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ],\n",
+ " [ 0. , 0. , 0. , 0.5 , 0. ,\n",
+ " 0. , 0. , 1. , 0. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ],\n",
+ " [ 0. , 0.125 , 0. , 0. , 0.5 ,\n",
+ " 0. , 0. , 0. , 1. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ],\n",
+ " [ 0.02083333, 0. , 0.125 , 0. , 0. ,\n",
+ " 0.5 , 0. , 0. , 0. , 1. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ],\n",
+ " [ 0. , 0. , 0. , 0. , 0. ,\n",
+ " 0. , 0. , 0. , 0. , 0. ,\n",
+ " 1. , 0. , 0. , 0. , 0. ],\n",
+ " [ 0. , 0. , 0. , 0. , 0. ,\n",
+ " 0. , 0.5 , 0. , 0. , 0. ,\n",
+ " 0. , 1. , 0. , 0. , 0. ],\n",
+ " [ 0. , 0. , 0. , 0.125 , 0. ,\n",
+ " 0. , 0. , 0.5 , 0. , 0. ,\n",
+ " 0. , 0. , 1. , 0. , 0. ],\n",
+ " [ 0. , 0.02083333, 0. , 0. , 0.125 ,\n",
+ " 0. , 0. , 0. , 0.5 , 0. ,\n",
+ " 0. , 0. , 0. , 1. , 0. ],\n",
+ " [ 0.00260417, 0. , 0.02083333, 0. , 0. ,\n",
+ " 0.125 , 0. , 0. , 0. , 0.5 ,\n",
+ " 0. , 0. , 0. , 0. , 1. ]])"
+ ]
+ },
+ "execution_count": 62,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "tmat"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 63,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "<matplotlib.image.AxesImage at 0x7f981a252748>"
+ ]
+ },
+ "execution_count": 63,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAADcBJREFUeJzt3X+s3XV9x/Hn29tS1gLSgijQZoVJ\nmoBZBmmw6uLMOn6OUP/wj5K5dWJCzOIGxkVLSGayv+ZcdFtmZhpwYxkBM4TZGBg0VbMsWaultoVa\naCtjUFspk4UKxpay9/4430su13N/fX/dc/08H8nN+fX59vvu99zX/Z7zPd/PeUdmIqk8b5vvAiTN\nD8MvFcrwS4Uy/FKhDL9UKMMvFcrwS4Uy/FKhDL9UqEV9ruz8FWO5etXiOS93cN/SDqqRfvn8nNc4\nlSdjNmN7Df/qVYv57mOr5rzcdRf9RgfVSL98dub2WY/1Zb9UKMMvFapR+CPi+oh4JiIOR8TmtoqS\n1L3a4Y+IMeDLwA3A5cAtEXF5W4VJ6laTPf/VwOHMfDYzTwEPABvaKUtS15qE/2LghQm3j1T3SVoA\nmoR/2GeJv/C1QBFxW0TsiohdL/3kjQark9SmJuE/Akz80H4lcHTyoMzckplrM3PtO84ba7A6SW1q\nEv7vAZdFxCURcQawEdjaTlmSulb7DL/MPB0RnwQeA8aAr2bm/tYqk9SpRqf3ZuYjwCMt1SKpR57h\nJxXK8EuF6nVW38F9S2vN0Hvs6J7a63RGoDSce36pUIZfKpThlwpl+KVCGX6pUIZfKpThlwpl+KVC\nGX6pUIZfKpThlwpl+KVCGX6pUL3O6ju5ahmHP71uzstdd1H9dTojUBrOPb9UKMMvFcrwS4Vq0qtv\nVUR8OyIORMT+iLi9zcIkdavJAb/TwKczc3dEnA08ERHbMvMHLdUmqUO19/yZeSwzd1fXfwocwF59\n0oLRynv+iFgNXAnsbOPfk9S9xuGPiLOArwN3ZOaJIY+/2ajzjVdfa7o6SS1pFP6IWMwg+Pdl5kPD\nxkxs1Dl21rImq5PUoiZH+wO4BziQmV9sryRJfWiy5/8A8PvAb0fEnurnxpbqktSxJl16/wOIFmuR\n1CPP8JMKZfilQkVm9rayc2JFvjfW97a+pupOB3YqsObLztzOiXx5Vm/H3fNLhTL8UqEMv1Qowy8V\nyvBLhTL8UqEMv1Qowy8VyvBLhTL8UqEMv1Qowy8VyvBLheq1UWddh7809+ae4979qR21l607O8/m\noFoI3PNLhTL8UqEMv1SoNpp2jEXE9yPim20UJKkfbez5b2fQp0/SAtK0Y89K4HeBu9spR1Jfmu75\n/xr4DPB/LdQiqUdN2nXdBBzPzCdmGPdmo87XOVl3dZJa1rRd180R8RzwAIO2Xf88edDERp2LWdJg\ndZLaVDv8mXlnZq7MzNXARuBbmfnR1iqT1Ck/55cK1cq5/Zn5HeA7bfxbkvrhnl8qlOGXCtXrlN44\ncwlj714z5+WaTMudj+nATabl2hxUfXHPLxXK8EuFMvxSoQy/VCjDLxXK8EuFMvxSoQy/VCjDLxXK\n8EuFMvxSoQy/VCjDLxUqMrO3lZ0TK/K9sb639ZXE5qAC2JnbOZEvx2zGuueXCmX4pUIZfqlQTdt1\nnRsRD0bE0xFxICLe11ZhkrrV9Gu8/gb4t8z8SEScASxtoSZJPagd/og4B/gg8IcAmXkKONVOWZK6\n1uRl/6XAS8A/RMT3I+LuiFjWUl2SOtYk/IuAq4C/z8wrgdeAzZMH2ahTGk1Nwn8EOJKZO6vbDzL4\nY/AWNuqURlOTRp0/Bl6IiPEv4l8P/KCVqiR1runR/j8G7quO9D8LfKx5SZL60Cj8mbkHWNtSLZJ6\n5Bl+UqEMv1SoXht1LjR1m3w2aSxal81BNVfu+aVCGX6pUIZfKpThlwpl+KVCGX6pUIZfKpThlwpl\n+KVCGX6pUIZfKpThlwpl+KVCLYhGnWNXrJl50BTe2P9M7WXrqjsbEOZnRmBdNgcdPTbqlDQjwy8V\nyvBLhWraqPNTEbE/Ip6KiPsj4sy2CpPUrdrhj4iLgT8B1mbme4AxYGNbhUnqVtOX/YuAX4mIRQw6\n9B5tXpKkPjTp2PMj4K+A54FjwCuZ+XhbhUnqVpOX/cuBDcAlwEXAsoj46JBxNuqURlCTl/2/A/xX\nZr6Uma8DDwHvnzzIRp3SaGoS/ueBdRGxNCKCQaPOA+2UJalrTd7z72TQlns38GT1b21pqS5JHWva\nqPNzwOdaqkVSjzzDTyqU4ZcK1e+U3redl+uW3DDn5fJk/Y8IS5kOvJCmAoPNQbvilF5JMzL8UqEM\nv1Qowy8VyvBLhTL8UqEMv1Qowy8VyvBLhTL8UqEMv1Qowy8VyvBLhVoQjTqlcTYHnZ6z+iTNyPBL\nhTL8UqFmDH9EfDUijkfEUxPuWxER2yLiUHW5vNsyJbVtNnv+fwSun3TfZmB7Zl4GbK9uS1pAZgx/\nZv478PKkuzcA91bX7wU+3HJdkjpW9z3/OzPzGEB1eUF7JUnqQ6OmHbMREbcBtwGcydKuVydpluru\n+V+MiAsBqsvjUw20Uac0muqGfyuwqbq+CfhGO+VI6stsPuq7H/hPYE1EHImIjwN/AVwTEYeAa6rb\nkhaQGd/zZ+YtUzzkSfrSAuYZflKhDL9UqM4/6lM/6jb4hIXV5LPJtFybg76Ve36pUIZfKpThlwpl\n+KVCGX6pUIZfKpThlwpl+KVCGX6pUIZfKpThlwpl+KVCGX6pUDbqnMbYFWtqLffG/mdarqRbdWcE\nLqTZgE0spOagNuqUNCPDLxXK8EuFqtuo8wsR8XRE7IuIhyPi3G7LlNS2uo06twHvycxfBw4Cd7Zc\nl6SO1WrUmZmPZ+bp6uYOYGUHtUnqUBvv+W8FHm3h35HUo0bf3hsRdwGngfumGWOjTmkE1Q5/RGwC\nbgLW5zRnCmXmFmALDE7yqbs+Se2qFf6IuB74LPBbmfmzdkuS1Ie6jTr/Djgb2BYReyLiKx3XKall\ndRt13tNBLZJ65Bl+UqEMv1SoBTGlN5Ysqb3OPHmy9rJ11Z0KDAtrOnApzUGb6Ls5qFN6Jc3I8EuF\nMvxSoQy/VCjDLxXK8EuFMvxSoQy/VCjDLxXK8EuFMvxSoQy/VCjDLxVqQczqa6KUGYELaTYg2Bx0\nJnVnA1593Qvs2vtzZ/VJmprhlwpl+KVC1WrUOeGxP42IjIjzuylPUlfqNuokIlYB1wDPt1yTpB7U\natRZ+RLwGcAuPNICVOs9f0TcDPwoM/e2XI+knsy5XVdELAXuAq6d5XgbdUojqM6e/9eAS4C9EfEc\nsBLYHRHvGjY4M7dk5trMXLuY+ifcSGrXnPf8mfkkcMH47eoPwNrM/J8W65LUsbqNOiUtcHUbdU58\nfHVr1UjqjWf4SYUy/FKhep3SGxEvAf89xcPnA6N00HDU6oHRq8l6pjcf9fxqZr5jNgN7Df90ImJX\nZq6d7zrGjVo9MHo1Wc/0Rq2eyXzZLxXK8EuFGqXwb5nvAiYZtXpg9GqynumNWj1vMTLv+SX1a5T2\n/JJ61Hv4I+L6iHgmIg5HxOYhjy+JiK9Vj++MiNUd1rIqIr4dEQciYn9E3D5kzIci4pWI2FP9/FlX\n9UxY53MR8WS1vl1DHo+I+NtqG+2LiKs6rGXNhP/7nog4ERF3TBrT6TYa9m1SEbEiIrZFxKHqcvkU\ny26qxhyKiE0d1vOFiHi6ej4ejohzp1h22ue2V5nZ2w8wBvwQuBQ4A9gLXD5pzB8BX6mubwS+1mE9\nFwJXVdfPBg4OqedDwDd73k7PAedP8/iNwKNAAOuAnT0+fz9m8Flyb9sI+CBwFfDUhPv+EthcXd8M\nfH7IciuAZ6vL5dX15R3Vcy2wqLr++WH1zOa57fOn7z3/1cDhzHw2M08BDwAbJo3ZANxbXX8QWB8R\ns/oe8rnKzGOZubu6/lPgAHBxF+tq2Qbgn3JgB3BuRFzYw3rXAz/MzKlO1OpEDv82qYm/J/cCHx6y\n6HXAtsx8OTP/F9jGkK+ka6OezHw8M09XN3cwmOo+0voO/8XACxNuH+EXw/bmmGpjvgKc13Vh1duL\nK4GdQx5+X0TsjYhHI+KKrmth8NVoj0fEE9WXoUw2m+3YhY3A/VM81vc2emdmHoPBH3EmTDOfYL62\n060MXpkNM9Nz25s5z+dvaNgefPLHDbMZ06qIOAv4OnBHZp6Y9PBuBi9zX42IG4F/BS7rsh7gA5l5\nNCIuALZFxNPV3ubNkocs0/U2OgO4GbhzyMPzsY1mYz62013AaeC+KYbM9Nz2pu89/xFg1YTbK4Gj\nU42JiEXA2xn+BaKtiIjFDIJ/X2Y+NPnxzDyRma9W1x8BFnf9VeWZebS6PA48zODt0kSz2Y5tuwHY\nnZkvTn5gPrYR8OL4W53q8viQMb1up+qA4k3A72X1Bn+yWTy3vek7/N8DLouIS6o9yUZg66QxW4Hx\no7IfAb411YZsqjqWcA9wIDO/OMWYd40fc4iIqxlss590UU+1jmURcfb4dQYHkib3TNgK/EF11H8d\n8Mr4S+AO3cIUL/n73kaVib8nm4BvDBnzGHBtRCyvPg24trqvdRFxPfBZ4ObM/NkUY2bz3Pan7yOM\nDI5UH2Rw1P+u6r4/Z7DRAM4E/gU4DHwXuLTDWn6TwcvAfcCe6udG4BPAJ6oxnwT2M/hkYgfw/o63\nz6XVuvZW6x3fRhNrCuDL1TZ8ksHXqHVZ01IGYX77hPt620YM/ugcA15nsDf/OIPjQNuBQ9Xlimrs\nWuDuCcveWv0uHQY+1mE9hxkcXxj/PRr/xOoi4JHpntv5+vEMP6lQnuEnFcrwS4Uy/FKhDL9UKMMv\nFcrwS4Uy/FKhDL9UqP8H8zJzDMuV3ZMAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ "<matplotlib.figure.Figure at 0x7f981a4e8780>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.imshow(tmat)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "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.6.4+"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
--
GitLab