This page was generated from unit-1.6-adaptivity/adaptivity.ipynb.

1.6 Error estimation & adaptive refinement

In this tutorial, we apply a Zienkiewicz-Zhu type error estimator and run an adaptive loop with these steps:

\[\text{Solve}\rightarrow \text{Estimate}\rightarrow \text{Mark}\rightarrow \text{Refine}\rightarrow \text{Solve} \rightarrow \ldots\]
[1]:
from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import SplineGeometry
import matplotlib.pyplot as plt

Geometry

The following geometry represents a heated chip embedded in another material that conducts away the heat.

[2]:
#   point numbers 0, 1, ... 11
#   sub-domain numbers (1), (2), (3)
#
#
#             7-------------6
#             |             |
#             |     (2)     |
#             |             |
#      3------4-------------5------2
#      |                           |
#      |             11            |
#      |           /   \           |
#      |         10 (3) 9          |
#      |           \   /     (1)   |
#      |             8             |
#      |                           |
#      0---------------------------1
#

def MakeGeometry():
    geometry = SplineGeometry()

    # point coordinates ...
    pnts = [ (0,0), (1,0), (1,0.6), (0,0.6), \
             (0.2,0.6), (0.8,0.6), (0.8,0.8), (0.2,0.8), \
             (0.5,0.15), (0.65,0.3), (0.5,0.45), (0.35,0.3) ]
    pnums = [geometry.AppendPoint(*p) for p in pnts]

    # start-point, end-point, boundary-condition, left-domain, right-domain:
    lines = [ (0,1,1,1,0), (1,2,2,1,0), (2,5,2,1,0), (5,4,2,1,2), (4,3,2,1,0), (3,0,2,1,0), \
              (5,6,2,2,0), (6,7,2,2,0), (7,4,2,2,0), \
              (8,9,2,3,1), (9,10,2,3,1), (10,11,2,3,1), (11,8,2,3,1) ]

    for p1,p2,bc,left,right in lines:
        geometry.Append(["line", pnums[p1], pnums[p2]], bc=bc, leftdomain=left, rightdomain=right)

    geometry.SetMaterial(1,"base")
    geometry.SetMaterial(2,"chip")
    geometry.SetMaterial(3,"top")

    return geometry

mesh = Mesh(MakeGeometry().GenerateMesh(maxh=0.2))

Spaces & forms

The problem is to find \(u\) in \(H_{0,D}^1\) satisfying

\[\int_\Omega \lambda \nabla u \cdot \nabla v = \int_\Omega f v\]

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\).

[3]:
fes = H1(mesh, order=3, dirichlet=[1])
u, v = fes.TnT()

# one heat conductivity coefficient per sub-domain
lam = CoefficientFunction([1, 1000, 10])
a = BilinearForm(fes)
a += lam*grad(u)*grad(v)*dx

# heat-source in inner subdomain
f = LinearForm(fes)
f += CoefficientFunction([0, 0, 1])*v * dx

c = Preconditioner(a, type="multigrid", inverse="sparsecholesky")

gfu = GridFunction(fes)
Draw (gfu)
[3]:
BaseWebGuiScene

Note that the linear system is not yet assembled above.

Solve

Since we must solve multiple times, we define a function to solve the boundary value problem, where assembly, update, and solve occurs.

[4]:
def SolveBVP():
    fes.Update()
    gfu.Update()
    a.Assemble()
    f.Assemble()
    inv = CGSolver(a.mat, c.mat)
    gfu.vec.data = inv * f.vec
[5]:
SolveBVP()
Draw(gfu)
[5]:
BaseWebGuiScene

Estimate

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.

[6]:
space_flux = HDiv(mesh, order=2)
gf_flux = GridFunction(space_flux, "flux")

flux = lam * grad(gfu)
gf_flux.Set(flux)

Element-wise error estimator: On each element \(T\), set

\[\eta_T^2 = \int_T \frac{1}{\lambda} |\lambda \nabla u_h - I_h(\lambda \nabla u_h) |^2\]

where \(u_h\) is the computed solution gfu and \(I_h\) is the interpolation performed by Set in NGSolve.

[7]:
err = 1/lam*(flux-gf_flux)*(flux-gf_flux)
Draw(err, mesh, 'error_representation')
[7]:
BaseWebGuiScene
[8]:
eta2 = Integrate(err, mesh, VOL, element_wise=True)
print(eta2)
 4.77713e-10
 3.62503e-08
 4.73264e-06
 1.39888e-10
 2.56218e-10
 1.68139e-08
 3.89785e-08
 1.7684e-07
 1.2484e-07
 2.62242e-07
 1.30604e-07
 6.66851e-09
 4.49293e-07
 3.21719e-06
 6.65763e-07
 2.49744e-08
 2.56941e-07
 4.0863e-06
 1.39265e-07
 5.13219e-07
 9.93691e-08
 1.22606e-06
 4.89833e-10
 1.05638e-06
 5.005e-07
 1.03982e-06
 9.39242e-10
 1.62172e-06
 3.01971e-08
 5.35648e-07
 4.09083e-09
 2.08152e-07
 1.12974e-10
 1.18657e-11
 1.58514e-10
 1.45496e-11
 1.91237e-11
 2.16696e-12
 8.31556e-07
 8.72613e-07

The above values, one per element, lead us to identify elements which might have large error.

Mark

We mark elements with large error estimator for refinement.

[9]:
maxerr = max(eta2)
print ("maxerr = ", maxerr)

for el in mesh.Elements():
    mesh.SetRefinementFlag(el, eta2[el.nr] > 0.25*maxerr)
maxerr =  4.73263819311393e-06

Refine & solve again

Refine marked elements:

[10]:
mesh.Refine()
SolveBVP()
Draw(gfu)
[10]:
BaseWebGuiScene

Automate the above steps

[11]:
l = []    # l = list of estimated total error

def CalcError():

    # compute the flux:
    space_flux.Update()
    gf_flux.Update()
    flux = lam * grad(gfu)
    gf_flux.Set(flux)

    # compute estimator:
    err = 1/lam*(flux-gf_flux)*(flux-gf_flux)
    eta2 = Integrate(err, mesh, VOL, element_wise=True)
    maxerr = max(eta2)
    l.append ((fes.ndof, sqrt(sum(eta2))))
    print("ndof =", fes.ndof, " maxerr =", maxerr)

    # mark for refinement:
    for el in mesh.Elements():
        mesh.SetRefinementFlag(el, eta2[el.nr] > 0.25*maxerr)
[12]:
CalcError()
mesh.Refine()
ndof = 454  maxerr = 2.7688817994871426e-06

Run the adaptive loop

[13]:
while fes.ndof < 50000:
    SolveBVP()
    Draw(gfu)
    CalcError()
    mesh.Refine()
ndof = 547  maxerr = 9.146694004280255e-07
ndof = 1168  maxerr = 3.3095640452852995e-07
ndof = 1984  maxerr = 1.1989055931850714e-07
ndof = 2890  maxerr = 4.3458640874366095e-08
ndof = 3760  maxerr = 1.8502847107724447e-08
ndof = 4567  maxerr = 9.244842649021228e-09
ndof = 5176  maxerr = 4.620482787246551e-09
ndof = 5821  maxerr = 2.3091632776959873e-09
ndof = 6493  maxerr = 1.1540222393413907e-09
ndof = 7120  maxerr = 5.767311562199303e-10
ndof = 7828  maxerr = 2.882221612084458e-10
ndof = 8623  maxerr = 1.4405135479307996e-10
ndof = 9478  maxerr = 7.199432111938219e-11
ndof = 10426  maxerr = 3.598117400058056e-11
ndof = 11533  maxerr = 1.7982714291509452e-11
ndof = 12796  maxerr = 8.987384603538555e-12
ndof = 14356  maxerr = 4.491721447675487e-12
ndof = 16321  maxerr = 2.2448916112611e-12
ndof = 18976  maxerr = 1.1219587713697661e-12
ndof = 21676  maxerr = 5.607214905263319e-13
ndof = 25231  maxerr = 2.8023556039347864e-13
ndof = 28792  maxerr = 1.40053937563383e-13
ndof = 33070  maxerr = 6.99947635297359e-14
ndof = 39082  maxerr = 3.4981400215296833e-14
ndof = 45721  maxerr = 1.7482171335208574e-14
ndof = 52705  maxerr = 8.736925640477952e-15

Plot history of adaptive convergence

[14]:
plt.yscale('log')
plt.xscale('log')
plt.xlabel("ndof")
plt.ylabel("H1 error-estimate")
ndof,err = zip(*l)
plt.plot(ndof,err, "-*")

plt.ion()
plt.show()
../../_images/i-tutorials_unit-1.6-adaptivity_adaptivity_24_0.png
[ ]: