{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 1.6 Error estimation & adaptive refinement\n", "\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" ] }, { "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", " base = Rectangle(1, 0.6).Face()\n", " chip = MoveTo(0.5,0.15).Line(0.15,0.15).Line(-0.15,0.15).Line(-0.15,-0.15).Close().Face()\n", " top = MoveTo(0.2,0.6).Rectangle(0.6,0.2).Face()\n", " base -= chip\n", "\n", " base.faces.name=\"base\"\n", " chip.faces.name=\"chip\"\n", " chip.faces.col=(1,0,0)\n", " top.faces.name=\"top\"\n", " geo = Glue([base,chip,top])\n", " geo.edges.name=\"default\"\n", " geo.edges.Min(Y).name=\"bot\"\n", " return OCCGeometry(geo, dim=2)\n", "\n", "mesh = Mesh(MakeGeometryOCC().GenerateMesh(maxh=0.2))\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 \\lambda \\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$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fes = H1(mesh, order=3, dirichlet=[1])\n", "u, v = fes.TnT()\n", "\n", "# one heat conductivity coefficient per sub-domain\n", "lam = CoefficientFunction([1, 1000, 10])\n", "a = BilinearForm(lam*grad(u)*grad(v)*dx)\n", "\n", "# heat-source in inner subdomain\n", "f = LinearForm(fes)\n", "f = LinearForm(1*v*dx(definedon=\"chip\"))\n", "\n", "c = Preconditioner(a, type=\"multigrid\", inverse=\"sparsecholesky\")\n", "\n", "gfu = GridFunction(fes)" ] }, { "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():\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" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "SolveBVP()\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", "gf_flux = GridFunction(space_flux, \"flux\")\n", "\n", "flux = lam * 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}{\\lambda} \n", "|\\lambda \\nabla u_h - I_h(\\lambda \\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/lam*(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", "print(eta2)" ] }, { "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", " # see below for vectorized alternative" ] }, { "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()\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 = lam * grad(gfu) \n", " gf_flux.Set(flux) \n", " \n", " # compute estimator:\n", " err = 1/lam*(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" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "CalcError()\n", "mesh.Refine()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run the adaptive loop" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "level = 0 \n", "while fes.ndof < 50000: \n", " SolveBVP()\n", " level = level + 1\n", " if level%5 == 0:\n", " print('adaptive step #', level)\n", " Draw(gfu)\n", " CalcError()\n", " mesh.Refine()" ] }, { "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.ion()\n", "plt.show()" ] } ], "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.11.3" } }, "nbformat": 4, "nbformat_minor": 4 }