demo_singular-poisson-rst.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Singular Poisson\n",
    "================\n",
    "\n",
    "This demo is implemented in a single Python file,\n",
    "demo\\_singular-poisson.py, which contains both the variational forms and\n",
    "the solver.\n",
    "\n",
    "This demo illustrates how to:\n",
    "\n",
    "-   Solve a linear partial differential equation\n",
    "-   Apply non-zero Neumann boundary conditions\n",
    "-   Define Expressions\n",
    "-   Define a FunctionSpace\n",
    "-   Use the Krylov solver\n",
    "-   Solve singular problems\n",
    "\n",
    "The solution for $u$ in this demo will look as follows:\n",
    "\n",
    "![image](singular_u.png)\n",
    "\n",
    "Equation and problem definition\n",
    "-------------------------------\n",
    "\n",
    "The Poisson equation is the canonical elliptic partial differential\n",
    "equation. For a domain $\\Omega \\in \\mathbb{R}^n$ with boundary\n",
    "$\\Gamma = \\partial \\Omega$, the Poisson equation with pure Neumann\n",
    "boundary conditions reads:\n",
    "\n",
    "$$-\\nabla^{2} u &= f \\quad {\\rm in} \\ \\Omega \\\\\n",
    "\\nabla u \\cdot n &= g \\quad {\\rm on} \\ \\Gamma \\\\$$\n",
    "\n",
    "Since only Neumann conditions are applied, $u$ is only determined up to\n",
    "a constant by the above equations. An addition constraint is thus\n",
    "required, for instance\n",
    "\n",
    "$$\\int u \\, {\\rm d} x = 0$$\n",
    "\n",
    "The most standard variational form of the Poisson equation reads: find\n",
    "$u \\in V$ such that\n",
    "\n",
    "$$a(u, v) = L(v) \\quad \\forall \\ v \\in V,$$\n",
    "\n",
    "where $V$ is a suitable function space and\n",
    "\n",
    "$$a(u, v) &= \\int_{\\Omega} \\nabla u \\cdot \\nabla v \\, {\\rm d} x, \\\\\n",
    "L(v)    &= \\int_{\\Omega} f v \\, {\\rm d} x\n",
    "+ \\int_{\\Gamma} g v \\, {\\rm d} s.$$\n",
    "\n",
    "The expression $a(u, v)$ is the bilinear form and $L(v)$ is the linear\n",
    "form.\n",
    "\n",
    "If we make the Ansatz that $u$ can be expressed as a linear combination\n",
    "of the basis functions of $V$, and discretize the equation, we can write\n",
    "our problem as a linear system:\n",
    "\n",
    "$$AU = b,$$\n",
    "\n",
    "where $U$ gives the coefficient for the basis functions expressing $u$.\n",
    "\n",
    "Since we have pure Neumann boundary conditions, the matrix $A$ is\n",
    "singular. There exists a non-trival vector $e$ such that\n",
    "\n",
    "$$Ae=0.$$\n",
    "\n",
    "span $\\{ e \\}$ is the null space of A. Consequently, the matrix $A$ is\n",
    "rank deficient and the right-hand side vector $b$ may fail to be in the\n",
    "column space of $A$. We therefore need to remove the components of $b$\n",
    "that do not lie in the column space to make the system solvable.\n",
    "\n",
    "In this demo, we shall consider the following definitions of the input\n",
    "functions, the domain, and the boundaries:\n",
    "\n",
    "-   $\\Omega = [0,1] \\times [0,1]$ (a unit square)\n",
    "-   $\\Gamma = \\partial \\Omega$ (boundary)\n",
    "-   $g = -\\sin(5x)$ (normal derivative)\n",
    "-   $f = 10\\exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02)$ (source term)\n",
    "\n",
    "Implementation\n",
    "--------------\n",
    "\n",
    "This description goes through the implementation (in\n",
    "demo\\_singular-poisson.py) of a solver for the above described Poisson\n",
    "equation step-by-step.\n",
    "\n",
    "First, the :pydolfin module is imported: :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specifying matplotlib as backend:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "parameters[\"plotting_backend\"]=\"matplotlib\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we check that dolfin is configured with the backend called PETSc,\n",
    "since it provides us with a wide range of methods used by\n",
    ":pyKrylovSolver &lt;dolfin.cpp.la.KrylovSolver&gt;. We set PETSc as our\n",
    "backend for linear algebra. :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test for PETSc\n",
    "if not has_linear_algebra_backend(\"PETSc\"):\n",
    "    info(\"DOLFIN has not been configured with PETSc. Exiting.\")\n",
    "    exit()\n",
    "\n",
    "parameters[\"linear_algebra_backend\"] = \"PETSc\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin by defining a mesh of the domain and a finite element function\n",
    "space $V$ relative to this mesh. We use a built-in mesh provided by the\n",
    "class :pyUnitSquareMesh\n",
    "&lt;dolfin.cpp.mesh.UnitSquareMesh&gt;. In order to create a mesh\n",
    "consisting of $64 \\times 64$ squares with each square divided into two\n",
    "triangles, we do as follows: :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create mesh and define function space\n",
    "mesh = UnitSquareMesh(64, 64)\n",
    "V = FunctionSpace(mesh, \"CG\", 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we need to specify the trial functions (the unknowns) and the test\n",
    "functions on the space $V$. This can be done using a\n",
    ":pyTrialFunction &lt;dolfin.functions.function.TrialFunction&gt; and a\n",
    ":pyTestFunction\n",
    "&lt;dolfin.functions.function.TrialFunction&gt; as follows: :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = TrialFunction(V)\n",
    "v = TestFunction(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Further, the source $f$ and the boundary normal derivative $g$ are\n",
    "involved in the variational forms, and hence we must specify these. Both\n",
    "$f$ and $g$ are given by simple mathematical formulas, and can be easily\n",
    "declared using the\n",
    ":pyExpression &lt;dolfin.functions.expression.Expression&gt; class. Note\n",
    "that the strings defining f and g use C++ syntax since, for efficiency,\n",
    "DOLFIN will generate and compile C++ code for these expressions at\n",
    "run-time. :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = Expression(\"10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)\", degree=2)\n",
    "g = Expression(\"-sin(5*x[0])\", degree=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With $u,v,f$ and $g$, we can write down the bilinear form $a$ and the\n",
    "linear form $L$ (using UFL operators). :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = inner(grad(u), grad(v))*dx\n",
    "L = f*v*dx + g*v*ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to transform our variational problem into a linear system we\n",
    "need to assemble the coefficient matrix `A` and the right-side vector\n",
    "`b`. We do this using the function :pyassemble\n",
    "&lt;dolfin.cpp.fem.Assembler.assemble&gt;: :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assemble system\n",
    "A = assemble(a)\n",
    "b = assemble(L)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We specify a Vector for storing the result by defining a\n",
    ":pyFunction &lt;dolfin.cpp.function.Function&gt;. :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution Function\n",
    "u = Function(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we specify the iterative solver we want to use, in this case a\n",
    ":pyPETScKrylovSolver &lt;dolfin.cpp.la.PETScKrylovSolver&gt; with the\n",
    "conjugate gradient (CG) method, and attach the matrix operator to the\n",
    "solver. :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Krylov solver\n",
    "solver = PETScKrylovSolver(\"cg\")\n",
    "solver.set_operator(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We impose our additional constraint by removing the null space component\n",
    "from the solution vector. In order to do this we need a basis for the\n",
    "null space. This is done by creating a vector that spans the null space,\n",
    "and then defining a basis from it. The basis is then attached to the\n",
    "matrix `A` as its null space. :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create vector that spans the null space and normalize\n",
    "null_vec = Vector(u.vector())\n",
    "V.dofmap().set(null_vec, 1.0)\n",
    "null_vec *= 1.0/null_vec.norm(\"l2\")\n",
    "\n",
    "# Create null space basis object and attach to PETSc matrix\n",
    "null_space = VectorSpaceBasis([null_vec])\n",
    "as_backend_type(A).set_nullspace(null_space)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Orthogonalization of `b` with respect to the null space makes sure that\n",
    "it doesn't contain any component in the null space. :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "null_space.orthogonalize(b);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we are able to solve our linear system :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solver.solve(u.vector(), b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and plot the solution :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plot(u, interactive=True)"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 1
}