Goal driven error estimates

Contents

28. Goal driven error estimates#

The above error estimators estimate the error in the energy norm \(V\). Some applications require to compute certain values (such as point values, average values, line integrals, fluxes through surfaces, …). These values are described by linear functionals \(b : V \rightarrow {\mathbb R}\). We want to design a method such that the error in this goal, i.e.,

\[ b(u) - b(u_h) \]

is small. The technique is to solve additionally the dual problem, where the right hand side is the goal functional:

\[ \mbox{Find } w \in V : \qquad A(v, w) = b(v) \qquad \forall \, v \in V. \]

Usually, one cannot solve the dual problem either, and one applies a Galerkin method also for the dual problem:

\[ \mbox{Find } w_h \in V_h : \qquad A(v_h, w_h) = b(v_h) \qquad \forall \, v_h \in V_h. \]

In the case of point values, the solution of the dual problem is the Green function (which is not in \(H^1\)). The error in the goal is

\[ b(u-u_h) = A(u-u_h, w) = A(u-u_h, w-w_h). \]

A rigorous upper bound for the error in the goal is obtained by using continuity of the bilinear-form, and energy error estimates \(\eta^1\) and \(\eta^2\) for the primal and dual problem, respectively:

\[ | b(u-u_h) | \preceq \| u - u_h \|_V \| w - w_h \|_V \preceq \eta^1(u_h, f) \, \eta^2(w_h, b). \]

A good heuristic is the following (unfortunately, not correct) estimate

\[\begin{align*} b(u-u_h) & = A(u-u_h, w-w_h) \\ & \preceq \sum_{T \in {\cal T} } \| u - u_h \|_{H^1(T)} \, \| w - w_h \|_{H^1(T)} \\ & \preceq \sum_{T} \eta_T^1 (u_h, f) \, \eta_T^2 (w_h, b) \end{align*}\]

The last step would require a local reliability estimate. But, this is not true.

We can interpret the last sum that way: The local estimators \(\eta^2_T(w_h)\) provide a way for weighting the primal local estimators according to the desired goal.

28.1. Example#

On \(\Omega = (0,1)^2-(0.4,0.4)^2\), consider the problem: find \(u \in H_0^1(\Omega)\) such that

\[ \int \lambda \nabla u \nabla v \, dx = \int_{(2,4)\times (6,8)} v \, dx \]

We are interested in the point-value \(u(0.7,0.3)\), i.e. the evaluation functional is

\[ b(v) = v(0.7,0.3) \]

Note that this point-evaluation is not a continuous linear-form on \(H^1\).

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw

rect = Rectangle(1,1).Face()
rect.edges.name="dirichlet"
inner = MoveTo(0.2,0.6).Rectangle(0.2,0.2).Face()
hole = Rectangle(0.4,0.4).Face()
outer = rect-inner-hole
inner.faces.name="source"
shape = Glue([outer,inner])
mesh = Mesh(OCCGeometry(shape,dim=2).GenerateMesh(maxh=0.1))
Draw (mesh);

fes = H1(mesh,order=3, autoupdate=True, dirichlet="dirichlet")
u,v = fes.TnT()
lam = mesh.MaterialCF( { "inner" : 5 }, default=1)

a = BilinearForm(lam*grad(u)*grad(v)*dx)
f = LinearForm(100*v*dx("source"))
b = LinearForm(fes)
b += v(0.7,0.3) 

gfu = GridFunction(fes, autoupdate=True)
gfudual = GridFunction(fes, autoupdate=True)
fesflux = HDiv(mesh,order=3, autoupdate=True)
gfflux = GridFunction(fesflux, autoupdate=True)

hist = []


def SolveEstimateMark():
    a.Assemble()
    f.Assemble()
    b.Assemble()
    gfu.vec.data = a.mat.Inverse(fes.FreeDofs())*f.vec
    gfudual.vec.data = a.mat.Inverse(fes.FreeDofs())*b.vec

    gfflux.Set(lam*grad(gfu))
    errest = Integrate ( lam*(1/lam*gfflux-grad(gfu))**2, mesh, element_wise=True)
    gfflux.Set(lam*grad(gfudual))
    errestdual = Integrate ( lam*(1/lam*gfflux-grad(gfudual))**2, mesh, element_wise=True)

    for i in range(len(errest)):
        errest[i] = sqrt(errest[i]*errestdual[i])

    errmax = max(errest)
    hist.append ( (fes.ndof, sum(errest), InnerProduct(f.vec, gfudual.vec )) )
    for i in range(mesh.GetNE(VOL)):
        mesh.SetRefinementFlag(ElementId(VOL,i), errest[i] > 0.25*errmax)

SolveEstimateMark()
while fes.ndof < 10000:
    mesh.Refine()
    SolveEstimateMark()

Draw (mesh);

copy mesh complete

Draw (gfu)
Draw (gfudual);

We obtain mesh-refinement which combines refinement due to singularities at corners, and taking the target functional into account !

Contents