Labs/Lab2_Tuesday/UnconstrainedMinimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Minimization of nonlinear energy functional\n",
    "\n",
    "In this example we solve the following nonlinear minimization problem\n",
    "\n",
    "*Find * $u^* \\in H^1_0(\\Omega)$ *such that*\n",
    "\n",
    "$$ u^* = \\operatorname*{argmin}_{u \\in H^1_0(\\Omega)} \\Pi(u). $$\n",
    "\n",
    "Here the energy functional $\\Pi(u)$ has the form\n",
    "\n",
    "$$ \\Pi(u) = \\frac{1}{2} \\int_\\Omega k(u) \\nabla u \\cdot \\nabla u \\, dx - \\int_\\Omega f\\,u \\, dx,$$\n",
    "\n",
    "where\n",
    "\n",
    "$$k(u) = k_1 + k_2 u^2. $$\n",
    "\n",
    "## Necessary optimality condition (Euler-Lagrange condition)\n",
    "\n",
    "Let $\\delta_u \\Pi(u, \\hat{u})$ denote the first variation of $\\Pi(u)$ in the *direction* $\\hat{u}$, i.e.\n",
    "\n",
    "$$\\delta_u \\Pi(u, \\hat{u}) := \\left. \\frac{d}{d \\varepsilon} \\Pi(u + \\varepsilon \\hat{u})\\right|_{\\varepsilon=0} = \\lim_{\\varepsilon \\rightarrow 0} \\frac{\\Pi(u + \\varepsilon \\hat{u}) - \\Pi(u)}{\\varepsilon}.$$\n",
    "\n",
    "The necessary condition is that the first variation of $\\Pi(u)$ equals to 0 for all directions $\\hat{u}$:\n",
    "\n",
    "$$ \\delta_u \\Pi = 0 \\Longleftrightarrow \\lim_{\\varepsilon \\rightarrow 0} \\frac{\\Pi(u + \\varepsilon \\hat{u}) - \\Pi(u)}{\\varepsilon} = 0 \\quad \\forall \\hat{u} \\in H_0^1(\\Omega).$$\n",
    "\n",
    "### Weak form:\n",
    "\n",
    "To obtain the weak form of the above necessary condition, we first expand the term $\\Pi(u + \\varepsilon \\hat{u})$ as\n",
    "\n",
    "$$\\Pi(u + \\varepsilon \\hat{u}) = \\frac{1}{2} \\int_\\Omega [k_1 + k_2(u + \\varepsilon\\hat{u})^2] (\\nabla u + \\varepsilon \\nabla \\hat{u})\\cdot(\\nabla u + \\varepsilon \\nabla \\hat{u}) dx - \\int_\\Omega f\\,(u+\\varepsilon\\hat{u}) dx.$$\n",
    "\n",
    "After some simplification, we obtain\n",
    "\n",
    "$$\\frac{\\Pi(u + \\varepsilon \\hat{u}) - \\Pi(u)}{\\epsilon} = \\int_\\Omega \\left[k_2 u \\hat{u} \\nabla u \\cdot \\nabla u + (k_1 + k_2u^2)\\nabla \\hat{u}\\cdot \\nabla u\\right] dx - \\int_\\Omega f \\hat{u} dx + \\mathcal{O}(\\epsilon).$$\n",
    "\n",
    "By neglecting the $\\mathcal{O}(\\epsilon)$ terms, we write the weak form of the necessary conditions as\n",
    "\n",
    "*Find *$u\\in H_0^1(\\Omega)$ *such that*\n",
    "\n",
    "$$ \\int_\\Omega \\left[k_2 u \\hat{u} \\nabla u \\cdot \\nabla u + (k_1 + k_2u^2)\\nabla \\hat{u}\\cdot \\nabla u\\right] dx = \\int_\\Omega f \\hat{u} dx \\quad \\forall \\hat{u} \\in H_0^1.$$\n",
    "\n",
    "### Strong form:\n",
    "To obtain the strong form, we invoke Green's first identity and write\n",
    "\n",
    "$$ \\int_\\Omega \\left[k_2 u \\nabla u \\cdot \\nabla u -  \\nabla \\cdot [(k_1 + k_2u^2) \\nabla u] \\right] \\hat{u} dx + \\int_{\\partial \\Omega} [(k_1 + k_2u^2) \\nabla u]\\cdot n \\hat{u} ds = \\int_\\Omega f \\hat{u} dx \\quad \\forall \\hat{u} \\in H_0^1.$$\n",
    "\n",
    "Since $\\hat{u}$ is arbitrary in $\\Omega$ and $\\hat{u} = 0$ on $\\partial \\Omega$, the strong form of the non-linear boundary problem reads\n",
    "\n",
    "$$ - \\nabla \\cdot [(k_1 + k_2u^2) \\nabla u + k_2 u \\nabla u \\cdot \\nabla u = f \\quad {\\rm in} \\; \\Omega; $$\n",
    "$$ u = 0 \\quad {\\rm on} \\; \\partial\\Omega.$$\n",
    "\n",
    "## Infinite-dimensional Newton's Method\n",
    "\n",
    "Consider the expansion of the first variation $\\delta_u \\Pi(u, \\hat{u})$ about $u$ in a *direction* $\\tilde{u}$\n",
    "\n",
    "$$\\delta_u \\Pi(u+\\tilde{u}, \\hat{u}) \\approx \\delta_u \\Pi(u, \\hat{u}) + \\delta_u^2\\Pi(u, \\hat{u}, \\tilde{u}),$$\n",
    "\n",
    "where\n",
    "\n",
    "$$ \\delta_u^2\\Pi(u, \\hat{u}, \\tilde{u}) = \\left. \\frac{d}{d\\varepsilon} \\delta_u \\Pi(u + \\varepsilon \\tilde{u}, \\hat{u}) \\right|_{\\varepsilon=0}.$$\n",
    "\n",
    "The infinite-dimensional Newton's method reads\n",
    "\n",
    "*Given the current solution *$u_k$, *find* $\\tilde{u} \\in H^1_0$ *such that*\n",
    "\n",
    "$$ \\delta_u^2 \\Pi(u_k, \\hat{u}, \\tilde{u}) = -\\delta_u \\Pi(u_k, \\hat{u}) \\quad \\forall \\, \\hat{u} \\in H_0^1.$$\n",
    "\n",
    "*Update the solution using the Newton direction* $\\tilde{u}$\n",
    "$$ u_{k+1} = u_k + \\tilde{u}.$$\n",
    "\n",
    "### Hessian\n",
    "\n",
    "To derive the weak form of the Hessian, we first expand the term $\\delta_u \\Pi(u +\\varepsilon \\tilde{u},\\hat{u})$ as\n",
    "\n",
    "$$\\delta_u \\Pi(u+\\varepsilon\\tilde{u}, \\hat{u}) = \\int_\\Omega \\left[k_2 (u+\\varepsilon\\tilde{u}) \\hat{u} \\nabla (u+\\varepsilon\\tilde{u}) \\cdot \\nabla (u+\\varepsilon\\tilde{u}) + (k_1 + k_2(u+\\varepsilon\\tilde{u})^2)\\nabla \\hat{u}\\cdot \\nabla (u+\\varepsilon\\tilde{u}) \\right] dx - \\int_\\Omega f \\hat{u} dx \\quad \\forall \\hat{u} \\in H_0^1.$$\n",
    "\n",
    "Then, after some simplification, we obtain\n",
    "\n",
    "$$\\delta^2 \\Pi(u, \\tilde{u}, \\hat{u}) := \\frac{d}{d\\varepsilon} \\delta_u \\Pi(u+\\varepsilon\\tilde{u}, \\hat{u}) = \n",
    "\\int_\\Omega \\left[k_2\\tilde{u}\\hat{u}\\nabla u \\cdot \\nabla u + 2k_2 u \\hat{u} \\nabla \\tilde{u} \\cdot \\nabla u + 2k_2 u \\tilde{u} \\nabla \\hat{u} \\cdot \\nabla u + (k_1 + k_2u^2) \\nabla \\hat{u} \\cdot \\nabla \\tilde{u} \\right] dx. $$\n",
    "\n",
    "### Weak form of Newton step:\n",
    "\n",
    "*Given *$u \\in H_0^1$, *find * $\\tilde{u} \\in H^1_0$ *such that*\n",
    "\n",
    "$$\\int_\\Omega \\left[k_2\\tilde{u}\\hat{u}\\nabla u \\cdot \\nabla u + 2k_2 u \\hat{u} \\nabla \\tilde{u} \\cdot \\nabla u + 2k_2 u \\tilde{u} \\nabla \\hat{u} \\cdot \\nabla u + (k_1 + k_2u^2) \\nabla \\hat{u} \\cdot \\nabla \\tilde{u} \\right] dx = - \\int_\\Omega \\left[k_2 u \\hat{u} \\nabla u \\cdot \\nabla u + (k_1 + k_2u^2)\\nabla \\hat{u}\\cdot \\nabla u -f \\hat{u} \\right] dx \\quad \\forall \\, \\hat{u} \\in H_0^1. $$\n",
    "\n",
    "The solution is then updated using the Newton direction $\\tilde{u}$\n",
    "\n",
    "$$ u^{\\rm new} = u + \\alpha \\tilde{u}.$$\n",
    "\n",
    "Here $\\alpha$ denotes a relaxation parameter (back-tracking/line-search) used to achieve global convergence of the Newton method.\n",
    "\n",
    "### Strong form of the Newton step\n",
    "\n",
    "$$ - \\nabla \\cdot \\left[ (k_1 + k_2 u^2) \\nabla \\tilde{u}\\right] + 2k_2u\\nabla\\tilde{u}\\cdot\\nabla u - \\nabla\\cdot(2k_2 u \\tilde{u} \\nabla u) + k_2 \\tilde{u} \\nabla u \\nabla u = \\nabla \\cdot\\left[(k_1 + k_2 u^2)\\nabla \\right]u - k_2 u \\nabla u\\cdot \\nabla u + f \\quad {\\rm in} \\, \\Omega.$$\n",
    "$$ \\tilde{u} = 0 \\quad {\\rm on} \\, \\partial \\Omega. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Load modules\n",
    "\n",
    "To start we load the following modules:\n",
    "\n",
    "- dolfin: the python/C++ interface to FEniCS\n",
    "\n",
    "- [math](https://docs.python.org/2/library/math.html): the python module for mathematical functions\n",
    "\n",
    "- [numpy](http://www.numpy.org/): a python package for linear algebra\n",
    "\n",
    "- [matplotlib](http://matplotlib.org/): a python package used for plotting the results\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division, absolute_import\n",
    "\n",
    "import dolfin as dl\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from hippylib import nb\n",
    "\n",
    "import math\n",
    "import numpy as np\n",
    "import logging\n",
    "\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 finite element spaces\n",
    "\n",
    "We construct a triangulation (mesh) $\\mathcal{T}_h$ of the computational domain $\\Omega := [0, 1]^2$ with `nx` elements in the *x*-axis direction and `ny` elements in the *y*-axis direction.\n",
    "\n",
    "On the mesh $\\mathcal{T}_h$, we then define the finite element space $V_h \\subset H^1(\\Omega)$ consisting of globally continuous piecewise linear functions and we create a function $u \\in V_h$.\n",
    "\n",
    "By denoting by $\\left[{\\phi_i(x)}\\right]_{i=1}^{{\\rm dim}(V_h)}$ the finite element basis for the space $V_h$ we have\n",
    "$$ u = \\sum_{i=1}^{{\\rm dim}(V_h)} {\\rm u}_i \\phi_i(x), $$\n",
    "where ${\\rm u}_i$ represents the coefficients in the finite element expansion of $u$.\n",
    "\n",
    "Finally we define two special types of functions: the `TestFunction` $\\hat{u}$ and the `TrialFunction` $\\tilde{u}$. These special types of functions are used by `FEniCS` to generate the finite element vectors and matrices which stem from the first and second variations of the energy functional $\\Pi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nx = 32\n",
    "ny = 32\n",
    "mesh = dl.UnitSquareMesh(nx,ny)\n",
    "Vh = dl.FunctionSpace(mesh, \"CG\", 1)\n",
    "\n",
    "uh = dl.Function(Vh)\n",
    "u_hat = dl.TestFunction(Vh)\n",
    "u_tilde = dl.TrialFunction(Vh)\n",
    "\n",
    "nb.plot(mesh, show_axis=\"on\")\n",
    "print( \"dim(Vh) = \", Vh.dim() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Define the energy functional\n",
    "\n",
    "We now define the energy functional\n",
    "$$ \\Pi(u) = \\frac{1}{2} \\int_\\Omega (k_1 + k_2u^2) \\nabla u \\cdot \\nabla u dx - \\int_\\Omega f\\,u dx.$$\n",
    "\n",
    "The parameters $k_1$, $k_2$ and the forcing term $f$ are defined in FEniCS using the keyword [`Constant`](https://fenicsproject.org/olddocs/dolfin/2017.2.0/python/programmers-reference/functions/constant/Constant.html). To define coefficients that are space dependent one should use the keyword [`Expression`](https://fenicsproject.org/olddocs/dolfin/2017.2.0/python/programmers-reference/functions/expression/Expression.html).\n",
    "\n",
    "The Dirichlet boundary condition\n",
    "$$ u = 0 \\quad {\\rm on} \\; \\partial\\Omega$$\n",
    "is imposed using the [`DirichletBC`](https://fenicsproject.org/olddocs/dolfin/2017.2.0/python/programmers-reference/fem/bcs/DirichletBC.html) class.\n",
    "\n",
    "To construct this object we need to provide\n",
    "\n",
    "- the finite element space `Vh`\n",
    "\n",
    "- the value `u_0` of the solution at the Dirichlet boundary. `u_0` can either be a `Constant` or an `Expression` object.\n",
    "\n",
    "- the object `Boundary` that defines on which part of $\\partial \\Omega$ we want to impose such condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = dl.Constant(1.)\n",
    "k1 = dl.Constant(0.05)\n",
    "k2 = dl.Constant(1.)\n",
    "\n",
    "Pi = dl.Constant(.5)*(k1 + k2*uh*uh)*dl.inner(dl.grad(uh), dl.grad(uh))*dl.dx - f*uh*dl.dx\n",
    "\n",
    "class Boundary(dl.SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        return on_boundary\n",
    "\n",
    "u_0 = dl.Constant(0.)    \n",
    "bc = dl.DirichletBC(Vh,u_0, Boundary() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. First variation\n",
    "\n",
    "The weak form of the first variation reads\n",
    "\n",
    "$$\\delta_u \\Pi(u, \\hat{u}) = \\int_\\Omega \\left[k_2 u \\hat{u} \\nabla u \\cdot \\nabla u + (k_1 + k_2u^2)\\nabla \\hat{u}\\cdot \\nabla u\\right] dx - \\int_\\Omega f \\hat{u} dx \\quad \\forall \\hat{u} \\in H_0^1.$$\n",
    "\n",
    "We use a **finite difference check** to verify that our derivation is correct.\n",
    "More specifically, we consider a function\n",
    "$$ u_0 = x(x-1)y(y-1) \\in H^1_0(\\Omega) $$\n",
    "and we verify that for a random direction $\\hat{u} \\in H^1_0(\\Omega)$ we have\n",
    "$$ r := \\left| \\frac{\\Pi(u_0 + \\varepsilon \\hat{u}) - \\Pi(u_0)}{\\varepsilon} - \\delta_u \\Pi(u, \\hat{u})\\right| = \\mathcal{O}(\\varepsilon).$$\n",
    "\n",
    "In the figure below we show in a loglog scale the value of $r$ as a function of $\\varepsilon$. We observe that $r$ decays linearly for a wide range of values of $\\varepsilon$, however we notice an increase in the error for extremely small values of $\\varepsilon$ due to numerical stability and finite precision arithmetic.\n",
    "\n",
    "**NOTE:** To compute the first variation we can also use the [automatic differentiation](https://fenicsproject.org/olddocs/dolfin/2017.2.0/python/programmers-reference/fem/formmanipulations/derivative.html#dolfin.fem.formmanipulations.derivative) of variational forms capabilities of FEniCS and write\n",
    "\n",
    "`grad = dl.derivative(Pi, u, u_hat)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grad = (k2*uh*u_hat)*dl.inner(dl.grad(uh), dl.grad(uh))*dl.dx + \\\n",
    "       (k1 + k2*uh*uh)*dl.inner(dl.grad(uh), dl.grad(u_hat))*dl.dx - f*u_hat*dl.dx\n",
    "\n",
    "u0 = dl.interpolate(dl.Expression(\"x[0]*(x[0]-1)*x[1]*(x[1]-1)\", degree=2), Vh)\n",
    "\n",
    "n_eps = 32\n",
    "eps = 1e-2*np.power(2., -np.arange(n_eps))\n",
    "err_grad = np.zeros(n_eps)\n",
    "\n",
    "uh.assign(u0)\n",
    "pi0 = dl.assemble(Pi)\n",
    "grad0 = dl.assemble(grad)\n",
    "\n",
    "uhat = dl.Function(Vh).vector()\n",
    "uhat.set_local(np.random.randn(Vh.dim()))\n",
    "uhat.apply(\"\")\n",
    "bc.apply(uhat)\n",
    "dir_grad0 = grad0.inner(uhat)\n",
    "\n",
    "for i in range(n_eps):\n",
    "    uh.assign(u0)\n",
    "    uh.vector().axpy(eps[i], uhat) #uh = uh + eps[i]*dir\n",
    "    piplus = dl.assemble(Pi)\n",
    "    err_grad[i] = abs( (piplus - pi0)/eps[i] - dir_grad0 )\n",
    "\n",
    "plt.figure()    \n",
    "plt.loglog(eps, err_grad, \"-ob\", label=\"Error Grad\")\n",
    "plt.loglog(eps, (.5*err_grad[0]/eps[0])*eps, \"-.k\", label=\"First Order\")\n",
    "plt.title(\"Finite difference check of the first variation (gradient)\")\n",
    "plt.xlabel(\"eps\")\n",
    "plt.ylabel(\"Error grad\")\n",
    "plt.legend(loc = \"upper left\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Second variation\n",
    "\n",
    "The weak form of the second variation reads\n",
    "\n",
    "$$\\delta_u^2 \\Pi(u, \\tilde{u}, \\hat{u}) := \\frac{d}{d\\varepsilon} \\delta_u \\Pi(u+\\varepsilon\\tilde{u}, \\hat{u}) = \n",
    "\\int_\\Omega \\left[k_2\\tilde{u}\\hat{u}\\nabla u \\cdot \\nabla u + 2k_2 u \\hat{u} \\nabla \\tilde{u} \\cdot \\nabla u + 2k_2 u \\tilde{u} \\nabla \\hat{u} \\cdot \\nabla u + (k_1 + k_2u^2) \\nabla \\hat{u} \\cdot \\nabla \\tilde{u} \\right] dx. $$\n",
    "\n",
    "As before, we verify that for a random direction $\\hat{u} \\in H^1_0(\\Omega)$ we have\n",
    "$$ r := \\left\\| \\frac{\\delta_u\\Pi(u_0 + \\varepsilon \\tilde{u}, \\hat{u}) - \\delta_u \\Pi(u_0, \\hat{u})}{\\varepsilon} - \\delta_u^2 \\Pi(u, \\tilde{u}, \\hat{u})\\right\\| = \\mathcal{O}(\\varepsilon).$$\n",
    "\n",
    "In the figure below we show in a loglog scale the value of $r$ as a function of $\\varepsilon$. As before, we observe that $r$ decays linearly for a wide range of values of $\\varepsilon$, however we notice an increase in the error for extremely small values of $\\varepsilon$ due to numerical stability and finite precision arithmetic.\n",
    "\n",
    "**NOTE:** To compute the second variation we can also use automatic differentiation and write\n",
    "\n",
    "`H = dl.derivative(grad, u, u_tilde)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = k2*u_tilde*u_hat*dl.inner(dl.grad(uh), dl.grad(uh))*dl.dx + \\\n",
    "     dl.Constant(2.)*(k2*uh*u_hat)*dl.inner(dl.grad(u_tilde), dl.grad(uh))*dl.dx + \\\n",
    "     dl.Constant(2.)*k2*u_tilde*uh*dl.inner(dl.grad(uh), dl.grad(u_hat))*dl.dx + \\\n",
    "     (k1 + k2*uh*uh)*dl.inner(dl.grad(u_tilde), dl.grad(u_hat))*dl.dx\n",
    "\n",
    "uh.assign(u0)\n",
    "H_0 = dl.assemble(H)\n",
    "H_0uhat = H_0 * uhat\n",
    "err_H = np.zeros(n_eps)\n",
    "\n",
    "for i in range(n_eps):\n",
    "    uh.assign(u0)\n",
    "    uh.vector().axpy(eps[i], uhat)\n",
    "    grad_plus = dl.assemble(grad)\n",
    "    diff_grad = (grad_plus - grad0)\n",
    "    diff_grad *= 1/eps[i]\n",
    "    err_H[i] = (diff_grad - H_0uhat).norm(\"l2\")\n",
    "    \n",
    "plt.figure()    \n",
    "plt.loglog(eps, err_H, \"-ob\", label=\"Error Hessian\")\n",
    "plt.loglog(eps, (.5*err_H[0]/eps[0])*eps, \"-.k\", label=\"First Order\")\n",
    "plt.title(\"Finite difference check of the second variation (Hessian)\")\n",
    "plt.xlabel(\"eps\")\n",
    "plt.ylabel(\"Error Hessian\")\n",
    "plt.legend(loc = \"upper left\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. The infinite dimensional Newton Method\n",
    "\n",
    "The infinite dimensional Newton step reads\n",
    "\n",
    "*Given *$u_n \\in H_0^1$, *find * $\\tilde{u} \\in H^1_0$ *such that*\n",
    "$$ \\delta_u^2 \\Pi(u_n, \\hat{u}, \\tilde{u}) =  - \\delta_u \\Pi(u_n, \\hat{u}) \\quad $$\n",
    "Update the solution $u_{n+1}$ using the Newton direction $\\tilde{u}$\n",
    "$$ u_{n+1} = u + \\alpha \\tilde{u}.$$\n",
    "\n",
    "Here, for simplicity, we choose $\\alpha$ equal to 1. In general, to guarantee global convergence of the Newton method the parameter $\\alpha$ should be appropriately chosen (e.g. *back-tracking* or *line search*).\n",
    "\n",
    "The linear systems to compute the Newton directions are solved using the conjugate gradient (CG) with algebraic multigrid preconditioner with a fixed tolerance. In practice, one should solve the Newton system inexactly by early termination of CG \n",
    "iterations via Eisenstat–Walker (to prevent oversolving) and Steihaug (to avoid negative curvature) criteria.\n",
    "\n",
    "In the output below, for each iteration we report the number of CG iterations, the value of the energy functional, the norm of the gradient, and the inner product between the gradient and the Newton direction $\\delta_u \\Pi(u_0, \\tilde{u})$.\n",
    "\n",
    "In the example, the stopping criterion is relative norm of the gradient $\\frac{\\delta_u \\Pi(u_n, \\hat{u})}{\\delta_u \\Pi(u_0, \\hat{u})} \\leq \\tau$. However robust implementation of the stopping criterion should monitor also the quantity $\\delta_u \\Pi(u_0, \\tilde{u})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uh.assign(dl.interpolate(dl.Constant(0.), Vh))\n",
    "\n",
    "rtol = 1e-9\n",
    "max_iter = 10\n",
    "\n",
    "pi0 = dl.assemble(Pi)\n",
    "g0 = dl.assemble(grad)\n",
    "bc.apply(g0)\n",
    "tol = g0.norm(\"l2\")*rtol\n",
    "\n",
    "du = dl.Function(Vh).vector()\n",
    "\n",
    "lin_it = 0\n",
    "print (\"{0:3} {1:3} {2:15} {3:15} {4:15}\".format(\n",
    "      \"It\", \"cg_it\", \"Energy\", \"(g,du)\", \"||g||l2\") )\n",
    "\n",
    "for i in range(max_iter):\n",
    "    [Hn, gn] = dl.assemble_system(H, grad, bc)\n",
    "    if gn.norm(\"l2\") < tol:\n",
    "        print (\"\\nConverged in \", i, \"Newton iterations and \", lin_it, \"linear iterations.\")\n",
    "        break\n",
    "    myit = dl.solve(Hn, du, gn, \"cg\", \"petsc_amg\")\n",
    "    lin_it = lin_it + myit\n",
    "    uh.vector().axpy(-1., du)\n",
    "    pi = dl.assemble(Pi)\n",
    "    print (\"{0:3d} {1:3d} {2:15e} {3:15e} {4:15e}\".format(\n",
    "      i, myit, pi, -gn.inner(du), gn.norm(\"l2\")) )\n",
    "    \n",
    "    plt.figure()\n",
    "    nb.plot(uh, mytitle=\"Iteration {0:1d}\".format(i))\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uh.assign(dl.interpolate(dl.Constant(0.), Vh))\n",
    "parameters={\"symmetric\": True, \"newton_solver\": {\"relative_tolerance\": 1e-9, \"report\": True, \\\n",
    "                                                 \"linear_solver\": \"cg\", \"preconditioner\": \"petsc_amg\"}}\n",
    "dl.solve(grad == 0, uh, bc, J=H, solver_parameters=parameters)\n",
    "final_g = dl.assemble(grad)\n",
    "bc.apply(final_g)\n",
    "\n",
    "print ( \"Norm of the gradient at converge\", final_g.norm(\"l2\") )\n",
    "print (\"Value of the energy functional at convergence\", dl.assemble(Pi) )\n",
    "nb.plot(uh)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hands on\n",
    "\n",
    "Consider the following nonlinear minimization problem\n",
    "\n",
    "*Find * $u^* \\in H^1(\\Omega)$ *such that*\n",
    "\n",
    "$$ u^* = \\operatorname*{argmin}_{u \\in H^1(\\Omega)} \\Pi(u), $$\n",
    "\n",
    "where\n",
    "\n",
    "$$ \\Pi(u) = \\frac{1}{2} \\int_\\Omega |\\nabla u|^2 \\, dx + \\int_\\Omega e^{-u} \\, dx + \\int_{\\partial\\Omega} u^2 \\, ds \\quad \\forall u \\in H^1(\\Omega). $$\n",
    "\n",
    "### Question 1\n",
    "\n",
    "Derive the first-order necessary condition for optimality using calculus of variations, in both weak and strong form.\n",
    "\n",
    "### Question 2\n",
    "\n",
    "Derive the infinite-dimensional Newton step, in both weak and strong form.\n",
    "\n",
    "### Question 3\n",
    "\n",
    "Discretize and solve the above nonlinear minimization problem using FEniCS.\n"
   ]
  }
 ],
 "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
}