demo_nonlinear-poisson.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nonlinear Poisson equation\n",
    "==========================\n",
    "\n",
    "This demo is implemented in a single Python file,\n",
    "demo\\_nonlinear-poisson.py, which contains both the variational form and\n",
    "the solver.\n",
    "\n",
    "This demo illustrates how to:\n",
    "\n",
    "-   Solve a nonlinear partial differential equation (in this case a\n",
    "    nonlinear variant of Poisson's equation)\n",
    "-   Create and apply Dirichlet boundary conditions\n",
    "-   Define an :pyExpression &lt;dolfin.cpp.function.Expression&gt;\n",
    "-   Define a :pyFunctionSpace &lt;dolfin.cpp.function.FunctionSpace&gt;\n",
    "-   Create a :pySubDomain &lt;dolfin.cpp.mesh.SubDomain&gt;\n",
    "\n",
    "The solution for $u$ in this demo will look as follows:\n",
    "\n",
    "![image](plot_u.png)\n",
    "\n",
    "and the gradient of $u$ will look like this:\n",
    "\n",
    "![image](plot_u_gradient.png)\n",
    "\n",
    "Equation and problem definition\n",
    "-------------------------------\n",
    "\n",
    "For a domain $\\Omega \\subset \\mathbb{R}^N$ with boundary\n",
    "$\\partial \\Omega = \\Gamma_{D} \\cup \\Gamma_{N}$, we consider the\n",
    "following nonlinear Poisson equation with given boundary conditions:\n",
    "\n",
    "$$n$$$$- \\nabla\\cdot((1 + u^2) \\nabla u) &= f \\quad {\\rm in}\\, \\Omega,\\\\\n",
    "u &= 1  \\quad  {\\rm on}\\, \\Gamma_D,\\\\\n",
    "\\nabla u\\cdot n &= 0 \\quad  {\\rm on}\\, \\Gamma_N.$$\n",
    "\n",
    "Here $f$ is input data and $n$ denotes the outward directed boundary\n",
    "normal. The nonlinear variational form can be written in the following\n",
    "canonical form: find $u \\in V$ such that\n",
    "\n",
    "$$F(u;v) = 0 \\quad \\forall \\, v \\in \\hat{V}$$\n",
    "\n",
    "Here $F:V\\times\\hat{V}\\rightarrow\\mathbb{R}$ is a semilinear form,\n",
    "linear in the argument subsequent to the semicolon, and $V$ is some\n",
    "suitable function space. The semilinear form is defined as follows:\n",
    "\n",
    "$$F(u;v) = \\int_\\Omega (1 + u^2)\\cdot\\nabla u \\cdot \\nabla v - f v \\,{\\rm dx} = 0.$$\n",
    "\n",
    "To solve the nonlinear system $b(U) = 0$ by Newton's method we compute\n",
    "the Jacobian $A = b'$, where $U$ is the coefficients of the linear\n",
    "combination in the finite element solution\n",
    "$u_h = \\sum_{j=1}^{N}U_j\\phi_j, \\;\n",
    "b:\\mathbb{R}^N\\rightarrow\\mathbb{R}^N$ and\n",
    "\n",
    "$$b_i(U) = F(u_h;\\hat{\\phi}_i),\\quad i = 1,2,\\dotsc,N.$$\n",
    "\n",
    "Linearizing the semilinear form $F$ around $u = u_h$, we obtain\n",
    "\n",
    "$$F'(u_h;\\delta u,v) = \\int_\\Omega [(2 \\delta u\\nabla u_h)\\cdot\\nabla v + ((1+u_h^2)\\nabla\\delta u)\\nabla v] \\,{\\rm dx}$$\n",
    "\n",
    "We note that for each fixed $u_h$, $a =\n",
    "F'(u_h;\\,\\cdot\\,,\\,\\cdot\\,)$ is a bilinear form and $L =\n",
    "F(u_h;\\,\\cdot\\,,\\,\\cdot\\,)$ is a linear form. In each Newton iteration,\n",
    "we thus solve a linear variational problem of the canonical form: find\n",
    "$\\delta u \\in V_{h,0}$ such that\n",
    "\n",
    "$$F'(u_h;\\delta u,v) = -F(u_h;v)\\quad\\forall\\,v\\in\\hat{V}_h.$$\n",
    "\n",
    "In this demo, we shall consider the following definitions of the input\n",
    "function, the domain, and the boundaries:\n",
    "\n",
    "-   $\\Omega = [0,1] \\times [0,1]\\,\\,\\,$ (a unit square)\n",
    "-   $\\Gamma_{D} = \\{(1, y) \\subset \\partial \\Omega\\}\\,\\,\\,$\n",
    "    (Dirichlet boundary)\n",
    "-   $\\Gamma_{N} = \\{(x, 0) \\cup (x, 1) \\cup (0, y) \\subset \\partial \\Omega\\}\\,\\,\\,$\n",
    "    (Neumann boundary)\n",
    "-   $f(x, y) = x\\sin(y)\\,\\,\\,$ (source term)\n",
    "\n",
    "Implementation\n",
    "--------------\n",
    "\n",
    "This description goes through the implementation (in\n",
    "demo\\_nonlinear-poisson.py) of a solver for the above described\n",
    "nonlinear Poisson 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": [
    "Next, we want to consider the Dirichlet boundary condition. A simple\n",
    "Python function, returning a boolean, can be used to define the\n",
    "subdomain for the Dirichlet boundary condition ($\\Gamma_D$). The\n",
    "function should return True for those points inside the subdomain and\n",
    "False for the points outside. In our case, we want to say that the\n",
    "points $(x, y)$ such that $x = 1$ are inside on the inside of\n",
    "$\\Gamma_D$. (Note that because of rounding-off errors, it is often wise\n",
    "to instead specify $|x - 1| < \\epsilon$, where $\\epsilon$ is a small\n",
    "number (such as machine precision).):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sub domain for Dirichlet boundary condition\n",
    "class DirichletBoundary(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then define a mesh of the domain and a finite element function space\n",
    "V relative to this mesh. We use the built-in mesh provided by the class\n",
    ":pyUnitSquareMesh\n",
    "&lt;dolfin.cpp.mesh.UnitSquareMesh&gt;. In order to create a mesh\n",
    "consisting of $32 \\times 32$ 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(32, 32)\n",
    "File(\"mesh.pvd\") << mesh\n",
    "\n",
    "V = FunctionSpace(mesh, \"CG\", 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second argument to :pyFunctionSpace\n",
    "&lt;dolfin.cpp.function.FunctionSpace&gt; is the finite element family,\n",
    "while the third argument specifies the polynomial degree. Thus, in this\n",
    "case, we use 'CG', for Continuous Galerkin, as a synonym for 'Lagrange'.\n",
    "With degree 1, we simply get the standard linear Lagrange element, which\n",
    "is a triangle with nodes at the three vertices (or in other words,\n",
    "continuous piecewise linear polynomials).\n",
    "\n",
    "The Dirichlet boundary condition can be created using the class\n",
    ":pyDirichletBC &lt;dolfin.cpp.fem.DirichletBC&gt;. A\n",
    ":pyDirichletBC &lt;dolfin.cpp.fem.DirichletBC&gt; takes three arguments:\n",
    "the function space the boundary condition applies to, the value of the\n",
    "boundary condition, and the part of the boundary on which the condition\n",
    "applies. In our example, the function space is V, the value of the\n",
    "boundary condition (1.0) can be represented using a Constant and the\n",
    "Dirichlet boundary is defined above. The definition of the Dirichlet\n",
    "boundary condition then looks as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define boundary condition\n",
    "g = Constant(1.0)\n",
    "bc = DirichletBC(V, g, DirichletBoundary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we want to express the variational problem. First, we need to\n",
    "specify the function u which represents the solution. Upon\n",
    "initialization, it is simply set to the zero function, which will\n",
    "represent the initial guess $u_0$. A :pyFunction\n",
    "&lt;dolfin.cpp.function.Function&gt; represents a function living in a\n",
    "finite element function space. The test function $v$ is specified, also\n",
    "living in the function space $V$. We do this by defining a\n",
    ":pyFunction &lt;dolfin.cpp.function.Function&gt; and a\n",
    ":pyTestFunction &lt;dolfin.functions.function.TestFunction&gt; on the\n",
    "previously defined :pyFunctionSpace\n",
    "&lt;dolfin.cpp.function.FunctionSpace&gt; V.\n",
    "\n",
    "Further, the source $f$ is involved in the variational forms, and hence\n",
    "we must specify this. We have $f$ given by a simple mathematical\n",
    "formula, which can be easily declared using the\n",
    ":pyExpression &lt;dolfin.cpp.function.Expression&gt; class. Note that\n",
    "the strings defining f use C++ syntax since, for efficiency, DOLFIN will\n",
    "generate and compile C++ code for this expression at run-time.\n",
    "\n",
    "By defining the function in this step and omitting the trial function we\n",
    "tell FEniCS that the problem is nonlinear. With these ingredients, we\n",
    "can write down the semilinear form F (using UFL operators). In summary,\n",
    "this reads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define variational problem\n",
    "u = Function(V)\n",
    "v = TestFunction(V)\n",
    "f = Expression(\"x[0]*sin(x[1])\", degree=2)\n",
    "F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we have specified the variational forms and can consider the\n",
    "solution of the variational problem. Next, we can call the solve\n",
    "function with the arguments F == 0, u, bc and solver parameters as\n",
    "follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute solution\n",
    "solve(F == 0, u, bc,\n",
    "      solver_parameters={\"newton_solver\":{\"relative_tolerance\":1e-6}})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Newton procedure is considered to have converged when the residual\n",
    "$r_n$ at iteration $n$ is less than the absolute tolerance or the\n",
    "relative residual $\\frac{r_n}{r_0}$ is less than the relative tolerance.\n",
    "\n",
    "A :pyFunction &lt;dolfin.cpp.function.Function&gt; can be manipulated in\n",
    "various ways, in particular, it can be plotted and saved to file. Here,\n",
    "we output the solution to a VTK file (using the suffix .pvd) for later\n",
    "visualization and also plot it using the plot command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot solution and solution gradient\n",
    "plt.figure()\n",
    "plot(u, title=\"Solution\")\n",
    "plt.figure()\n",
    "plot(grad(u), title=\"Solution gradient\")\n",
    "interactive()\n",
    "\n",
    "# Save solution in VTK format\n",
    "file = File(\"nonlinear_poisson.pvd\")\n",
    "file << u"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 1
}