6-NewtonRaphson.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "21cf672e-ba27-480c-a1e0-39c6dddd09f8",
   "metadata": {},
   "source": [
    "# Your first nonlinear solver- newton-raphson iteration\n",
    "\n",
    "- [x] newton-raphson method\n",
    "- [x] weak form\n",
    "- [x] residual and jacobian\n",
    "- [x] nonlinear iteration\n",
    "- [x] convergence criterion\n",
    "- [x] nonlinear poisson equation\n",
    "\n",
    "Author: Yang Bai @ M3 Group\n",
    "\n",
    "Date  : 2022.04.17\n",
    "\n",
    "QQ group: 628204857"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "279b0813-25ab-4118-b26b-f08ef9a516e0",
   "metadata": {},
   "source": [
    "# Newton-raphson method\n",
    "\n",
    "How to find the root of $f(x)=0$?\n",
    "\n",
    "Find $x_{0}$ which s.a. $f(x_{0})=0$.\n",
    "\n",
    "## Taylor expansion\n",
    "$f(x^{})=f(x^{k})+\\frac{\\partial f(x^{k})}{\\partial x}(x^{}-x^{k})+\\frac{\\partial^{2} f(x^{k})}{\\partial x^{2}}\\frac{1}{2!}(x^{}-x^{k})^{2}+\\dots$\n",
    "\n",
    "$f(x_{0})=f(x^{k})+\\frac{\\partial f(x^{k})}{\\partial x}(x_{0}-x^{k})=0$\n",
    "\n",
    "$x_{0}=x^{k}-\\frac{f(x^{k})}{\\frac{\\partial f(x^{k})}{\\partial x}}$\n",
    "\n",
    "$f(x^{k})=-\\frac{\\partial f(x^{k})}{\\partial x}(x_{0}-x^{k})=K\\Delta x$, with $K=-\\frac{\\partial f(x^{k})}{\\partial x}$, $\\Delta x=x_{0}-x^{k}$\n",
    "\n",
    "## Example\n",
    "\n",
    "### Ex. 1 \n",
    "$f(x)=x-1=0$\n",
    "\n",
    "$K=-\\frac{\\partial f}{\\partial x}=-1$\n",
    "\n",
    "1st iteration: set x=0, f(x)=-1 ===> dx=K\\f=1\n",
    "\n",
    "--- update x=x+dx=1  check f(x)=0, done!\n",
    "\n",
    "### Ex 2\n",
    "$f(x)=1.0+sin(x)-x$\n",
    "\n",
    "$K=-\\frac{\\partial f}{\\partial x}=-cos(x)+1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9fb9ad30-27e4-4371-afe9-9d4ba81c3323",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iters=  1, |R|=  9.794255e-01, |dx|=  8.000703e+00\n",
      "iters=  2, |R|=  6.702639e+00, |dx|=  4.182424e+00\n",
      "iters=  3, |R|=  4.241618e+00, |dx|=  3.064783e+00\n",
      "iters=  4, |R|=  6.965849e-01, |dx|=  1.012482e+00\n",
      "iters=  5, |R|=  4.980408e-01, |dx|=  3.035862e-01\n",
      "iters=  6, |R|=  3.809059e-02, |dx|=  2.756864e-02\n",
      "iters=  7, |R|=  3.525588e-04, |dx|=  2.599914e-04\n",
      "iters=  8, |R|=  3.158407e-08, |dx|=  2.329557e-08\n",
      "iters=  9, |R|=  4.440892e-16, |dx|=  3.275484e-16\n",
      "the converged solution is: x=  1.934563e+00, f(x)=-2.22044605e-16\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "def f(x):\n",
    "    return 1.0+np.sin(x)-x\n",
    "def df(x):\n",
    "    return np.cos(x)-1\n",
    "\n",
    "x=0.5;rnorm=1.0;iters=0;tolerance=1.0e-10\n",
    "while iters<200 and rnorm>tolerance:\n",
    "    R=f(x);K=-df(x)\n",
    "    dx=R/K\n",
    "    x+=dx\n",
    "    rnorm=np.abs(R)\n",
    "    iters+=1\n",
    "    print('iters=%3d, |R|=%14.6e, |dx|=%14.6e'%(iters,rnorm,np.abs(dx)))\n",
    "if(not rnorm>tolerance):\n",
    "    print('the converged solution is: x=%14.6e, f(x)=%15.8e'%(x,f(x)))\n",
    "else:\n",
    "    print('oops, your newton-raphson failed!')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f86fd16b-92e3-4fa0-84c9-2dea6eaabb84",
   "metadata": {},
   "source": [
    "# Newton-raphson in FEM\n",
    "\n",
    "For a FEM problem $\\Pi\\sim\\Pi(u,\\nabla u,\\dot{u},\\dots)$, one need to find $u$ s.a. $\\delta\\Pi=R(u,\\nabla u,\\dot{u},\\dots)\\delta u=0$ with $R=[K(u,\\nabla u,\\dot{u})]\\{u\\}-f(u,\\nabla u,\\dot{u})$\n",
    "\n",
    "## Nonlinear poisson equation\n",
    "\n",
    "$\\nabla\\cdot(\\sigma\\nabla\\phi)=f$ with $\\sigma\\sim\\sigma(\\phi)$ and $f\\sim f(\\phi)$, $-\\sigma\\nabla\\phi\\cdot\\vec{n}=0$\n",
    "\n",
    "### weak form\n",
    "$\\int_{\\partial\\Omega}\\sigma\\nabla\\phi\\cdot\\vec{n}\\delta\\phi dS-\\int_{\\Omega}\\sigma\\nabla\\phi\\nabla\\delta\\phi dV=\\int_{\\Omega}f\\delta\\phi dV$\n",
    "\n",
    "then the residual and jacobian matrix can be read as follows:\n",
    "\n",
    "$R^{I}=\\int_{\\Omega}\\sigma\\nabla\\phi\\nabla N^{I}dV+\\int_{\\Omega}fN^{I}dV$\n",
    "\n",
    "$K^{IJ}=-\\frac{\\partial R^{I}}{\\partial\\phi^{J}}=-\\int_{\\Omega}\\frac{\\partial\\sigma}{\\partial\\phi}N^{J}\\nabla\\phi\\nabla N^{I}dV-\\int_{\\Omega}\\sigma\\nabla N^{J}\\nabla N^{I}dV-\\int_{\\Omega}\\frac{\\partial f}{\\partial\\phi}N^{J}N^{I}dV$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f94144e7-81f3-4fc6-b65f-21395e167413",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iters=  1, |R0 |=   1.93746e-01, |R |=   1.93746e-01\n",
      "           |dU0|=   1.05849e+01, |dU|=   1.05849e+01\n",
      "           |E0 |=   2.05079e+00, |E |=   2.05079e+00\n",
      "iters=  2, |R0 |=   1.93746e-01, |R |=   4.72237e-02\n",
      "           |dU0|=   1.05849e+01, |dU|=   1.55848e+00\n",
      "           |E0 |=   2.05079e+00, |E |=   7.35972e-02\n",
      "iters=  3, |R0 |=   1.93746e-01, |R |=   9.57083e-04\n",
      "           |dU0|=   1.05849e+01, |dU|=   2.46349e-02\n",
      "           |E0 |=   2.05079e+00, |E |=   2.35776e-05\n",
      "iters=  4, |R0 |=   1.93746e-01, |R |=   4.03722e-07\n",
      "           |dU0|=   1.05849e+01, |dU|=   6.51174e-06\n",
      "           |E0 |=   2.05079e+00, |E |=   2.62893e-12\n",
      "iters=  5, |R0 |=   1.93746e-01, |R |=   5.72883e-14\n",
      "           |dU0|=   1.05849e+01, |dU|=   5.33811e-13\n",
      "           |E0 |=   2.05079e+00, |E |=   3.05811e-26\n",
      "write result to nlpoisson-000000.vtu\n"
     ]
    }
   ],
   "source": [
    "from FEToy.mesh.lagrange2dmesh import mesh2d\n",
    "from FEToy.fe.shapefun import shape2d\n",
    "from FEToy.fe.gaussrule import gausspoint2d\n",
    "from FEToy.postprocess.Result import ResultIO\n",
    "\n",
    "# define your mesh\n",
    "mymesh=mesh2d(xmax=2.0,ymax=2.0,nx=20,ny=20,meshtype='quad4')\n",
    "mymesh.createmesh()\n",
    "\n",
    "# define your shape function\n",
    "shp=shape2d(meshtype='quad4')\n",
    "shp.update()\n",
    "\n",
    "# define your gauss points\n",
    "qpoints=gausspoint2d(ngp=3)\n",
    "qpoints.creategausspoint()\n",
    "\n",
    "\n",
    "### defines your physical parameters here\n",
    "# since sigma and f both are nonlinear, we can offer functions for them\n",
    "def sigma(x):\n",
    "    return 1+np.exp(x)\n",
    "    # return 1.0+x*x\n",
    "def dsigma(x):\n",
    "    return np.exp(x)\n",
    "    # return 2*x\n",
    "### for f\n",
    "def f(x):\n",
    "    return 1.0+0.1*x\n",
    "def df(x):\n",
    "    return 0.1\n",
    "\n",
    "# define your K and rhs, and other vectors\n",
    "nDofs=mymesh.nodes*1\n",
    "K=np.zeros((nDofs,nDofs))\n",
    "R=np.zeros(nDofs)\n",
    "U=np.zeros(nDofs)\n",
    "penalty=1.0e16 # for dirichlet bc\n",
    "\n",
    "# for parameters used in NR iteration\n",
    "tolerance=1.0e-10;reltol=1.0e-15 # absolute tolerenace and relative tolerance\n",
    "maxiters=50\n",
    "rnorm=1.0;rnorm0=1.0;iters=0 # init\n",
    "dunorm=1.0;dunorm0=1.0\n",
    "enorm=1.0;enorm0=1.0\n",
    "converged=False\n",
    "\n",
    "gradphi=np.zeros(2)\n",
    "# now we can do the newton-raphson iteration\n",
    "while iters<maxiters and (not converged):\n",
    "    K[:,:]=0.0;R[:]=0.0 # init\n",
    "    \n",
    "    # do the element loop\n",
    "    for e in range(mymesh.elements):\n",
    "        elconn=mymesh.elementconn[e,:]\n",
    "        nodes=mymesh.nodecoords[elconn]\n",
    "        for gp in range(qpoints.ngp2):\n",
    "            xi =qpoints.gpcoords[gp,1]\n",
    "            eta=qpoints.gpcoords[gp,2]\n",
    "            w  =qpoints.gpcoords[gp,0]\n",
    "            shp_val,shp_grad,j=shp.calc(xi,eta,nodes[:,0],nodes[:,1]) # xi,eta,x,y\n",
    "            JxW=j*w # the weight times the determinate of jacobian transformation\n",
    "            \n",
    "            phi=0.0 # for the phi on current gauss point\n",
    "            gradphi[0]=0.0;gradphi[1]=0.0\n",
    "            x=0.0;y=0.0 # coordinates of current gauss point\n",
    "            for i in range(mymesh.nodesperelement):\n",
    "                iInd=elconn[i] # global id\n",
    "                phi+=shp_val[i]*U[iInd]\n",
    "                gradphi[0]+=shp_grad[i,0]*U[iInd]\n",
    "                gradphi[1]+=shp_grad[i,1]*U[iInd]\n",
    "                x+=shp_val[i]*nodes[i,0]\n",
    "                y+=shp_val[i]*nodes[i,1]\n",
    "            scale=1.0+np.sin(x*y)\n",
    "            # scale=1.0\n",
    "            for i in range(mymesh.nodesperelement):\n",
    "                iInd=elconn[i]\n",
    "                # assemble to global F\n",
    "                R[iInd]+=scale*sigma(phi)*(gradphi[0]*shp_grad[i,0]+gradphi[1]*shp_grad[i,1])*JxW\\\n",
    "                        +f(phi)*shp_val[i]*JxW\n",
    "                # now we can calculate the K matrix\n",
    "                for j in range(mymesh.nodesperelement):\n",
    "                    jInd=elconn[j]\n",
    "                    K[iInd,jInd]+=(-1.0)*scale*dsigma(phi)*shp_val[j]*(gradphi[0]*shp_grad[i,0]+gradphi[1]*shp_grad[i,1])*JxW\\\n",
    "                                 +(-1.0)*scale*sigma(phi)*(shp_grad[j,0]*shp_grad[i,0]+shp_grad[j,1]*shp_grad[i,1])*JxW\\\n",
    "                                 +(-1.0)*df(phi)*shp_val[j]*shp_val[i]*JxW\n",
    "    ### end-of-element-loop\n",
    "    \n",
    "    ###################################################################\n",
    "    ### now we apply the boundary condition\n",
    "    ### for dirichlet bc, please keep in mind, you are within a\n",
    "    ### NR iteration, your dx should remain zero when the \n",
    "    ### dof is a 'dirichlet boundary dof'(because x=x+dx) !\n",
    "    ###################################################################\n",
    "    # here we fix all the boundary 'phi' to a constant value\n",
    "    iInd=mymesh.bcnodeids['left']\n",
    "    K[iInd,iInd]+=penalty\n",
    "    R[iInd]=0.0\n",
    "    # iInd=mymesh.bcnodeids['right']\n",
    "    # K[iInd,iInd]+=penalty\n",
    "    # R[iInd]=0.0\n",
    "    # iInd=mymesh.bcnodeids['bottom']\n",
    "    # K[iInd,iInd]+=penalty\n",
    "    # R[iInd]=0.0\n",
    "    # iInd=mymesh.bcnodeids['top']\n",
    "    # K[iInd,iInd]+=penalty\n",
    "    # R[iInd]=0.0\n",
    "    \n",
    "    # now all your K and R are assembled, you can solve them\n",
    "    dU=np.linalg.solve(K,R)\n",
    "    \n",
    "    # update info\n",
    "    iters+=1\n",
    "    U+=dU\n",
    "    rnorm=np.linalg.norm(R)\n",
    "    dunorm=np.linalg.norm(dU)\n",
    "    enorm=rnorm*dunorm\n",
    "    if iters==1:\n",
    "        rnorm0=rnorm\n",
    "        dunorm0=dunorm\n",
    "        enorm0=enorm\n",
    "    print('iters=%3d, |R0 |=%14.5e, |R |=%14.5e'%(iters,rnorm0,rnorm))\n",
    "    print('           |dU0|=%14.5e, |dU|=%14.5e'%(dunorm0,dunorm))\n",
    "    print('           |E0 |=%14.5e, |E |=%14.5e'%(enorm0,enorm))\n",
    "    \n",
    "    # convergence check\n",
    "    if rnorm<tolerance or rnorm<reltol*rnorm0:\n",
    "        converged=True\n",
    "        break\n",
    "### end-of-NR-loop\n",
    "result=ResultIO('nlpoisson')\n",
    "result.save2vtu(mesh=mymesh,solution=U,varnamelist=['phi'],step=0)"
   ]
  }
 ],
 "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}