3.2. 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{ESIMATE}\rightarrow \text{MARK}\rightarrow \text{REFINE}\rightarrow \text{SOLVE} \rightarrow \ldots \)$

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
import matplotlib.pyplot as plt
import numpy as np

3.2.1. Geometry

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

def MakeGeometryOCC():
    square = Circle( (0,0), 2).Face() # create a rectangle with boundary conditions
    square.edges.name = "outer"
    square.faces.name = "air"

    el1 = MoveTo(-0.4, 0.2).Rectangle( 0.8, 0.1).Face() # create a rectangle with boundary conditions
    el1.edges.name = "el1"
    el1.vertices.name = "el1"

    el2 = MoveTo(-0.4, -0.2).Rectangle( 0.8, 0.1).Face() # create a rectangle with boundary conditions
    el2.edges.name = "el2"
    el2.vertices.name = "el2"

    dielec = MoveTo(-0.9, -0.1).Rectangle( 1.8, 0.3).Face()
    dielec.faces.name = "dielec"

    air = square  # subtract the rectangles from the air rectangle
    shape = Glue([air - dielec, dielec])
    shape =shape - el1 - el2

    ### adding extra specifications on the shape
    #predefined mesh size for the shape


    return OCCGeometry(shape, dim=2)

mesh = Mesh(MakeGeometryOCC().GenerateMesh(maxh=0.5))
mesh.Curve(3)
Draw(mesh);

3.2.2. Spaces & forms

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

\[ \int_\Omega \varepsilon \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\).

On the notes adaptivity is explained using \(\lambda\) instead of \(\varepsilon\).

fes = H1(mesh, order=3, dirichlet="el.*")
u, v = fes.TnT()

# one heat conductivity coefficient per sub-domain
eps = CoefficientFunction([1, 2])
a = BilinearForm(eps*grad(u)*grad(v)*dx)

# heat-source in inner subdomain
f = LinearForm(fes)

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

electrode = mesh.BoundaryCF({"el1":1, "el2":-1}, default = 0 ) # define the boundary conditions
gfu = GridFunction(fes)
gfu.Set(electrode, definedon=mesh.Boundaries("el.*")) # set the boundary conditions

Note that the linear system is not yet assembled above.

3.2.3. Solve

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

def SolveBVP(gfu):
    fes.Update()
    gfu.Update()
    a.Assemble()
    f.Assemble()
    inv = CGSolver(a.mat, c.mat)
    gfu.vec.data +=inv* (f.vec - a.mat * gfu.vec) # solve the system

SolveBVP(gfu)
Draw(gfu);

3.2.4. 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.

space_flux = HDiv(mesh, order=2)
#space_flux = VectorH1(mesh, order=2)

gf_flux = GridFunction(space_flux, "flux")

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

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

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

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

err = 1/eps*(flux-gf_flux)*(flux-gf_flux)
Draw(err, mesh, 'error_representation')

BaseWebGuiScene

eta2 = Integrate(err, mesh, VOL, element_wise=True)

# print only the first 10 values
print(np.array(eta2)[:10], "...")

[4.21842882e-08 3.51994462e-08 1.15858338e-07 9.16675035e-07
 8.56871323e-08 4.69599389e-07 1.35455371e-07 1.08785321e-06
 3.81165502e-08 2.91867924e-07] ...

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

3.2.5. Mark

We mark elements with large error estimator for refinement.

maxerr = max(eta2)
print ("maxerr = ", maxerr)

for el in mesh.Elements():
    mesh.SetRefinementFlag(el, eta2[el.nr] > 0.25*maxerr)

maxerr =  0.17618067290069273

3.2.6. Refine & solve again

Refine marked elements:

mesh.Refine()
SolveBVP(gfu)
Draw(gfu);

3.2.7. Automate the above steps

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

def CalcError():

    # compute the flux:
    space_flux.Update()      
    gf_flux.Update()
    flux = eps * grad(gfu)        
    gf_flux.Set(flux) 
    
    # compute estimator:
    err = 1/eps*(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 (vectorized alternative)
    mesh.ngmesh.Elements2D().NumPy()["refine"] = eta2.NumPy() > 0.25*maxerr
    

CalcError()
mesh.Refine()
mesh.Curve(4);

ndof = 1652  maxerr = 0.09112441450626464

3.2.8. Run the adaptive loop

level = 0 
while fes.ndof < 100000:  
    SolveBVP(gfu)
    level = level + 1
    if level%5 == 0:
        print('adaptive step #', level)
        Draw(gfu)
    CalcError()
    mesh.Refine()
    mesh.Curve(4)

SolveBVP(gfu)

ndof = 2261  maxerr = 0.04002281240824464
ndof = 3026  maxerr = 0.016862384250682788
ndof = 3884  maxerr = 0.007164702272429664
ndof = 4763  maxerr = 0.0030597661875181887
adaptive step # 5

ndof = 5519  maxerr = 0.0013107177167427295

ndof = 6428  maxerr = 0.0005624494605382381

ndof = 7277  maxerr = 0.0002416603668252826
ndof = 8507  maxerr = 0.00010387514478162288

ndof = 9851  maxerr = 4.469182775061902e-05

adaptive step # 10

ndof = 11477  maxerr = 1.9224574263818197e-05
ndof = 12926  maxerr = 8.273945688343705e-06

ndof = 14885  maxerr = 3.559813510266537e-06
ndof = 17162  maxerr = 1.5322545652916487e-06

ndof = 20009  maxerr = 9.539191899192556e-07
adaptive step # 15

ndof = 21641  maxerr = 3.8866562381942105e-07

ndof = 26072  maxerr = 1.673313129511147e-07

ndof = 31460  maxerr = 6.298646211711452e-08

ndof = 38234  maxerr = 2.5169140936195947e-08

ndof = 45122  maxerr = 1.049478136859209e-08

adaptive step # 20

ndof = 53699  maxerr = 4.3900170939369345e-09

ndof = 64940  maxerr = 1.8902728276788e-09

ndof = 78116  maxerr = 8.385455776386515e-10

ndof = 93812  maxerr = 4.85783536246917e-10

ndof = 105272  maxerr = 1.7098027937699436e-10

3.2.9. Plot history of adaptive convergence

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

plt.plot(ndof, [(ndof[0]/n)**(1.5)for n in ndof], "--", label="3/2-th order")

# plot the convergence rate line

plt.legend()
plt.ion()
plt.show()

N = 10
# fac = 0 if mesh.dim == 2 else 1
# 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)]
p = [(i / N, -1 + (2*j) / N, 0) for i in range(N+1) for j in range(N+1)]

fieldlines = flux._BuildFieldLines(
    mesh, p, num_fieldlines=N**3//5, randomized=False, length=0.3)
clipping = {"clipping": {"y": 0, "z": -1}}


Draw(-eps*grad(gfu), mesh, "X" , min = 0, max = 10);
Draw(-eps*grad(gfu), mesh,  "X", draw_vol=True, draw_surf=True, objects=[fieldlines],
     autoscale=True, min=0, max=10, settings={"Objects": {"Surface": False}}, **clipping);