BEM-kernelgen

BEM-kernelgen.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kernel transformation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Input: kernel function $K(r)$.  Output \"first integral\" of $K$: \n",
    "$$\n",
    "\\mathcal{K}_n(X) = \\int_0^1 w^n K(wX) dw\n",
    "$$\n",
    "for $X \\in (0,\\infty)$.   Some small set of $n$ values, to be determined, will be supplied at compile time.\n",
    "\n",
    "$K$ will be specified as a subtype of `AbstractKernel`.\n",
    "\n",
    "* For some kernel types, $\\mathcal{K}_n$ is known analytically.  (See Table I in [Homer's paper](http://math.mit.edu/~stevenj/papers/ReidWhJo14.pdf).)\n",
    "\n",
    "* $K$ may depend on some compile-time numeric parameters, specified as type parameters.\n",
    "\n",
    "* If the $\\mathcal{K}_n$ integral is *not* known analytically, we may have to compute it numerically.  This will be done at *compile time* by staged functions, with the numeric integration results used to compute the coefficients of a Chebyshev polynomial fit, which can then be compiled into an efficient polynomial approximation of $\\mathcal{K}_n$.   However, to compute a Chebyshev approximation for a function defined on $(0,\\infty)$, we will have to perform a coordinate transform from $X \\to (0,1)$, and the type of coordinate transformation will depend on how fast $\\mathcal{K}_n$ decays asymptotically as $X\\to\\infty$.  This decay rate can be *specified via the type* of `K`.\n",
    "\n",
    "The $\\mathcal{K}_n$ function will be parameterized by a `FirstIntegral{K,n}` type parameterized by an `AbstractKernel` type `K` and an integer `n`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "abstract AbstractKernel\n",
    "\n",
    "# any kernel ~ X^p for X ≪ s and ~ X^q for X ≫ s\n",
    "abstract PowerLawScaling{p,q,s} <: AbstractKernel\n",
    "\n",
    "type FirstIntegral{K<:AbstractKernel,n}; end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytically known integrals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "type PowerLaw{p} <: PowerLawScaling{p,p}; end # rᵖ power law"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#function Base.call{p,n}(::FirstIntegral{PowerLaw{p},n}, X::Number)\n",
    "#    return p >= 0 ? X^p / (1 + n + p) : inv(X^(-p) * (1 + n + p))\n",
    "#end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#F = FirstIntegral{PowerLaw{-1}, 3}()\n",
    "#F(3.7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#@code_llvm call(F, 3.7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerically computed integrals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a general multi-physics BEM package, one might conceivably have a user-specified kernel $K$ for which the first integral $\\mathcal{K}_n$ is *not* known analytically.   Performing the integral numerically at runtime would be too expensive, however, especially since the integrand may have an integrable singularity.\n",
    "\n",
    "Instead, we will perform the integral $\\mathcal{K}_n$ at *compile-time* (or, at least, outside the innermost BEM loops) for various $X$ and use these values to compute a *Chebyshev interpolating polynomial* $C(\\xi)$.  This polynomial will then be used to generate an efficient $\\mathcal{K}_n(X)$.\n",
    "\n",
    "There are three tricky points:\n",
    "\n",
    "* $\\mathcal{K}_n(X)$ will probably be singular as $X\\to 0$, which means we can't fit it directly to a polynomial.  (This is not a problem for how $\\mathcal{K}_n$ is *used*, because $\\mathcal{K}_n$ is always used for integration over domains that do not include $X=0$.)   This will be dealt with by requiring the user to specify the degree $p$ of the singularity as $X \\to 0$: i.e. $\\mathcal{K}_n(X) = O(X^p)$ for $X\\to 0$.  We will factor out this singularity from $\\mathcal{K}_n$ and fit the remaining (non-singular) function to a polynomial.\n",
    "\n",
    "* We need $\\mathcal{K}_n(X)$ for $X \\in (0,\\infty)$, whereas polynomial interpolation requires a finite interval, typically $(-1,1)$.   We will handle this by choosing a coordinate mapping $\\xi(X) \\in (-1,1)$.  Because such a coordinate mapping is necessarily singular, however, it will screw up the convergence of polynomial interpolation if we choose the wrong degree of singularity — we want a mapping such that $\\mathcal{K}_n(\\xi(X))$ is nonsingular (e.g. a low-degree polynomial in $\\xi$) as $X\\to\\infty$.   To choose this, the user will specify a degree $q$ of the decay rate as $X\\to\\infty$, i.e. $\\mathcal{K}_n(X) = O(X^q)$ for $X\\to \\infty$.\n",
    "\n",
    "* $X$ is dimensionful (it is a physical distance within one of the triangles or other geometric elements of the BEM basis).  Mapping it to a dimensionless $\\xi$ inevitably involves choosing a scale $s$ of $X$.   This $s$ should be user-specified (e.g. it can be the median diameter of the BEM elements).  That is, $\\mathcal{K}_n(X) \\sim X^p$ for $X \\ll s$ and $\\mathcal{K}_n(X) \\sim X^q$ for $X \\gg s$.\n",
    "\n",
    "In summary, the user will specify a `PowerLawScaling{p,q,s}` kernel type parameterized by `p`, `q`, and `s`.  They will also define a `call(::PowerLawScaling{p,q,s}, r::Number)` method that computes $K(r)$.\n",
    "\n",
    "The polynomial fit will be performed as follows, assuming $p \\le 0$ and $q \\le 0$.  First, we let $L_n(X) = \\mathcal{K}_n(X) / (s^p + X^p)$, which eliminates the $x\\to 0$ singularity while still having $L_n \\sim X^q$ for $X \\gg s$.   Second, we let $X = (1-\\xi)^{1/q} - 2^{1/q}$ [or equivalently $\\xi = 1 - (x+2^{1/q})^q$], which maps $\\xi \\in (-1,1)$ to $X \\in (0,\\infty)$, and has the property that $X^q \\approx 1-\\xi$ as $\\xi\\to 1$, so the coordinate remapping produces a nice degree-1 polynomial.   Finally, we fit $L_n(\\xi(X)) = C(\\xi)$ to a Chebyshev polynomial $C$, and compute $\\mathcal{K}_n(X)$ via $\\mathcal{K}_n(X) = (s^p + X^p) \\times L_n(\\xi(X))$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Chebyshev interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following routines compute the coefficients $c_n$ of a Chebyshev interpolating polynomial $C(x) = \\sum_{n=0}^{N-1} c_n T_n(x)$ for a function $f(x)$ on $(-1,1)$, where $T_n(x) = \\cos(n \\cos^{-1}x)$ are the first-kind Chebyshev polynomials.\n",
    "\n",
    "We compute these coefficients $c_n$ by first evaluating $f(x)$ at the Chebyshev points $\\cos\\left(\\pi\\frac{n+1/2}{N}\\right)$ for $n=0,\\ldots,N-1$, for which the Chebyshev sum is equivalent to a type-III discrete cosine transform (DCT-III), so that the coefficients $c_n$ are computed by a DCT-II.   These are *not* quite the typical Chebyshev points, which correspond to a DCT-I: the difference is that the DCT-I corresponds to the closed interval $[-1,1]$, i.e. it includes the endpoints, whereas our function may involve terms that blow up at the endpoints (althought the overall function should be okay) so we don't want to evaluate it there.\n",
    "\n",
    "We also provide a function `evalcheb` to evaluate $C(x)$ for any $x\\in(-1,1)$ by a Clenshaw recurrence, and a macro version `@evalcheb` (analogous to `Base.@horner`) that generates a completely inlined version of this recurrence for the case where $c$ is fixed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# N chebyshev points (order N) on the interval (-1,1)\n",
    "chebx(N) = [cos(π*(n+0.5)/N) for n in 0:N-1]\n",
    "\n",
    "# N chebyshev coefficients for vector of f(x) values on chebx points x\n",
    "function chebcoef(f::AbstractVector)\n",
    "    a = FFTW.r2r(f, FFTW.REDFT10) / length(f)\n",
    "    a[1] /= 2\n",
    "    return a\n",
    "end\n",
    "\n",
    "# given a function f and a tolerance, return enough Chebyshev coefficients to\n",
    "# reconstruct f to that tolerance on (-1,1)\n",
    "function chebcoef(f, tol=1e-13)\n",
    "    N = 10\n",
    "    local c\n",
    "    while true\n",
    "        x = chebx(N)\n",
    "        c = chebcoef(Float64[f(y) for y in x])\n",
    "        # look at last 3 coefs, since individual c's might be zero\n",
    "        if max(abs(c[end]),abs(c[end-1]),abs(c[end-2])) < tol * maxabs(c)\n",
    "            break\n",
    "        end\n",
    "        N *= 2\n",
    "    end\n",
    "    v₀ = maxabs(c) * tol\n",
    "    return c[1:findlast(v -> abs(v) > tol, c)] # shrink to minimum length\n",
    "end\n",
    "\n",
    "# given cheb coefficients a, evaluate them for x in (-1,1) by Clenshaw recurrence\n",
    "function evalcheb(x, a)\n",
    "    isempty(a) && throw(BoundsError())\n",
    "    -1 ≤ x ≤ 1 || throw(DomainError())\n",
    "    bₖ₊₁ = bₖ₊₂ = zero(x)\n",
    "    for k = length(a):-1:2\n",
    "        bₖ = a[k] + 2x*bₖ₊₁ - bₖ₊₂\n",
    "        bₖ₊₂ = bₖ₊₁\n",
    "        bₖ₊₁ = bₖ\n",
    "    end\n",
    "    return a[1] + x*bₖ₊₁ - bₖ₊₂\n",
    "end\n",
    "\n",
    "# inlined version of evalcheb given coefficents a, and x in (-1,1)\n",
    "macro evalcheb(x, a...)\n",
    "    isempty(a) && throw(BoundsError())\n",
    "    # Clenshaw recurrence, evaluated symbolically:\n",
    "    bₖ₊₁ = bₖ₊₂ = 0\n",
    "    for k = length(a):-1:2\n",
    "        bₖ = esc(a[k])\n",
    "        if bₖ₊₁ != 0\n",
    "            bₖ = :(muladd(t2, $bₖ₊₁, $bₖ))\n",
    "        end\n",
    "        if bₖ₊₂ != 0\n",
    "            bₖ = :($bₖ - $bₖ₊₂)\n",
    "        end\n",
    "        bₖ₊₂ = bₖ₊₁\n",
    "        bₖ₊₁ = bₖ\n",
    "    end\n",
    "    ex = esc(a[1])\n",
    "    if bₖ₊₁ != 0\n",
    "        ex = :(muladd(t, $bₖ₊₁, $ex))\n",
    "    end\n",
    "    if bₖ₊₂ != 0\n",
    "        ex = :($ex - $bₖ₊₂)\n",
    "    end\n",
    "    Expr(:block, :(t = $(esc(x))), :(t2 = 2t), ex)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try a simple test case: performing Chebyshev interpolation of $\\exp(x)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "c = chebcoef(exp)\n",
    "x = linspace(-1,1,100)\n",
    "maximum(abs(Float64[evalcheb(y,c) for y in x] - exp(x))) # the maximum error on [-1,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# check that the evalcheb macro works\n",
    "evalcheb(0.1234, c[1:4]) - @evalcheb(0.1234, c[1],c[2],c[3],c[4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First-integral generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# extract parameters from PowerLawScaling type\n",
    "pqsPowerLawScaling{p,q,s}(::PowerLawScaling{p,q,s}) = (p,q,s)\n",
    "\n",
    "@generated function Base.call{P<:PowerLawScaling,n}(::FirstIntegral{P,n}, X::Real)\n",
    "    # compute the Chebyshev coefficients (of the rescaled 𝒦ₙ as described above)\n",
    "    K = P()\n",
    "    p,q,s = pqsPowerLawScaling(K)\n",
    "    \n",
    "    𝒦ₙ = X -> quadgk(w -> w^n * K(w*X), 0,1, abstol=1e-12, reltol=1e-10)[1]\n",
    "    Lₙ = p < 0 ? X -> 𝒦ₙ(X) / (s^p + X^p) : 𝒦ₙ # scale out X ≪ s singularity\n",
    "    q > 0 && throw(DomainError()) # don't know how to deal with growing kernels\n",
    "    qinv = 1/q\n",
    "    c = chebcoef(ξ -> Lₙ((1-ξ)^qinv - 2^qinv), 1e-9)\n",
    "    \n",
    "    # return an expression that inlines the evaluation of 𝒦ₙ via C(ξ)\n",
    "    quote\n",
    "        X <= 0 && throw(DomainError())\n",
    "        ξ = 1 - (X + $(2^qinv))^$q\n",
    "        C = @evalcheb ξ $(c...)\n",
    "        return $p < 0 ? C * (X^$p + $(s^p)) : C\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# A simple example where the result is known analytically:\n",
    "\n",
    "type DumbPowerLaw{p,s} <: PowerLawScaling{p,p,s}; end # rᵖ power law\n",
    "Base.call{p}(::DumbPowerLaw{p}, r) = r^p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "F = FirstIntegral{DumbPowerLaw{-1,1.0}, 3}()\n",
    "F(3.7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@code_llvm call(F,3.7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using PyPlot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = [0.01:.01:1.0;]; plot(x, map(FirstIntegral{DumbPowerLaw{-1,1.}, 3}(),x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.0-dev",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}