Error estimation for SANS

The research projects I work on are all centered around a code called SANS and methodology for simulating field equations. SANS stands for Solution-Adaptive Numerical Simulator. The crux of the solution-adaptivity part is an algorithm called MOESS which was developed by Prof. Masayuki Yano of University of Toronto, who worked with my advisor back in the day. MOESS was developed to well approximate the change in error that would be induced by a change in a mesh, without costing an intractible amount of computational effort. If you're following, this means you can try to find the mesh that minimizes the error-- instead of throwing money (or computational cores) at a PDE, you can try to finesse your way to a better answer for what you're interested in. At least, that's the theory.

Discretization of conservation laws and engineering quantities of interest

The point of this post is to freshen up my math on how the error is formulated to hand into MOESS. The solution methodolgy in SANS is all about solving the so-called weak formulation of the conservation laws that dictate physical systems of the type we are interested in– typically fluid (and particularly aerodynamic) flows. The conservation law can be written in a residual form:
$$
\mathcal{R}(u, v)= 0 \quad \quad \forall v \in \mathcal{V}
$$
which holds on a domain $\Omega$ and where $\mathcal{V}$ is the sufficiently smooth space of possible solutions^[1] to the conservation law, and $u \in \mathcal{V}$.

When we do simulations of these type of systems, we typically are looking for a specific engineering quantity of interest, like the drag around an airfoil. We call the generic engineering quantity of interest an output:
$$
J= \mathcal{J}(u)
$$
with the form shown here that indicates the dependence of the output on the whole solution that results. Of course the complete solution of the simulation is good to have, but strictly speaking we have a goal in mind: calculating the value of the engineering output.

In order to solve the weak form residual of the governing equation, we seek to solve a discretized form:
$$
\mathcal{R}_{h, p}(u_{h, p}, v_{h, p})= 0 \quad \quad \forall v_{h, p} \in \mathcal{V}_{h, p}
$$
where here, $\mathcal{V}_{h, p}$ is a finite element space, and $\mathcal{R}_{h, p}$ is a discretized form of the governing equations, which includes a choice of discretizion method (like discontinuous Galerkin with a particular type of stabilization) that is satisfactorially designed^[2]. This of course implies an approximation to the engineering quantity of interest:
$$
J_{h, p}= \mathcal{J}_{h, p}(u_{h, p})
$$
where $\mathcal{J}_{h, p}:\mathcal{V}_{h, p} \to \mathbb{R}$ and is appropriately formulated^[2:1].

Error estimation and the DWR method

We want to estimate the error in the discretized output vs. the true output:
$$
\varepsilon_\mathrm{true}= J - J_{h, p}= \mathcal{J}(u) - \mathcal{J}_{h, p}(u_{h, p})
$$
The dual-weighted residual (DWR) method^[3] allows you to express that quantity by:
$$
\varepsilon_\mathrm{true}= \mathcal{J}(u) - \mathcal{J}_{h, p}(u_{h, p})= -\mathcal{R}_{h, p}(u_{h, p}, \psi)
$$
where $\psi \in \mathcal{W} \equiv \mathcal{V} \oplus \mathcal{V}_{h, p}$ is the "adjoint" solution satisfying a transpose equation that describes the sensitivities. By Galerkin orthoganality, we can write:
$$
\varepsilon_\mathrm{true}= -\mathcal{R}_{h, p}(u_{h, p}, \psi - v_{h, p}) \quad \quad \forall v_{h, p}
$$
Of course, $\psi$ is not in a computable space, so we approximate using
$$
\varepsilon_\mathrm{true} \approx -\mathcal{R}_{h, p}(u_{h, p}, \psi_{h, \hat{p}})
$$
with an enriched space $\hat{p}= p + 1$, which puts $\psi_{h, p} \in V_{h, \hat{p}} \subset W$. Masa^[1:1] says it better than I could:

The linearization and the adjoint-approximation errors may be significant,
especially on coarse meshes; however, in practice, the error estimate is
sufficiently accurate for the purpose of mesh adaptation.

From here, we can describe the local error estimate-- the error associated with an individual element-- on element $\kappa$ by:
$$
\eta_\kappa= |\mathcal{R}_{h, p}(u_{h, p}, \psi_{h, \hat{p}}|_\kappa)|
$$
Neato, from here, we hand it off to Masa and MOESS to optimize and Philip's avro remeshing code to get a new grid.

I can talk about those some other time.