kifmm-laplacian.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Pkg.add(\"PyPlot\")\n",
    "Pkg.add(\"DataStructures\")\n",
    "Pkg.update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "getInversion (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "include(\"kernel.jl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# InitialIze KIFMM kernel \n",
    "\n",
    "The aim of this notebook is to illustrate usage of ``KIFMM`` to calculate ``p = K * q``.\n",
    "\n",
    "necessary notation:\n",
    "\n",
    "- ``k``: kernel, a struct to hold everything\n",
    "- ``source``: source points in 2D\n",
    "- ``normal``: normal vectors(optional) at each source points in 2D\n",
    "- ``target``: target points in 2D\n",
    "\n",
    "- ``f``: self defined kernel function, signature is\n",
    "\n",
    "    ```\n",
    "    function f(source, normal, target)\n",
    "        ...\n",
    "        return someNumber\n",
    "    end\n",
    "    ```\n",
    "\n",
    "    example usage:\n",
    "    ```\n",
    "    output = f([0.5, 0.2], [0.6, 0.8], [0.7, 0.4])\n",
    "    ```\n",
    "    \n",
    "### Calling ``initialize! `` \n",
    "\n",
    "When to initialize the kernel, one can call ``initialize!``.\n",
    "\n",
    "```\n",
    "initialize!(kernel, np, source, normal, target, charge, N, M, rank, maxLevel, f)\n",
    "```\n",
    "\n",
    "where\n",
    "\n",
    "- ``np`` is the number of surface points to approximate boundary integral(potential theory).\n",
    "- ``charge`` is the vector ``q`` to be applied.\n",
    "- ``N``, number of source\n",
    "- ``M``, number of target\n",
    "- ``rank``, max source points in each leaf box\n",
    "- ``maxLevel``, max level of tree\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.5177, 0.0001) is outside of source box.\n",
      "(0.4192, 0.0001) is outside of source box.\n",
      "(0.0001, 0.0707) is outside of source box.\n",
      "(0.0000, 0.6179) is outside of source box.\n",
      "(0.4192, 0.0001) is outside of source box.\n",
      "(0.0001, 0.0707) is outside of source box.\n",
      "(0.5177, 0.0001) is outside of source box.\n",
      "(0.0000, 0.6179) is outside of source box.\n",
      "(0.0001, 0.0707) is outside of source box.\n",
      "(0.4192, 0.0001) is outside of source box.\n",
      "(0.0001, 0.0707) is outside of source box.\n",
      "(0.0001, 0.0707) is outside of source box.\n",
      "(0.4192, 0.0001) is outside of source box.\n",
      "(0.4192, 0.0001) is outside of source box.\n",
      "(0.5177, 0.0001) is outside of source box.\n",
      "(0.5177, 0.0001) is outside of source box.\n",
      "(0.5177, 0.0001) is outside of source box.\n",
      "(0.0000, 0.6179) is outside of source box.\n",
      "(0.0000, 0.6179) is outside of source box.\n",
      "(0.0000, 0.6179) is outside of source box.\n",
      "(0.0000, 0.6179) is outside of source box.\n",
      "elapsed time: 1.676239151 seconds\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.676239151"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k = KIFMM()\n",
    "N = 9000\n",
    "source = rand(N, 2)\n",
    "target = rand(N, 2)\n",
    "charge = rand(N)\n",
    "\n",
    "### kernel function for Laplacian, this kernel does not require normal vectors\n",
    "function f(s, n, t)\n",
    "    r =  norm(s-t)\n",
    "    if (r < 1e-4)\n",
    "        return 0.\n",
    "    else\n",
    "        return log(norm(s-t))\n",
    "    end\n",
    "end\n",
    "\n",
    "tic()\n",
    "initialize!(k, 16, source, source, target, charge, N, N, 40, 10, f);\n",
    "toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### display the tree\n",
    "\n",
    "To display the tree and all interaction lists (UList, VList, WList, XList), use\n",
    "```\n",
    "display(tree, [index of nodes])\n",
    "```\n",
    "\n",
    "- black node is the central node\n",
    "- red is Ulist\n",
    "- blue is Wlist\n",
    "- green is Vlist\n",
    "- grey is Xlist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "display(k.tree, [89, 180])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FMM!\n",
    "simply call \n",
    "```\n",
    "p = FMM!(kernel)\n",
    "```\n",
    "\n",
    "to get the result ``p = K * q`` using ``KIFMM``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tic()\n",
    "potential = FMM!(k)\n",
    "toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### check result using pure interaction (without any matrix operations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tic()\n",
    "checkPotential = zeros(N)\n",
    "for i = 1:N\n",
    "    for j = 1:N\n",
    "        checkPotential[i] += k.kernel(k.tree.sourceTree[j,:], k.tree.normalTree[j,:], k.tree.targetTree[i,:]) * k.chargeTree[j]\n",
    "    end\n",
    "end\n",
    "toc()\n",
    "norm(potential - checkPotential)/norm(checkPotential)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}