{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "center",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.3 Nonlinear problems\n",
    "\n",
    "In this unit we turn our attention to nonlinear PDE problems and the tools that `NGSolve` provides to simplify it."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A simple scalar PDE\n",
    "Let us start with a simple PDE with a nonlinearity: \n",
    "\n",
    "$$\n",
    "- \\Delta u + \\frac13 u^3 = 10 \\text{ in } \\Omega\n",
    "$$\n",
    "\n",
    "on the unit square $\\Omega = (0,1)^2$. \n",
    "\n",
    "We note that this PDE can also be formulated as a nonlinear minimization problem (cf. [3.4](../unit-3.4-nonlmin/nonlmin.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# define geometry and generate mesh\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import *\n",
    "from netgen.occ import *\n",
    "shape = Rectangle(1,1).Face()\n",
    "shape.edges.Min(X).name=\"left\"\n",
    "shape.edges.Max(X).name=\"right\"\n",
    "shape.edges.Min(Y).name=\"bottom\"\n",
    "shape.edges.Max(Y).name=\"top\"\n",
    "geom = OCCGeometry(shape, dim=2)\n",
    "mesh = Mesh(geom.GenerateMesh(maxh=0.3))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In NGSolve we can solve the PDE conveniently using the *linearization* feature of `SymbolicBFI`.\n",
    "\n",
    "The `BilinearForm` (**which is not bilinear!**) needed in the weak formulation is\n",
    "$$\n",
    "  A(u,v) = \\int_{\\Omega} \\nabla u \\nabla v + 1/3 u^3 v - 10 v ~ dx \\quad ( = 0 ~ \\forall~v \\in H^1_0)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "V = H1(mesh, order=3, dirichlet=[1,2,3,4])\n",
    "u,v = V.TnT()\n",
    "a = BilinearForm(V)\n",
    "a += (grad(u) * grad(v) + 1/3*u**3*v- 10 * v)*dx"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Newton's method\n",
    "\n",
    "In the preparation of the formulation of Newton's method we identify the semilinear form $A(\\cdot,\\cdot)$ with the corresponding operator $A: V \\to V'$ and write (assuming sufficient regularity)\n",
    "\n",
    " * $A(u)$ for the linear form $[A(u)](v) = A(u,v)$ and \n",
    " * the derivative of $A$ at an evaluation point $u$ is the bilinear form $\\delta A(u)$ with\n",
    " $$[\\delta A(u)] (w,v) = \\lim_{h \\to 0} \\frac{ A(u+h\\cdot w,v) - A(u,v)}{h}$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "You could compute these linear and bilinear forms manually. However, `NGSolve` provides functions for both operations: \n",
    "\n",
    "* The linear form $A(u)$ for a given (vector) u is obtained from `a.Apply(vec,res)` (resulting in a vector `res` representing the linear form)\n",
    "* The bilinear form (represented by the corresponding matrix) is stored in `a.mat` after a call of `a.AssembleLinearization(vec)`"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Under the hood it uses its functionality to derive `CoefficientFunction`s symbolically. \n",
    "\n",
    "For example, `NGSolve` derives the bilinear form integrand $ \\frac13 u^3 v $ w.r.t. $u$ at $\\tilde{u}$ in direction $w$ resulting in the integrand $\\tilde{u}^2 w v$.\n",
    "\n",
    "This allows to form the corresponding bilinear form integrals automatically for you."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To obtain a Newton algorithm we hence only need to translate the following pseudo-code formulation of Newton's method to `NGSolve`. The pseudo code is:\n",
    "\n",
    "`NewtonSolve`(Pseudo code):\n",
    "\n",
    "* Given an initial guess $u^0$\n",
    "* loop over $i=0,..$ until convergence:\n",
    "  * Compute linearization: $A (u^i) + \\delta A(u^i) \\Delta u^{i} = 0$:\n",
    "    * $f^i = A (u^i)$ \n",
    "    * $B^i = \\delta A(u^i)$ \n",
    "    * Solve $B^i \\Delta u^i = -f^i$\n",
    "  * Update $u^{i+1} = u^i + \\Delta u^{i}$\n",
    "  * Evaluate stopping criteria\n",
    "\n",
    "As a stopping criteria we take $\\langle A u^i,\\Delta u^i \\rangle = \\langle A u^i, A u^i \\rangle_{(B^i)^{-1}}< \\varepsilon$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now, here comes the same thing in `NGSolve` syntax:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def SimpleNewtonSolve(gfu,a,tol=1e-13,maxits=10, callback=lambda gfu: None):\n",
    "    res = gfu.vec.CreateVector()\n",
    "    du = gfu.vec.CreateVector()\n",
    "    fes = gfu.space\n",
    "    callback(gfu)\n",
    "    for it in range(maxits):\n",
    "        print (\"Iteration {:3}  \".format(it),end=\"\")\n",
    "        a.Apply(gfu.vec, res)\n",
    "        a.AssembleLinearization(gfu.vec)\n",
    "        du.data = a.mat.Inverse(fes.FreeDofs()) * res\n",
    "        gfu.vec.data -= du\n",
    "        callback(gfu)\n",
    "        #stopping criteria\n",
    "        stopcritval = sqrt(abs(InnerProduct(du,res)))\n",
    "        print (\"<A u\",it,\", A u\",it,\">_{-1}^0.5 = \", stopcritval)\n",
    "        if stopcritval < tol:\n",
    "            break"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "center",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's apply this to the previous PDE problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(V)\n",
    "gfu.Set((x*(1-x))**4*(y*(1-y))**4) # initial guess\n",
    "gfu_it = GridFunction(gfu.space,multidim=0)\n",
    "cb = lambda gfu : gfu_it.AddMultiDimComponent(gfu.vec) # store current state\n",
    "SimpleNewtonSolve(gfu, a, callback = cb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_style": "center",
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Draw(gfu,mesh,\"u\", deformation = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_style": "center",
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "Draw(gfu_it,mesh,\"u\", deformation = True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "center",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Use the `multidim`-Slider to inspect the results after the iterations."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "There are also some solvers shipped with NGSolve now. The `ngsolve.solvers.Newton` method allows you to use static condensation and is overall more refined (e.g. with damping options). For this tutorial we will mostly stay with the simple hand-crafted `SimpleNewtonSolve`. \n",
    "Here is the alternative demonstrated once, nevertheless:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from ngsolve.solvers import *\n",
    "help(Newton)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We call this Newton method as well here. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu.Set((x*(1-x))**4*(y*(1-y))**4) # initial guess\n",
    "Newton(a,gfu,freedofs=gfu.space.FreeDofs(),maxit=100,maxerr=1e-11,inverse=\"umfpack\",dampfactor=1,printing=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A trivial problem:\n",
    "As a second simple problem, let us consider a trivial scalar problem:\n",
    "$$\n",
    "  5 u^2 = 1, \\qquad u \\in \\mathbb{R}.\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "We chose this problem and put this in the (somewhat artificially) previous PDE framework and solve with Newton's method in a setting that you could easily follow with pen and paper:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "V = NumberSpace(mesh)\n",
    "u,v = V.TnT()\n",
    "a = BilinearForm(V)\n",
    "a += ( 5*u*u*v - 1 * v)*dx\n",
    "gfu = GridFunction(V)\n",
    "gfu.vec[:] = 1\n",
    "SimpleNewtonSolve(gfu,a, callback = lambda gfu : print(f\"u^k = {gfu.vec[0]}, u^k**2 = {gfu.vec[0]**2}\"))\n",
    "print(f\"\\nscalar solution, {gfu.vec[0]}, exact: {sqrt(0.2)}, error: {abs(sqrt(0.2)-gfu.vec[0])}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Another example: Stationary Navier-Stokes:\n",
    "\n",
    "Next, we consider incompressible Navier-Stokes equations again. This time however stationary.\n",
    "\n",
    "Find $\\mathbf{u} \\in \\mathbf{V}$, $p \\in Q$, $\\lambda \\in \\mathbb{R}$ so that\n",
    "\n",
    "\\begin{align*}\n",
    "\\int_{\\Omega} \\nu \\nabla \\mathbf{u} : \\nabla \\mathbf{v} + (\\mathbf{u} \\cdot \\nabla) \\mathbf{u} \\cdot \\mathbf{v}& - \\int_{\\Omega} \\operatorname{div}(\\mathbf{v}) p & &= \\int \\mathbf{f}  \\cdot \\mathbf{v}  && \\forall \\mathbf{v} \\in \\mathbf{V}, \\\\ \n",
    "- \\int_{\\Omega} \\operatorname{div}(\\mathbf{u}) q & & \n",
    "+ \\int_{\\Omega} \\lambda q\n",
    "&= 0 && \\forall q \\in Q, \\\\\n",
    "& \\int_{\\Omega} \\mu p & &= 0 && \\forall \\mu \\in \\mathbb{R}.\n",
    "\\end{align*}\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The domain $\\Omega$ is still $(0,1)^2$ and we prescribe homogenuous Dirichlet bnd. conditions for the velocity, except for the top boundary where we prescribe a tangential velocity. This setup is known as \"driven cavity\".\n",
    "\n",
    "Note that we use a scalar constraint to fix the pressure level that is otherwise not controlled in the presence of pure Dirichlet conditions."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We use a higher order Taylor-Hood discretization again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mesh = Mesh (geom.GenerateMesh(maxh=0.05)); nu = Parameter(1)\n",
    "V = VectorH1(mesh,order=3,dirichlet=\"bottom|right|top|left\")\n",
    "Q = H1(mesh,order=2); \n",
    "N = NumberSpace(mesh); \n",
    "X = V*Q*N\n",
    "(u,p,lam), (v,q,mu) = X.TnT()\n",
    "a = BilinearForm(X)\n",
    "a += (nu*InnerProduct(grad(u),grad(v))+InnerProduct(grad(u)*u,v)\n",
    "      -div(u)*q-div(v)*p-lam*q-mu*p)*dx"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The boundary condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(X)\n",
    "gfu.components[0].Set(CF((4*x*(1-x),0)),\n",
    "                      definedon=mesh.Boundaries(\"top\"))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now, let's apply the Newton:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     3,
     11
    ],
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def SolveAndVisualize(multidim=True):\n",
    "    gfu.components[0].Set(CF((4*x*(1-x),0)),\n",
    "                      definedon=mesh.Boundaries(\"top\"))\n",
    "    if multidim:\n",
    "        gfu_it = GridFunction(gfu.space,multidim=0)\n",
    "        cb = lambda gfu : gfu_it.AddMultiDimComponent(gfu.vec) # store current state\n",
    "        SimpleNewtonSolve(gfu, a, callback = cb)\n",
    "    else:\n",
    "        SimpleNewtonSolve(gfu, a)\n",
    "    Draw(gfu.components[0],mesh, vectors = {\"grid_size\" : 25})\n",
    "    print(\"above you see the solution after the Newton solve.\")\n",
    "    if multidim:\n",
    "        Draw(gfu_it.components[0], mesh, vectors = {\"grid_size\" : 25})\n",
    "        print(\"above you can inspect the results after each iteration of the Newton solve (use multidim-slider).\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "SolveAndVisualize()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The problem becomes more interesting if we decrease the viscosity $\\nu$ (`nu`), i.e. if we increase the Reynolds number *Re*. Let's consider the previous setup for decreasing values of `nu`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "nu.Set(0.01)\n",
    "SolveAndVisualize(multidim=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "nu.Set(0.01)\n",
    "SolveAndVisualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "nu.Set(0.001)\n",
    "SolveAndVisualize()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now, viscosity is so small that Newton does not converge anymore.\n",
    "In this case we can fix it using a small damping parameter with the `ngsolve.solvers.Newton`-solver:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "nu.Set(0.001)\n",
    "gfu.components[0].Set(CF((4*x*(1-x),0)),definedon=mesh.Boundaries(\"top\"))\n",
    "Newton(a,gfu,maxit=20,dampfactor=0.1)\n",
    "Draw(gfu.components[0],mesh, vectors = {\"grid_size\" : 25})"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Tasks** \n",
    "\n",
    "After these explanations, here are a few suggestions for simple play-around-tasks:\n",
    "\n",
    "* Take the first PDE example and set up the linearization linear and bilinear forms by hand and implement a Newton solver without exploiting the convenience functions `NGSolve` provides to you.\n",
    "* Combine [unit 3.1](../unit-3.1-parabolic/parabolic.ipynb) and this unit and write an implicit time integration solver for the unsteady Navier-Stokes equations (start with an implicit Euler)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}