biharmonic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Warning!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The current (Sep 2017) implementation of Hermite elements for $d \\ge 2$ is bogus (see [Poisson2D](Poisson2D.ipynb) too): The Hermite trafo is applied to the dof at the barycenter, which is just an evaluation dof. This means that we try to access derivatives at this node which don't exist, resulting in access beyond the end of the array of basis values at integration points in the generated C++ code in `tabulate_tensor()`. We need to modify the intermediate representation of the terms in the integral to account for this dof, by looping first through the ones at the vertices then considering the last one without multiplication by the $H$ matrix, etc. An example value for `trans_integrals` for a form `v*dx` follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "trans_integrals = [(6, {(('j', 0, 9),): [[{'H', 'det'}, {6}, {'FE0'}, set(), {}],\n",
    "                             [('j', '((FE0[ip][3*(j/3)+0])*(H[3 * (j % 3) + 0]) + ' +\n",
    "                                    '(FE0[ip][3*(j/3)+1])*(H[3 * (j % 3) + 1]) + ' +\n",
    "                                    '(FE0[ip][3*(j/3)+2])*(H[3 * (j % 3) + 2]))*W6[ip]*det', 7)]]},\n",
    "                   {}, {}, None, {}),\n",
    "                  (6, {(('j', 9, 10),): [[{'H', 'det'}, {6}, {'FE0'}, set(), {}],\n",
    "                             [('j', '(FE0[ip][3*(j/3)+0])*W6[ip]*det', 7)]]},\n",
    "                   {}, {}, None, {})]\n",
    "#points, terms, functions, ip_consts, coordinate, conditionals = trans_integrals[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "\n",
    "parameters['form_compiler']['representation'] = 'quadrature'\n",
    "# Optimization options for the form compiler\n",
    "# MBD: Disable these until I've finished optimizedquadraturetransformer.py\n",
    "parameters[\"form_compiler\"][\"cpp_optimize\"] = False\n",
    "parameters[\"form_compiler\"][\"optimize\"] = False\n",
    "\n",
    "\n",
    "# Make mesh ghosted for evaluation of DG terms\n",
    "parameters[\"ghost_mode\"] = \"shared_facet\"\n",
    "def domain(nx, ny):\n",
    "    return RectangleMesh(Point(0,-pi/2), Point(pi, pi/2), nx, ny, 'crossed')\n",
    "\n",
    "W = FunctionSpace(domain(10,10), \"Hermite\", 3)\n",
    "p = TrialFunction(W)\n",
    "q = TestFunction(W)\n",
    "a = inner(nabla_grad(p), nabla_grad(q))*dx\n",
    "L = 10*q*dx\n",
    "u = Function(W)\n",
    "solve(a==L, u)\n",
    "#b = assemble(L)\n",
    "plot(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook solves the Biharmonic equation,\n",
    "\n",
    "$$\\nabla^4 u(x, y) = f(x, y)$$\n",
    "\n",
    "on the square $[0, \\pi] \\times [-\\pi/2,\\pi/2]$ with source f given by\n",
    "\n",
    "$$f(x, y) = 4 \\sin(x) \\cos(y)$$\n",
    "\n",
    "and boundary conditions given by\n",
    "\n",
    "$$u(x, y)         = 0$$\n",
    "$$\\nabla^2 u(x, y) = 0$$\n",
    "\n",
    "using a discontinuous Galerkin formulation (interior penalty method). The analytic solution is\n",
    "\n",
    "$$u(x,y) = \\sin(x) \\cos(y).$$\n",
    "\n",
    "**TODO:** I should check whether the formulation makes sense for Hermite elements and what guarantees they provide."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "%matplotlib inline\n",
    "import nbimporter\n",
    "from boundary import apply_dirichlet_hermite\n",
    "from dofs import list_hermite_dofs\n",
    "from utils import ExpressionAD, Msg\n",
    "from interpolation import interpolate_hermite\n",
    "import matplotlib.pyplot as pl\n",
    "import numpy as np\n",
    "\n",
    "parameters['form_compiler']['representation'] = 'quadrature'\n",
    "# Optimization options for the form compiler\n",
    "# MBD: Disable these until I've finished optimizedquadraturetransformer.py\n",
    "parameters[\"form_compiler\"][\"cpp_optimize\"] = False\n",
    "parameters[\"form_compiler\"][\"optimize\"] = False\n",
    "\n",
    "\n",
    "# Make mesh ghosted for evaluation of DG terms\n",
    "parameters[\"ghost_mode\"] = \"shared_facet\"\n",
    "\n",
    "\n",
    "def sol(x, y):\n",
    "    \"\"\" Analytic solution. \"\"\"\n",
    "    import autograd.numpy as np\n",
    "    return np.sin(x)*np.cos(y)\n",
    "\n",
    "def biharmonic(V:FunctionSpace) -> Function:\n",
    "    \"\"\" Solves the biharmonic equation with a non-conforming discretisation.\n",
    "    Accepts either Lagrange or Hermite FunctionSpaces.\n",
    "\n",
    "    This code is based on DOLFIN's homonymous demo (C) 2009 Kristian B. Oelgaard\n",
    "    \"\"\"\n",
    "    class DirichletBoundary(SubDomain):\n",
    "        def inside(self, x, on_boundary):\n",
    "            return on_boundary\n",
    "\n",
    "    u0 = project(Constant(0.0), V)  # MBD: Need this for Hermite (interpolation doesn't work)\n",
    "    bc = DirichletBC(V, u0, DirichletBoundary())\n",
    "\n",
    "    u = TrialFunction(V)\n",
    "    v = TestFunction(V)\n",
    "\n",
    "    # Define normal component, mesh size and right-hand side\n",
    "    h = CellSize(V.mesh())\n",
    "    h_avg = (h('+') + h('-'))/2.0\n",
    "    n = FacetNormal(V.mesh())\n",
    "    \n",
    "    def source(x,y):\n",
    "        #import autograd.numpy as np\n",
    "        #return 4.0*np.pi**4*np.sin(np.pi*x)*np.sin(np.pi*y)\n",
    "        return 4*sol(x,y)\n",
    "    f = ExpressionAD(fun=source, degree=3)\n",
    "\n",
    "    # Penalty parameter\n",
    "    alpha = Constant(6.0)\n",
    "\n",
    "    # Define bilinear form\n",
    "    a = inner(div(grad(u)), div(grad(v)))*dx \\\n",
    "        - inner(avg(div(grad(u))), jump(grad(v), n))*dS \\\n",
    "        - inner(jump(grad(u), n), avg(div(grad(v))))*dS \\\n",
    "        + alpha/h_avg*inner(jump(grad(u),n), jump(grad(v),n))*dS\n",
    "    \n",
    "    # Solve variational problem\n",
    "    u = Function(V)\n",
    "    if V.ufl_element().family().lower() == 'hermite':\n",
    "        #L = interpolate_hermite(f, V)*v*dx  # nodal interpolation is much worse than projection\n",
    "        with Msg(\"Assembling\"):\n",
    "            L = project(f, V)*v*dx\n",
    "            A = assemble(a)\n",
    "            b = assemble(L)\n",
    "        with Msg(\"Applying Hermite BCs\"):\n",
    "            apply_dirichlet_hermite(A, b, bc)\n",
    "        with Msg(\"Solving\"):\n",
    "            solve(A, u.vector(), b)\n",
    "    else:\n",
    "        #%debug -b /home/fenics/local/lib/python2.7/site-packages/ffc/quadrature/quadraturetransformerbase.py:347 solve(a == L, u, bc)\n",
    "        L = f*v*dx\n",
    "        with Msg(\"Solving\"):\n",
    "            solve(a == L, u, bc)\n",
    "    \n",
    "    return u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def domain(nx, ny):\n",
    "    return RectangleMesh(Point(0,-pi/2), Point(pi, pi/2), nx, ny, 'crossed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = FunctionSpace(domain(10,10), \"Lagrange\", 3)\n",
    "solution = project(ExpressionAD(fun=sol, degree=3), W)\n",
    "lag = biharmonic(W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "V = FunctionSpace(W.mesh(), \"Hermite\", 3)\n",
    "her = biharmonic(V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(14,4))\n",
    "pl.subplot(1,3,1)\n",
    "plot(lag, title=\"Lagrange elements\", cmap='hot')\n",
    "pl.subplot(1,3,2)\n",
    "plot(her, title=\"Hermite elements\", cmap='hot')\n",
    "pl.subplot(1,3,3)\n",
    "plot(solution, title=\"Analytic solution\", cmap='hot')\n",
    "#pl.savefig(\"img/biharmonic-solutions.eps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute now the pointwise relative difference, defined as the quotient \n",
    "\n",
    "$$\\frac{u_h - u_l}{||u||_\\infty}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# We evaluate the solutions only at the vertices of the mesh to avoid any weird effects\n",
    "cc = W.mesh().coordinates()\n",
    "xx = np.array(sorted(list(set(cc[:,0]))))\n",
    "yy = np.array(sorted(list(set(cc[:,1]))))\n",
    "\n",
    "infnorm = sol(xx,yy).max()\n",
    "\n",
    "pl.figure(figsize=(16,6))\n",
    "pl.subplot(1,2,1)\n",
    "for y in yy[::7]:\n",
    "    pl.plot(xx, [np.abs((lag(x,y) - sol(x,y)))/infnorm for x in xx],\n",
    "            label='@%.1f' % y)\n",
    "pl.title(\"Relative error at several ordinates (Lagrange)\")\n",
    "_ = pl.legend()\n",
    "pl.subplot(1,2,2)\n",
    "for y in yy[::7]:\n",
    "    pl.plot(xx, [np.abs((her(x,y) - sol(x,y)))/infnorm for x in xx],\n",
    "            label='@%.1f' % y)\n",
    "pl.title(\"Relative error at several ordinates (Hermite)\")\n",
    "_ = pl.legend()\n",
    "\n",
    "pl.savefig(\"img/biharmonic-approx-error.eps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The difference between both approximations is strange:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(12,8))\n",
    "for y in yy[::7]:\n",
    "    pl.plot(xx, [abs(lag(x, y) - her(x, y)) for x in xx], label='@%.1f' % y)\n",
    "pl.title(\"$u_h - u_l$\")\n",
    "_ = pl.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The maximal difference is below 0.1%:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def relative_error(lag, her):\n",
    "    lv = lag.vector().array()\n",
    "    hv = her.vector().array()\n",
    "    lagrange_dofs = list(set(range(V.dim())) - set(list_hermite_dofs(V)))\n",
    "    lmin, lmax, hmin, hmax = lv.min(), lv.max(), hv[lagrange_dofs].min(), hv[lagrange_dofs].max()\n",
    "    return 100*abs(hmax - lmax)/lmax\n",
    "\n",
    "relative_error(lag, her)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "errors = []\n",
    "numdofs_her = []\n",
    "numdofs_lag = []\n",
    "times_lag = []\n",
    "times_her = []\n",
    "sizes = range(5, 35, 5)\n",
    "dolfin.DEBUG = Msg.output_level - 1  # Disable Msg's default output\n",
    "for s in sizes:\n",
    "    W = FunctionSpace(domain(s,s), \"Lagrange\", 3)\n",
    "    with Msg(\"Computing Lagrange solution with s=%d\" % s, level=0):\n",
    "        lag = biharmonic(W)\n",
    "    times_lag.append(Msg.last)\n",
    "    numdofs_lag.append(W.dim())\n",
    "    \n",
    "    V = FunctionSpace(W.mesh(), \"Hermite\", 3)\n",
    "    with Msg(\"Computing Hermite solution with s=%d\" % s, level=0):\n",
    "        her = biharmonic(V)\n",
    "    times_her.append(Msg.last)\n",
    "    numdofs_her.append(V.dim())\n",
    "\n",
    "    errors.append(relative_error(lag, her))\n",
    "\n",
    "numdofs_lag = np.array(numdofs_lag)\n",
    "numdofs_her = np.array(numdofs_her)\n",
    "times_lag = np.array(times_lag)\n",
    "times_her = np.array(times_her)\n",
    "errors = np.array(errors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Degree 3, sizes: 10,15,20,25,30,35\n",
    "```python\n",
    ">>> errors, numdofs_lag, numdofs_her\n",
    "\n",
    "   (array([ 0.63554695,  0.44286782,  0.35878701,  0.29271231,  0.25026288,\n",
    "         0.21675917]),\n",
    "    array([ 1861,  4141,  7321, 11401, 16381, 22261]),\n",
    "    array([ 1063,  2343,  4123,  6403,  9183, 12463]))\n",
    "```\n",
    "and \n",
    "```python\n",
    "# Manually copied\n",
    "times_lag = [0.221, 0.453, 0.800, 1.253, 1.809, 2.483]\n",
    "times_her = [0.858, 12.935, 70.857, 240.694, 1012.162, 3414.528]\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(10,5))\n",
    "pl.plot(sizes, errors/100, label='diff')\n",
    "_ = pl.title(\"Relative maximal difference\")\n",
    "#pl.savefig(\"img/biharmonic-diffs.eps\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(20,8))\n",
    "pl.subplot(1,2,1)\n",
    "pl.tight_layout(w_pad=2.5, rect=(0.05, 0.05, 0.95, 0.95))\n",
    "pl.plot(numdofs_her, times_her, label='her')\n",
    "pl.ylabel('seconds')\n",
    "pl.xlabel('dofs')\n",
    "_ = pl.title(\"Hermite\")\n",
    "pl.subplot(1,2,2)\n",
    "pl.plot(numdofs_lag, times_lag, label='lag')\n",
    "pl.xlabel('dofs')\n",
    "_ = pl.title(\"Lagrange\")\n",
    "#pl.savefig(\"img/biharmonic-times.eps\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "(max(times_lag) - min(times_lag)) / (numdofs_lag.max() - numdofs_lag.min())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following plot is weird, but this is likely to be because of dolfin being confused by Hermite dofs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "z = project(lag-her, V)\n",
    "print(norm(z))\n",
    "_ = plot(z, title=\"Difference\", cmap='bone')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the solutions to files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "File(\"biharmonic-hermite.pvd\") << her\n",
    "File(\"biharmonic-lagrange.pvd\") << lag"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.5.2"
  },
  "toc": {
   "nav_menu": {
    "height": "70px",
    "width": "254px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}