NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(complementary:amr)=\n",
    "# Error estimation & adaptive refinement\n",
    "\n",
    "In this tutorial, we apply a Zienkiewicz-Zhu type error estimator and run an adaptive loop with these steps:\n",
    "$$\n",
    "\\text{SOLVE}\\rightarrow\n",
    "\\text{ESIMATE}\\rightarrow\n",
    "\\text{MARK}\\rightarrow\n",
    "\\text{REFINE}\\rightarrow\n",
    "\\text{SOLVE} \\rightarrow \\ldots\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometry\n",
    "\n",
    "The following geometry represents a heated chip embedded in another material that conducts away the heat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MakeGeometryOCC():\n",
    "    square = Circle( (0,0), 2).Face() # create a rectangle with boundary conditions\n",
    "    square.edges.name = \"outer\"\n",
    "    square.faces.name = \"air\"\n",
    "\n",
    "    el1 = MoveTo(-0.4, 0.2).Rectangle( 0.8, 0.1).Face() # create a rectangle with boundary conditions\n",
    "    el1.edges.name = \"el1\"\n",
    "    el1.vertices.name = \"el1\"\n",
    "\n",
    "    el2 = MoveTo(-0.4, -0.2).Rectangle( 0.8, 0.1).Face() # create a rectangle with boundary conditions\n",
    "    el2.edges.name = \"el2\"\n",
    "    el2.vertices.name = \"el2\"\n",
    "\n",
    "    dielec = MoveTo(-0.9, -0.1).Rectangle( 1.8, 0.3).Face()\n",
    "    dielec.faces.name = \"dielec\"\n",
    "\n",
    "    air = square  # subtract the rectangles from the air rectangle\n",
    "    shape = Glue([air - dielec, dielec])\n",
    "    shape =shape - el1 - el2\n",
    "\n",
    "    ### adding extra specifications on the shape\n",
    "    #predefined mesh size for the shape\n",
    "\n",
    "\n",
    "    return OCCGeometry(shape, dim=2)\n",
    "\n",
    "mesh = Mesh(MakeGeometryOCC().GenerateMesh(maxh=0.5))\n",
    "mesh.Curve(3)\n",
    "Draw(mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spaces & forms\n",
    "\n",
    "The problem is to find $u$ in $H_{0,D}^1$ satisfying \n",
    "\n",
    "$$\n",
    "\\int_\\Omega \\varepsilon \\nabla u \\cdot \\nabla v = \\int_\\Omega f v \n",
    "$$\n",
    "\n",
    "for all $v$ in $H_{0,D}^1$. We expect the solution to have singularities due to the nonconvex re-enrant angles and discontinuities in $\\lambda$.\n",
    "\n",
    "**On the notes adaptivity is explained using $\\lambda$ instead of $\\varepsilon$.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=3, dirichlet=\"el.*\")\n",
    "u, v = fes.TnT()\n",
    "\n",
    "# one heat conductivity coefficient per sub-domain\n",
    "eps = CoefficientFunction([1, 2])\n",
    "a = BilinearForm(eps*grad(u)*grad(v)*dx)\n",
    "\n",
    "# heat-source in inner subdomain\n",
    "f = LinearForm(fes)\n",
    "\n",
    "c = Preconditioner(a, type=\"multigrid\", inverse=\"sparsecholesky\")\n",
    "\n",
    "electrode = mesh.BoundaryCF({\"el1\":1, \"el2\":-1}, default = 0 ) # define the boundary conditions\n",
    "gfu = GridFunction(fes)\n",
    "gfu.Set(electrode, definedon=mesh.Boundaries(\"el.*\")) # set the boundary conditions\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the linear system is not yet assembled above.\n",
    "\n",
    "### Solve \n",
    "\n",
    "Since we must solve multiple times, we define a function to solve the boundary value problem, where assembly, update, and solve occurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveBVP(gfu):\n",
    "    fes.Update()\n",
    "    gfu.Update()\n",
    "    a.Assemble()\n",
    "    f.Assemble()\n",
    "    inv = CGSolver(a.mat, c.mat)\n",
    "    gfu.vec.data +=inv* (f.vec - a.mat * gfu.vec) # solve the system\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SolveBVP(gfu)\n",
    "Draw(gfu);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate\n",
    "\n",
    "We implement a gradient-recovery-type error estimator. For this, we need an H(div) space for flux recovery. We must compute the flux  of the computed solution and interpolate it into this H(div) space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "space_flux = HDiv(mesh, order=2)\n",
    "#space_flux = VectorH1(mesh, order=2)\n",
    "\n",
    "gf_flux = GridFunction(space_flux, \"flux\")\n",
    "\n",
    "flux = eps * grad(gfu)\n",
    "gf_flux.Set(flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Element-wise error estimator:** On each element $T$, set \n",
    "\n",
    "$$\n",
    "\\eta_T^2 = \\int_T \\frac{1}{\\varepsilon} \n",
    "|\\varepsilon \\nabla u_h - I_h(\\varepsilon \\nabla u_h) |^2\n",
    "$$\n",
    "\n",
    "where $u_h$ is the computed solution `gfu` and $I_h$ is the interpolation performed by `Set` in NGSolve.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "err = 1/eps*(flux-gf_flux)*(flux-gf_flux)\n",
    "Draw(err, mesh, 'error_representation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eta2 = Integrate(err, mesh, VOL, element_wise=True)\n",
    "\n",
    "# print only the first 10 values\n",
    "print(np.array(eta2)[:10], \"...\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above values, one per element, lead us to identify elements which might have large error.\n",
    "\n",
    "\n",
    "### Mark \n",
    "\n",
    "We mark elements with large error estimator for refinement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maxerr = max(eta2)\n",
    "print (\"maxerr = \", maxerr)\n",
    "\n",
    "for el in mesh.Elements():\n",
    "    mesh.SetRefinementFlag(el, eta2[el.nr] > 0.25*maxerr)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Refine & solve again \n",
    "\n",
    "Refine marked elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.Refine()\n",
    "SolveBVP(gfu)\n",
    "Draw(gfu);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automate the above steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = []    # l = list of estimated total error\n",
    "\n",
    "def CalcError():\n",
    "\n",
    "    # compute the flux:\n",
    "    space_flux.Update()      \n",
    "    gf_flux.Update()\n",
    "    flux = eps * grad(gfu)        \n",
    "    gf_flux.Set(flux) \n",
    "    \n",
    "    # compute estimator:\n",
    "    err = 1/eps*(flux-gf_flux)*(flux-gf_flux)\n",
    "    eta2 = Integrate(err, mesh, VOL, element_wise=True)\n",
    "    maxerr = max(eta2)\n",
    "    l.append ((fes.ndof, sqrt(sum(eta2))))\n",
    "    print(\"ndof =\", fes.ndof, \" maxerr =\", maxerr)\n",
    "    \n",
    "    # mark for refinement (vectorized alternative)\n",
    "    mesh.ngmesh.Elements2D().NumPy()[\"refine\"] = eta2.NumPy() > 0.25*maxerr\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CalcError()\n",
    "mesh.Refine()\n",
    "mesh.Curve(4);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the adaptive loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "level = 0 \n",
    "while fes.ndof < 100000:  \n",
    "    SolveBVP(gfu)\n",
    "    level = level + 1\n",
    "    if level%5 == 0:\n",
    "        print('adaptive step #', level)\n",
    "        Draw(gfu)\n",
    "    CalcError()\n",
    "    mesh.Refine()\n",
    "    mesh.Curve(4)\n",
    "\n",
    "SolveBVP(gfu)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot history of adaptive convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.xlabel(\"ndof\")\n",
    "plt.ylabel(\"H1 error-estimate\")\n",
    "ndof,err = zip(*l)\n",
    "plt.plot(ndof,err, \"-*\")\n",
    "\n",
    "plt.plot(ndof, [(ndof[0]/n)**(1.5)for n in ndof], \"--\", label=\"3/2-th order\")\n",
    "\n",
    "# plot the convergence rate line\n",
    "\n",
    "plt.legend()\n",
    "plt.ion()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10\n",
    "# fac = 0 if mesh.dim == 2 else 1\n",
    "# p = [(-1+4*i/N,-2+4*j/N,fac * 2*k/N) for i in range(1,2*N) for j in range(1,2*N) for k in range(1,N)]\n",
    "p = [(i / N, -1 + (2*j) / N, 0) for i in range(N+1) for j in range(N+1)]\n",
    "\n",
    "fieldlines = flux._BuildFieldLines(\n",
    "    mesh, p, num_fieldlines=N**3//5, randomized=False, length=0.3)\n",
    "clipping = {\"clipping\": {\"y\": 0, \"z\": -1}}\n",
    "\n",
    "\n",
    "Draw(-eps*grad(gfu), mesh, \"X\" , min = 0, max = 10);\n",
    "Draw(-eps*grad(gfu), mesh,  \"X\", draw_vol=True, draw_surf=True, objects=[fieldlines],\n",
    "     autoscale=True, min=0, max=10, settings={\"Objects\": {\"Surface\": False}}, **clipping);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}