Labs/Lab2_Tuesday/Poisson2D.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Poisson Equation in 2D\n",
    "\n",
    "In this example we solve the Poisson equation in two space dimensions.\n",
    "\n",
    "For a domain $\\Omega \\subset \\mathbb{R}^2$ with boundary $\\partial \\Omega = \\Gamma_D \\cup \\Gamma_N$, we write the boundary value problem (BVP):\n",
    "\n",
    "$$ \n",
    "\\left\\{\n",
    "\\begin{array}{ll}\n",
    "- \\Delta u = f  & \\text{in} \\; \\Omega, \\\\\n",
    "u = u_D & \\text{on} \\; \\Gamma_D, \\\\\n",
    "\\nabla u \\cdot \\boldsymbol{n} = g & \\text{on} \\; \\Gamma_N.\n",
    "\\end{array}\n",
    "\\right.$$\n",
    "\n",
    "Here, $\\Gamma_D \\subset \\Omega$ denotes the part of the boundary where we prescribe Dirichlet boundary conditions, and $\\Gamma_N \\subset \\Omega$ denotes the part of the boundary where we prescribe Neumann boundary conditions. $\\boldsymbol{n}$ denotes the unit normal of $\\partial \\Omega$ pointing outside $\\Omega$. \n",
    "\n",
    "To obtain the weak form we define the functional spaces $V_{u_D} := \\left\\{ u \\in H^1(\\Omega) \\, |\\, u = u_D \\text{ on } \\Gamma_D \\right\\}$ and $V_{0} := \\left\\{ u \\in H^1(\\Omega) \\, |\\, u = 0 \\text{ on } \\Gamma_D \\right\\}$.\n",
    "Then we multiply the strong form by an arbitrary function $v \\in V_0$ and integrate over $\\Omega$:\n",
    "\n",
    "$$ - \\int_\\Omega \\Delta u \\, v \\, dx = \\int_\\Omega f\\,v \\, dx, \\quad \\forall v \\in V_0. $$\n",
    "\n",
    "Integration by parts of the non-conforming term gives\n",
    "\n",
    "$$ - \\int_\\Omega \\Delta u \\, v \\, dx = \\int_\\Omega \\nabla u \\cdot \\nabla v \\, dx - \\int_{\\partial\\Omega} (\\nabla u \\cdot \\boldsymbol{n}) \\,v\\, ds $$\n",
    "\n",
    "Recalling that $v = 0$ on $\\Gamma_D$ and that $\\nabla u \\cdot \\boldsymbol{n} = g$ on $\\Gamma_N$, the weak form of the BVP is the following.\n",
    "\n",
    "Find $u \\in V_{u_D}$:\n",
    "$$ \\int_\\Omega \\nabla u \\cdot \\nabla v \\, dx = \\int_\\Omega f\\,v \\, dx + \\int_{\\Gamma_N} g\\,v\\,ds, \\quad \\forall v \\in V_0. $$\n",
    "\n",
    "To obtain the finite element discretization we then introduce a triangulation (mesh) $\\mathcal{T}_h$ of the domain $\\Omega$ and we define a finite dimensional subspace $V_h \\subset H^1(\\Omega)$ consisting of globally continuous functions that are piecewise polynomial on each element of $\\mathcal{T}_h$.\n",
    "\n",
    "By letting $V_{h, u_D} := \\{ v_h \\in V_h \\,|\\, v_h = u_D \\text{ on } \\Gamma_D\\}$ and $V_{h, 0} := \\{ v_h \\in V_h \\,|\\, v_h = 0 \\text{ on } \\Gamma_D\\}$, the finite element method then reads:\n",
    "\n",
    "Find $u_h \\in V_{h, u_D}$ such that:\n",
    "$$ \\int_\\Omega \\nabla u_h \\cdot \\nabla v_h \\, dx = \\int_\\Omega f\\,v_h \\, dx + \\int_{\\Gamma_N} g\\,v_h\\,ds, \\quad \\forall v_h \\in V_{h,0}. $$\n",
    "\n",
    "In what follow, we will let $\\Omega := [0,1]\\times[0,1]$ be the unit square, $\\Gamma_N := \\{ (x,y) \\in \\partial\\Omega \\,|\\, y = 1\\}$ be the top boundary, and $\\Gamma_D := \\partial\\Omega \\setminus \\Gamma_N$ be the union of the left, bottom, and right boundaries.\n",
    "\n",
    "The coefficient $f$, $g$, $u_D$ are chosen such that the analytical solution is $u_{ex} = e^{\\pi y} \\sin(\\pi x)$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Imports\n",
    "We import the following Python packages:\n",
    "\n",
    "- `dolfin` is the python interface to FEniCS.\n",
    "- `matplotlib` is a plotting library that produces figure similar to the Matlab ones.\n",
    "- `math` is the python built-in library of mathematical functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, absolute_import, division\n",
    "\n",
    "# Import FEniCS\n",
    "import dolfin as dl\n",
    "import math\n",
    "\n",
    "# Enable plotting inside the notebook\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from hippylib import nb\n",
    "\n",
    "import logging\n",
    "logging.getLogger('FFC').setLevel(logging.WARNING)\n",
    "logging.getLogger('UFL').setLevel(logging.WARNING)\n",
    "dl.set_log_active(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Define the mesh and the finite element space\n",
    "We define a triangulation (mesh) of the unit square $\\Omega = [0,1]\\times[0,1]$ with `n` elements in each direction. The mesh size $h$ is $\\frac{1}{n}$.\n",
    "\n",
    "We also define the finite element space $V_h$ as the space of globally continuos functions that are piecewise polinomial (of degree $d$) on the elements of the mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 64\n",
    "d = 1\n",
    "mesh = dl.UnitSquareMesh(n, n)\n",
    "Vh = dl.FunctionSpace(mesh, \"Lagrange\", d)\n",
    "print(\"Number of dofs\", Vh.dim())\n",
    "nb.plot(mesh, mytitle=\"Finite Element Mesh\", show_axis='on')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define the Dirichlet boundary condition\n",
    "\n",
    "We define the Dirichlet boundary condition $u = u_d := \\sin(\\pi x)$ on $\\Gamma_D$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def boundary_d(x, on_boundary):\n",
    "    return (x[1] < dl.DOLFIN_EPS or x[0] < dl.DOLFIN_EPS or x[0] > 1.0 - dl.DOLFIN_EPS) and on_boundary\n",
    "\n",
    "u_d  = dl.Expression(\"sin(DOLFIN_PI*x[0])\", degree = d+2)\n",
    "bcs = [dl.DirichletBC(Vh, u_d, boundary_d)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Define the variational problem\n",
    "\n",
    "We write the variational problem $a(u_h, v_h) = F(v_h)$. Here, the bilinear form $a$ and the linear form $L$ are defined as\n",
    "\n",
    "- $a(u_h, v_h) := \\int_\\Omega \\nabla u_h \\cdot \\nabla v_h \\, dx$\n",
    "- $L(v_h) := \\int_\\Omega f v_h \\, dx + \\int_{\\Gamma_N} g \\, v_h \\, dx$.\n",
    "\n",
    "$u_h$ denotes the trial function and $v_h$ denotes the test function.  The coefficients $f = 0$ and $g = \\pi\\, e^{\\pi y} \\sin( \\pi x) $ are also given."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uh = dl.TrialFunction(Vh)\n",
    "vh = dl.TestFunction(Vh)\n",
    "\n",
    "f = dl.Constant(0.)\n",
    "g = dl.Expression(\"DOLFIN_PI*exp(DOLFIN_PI*x[1])*sin(DOLFIN_PI*x[0])\", degree=d+2)\n",
    "a = dl.inner(dl.grad(uh), dl.grad(vh))*dl.dx\n",
    "L = f*vh*dl.dx + g*vh*dl.ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Assemble and solve the finite element discrete problem\n",
    "\n",
    "We now assemble the finite element stiffness matrix $A$ and the right hand side vector $b$. Dirichlet boundary conditions are applied at the end of the finite element assembly procedure and before solving the resulting linear system of equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A, b = dl.assemble_system(a, L, bcs)\n",
    "uh = dl.Function(Vh)\n",
    "dl.solve(A, uh.vector(), b)\n",
    "\n",
    "nb.plot(uh, mytitle=\"Finite Element Solution\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Compute error norms\n",
    "\n",
    "We then compute the $L^2(\\Omega)$ and the energy norm of the difference between the exact solution and the finite element approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_ex = dl.Expression(\"exp(DOLFIN_PI*x[1])*sin(DOLFIN_PI*x[0])\", degree = d+2, domain=mesh)\n",
    "grad_u_ex = dl.Expression( (\"DOLFIN_PI*exp(DOLFIN_PI*x[1])*cos(DOLFIN_PI*x[0])\",\n",
    "                         \"DOLFIN_PI*exp(DOLFIN_PI*x[1])*sin(DOLFIN_PI*x[0])\"), degree = d+2, domain=mesh )\n",
    "\n",
    "norm_u_ex   = math.sqrt(dl.assemble(u_ex**2*dl.dx))\n",
    "norm_grad_ex = math.sqrt(dl.assemble(dl.inner(grad_u_ex, grad_u_ex)*dl.dx))\n",
    "\n",
    "err_L2   = math.sqrt(dl.assemble((uh - u_ex)**2*dl.dx))\n",
    "err_grad = math.sqrt(dl.assemble(dl.inner(dl.grad(uh) - grad_u_ex, dl.grad(uh) - grad_u_ex)*dl.dx))\n",
    "\n",
    "print (\"|| u_ex - u_h ||_L2 / || u_ex ||_L2 = \", err_L2/norm_u_ex)\n",
    "print (\"|| grad(u_ex - u_h)||_L2 / = || grad(u_ex)||_L2 \", err_grad/norm_grad_ex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hands on: Advection-Diffusion-Reaction PDEs\n",
    "\n",
    "For a bounded domain $\\Omega \\subset \\mathbb{R}^2$ with boundary $\\partial \\Omega = \\Gamma_D \\cup \\Gamma_N$, consider the boundary value problem\n",
    "\n",
    "$$ \n",
    "\\left\\{\n",
    "\\begin{array}{ll}\n",
    "- \\nabla \\cdot(k \\nabla u) + \\boldsymbol{w} \\cdot \\nabla u + c u = f  & \\text{in} \\; \\Omega, \\\\\n",
    "u = u_D & \\text{on} \\; \\Gamma_D, \\\\\n",
    "k \\nabla u \\cdot \\boldsymbol{n} = g & \\text{on} \\; \\Gamma_N.\n",
    "\\end{array}\n",
    "\\right.$$\n",
    "\n",
    "Here, $\\Gamma_D \\subset \\Omega$ denotes the part of the boundary where we prescribe Dirichlet boundary conditions, and $\\Gamma_N \\subset \\Omega$ denotes the part of the boundary where we prescribe Neumann boundary conditions. $\\boldsymbol{n}$ denotes the unit normal of $\\partial \\Omega$ pointing outside $\\Omega$.\n",
    "\n",
    "### Question 1\n",
    "> Derive the weak formulation of the above boundary value problem.\n",
    "\n",
    "### Question 2\n",
    "> Compute the finite element solution of the above problem using FEniCS with $k = k(\\boldsymbol{x}) = e^{-\\|\\boldsymbol{x}\\|^2}$, $\\boldsymbol{w} = [0,1]^t$, $c=1$, $f = 1$, $u_D = 0$, $g=0$. Choose $\\Omega$, $\\Gamma_N$, and $\\Gamma_D$ as in the example above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copyright &copy; 2018, The University of Texas at Austin & University of California, Merced. All Rights reserved. See file COPYRIGHT for details.\n",
    "\n",
    "This file is part of the hIPPYlib library. For more information and source code availability see https://hippylib.github.io.\n",
    "\n",
    "hIPPYlib is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License (as published by the Free Software Foundation) version 2.0 dated June 1991."
   ]
  }
 ],
 "metadata": {
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}