# p-adaptation techniques for unstructured grids: a finite volume approach

Filipe Forte Tenreiro, Duarte M. S. Albuquerque

<figure><img src="image_1.png" /></figure>

## To cite this version:

Filipe Forte Tenreiro, Duarte M. S. Albuquerque. p-adaptation techniques for unstructured grids: a finite volume approach. 2023. hal-05017403

HAL Id: hal-05017403

https://inria.hal.science/hal-05017403v1

Preprint submitted on 3 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# *p*-adaptation techniques for unstructured grids: a finite volume approach

Filipe F. Tenreiro, Duarte M. S. Albuquerque*

IDMEC, Technical University of Lisbon, Av. Rovisco Pais, Nº1, 1049-001, Lisbon, Portugal

## Abstract

With the ever-increasing power of computing resources, the development of faster algorithms to keep up with this growth is of the utmost interest to the research community. However, many computational fluid dynamics software packages still lack efficiency when high-accuracy simulations are required. Through numerical error estimation and solution-based order adaptation (or *p*-adaptation), we strive to dramatically improve their efficiency and robustness.

In order to construct an estimate for the criterion used to drive the (adaptive) process, it is first necessary to examine the sources of the error. Irrespective of the discretisation method employed, the fact that profile assumptions are necessary to derive the discretised equations means that the discretisation process inevitably introduces an error unless the underlying solution has a profile that is exactly represented by the assumed one. Provided that both iterative and round-off errors are small enough to where they do not pollute the solution, the primary source of numerical error should be associated with discretisation.

Even though it is the discretisation error that one wishes to reduce through adaptation, we should instead adapt based on the local contributions to this error, i.e., the truncation error (or perhaps an estimate of it). The two proposed methods for estimating the truncation error use a higher-order flux reconstruction about the face. While the higher-order polynomial regression estimator uses a solution that is of higher order of accuracy to approximate the leading term of the truncation error, the tau-extrapolation estimator simply attempts to extrapolate it from the underlying discrete one.

Finally, the automatic error-controlled (adaptive) algorithm uses the previous estimates to refine elements (or faces) in restricted portions of the domain deemed inaccurate. The process is applied to the steady-state form of a general scalar transport equation by choosing a solution *a priori* and then operating the governing equation onto the chosen solution, thereby generating an additional source term that requires no discretisation. Even though this method proved to be tailor-made to handle smooth problems with rapid solution changes, the presence of strong singularities appeared to significantly reduce its effectiveness.

**Keywords:** Finite volume method, Unstructured grids, Higher-order methods, Numerical error estimation, *p*-adaptation

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: filipe.tenreiro@tecnico.ulisboa.pt (Filipe F. Tenreiro), duartealbuquerque@tecnico.ulisboa.pt (Duarte M. S. Albuquerque)</footnote>
1

# 1. Literature review

2

3

4

5

6

7

8

9

10

The predictive capability of computational fluid dynamics (CFD) does not only depend on the validity of the sub-models employed but also on the ability to accurately and reliably estimate and reduce numerical errors. Most verification efforts that are documented for commercial codes seem to be limited to simple benchmark examples that demonstrate “engineering accuracy”, rather than verifying the order of accuracy of the code [1]. Abanto et al. [2] performed order verification studies on three different CFD codes, which were formally said to be second-order accurate. Most tests resulted in either first-order accuracy or non-convergent behaviour. It is our opinion that code users should be aware that commercial software companies are unlikely to perform rigorous code verification studies unless users request it.

11

## 1.1. Adaptive mesh refinement

12

13

14

15

16

17

18

Adaptive mesh refinement (or simply mesh adaptation) refers to the modification of the existing mesh to more accurately capture solution features without excessively increasing the computational cost, i.e., it aims to close the gap between high accuracy and low computational load. Central to any mesh adaptation technique are two essential requirements: a means of recognising where further resolution is needed and a mechanism to actually perform adaptation. However, as Hawken et al. [3] rightfully point out, even though a number of monitors are used in the literature, most of them are often chosen “heuristically”, with “little to no justification” for them.

19

### 1.1.1. Adaptation drivers

20

21

22

23

24

One of the most challenging aspects of mesh adaptation is finding a good criterion to drive the (adaptive) process. The least defensible methods for adaptation are based on either solution features or estimates of the discretisation error (feature-based adaptation). Conversely, the most rigorous ones are based on evaluating or estimating the local truncation error or, equivalently, the residual (truncation error-based adaptation).

25

### 1.1.2. Adaptation methods

26

27

Adaptive mesh refinement methods are generally classified as one of the following: h-, p-, or r-adaptation.

28

29

30

31

32

33

34

35

**Order-based adaptation (p-adaptation).** Order-based adaptation (or p-adaptation) consists of modifying (usually increasing) the local formal order of accuracy of the discretisation scheme without changing the geometry of the grid (see Fig. 1b). Even though this technique has proven accurate for smooth problems, particularly if the partial differential equations (PDEs) are elliptic or parabolic, it is typically not as successful when applied in the presence of discontinuities and/or singularities, where higher-order methods become less effective [4]. Despite its popularity in the finite element (FE) community, not much work has been reported in the finite difference (FD) and finite volume (FV) frameworks.

36

37

38

39

40

41

**Mesh-based adaptation (h- and r-adaptation).** Local grid coarsening and/or refinement (or h-adaptation) alters grid resolution without distortion by introducing and/or removing points hierarchically (see Fig. 1c) [5]. One drawback of using pure h-adaptation, though, is that the quality of the grid is adversely affected at the interface between adapted and unadapted regions. The algorithm proposed by Albuquerque et al. [6] uses a grid interface correction step to provide a smoother and higher-quality grid.
Grid movement (or *r*-adaptation) involves stretching the grid without changing the number of nodes or the connectivity list (see Fig. 1d). Although this technique comes at little to no added cost, care must be taken to ensure that problems due to excessive skewness do not arise. The movement of the nodes can be controlled in various ways [7]; one common technique is to treat the grid as if it were an elastic solid and solve a system of equations (subject to some forcing) that deforms the original grid.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Adaptation methods. Adapted from [8].</figcaption></figure>

A number of researchers have found that a combination of different adaptation methods provides the best approach to reducing the discretisation error. In the finite element method (FEM), *h*- and *p*-adaptation are often combined in what is known as *hp*-adaptive FEs [9].

## 1.2. Higher-order methods

Even though high-order does not have a brick-wall definition, CFD practitioners consent that a given numerical method is said to be of higher order when its order of accuracy is higher than two [10]. High-order is not a new concept, and the number of higher-order methods in the literature is considerable; these are reviewed in [11] and references therein. Although higher-order methods have been around for quite some time, not much attention was paid to them, as in the earlier stages of development, researchers were keen to make first- and second-order accurate methods more efficient. Now, having hit the saturation point in their development, the focus has turned to higher-order methods.

Although it is well known that a second-order accurate method is computationally more expensive than a first-order one, we can only compare method efficiency on the basis of the computational cost required to achieve a given error threshold. Broadly speaking, method efficiency has to do with the balance between accuracy and cost:

$$ \text{Efficiency} = \frac{\text{Accuracy}}{\text{Cost}} = \underbrace{\frac{\text{Accuracy}}{\text{Unknowns} \times \text{Time steps}}}_{\text{Discretisation}} \times \underbrace{\frac{\text{Unknowns} \times \text{Time steps}}{\text{Cost}}}_{\text{Implementation}}. \quad (1) $$

Even though higher-order methods have superior accuracy per unknown, i.e., per degree of freedom (DoF), discretisation efficiency is only half of the picture; implementation efficiency is just as important. The throughput (number of unknowns at a given cost) of higher-order methods can never be higher than that of lower-order ones [12]. Nevertheless, because they possess a higher degree of data reuse, higher-order methods tend to have a higher arithmetic intensity than well-optimised lower-order ones [13], i.e., for every byte of information transferred from main memory to cache, more floating-point operations are made.
## 2. Contributions

The main contributions made during this paper are as follows:

* Proposal of a numerical reconstruction technique whose point-value interpolation is done about (the centroid of) the face.

* Derivation of two distinct a posteriori error estimators.

- The consistency (i.e., asymptotic exactness) of these estimates is verified by evaluating the rate of reduction of the truncation error (or discrete residual) on successively finer grids.

- To relate the discretisation and truncation errors, we use the discrete form of the (discretisation) error transport equation (28).

* Proposal of a novel order-based adaptation (or p-adaptation) method.

- Since the proposed (adaptive) process does not involve changing the geometry of the grid, the error estimators derived in Sections 5.3.1 and 5.3.2 use a higher-order flux reconstruction technique.

- The proposed refinement rules were inspired by those of the same names in Red-Green (h-) refinement [14].

## 3. The governing equation

The conservation equations that govern the transport of mass, momentum, energy, and other specific quantities have a common form embodied in the general scalar transport equation. This equation simply resembles a balance between accumulation, transport, and generation of a generic scalar property whose transport can be dominated by convection (associated with a macroscopic movement of the medium) or by diffusion (associated with molecular diffusion). The changes in this property are described in terms of $\phi$: the unknown of the equation.

The standard form of the transport equation reads

$$ \frac{\partial \phi}{\partial t} + \nabla \cdot (\mathbf{u}\phi) - \nabla \cdot (\Gamma \nabla \phi) = \varphi(\phi), \qquad (2) $$

where $\mathbf{u}$ and $\Gamma$ are the velocity and diffusivity fields, respectively, and $\varphi(\phi)$ is the volumetric source (or sink) contribution to the balance of $\phi$.

By integrating the steady-state form of Eq. (2) over a general fixed control volume (CV), this equation is transformed into

$$ \int_V \nabla \cdot (\mathbf{J}^{\phi,C} - \mathbf{J}^{\phi,D}) dV = \int_V \varphi(\phi) dV, \qquad (3) $$

where $\mathbf{J}^{\phi,C} = \mathbf{u}\phi$ and $\mathbf{J}^{\phi,D} = \Gamma\nabla\phi$ are the convective and diffusive components of the total flux vector $\mathbf{J}^{\phi}$, respectively.

Using the divergence theorem to replace the volume integral with a surface integral, Eq. (3) becomes

$$ \oint_{\partial V} (\mathbf{J}^{\phi,C} - \mathbf{J}^{\phi,D}) \cdot d\mathbf{S} = \int_V \varphi(\phi) dV, \qquad (4) $$
where bold letters denote vectors, (·) the dot product operator, ∂V a closed surface bounding volume V, ∮$_{∂V}$ the surface integral over volume V and S the outward pointing surface vector with associated normal n, such that n dS = dS.

4. Finite volume (FV) discretisation

The standard form of the transport equation and the formal CV integration were described in Section 3. Here, we develop the numerical method based on this integration, the finite volume method (FVM). By working with the two-dimensional (2D) form of Eq. (4), the approximation techniques that are needed to obtain the so-called discretised equations are introduced.

For the sake of brevity, we will only resort to a one-dimensional (1D) (direct) approach to
derive the general form of the two a posteriori error estimators introduced in Section 5.3.

4.1. Flux integration over element faces

Replacing the surface integral by a summation of fluxes over a generic cell C, the left-hand side (LHS) of Eq. (4) becomes

$$
\oint_{\partial C} \mathbf{J}^{\phi} \cdot d\mathbf{S} = \sum_{f \sim \text{faces}(C)} \left( \int_{f} \mathbf{J}^{\phi} \cdot d\mathbf{S} \right). \qquad (5)
$$

No approximations have been made so far, and although Eq. (5) is still formally equivalent to the set of integral equations over all the CVs, the surface fluxes are evaluated at the faces of the cell rather than integrated within it. Using a Gaussian quadrature, the integral about a generic face *f* becomes

$$
\int_f \mathbf{J}^\phi \cdot d\mathbf{S} \approx \sum_{g \sim \mathcal{G}(f)} [\omega_g \mathbf{J}^\phi(\mathbf{x}_g)] \cdot \mathbf{S}_f, \quad (6)
$$
where $g$ refers to an integration point and $\mathcal{G}(f)$ to the subset of integration points along face $f$.

4.2. Source term volume integration

Volume integration is used for the source term. Since the integration domain is a subregion of $\mathbb{R}^2$, the right-hand side (RHS) of Eq. (4) becomes

$$
\begin{align*}
\iint_C \varphi \, dS &= \iint_C \left( \frac{\partial \varphi_x}{\partial x} + \frac{\partial \varphi_y}{\partial y} \right) \, dS \\
&= \oint_{\partial C} (-\varphi_y \, dx + \varphi_x \, dy) \\
&= \sum_{f \sim \text{faces}(C)} \left( \int_f \varphi \cdot d\mathbf{s} \right) \\
&= \sum_{f \sim \text{faces}(C)} \left( \int_f \varphi(\mathbf{r}(t)) \cdot \mathbf{r}'(t) \, dt \right), \tag{7}
\end{align*}
$$
where the vector field and the parametrisation of a generic line segment on the contour of C are respectively given by

$$
\begin{align}
\varphi(x, y) &= -\varphi_y(x, y)\mathbf{i} + \varphi_x(x, y)\mathbf{j} \nonumber \\
&= -\left(v\Phi - \Gamma_y \frac{\partial\Phi}{\partial y}\right)\mathbf{i} + \left(u\Phi - \Gamma_x \frac{\partial\Phi}{\partial x}\right)\mathbf{j}, \quad \text{with } \Phi \equiv \Phi(x, y) \tag{8a}
\end{align}
$$

$$
\mathbf{r}(t) = x(t)\mathbf{i} + y(t)\mathbf{j}. \tag{8b}
$$

4.3. Numerical solution reconstruction

The approximations to the integral in Eq. (6) require the values of $\phi$ and its gradient at loca-
tions other than the computational nodes (CV centroids). Here, since the point-value interpola-
tion is done about the centroid of the face, there is no need to interpolate the cell/vertex-centred
reconstruction to obtain the surface fluxes.

The local approximation of the underlying solution about the centroid of f can be modelled
as a nth-degree polynomial, yielding the general polynomial regression model

$$
\begin{align}
\phi_f(\mathbf{x}) &= \mathbf{c}_f^\top \mathbf{p}_n (\mathbf{x} - \mathbf{x}_f) + \mathbf{r}_f \nonumber \\
&= \mathbf{c}_f^\top \sum_{i=0}^{n} \sum_{j=0}^{n-i} (x - x_f)^i (y - y_f)^j + \mathbf{r}_f, \tag{9}
\end{align}
$$

where $\mathbf{p}_n(\mathbf{x})$ is a 2D polynomial basis function vector of degree $n$ and $\mathbf{r}_f$ is a residual vector.

The polynomial regression model

$$
\phi_f(x_i, y_i) = c_0 + c_1(x_i - x_f) + c_2(y_i - y_f) + \dots + r_{fi}, \text{ with } i = 1, 2, \dots \quad (10)
$$

can be expressed in matrix form (for each face) in terms of a design matrix **D**<sub>f</sub> ∈ ℝ<sup>m×n</sup>, a param-
eter vector **c**<sub>f</sub> ∈ ℝ<sup>n</sup>, a (fitted) response vector **φ**<sub>s</sub> ∈ ℝ<sup>m</sup>, and a residual vector **r**<sub>f</sub> ∈ ℝ<sup>m</sup> as

$$
\left[
\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}
1 & (x_1 - x_f) & (y_1 - y_f) & \cdots & \multicolumn{1}{c}{\left[ \begin{array}{c} c_1 \\ \vdots \\ c_b \end{array} \right]} & \multicolumn{1}{c}{\left[ \begin{array}{c} \phi_1 \\ \vdots \\ \phi_b \end{array} \right]} & \multicolumn{1}{c}{\left[ \begin{array}{c} r_1 \\ \vdots \\ r_b \end{array} \right]} \\
\vdots & \vdots & \vdots & \ddots & & & \\
1 & (x_b - x_f) & (y_b - y_f) & \cdots & \multicolumn{1}{c}{\left[ \begin{array}{c} c_1 \\ \vdots \\ c_b \end{array} \right]} & \multicolumn{1}{c}{\left[ \begin{array}{c} \phi_1 \\ \vdots \\ \phi_b \end{array} \right]} & \multicolumn{1}{c}{\left[ \begin{array}{c} r_1 \\ \vdots \\ r_b \end{array} \right]} \\
0 & n_{b_x} & n_{b_y} & \cdots & \multicolumn{1}{c}{\left[ \begin{array}{c} c_1 \\ \vdots \\ c_b \end{array} \right]} & \multicolumn{1}{c}{\left[ \begin{array}{c} \nabla \phi_b \cdot \mathbf{n}_b \\ \vdots \\ \nabla \phi_b \cdot \mathbf{n}_b \end{array} \right]} & \multicolumn{1}{c}{\left[ \begin{array}{c} r_1 \\ \vdots \\ r_b \end{array} \right]} \\
\multicolumn{2}{c}{\underbrace{n_{b_x} + n_{b_y} \left[ (x_b - x_f) + d_{b_x} \right] \cdot n_{b_x} \quad \left[ (y_b - y_f) + d_{b_y} \right] \cdot n_{b_y} \quad \dots}_{\mathbf{D}_f}} & & & & & \\
\multicolumn{2}{c}{\multicolumn{5}{c}{\mathbf{c}_f}} & \multicolumn{1}{c}{\multicolumn{1}{c}{\mathbf{r}_f}}
\end{array}
\right], \quad (11)
$$

which when using pure matrix notation can be written as

$$
\mathbf{D}_f \mathbf{c}_f = \hat{\mathbf{\Phi}}_s, \tag{12}
$$

with $\hat{\mathbf{\Phi}}_s = \mathbf{\Phi}_s - \mathbf{r}_f$ and

$$
\mathbf{d}_b \equiv \mathbf{d}(\mathbf{x}_b) = \frac{\Gamma(\mathbf{x}_b)}{\mathbf{u}(\mathbf{x}_b)}. \qquad (13)
$$

The last three rows of this system correspond to boundary condition equations. The one
below the red line corresponds to a Dirichlet condition, whereas the ones below the green and
blue lines correspond to Neumann and Robin conditions, respectively. If the discretisation sten-
cil contains a boundary face for which either of these conditions is prescribed, the respective
equation should be written as in Eq. (11).
### 4.3.1. Unconstrained polynomial reconstruction

Consider the previous “overdetermined” (linear) system of equations, where $\mathbf{D}_f$ is assumed to have full (column) rank. If there is no solution to this system because $\boldsymbol{\phi}_s$ does not lie in the range of $\mathbf{D}_f$, then the best approximation is one that minimises the “sum of squares” of the residuals.

In most situations, however, it is desirable to minimise the weighted square residual instead:

$$ \min_{\mathbf{c}_f} \| \sqrt{\mathbf{W}_f} \mathbf{r}_f \|_2^2 = \min_{\mathbf{c}_f} \| \sqrt{\mathbf{W}_f} (\boldsymbol{\phi}_s - \mathbf{D}_f \mathbf{c}_f) \|_2^2, \quad (14) $$

where $\mathbf{W}_f$ is a diagonal matrix with positive weights $w_i$<fnref n="1" />.

As detailed in Appendix A.1.1, practical solution methods do not explicitly form the product $\mathbf{D}_f^\mathrm{T} \mathbf{W}_f \mathbf{D}_f$ to solve the (weighted) unconstrained least-squares (ULS) problem in Eq. (14) for $\mathbf{c}_f$ (recall that the behaviour of the normal equations is governed by $\kappa^2$, not $\kappa$). Instead, we use the QR decomposition of the (weighted) design matrix $\sqrt{\mathbf{W}_f} \mathbf{D}_f^{-2}$ to solve the triangular system $\tilde{\mathbf{R}}_f \mathbf{c}_f^* = \tilde{\mathbf{Q}}_f^\mathrm{T} \sqrt{\mathbf{W}_f} \boldsymbol{\phi}_s$ for a unique solution $\mathbf{c}_f^* = \mathbf{P}_f \boldsymbol{\phi}_s$, where $\mathbf{P}_f = \tilde{\mathbf{R}}_f^{-1} \tilde{\mathbf{Q}}_f^\mathrm{T} \sqrt{\mathbf{W}_f}$ (see Appendix A.1.2).

Using these expressions, the general discretised relations for the convective and diffusive terms of the integral about face $f$ in Eq. (6) can now be written as

$$ \int_f \mathbf{J}^{\phi,C} \cdot d\mathbf{S} \approx \begin{bmatrix} u(\mathbf{x}_{g1}) \omega_{g1} & u(\mathbf{x}_{g2}) \omega_{g2} & \dots & \mathbf{d}_f^C(\mathbf{x}_g) \\ v(\mathbf{x}_{g1}) \omega_{g1} & v(\mathbf{x}_{g2}) \omega_{g2} & \dots & \mathbf{d}_f^C(\mathbf{x}_g) \end{bmatrix} \begin{bmatrix} \mathbf{P}_f \boldsymbol{\phi}_s \end{bmatrix} \cdot \mathbf{S}_f \quad (16a) \\ = \begin{bmatrix} \mathbf{t}_f^C \boldsymbol{\phi}_s \end{bmatrix} \cdot \mathbf{S}_f $$

$$ \int_f \mathbf{J}^{\phi,D} \cdot d\mathbf{S} \approx \begin{bmatrix} \Gamma_x(\mathbf{x}_{g1}) \omega_{g1} & \Gamma_x(\mathbf{x}_{g2}) \omega_{g2} & \dots & \mathbf{d}_f^{D_x}(\mathbf{x}_g) \\ \Gamma_y(\mathbf{x}_{g1}) \omega_{g1} & \Gamma_y(\mathbf{x}_{g2}) \omega_{g2} & \dots & \mathbf{d}_f^{D_y}(\mathbf{x}_g) \end{bmatrix} \begin{bmatrix} \mathbf{P}_f \boldsymbol{\phi}_s \end{bmatrix} \cdot \mathbf{S}_f \quad (16b) \\ = \begin{bmatrix} \mathbf{t}_f^D \boldsymbol{\phi}_s \end{bmatrix} \cdot \mathbf{S}_f, $$

where

$$ \mathbf{d}_f^C(\mathbf{x}) = \begin{bmatrix} 1 & (x-x_f) & (y-y_f) & (x-x_f)^2 & (x-x_f)(y-y_f) & (y-y_f)^2 & \dots \end{bmatrix}^\mathrm{T} \quad (17a) $$

$$ \mathbf{d}_f^{D_x}(\mathbf{x}) = \begin{bmatrix} 0 & 1 & 0 & 2(x-x_f) & (y-y_f) & 0 & \dots \end{bmatrix}^\mathrm{T} \quad (17b) $$

$$ \mathbf{d}_f^{D_y}(\mathbf{x}) = \begin{bmatrix} 0 & 0 & 1 & 0 & (x-x_f) & 2(y-y_f) & \dots \end{bmatrix}^\mathrm{T} \quad (17c) $$

and $\mathbf{t}_f$ is the (discrete) operator vector.

$$ w_i = \frac{1}{d_i^m}, \quad (15) $$

where $n = 0, 1, 2, 3$, etc., and $d_i$ is the distance between the centroid of face $f$ and the reference location of the $i$th element of the discretisation stencil. Here, we simply set $n$ to 2.

<footnote n="1">¹The following inverse-distance weighting function is widely used in FV methods [15]:</footnote>

<footnote n="2">²Once assembled, the (weighted) design matrix should be normalised using the procedure described in Appendix B.</footnote>
### 4.3.2. Constrained polynomial reconstruction

For the method described here to satisfy the imposed boundary condition exactly, the reconstruction of boundary faces must be constrained. Rather than enforcing the boundary condition all along the edge of the respective boundary CV, we enforce one constraint per integration (or Gauss) point.

The (linearly) constrained least-squares (CLS) problem refers to the problem of finding a least-squares solution that exactly satisfies $q$ additional scalar constraints of the form:

$$ \underbrace{\begin{bmatrix} 1 & (x_g - x_f) & (y_g - y_f) & \dots \\ \vdots & \vdots & \vdots & \vdots \\ 0 & n_{g_x} & n_{g_y} & \dots \\ \vdots & \vdots & \vdots & \vdots \\ n_{g_x} + n_{g_y} & [(x_g - x_f) + d_{g_x}] \cdot n_{g_x} & [(y_g - y_f) + d_{g_y}] \cdot n_{g_y} & \dots \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix}}_{\mathbf{H}_f} \underbrace{\begin{bmatrix} c_1 \\ \vdots \\ c_1 \\ \vdots \\ c_1 \end{bmatrix}}_{\mathbf{c}_f} = \underbrace{\begin{bmatrix} \phi_g \\ \vdots \\ \nabla \phi_g \cdot \mathbf{n}_g \\ \vdots \\ [\phi_g + \mathbf{d}_g \nabla \phi_g] \cdot \mathbf{n}_g \\ \vdots \end{bmatrix}}_{\mathbf{\Phi}_g}, \quad (18) $$

which when using pure matrix notation can be written as

$$ \mathbf{H}_f \mathbf{c}_f = \mathbf{\Phi}_g, \quad (19) $$

where $\mathbf{H}_f \in \mathbb{R}^{q \times n}$ and $\mathbf{\Phi}_g \in \mathbb{R}^q$, and

$$ \mathbf{d}_g = \mathbf{d}(\mathbf{x}_g) = \frac{\Gamma(\mathbf{x}_g)}{\mathbf{u}(\mathbf{x}_g)}. \quad (20) $$

The three blocks of equations in Eq. (18) correspond to boundary condition constraints. While the first block denotes a Dirichlet constraint, the second (below the green line) and third (below the blue line) ones denote Neumann and Robin constraints, respectively.

Following the procedure detailed in Appendix A.2, we can now write the optimality conditions for the (weighted) CLS problem as one set of $n + q$ (linear) equations as

$$ \left[ \frac{2\mathbf{D}_f^\mathrm{T}\mathbf{W}_f\mathbf{D}_f}{\mathbf{H}_f} \frac{\mathbf{H}_f^\mathrm{T}}{\mathbf{0}} \right] \left[ \frac{\mathbf{c}_f}{\mu} \right] = \left[ \frac{2\mathbf{D}_f \sqrt{\mathbf{W}_f} \mathbf{\Phi}_s}{\mathbf{\Phi}_g} \right], \quad (21) $$

where $\mu \in \mathbb{R}^q$ is a set of Lagrange multipliers.

These equations are called the Karush-Kuhn-Tucker (KKT) equations for the CLS problem (KKT are the initials of the last names of William Karush, Harold Kuhn, and Albert Tucker, the three researchers who derived the optimality conditions for a more general form of the constrained minimisation problem).

Assuming that $\mathbf{H}_f$ has linearly independent rows and the stacked matrix $\begin{bmatrix} \sqrt{\mathbf{W}_f} \mathbf{D}_f \\ \mathbf{H}_f \end{bmatrix}$ has linearly independent columns, we can solve this set of equations for a unique solution $\mathbf{c}_f^* = \mathbf{P}_f \mathbf{\Phi}_s + \mathbf{k}_f$, where $\mathbf{P}_f = \tilde{\mathbf{R}}_f^{-1} \{[\tilde{\mathbf{Q}}_{1f}^\mathrm{T} - \tilde{\mathbf{Q}}_{2f}^\mathrm{T}(\tilde{\mathbf{U}}_f^{-1}\tilde{\mathbf{T}}_f^\mathrm{T})\tilde{\mathbf{Q}}_{1f}^\mathrm{T}] \sqrt{\mathbf{W}_f}\} \text{ and } \mathbf{k}_f = \tilde{\mathbf{R}}_f^{-1} \{(\tilde{\mathbf{U}}_f^{-1}\tilde{\mathbf{T}}_f^\mathrm{T})^\mathrm{T} \mathbf{\Phi}_g\} \text{ (see Appendix A.2.1).}$
Similarly to the unconstrained formulation, the general discretised relations for the convec-
tive and diffusive terms of the integral about face f in Eq. (6) can now be written as

$$
\begin{align}
\int_f \mathbf{J}^{\phi,C} \cdot d\mathbf{S} &\approx \left[ \begin{bmatrix} u(\mathbf{x}_{g1}) \omega_{g1} & u(\mathbf{x}_{g2}) \omega_{g2} & \dots \\ v(\mathbf{x}_{g1}) \omega_{g1} & v(\mathbf{x}_{g2}) \omega_{g2} & \dots \end{bmatrix} \mathbf{d}_f^C(\mathbf{x}_g) \right] \left[ \mathbf{P}_f \boldsymbol{\phi}_s + \mathbf{k}_f \right] \cdot \mathbf{S}_f \tag{22a} \\
&= \left[ \mathbf{t}_f^C \boldsymbol{\phi}_s + \mathbf{w}_f^C \right] \cdot \mathbf{S}_f \nonumber
\end{align}
$$

$$
\begin{align}
\int_f \mathbf{J}^{\phi,D} \cdot d\mathbf{S} &\approx \left[ \begin{bmatrix} \Gamma_x(\mathbf{x}_{g1}) \omega_{g1} & \Gamma_x(\mathbf{x}_{g2}) \omega_{g2} & \dots \\ \Gamma_y(\mathbf{x}_{g1}) \omega_{g1} & \Gamma_y(\mathbf{x}_{g2}) \omega_{g2} & \dots \end{bmatrix} \mathbf{d}_f^{D_x}(\mathbf{x}_g) \right] \left[ \mathbf{P}_f \boldsymbol{\phi}_s + \mathbf{k}_f \right] \cdot \mathbf{S}_f \tag{22b} \\
&= \left[ \mathbf{t}_f^D \boldsymbol{\phi}_s + \mathbf{w}_f^D \right] \cdot \mathbf{S}_f, \nonumber
\end{align}
$$

where $\mathbf{t}_f$ and $\mathbf{w}_f$ are the (discrete) operator and constraint vectors, respectively.

4.3.3. Discretisation stencil

Due to the predominantly elliptic nature of the problem at hand, we will employ a central-
biased stencil for faces on the bulk of the domain and a non-central (e.g., upwind/downwind) one
for the ones near/on the boundary of the domain, where we cannot use central discretisation.

The number of elements that make up the discretisation stencil (m) is inherently determined
by the number of monomial terms in the 2D polynomial basis function (n) and by the number
of constraints (q) applied to the respective polynomial fit. The more elements (or equations)
used, the more computationally expensive the scheme becomes. It is, therefore, crucial to reduce
computational cost<fnref n="3" /> without sacrificing the robustness of the scheme.

The most common approaches adopted in [16–23] and references therein are based on the
recursive addition of elements (one by one [18, 19] or layer by layer [21]) until a (minimum)
number of elements determined beforehand is reached or a (maximum) condition number is
attained. However, one must ensure that some stopping criteria are introduced in case the target
metric cannot be reached [23].

Here, we simply add elements (one by one) based on the distance to the centroid of the face about which the polynomial is reconstructed (or fitted) until the problem becomes “sufficiently” well-conditioned (see Fig. 2). Since, in this work, we only deal with regular grids, it suffices to guarantee that the problem is both locally and globally non-singular. If irregular (or very skewed) grids were to be used, then an upper bound should be placed on some $\ell_k$-norm of the condition number of matrices $\mathbf{D}_f$ (ULS formulation) or $\left[ \frac{\mathbf{D}_f}{\mathbf{H}_f} \right]$ (CLS formulation), so as to not compromise the robustness of the scheme.

4.3.4. General approach

The result of the discretisation process is a (linear) system of equations of the form

$$
A\boldsymbol{\phi} = \mathbf{b}, \tag{23}
$$

where $\boldsymbol{\phi}$ is the discrete solution vector. In this system, the coefficients of the unknown variables that constitute matrix $\mathbf{A}$ are a result of the linearisation procedure and the grid geometry, while vector $\mathbf{b}$ contains all sources, constants, and boundary conditions.

<footnote n="3">³Here, we use the density (or the number of non-zero entries) of the global coefficient matrix (A) in Eq. (23) as an alias for the computational cost (or the number of DoFs).</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 2. Central-biased stencil about the face marked in red (elements are numbered in the order that they are added).</figcaption></figure>

Denoting the cell- and face-dependent entries of the operator (and constraint) vectors by $s_c$ and $s_f$, respectively, **A** and **b** can be assembled in the following way:

$$ \mathbf{A}_{ij} = \sum_{f \sim \text{faces}(i)} [\mathbf{t}_f^C(s_c) - \mathbf{t}_f^D(s_c)] \cdot \mathbf{S}_f \quad (24a) $$

$$ \mathbf{b}_i = \iint_i \varphi dS - \sum_{f \sim \text{faces}(i)} \left\{ \left[ \mathbf{t}_f^C(s_f) \mathbf{\phi}_s(s_f) - \mathbf{t}_f^D(s_f) \mathbf{\phi}_s(s_f) \right] + \underbrace{(\mathbf{w}_f^C - \mathbf{w}_f^D)}_{\text{CLS formulation (only)}} \right\} \cdot \mathbf{S}_f. \quad (24b) $$

## 5. Numerical verification

The validation of CFD codes for relevant applications is an extremely complex process, as it is hardly possible to generate experimental data at the same level of detail as that provided by the simulation. In addition, since both experimental and numerical data are constantly subjected to many error and uncertainty sources, an absolute validation of those results remains an unreachable objective. These uncertainties make the verification phase even more essential to ensure that, at least, the numerical discretisation is foolproof. Consequently, we will restrict our analysis to this last phase.

### 5.1. Exact solutions

For rigorous code verification, we must have access to the exact solution of the underlying governing equations. However, whenever exact solutions are found for complex PDEs, they often rely on significant simplifications in dimensionality, geometry, physics, etc. (e.g., Couette flow). For this reason, we will employ the method of manufactured solutions (MMS).

### 5.2. Error sources

Suppose that the (possibly non-linear) differential equation is given by

$$ \mathcal{L}(\Phi) = 0 \quad (25a) $$
and the difference approximation by

$$ \mathcal{L}_{(p)}(\phi) = 0, \qquad (25b) $$

where $\mathcal{L}_{(p)}(\cdot)$ is the pth-order accurate discrete form of the symbolic differential operator $\mathcal{L}(\cdot)$.

In general, the (sufficiently smooth) exact solution to the governing equations will not satisfy the difference approximation in Eq. (25b). The difference (or remainder),

$$ \bar{\tau}^{\Phi} \equiv \bar{\tau}(\Phi) = \mathcal{L}_{(p)}(\Phi) - \mathcal{L}_{(p)}(\phi) \qquad (26) $$

is called truncation error.

The solution to the discretised equations differs from the exact solution to the governing equations by the discretisation error, i.e.,

$$ e^{\Phi} \equiv e(\Phi) = \phi - \Phi. \qquad (27) $$

If the symbolic operator $\mathcal{L}_{(p)}(\cdot)$ is linear, then $\mathcal{L}_{(p)}(\Phi) - \mathcal{L}_{(p)}(\phi) = \mathcal{L}_{(p)}(\Phi - \phi)$, and thus Eq. (26) simplifies to

$$ \bar{\tau}^{\Phi} = -\mathcal{L}_{(p)}(e^{\Phi}). \qquad (28a) $$

Otherwise, for small enough $e^{\Phi}$, it becomes

$$ \bar{\tau}^{\Phi} \approx -\mathcal{L}'_{(p)}(e^{\Phi}), \qquad (28b) $$

where $\mathcal{L}'_{(p)}$ is the Jacobian matrix of $\mathcal{L}_{(p)}$. This equation governs the transport of the discretisation error and is called the discretisation error transport equation since it employs the discrete symbolic differential operator [24].

### 5.2.1. General approach

To derive the analytical expression for the local truncation error about a generic face $e$ in the bulk of the domain, consider the uniformly spaced 1D grid in Fig. 3.

<figure><img src="image_5.png" /><figcaption>Figure 3. Uniformly spaced 1D grid (spacing $h$).</figcaption></figure>

Firstly, we start by writing the facial values of $\phi(x)$ and $\phi(x)'$ about $e$ at $x = x_e$ in terms of the parameter vector of the polynomial regression model as

$$ \phi_e(x_e) = \frac{1}{2}\phi_C + \frac{1}{2}\phi_E \qquad (29a) $$

$$ \phi'_e(x_e) = -\frac{1}{h}\phi_C + \frac{1}{h}\phi_E. \qquad (29b) $$
Then, using the general expression in Eq. (26), we may write the local truncation error about face e at $x = x_e$ as

$$ \begin{aligned} \bar{\tau}_e^\Phi (x_e) &= \bar{\tau}_e^{\Phi,C} (x_e) - \bar{\tau}_e^{\Phi,D} (x_e) \\ &= u_e \underbrace{\left( \Phi_e + \frac{h^2}{8} \Phi_e'' + \dots \right)}_{\Phi_e^{(p)}(x_e)} - \Gamma_e \underbrace{\left( \Phi_e' + \frac{h^2}{24} \Phi_e''' + \dots \right)}_{\Phi_e'^{(p)}(x_e)} - \underbrace{\left[ u_e \underbrace{\left( \frac{\phi_C + \phi_E}{2} \right)}_{\phi_e^{(p)}(x_e)} - \Gamma_e \underbrace{\left( \frac{\phi_E - \phi_C}{h} \right)}_{\phi_e'^{(p)}(x_e)} \right]}_{\mathbf{t}_e \Phi_s} \quad (30) \\ &= u_e \underbrace{\left[ \Phi_e - \left( \frac{\phi_C + \phi_E}{2} \right) \right]}_{\bar{\tau}_0^\Phi(x_e)} - \Gamma_e \underbrace{\left[ \Phi_e' - \left( \frac{\phi_E - \phi_C}{h} \right) \right]}_{a_0^\Phi(x_e)} + \underbrace{u_e \left( \frac{h^2}{8} \Phi_e'' + \dots \right)}_{\bar{\tau}_1^\Phi(x_e) + \dots} - \Gamma_e \underbrace{\left( \frac{h^2}{24} \Phi_e''' + \dots \right)}_{a_1^\Phi(x_e) h^2 + \dots} \\ &= \underbrace{a_0^\Phi(x_e)}_{\bar{\tau}_0^\Phi(x_e)} \end{aligned} $$

where

$$ \begin{aligned} \bar{\tau}_e^{\Phi,C} (x_e) &= \mathbf{t}_e^C \Phi_s - \mathbf{t}_e^C \phi_s \\ &= u_e \left[ \Phi_e^{(p)}(x_e) - \phi_e(x_e) \right] \\ &= u_e \left[ \left( \frac{\Phi_C + \Phi_E}{2} \right) - \left( \frac{\phi_C + \phi_E}{2} \right) \right] \quad (31a) \\ &= u_e \left[ \left( \Phi_e + \frac{h^2}{8} \Phi_e'' + \dots \right) - \left( \frac{\phi_C + \phi_E}{2} \right) \right] \end{aligned} $$

and

$$ \begin{aligned} \bar{\tau}_e^{\Phi,D} (x_e) &= \mathbf{t}_e^D \Phi_s - \mathbf{t}_e^D \phi_s \\ &= \Gamma_e \left[ \Phi_e'^{(p)}(x_e) - \phi_e'(x_e) \right] \\ &= \Gamma_e \left[ \left( \frac{\Phi_E - \Phi_C}{h} \right) - \left( \frac{\phi_E - \phi_C}{h} \right) \right] \quad (31b) \\ &\stackrel{4}{=} \Gamma_e \left[ \left( \Phi_e' + \frac{h^2}{24} \Phi_e'' + \dots \right) - \left( \frac{\phi_E - \phi_C}{h} \right) \right] \end{aligned} $$

236 are its (weighted) convective and diffusive components at $x = x_e$, respectively.

237 Rewriting the face truncation (or tau) error in Eq. (30) in a more compact form yields

$$ \begin{aligned} \underbrace{\bar{\tau}_e^\Phi (x_e)}_{\text{Tau}} &= \underbrace{\bar{\tau}_0^\Phi |_{e} (x_e)}_{\text{Tau-0}} + \underbrace{\bar{\tau}_{\Sigma k}^\Phi |_{e} (x_e)}_{\text{Tau-k}} \\ &\equiv \bar{\tau}_0^\Phi |_{e} (x_e) + \sum_{k=m}^{\infty} \bar{\tau}_k^\Phi |_{e} (x_e), \end{aligned} \quad (33) $$

238 where $m$ refers to the leading term of its tau-k component.

$$ \Phi_C = \Phi_e - \frac{h}{2}\Phi_e' + \frac{h^2}{8}\Phi_e'' - \frac{h^3}{48}\Phi_e''' + \frac{h^4}{384}\Phi_e^{iv} - \frac{h^5}{3840}\Phi_e^v + \dots \quad (32a) $$

$$ \Phi_E = \Phi_e + \frac{h}{2}\Phi_e' + \frac{h^2}{8}\Phi_e'' + \frac{h^3}{48}\Phi_e''' + \frac{h^4}{384}\Phi_e^{iv} + \frac{h^5}{3840}\Phi_e^v + \dots, \quad (32b) $$

where $\Phi_e'$, $\Phi_e''$, $\Phi_e'''$, $\Phi_e^{iv}$ and $\Phi_e^v$ denote the first, second, third, fourth, and fifth derivatives of $\Phi(x)$ at $x = x_e$.

<footnote n="4">4 The 1D Taylor series expansions of $\Phi_C$ and $\Phi_E$ about face e are given by</footnote>
### 5.2.2. Closure

Using the 1D expressions derived in Section 5.2.1, we can now write the 2D form of the two components of the truncation error of a generic cell C in a more general form as

$$ 
\begin{aligned}
\bar{\tau}_{0}^{\Phi}|_{C} &= \iint_{C} \varphi \, dS - \sum_{f \sim \text{faces}(C)} \left( \int_{f} \mathbf{J}^{\Phi} \cdot d\mathbf{S} \right) \\
&\approx \sum_{f \sim \text{faces}(C)} \left( \int_{f} \varphi \cdot d\mathbf{s} - \sum_{g \sim \mathcal{G}(f)} [\omega_{g} \mathbf{J}^{\Phi}(\mathbf{x}_{g})] \cdot \mathbf{S}_{f} \right) \\
&= \sum_{f \sim \text{faces}(C)} \left( \int_{f} \varphi \cdot d\mathbf{s} - [\mathbf{t}_{f} \Phi_{s} + \mathbf{w}_{f}] \cdot \mathbf{S}_{f} \right) \\
&= 0
\end{aligned}
\quad (34a) $$

$$ 
\begin{aligned}
\sum_{k=m}^{\infty} \bar{\tau}_{k}^{\Phi}|_{C} &= \sum_{f \sim \text{faces}(C)} \left( \int_{f} \mathbf{J}^{\Phi} \cdot d\mathbf{S} \right) - \iint_{C} \varphi \, dS \\
&= \sum_{f \sim \text{faces}(C)} \left( [\mathbf{t}_{f} \Phi_{s} + \mathbf{w}_{f}] \cdot \mathbf{S}_{f} - \int_{f} \varphi \cdot d\mathbf{s} \right) \\
&= \sum_{f \sim \text{faces}(C)} \left( \sum_{k=m}^{\infty} \bar{\tau}_{k}^{\Phi}|_{f}(x_{f}) \cdot \mathbf{S}_{f} \right),
\end{aligned}
\quad (34b) $$

where m refers to the leading term of its tau-k component.

### 5.3. Error estimation

In verification, the exact solution of the mathematical model is generally not known, and thus errors must be estimated rather than simply evaluated. Current efforts in estimating the discretisation error are focused on a posteriori methods. These methods use the solution of the discretised equations to estimate the discretisation error relative to the exact solution of the mathematical model.

In general, the level of reliability of a posteriori methods is strongly problem-dependent. As detailed in Section 4.3.3, all the methods developed here are predominantly suited for elliptic problems, and as a result, they are not as well-developed for mathematical models that are parabolic or hyperbolic in nature. The level of complexity of the problem is also an important issue. These error estimators work well for smooth, linear problems with simple physics and geometries; however, strong non-linearities, discontinuities, singularities, and physical complexity can significantly reduce their applicability.

**Example 1.** Consider the 1D steady-state form of Eq. (2) discretised with a second-order accurate scheme over the domain $[-1, 1]$, where $u = 1$ [$m/s$] and $\Gamma = 1$ [$m^2/s$]. Dirichlet conditions were prescribed for both boundaries, and the (manufactured) solution to this equation was assumed to be:

$$ \Phi(x) = x^5. \quad (35) $$

<footnote>$^5\bar{\tau}_{0}^{\Phi}|_{C}$ corresponds to the difference between the integral of the source term and the (discrete) surface flux integral. Since the source term is computed exactly, then $\bar{\tau}_{0}^{\Phi}|_{C} = 0$. If a Gaussian quadrature integration was instead used to discretise the RHS of Eq. (4), this would not be the case.</footnote>
The 1D distribution of the analytical scalar property and respective volumetric source (or sink) term is shown in Figs. 4a and 4b, respectively.

<figure><img src="image_6.png" /><figcaption>Figure 4. 1D distribution of the analytical scalar property and respective volumetric source (or sink) term.</figcaption></figure>

In the next two Sections (5.3.1 and 5.3.2), we present two distinct forms of estimating the numerical error, both of which explicitly require the formulation of a higher-order problem in the original grid.

### 5.3.1. Higher-order polynomial regression

The estimate presented here seeks to find an approximation of the exact solution that is of higher order of accuracy than the underlying discrete one. Examples of this approach can be found in the works of Hay and Visonneau [25] and Hay et al. [26].

Replacing the analytical response vector with a discrete one of higher order (e.g., fourth-order accuracy), Eq. (30) becomes

$$ \begin{aligned} \bar{\tau}_e^{\phi(h)}(x_e) &= \underbrace{u_e\left(\phi_e + \frac{h^2}{8}\phi_e''\right) - \Gamma_e\left(\phi_e' + \frac{h^2}{24}\phi_e'''\right)}_{\mathbf{t}_{e(l)}\mathbf{\Phi}_{s(l)}} - \underbrace{\left[u_e\left(\frac{\phi_C + \phi_E}{2}\right) - \Gamma_e\left(\frac{\phi_E - \phi_C}{h}\right)\right]}_{\mathbf{t}_{e(l)}\mathbf{\Phi}_{s(l)}} \\ &= \underbrace{u_e\left[\phi_e - \left(\frac{\phi_C + \phi_E}{2}\right)\right] - \Gamma_e\left[\phi_e' - \left(\frac{\phi_E - \phi_C}{h}\right)\right]}_{\bar{\tau}_0^{\phi(h)}|_e(x_e)} + \underbrace{u_e\left(\frac{h^2}{8}\phi_e''\right) - \Gamma_e\left(\frac{h^2}{24}\phi_e'''\right)}_{\bar{\tau}_1^{\phi(h)}|_e(x_e)} = a_0^{\phi(h)}|_e(x_e)h^2 \end{aligned} \quad (36) $$

where the subscripts (l) and (h) refer to the lower- and higher-order accurate discretisation schemes, respectively (note that the facial values $\phi_e, \phi_e', \phi_e''$, and $\phi_e'''$ are approximated by a cubic variation of $\phi(x)$ about the face).

If the difference between the analytical expression in Eq. (30) and the finite estimate in Eq. (36) is “small” in comparison to the latter, then this estimate is a “good” approximation of the local face truncation error (see Fig. 5).

As shown in Fig. 6, we can see that (in the asymptotic range) this difference is indeed negligible, i.e.,

$$ \mathcal{L}_{(2)}(\Phi - \Phi_{(4)}) \approx O(h^4). \quad (37) $$
<figure><img src="image_7.png" /><figcaption>**Figure 5.** 1D distribution of the magnitude of the error for a uniformly spaced grid with 100 CVs (2nd-order accuracy).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Error norm decay for increasingly refined uniformly spaced grids, with $h_0 = 1e^{-02}$ (2nd-order accuracy).</figcaption></figure>
### 5.3.2. Tau-extrapolation

The estimate presented here uses an extrapolation of the lower-order underlying discrete solution to obtain a higher-order variation of $\phi(x)$ about the face. In other words, this variation simply corresponds to substituting the qth-order accurate discrete solution in a pth-order accurate discrete operator $\mathcal{L}_{(p)}(\cdot)$, with $p > q$. This approach is referred to as the defect correction method in the works of Ervin and Layton [27], and Pierce and Giles [28].

Replacing the lower-order discrete operator vector with one of higher order (e.g., fourth-order accuracy), Eq. (30) becomes

$$ 
\begin{aligned}
\bar{\tau}_e^{\phi_{(l)}}(x_e) &= (\mathbf{t}_{e(l)} - \mathbf{t}_{e(h)}) \phi_{s(l)} \\
&= u_e \left(\frac{h^2}{8} \phi_e''\right) - \Gamma_e \left(\frac{h^2}{24} \phi_e'''\right) \\
&= a_1^{\phi_{(l)}}|_e (x_e) h^2,
\end{aligned}
\quad (38) $$

where

$$ 
\begin{aligned}
(\mathbf{t}_{e(l)}^C - \mathbf{t}_{e(h)}^C) \phi_{s(l)} &= u_e \left[ \phi_e^{(p)}(x_e) - \phi_e(x_e) \right] \\
&= u_e \left[ \left( \frac{\phi_C + \phi_E}{2} \right) - \phi_e(x_e) \right] \\
&= u_e \left[ \left( \frac{\phi_C + \phi_E}{2} \right) - \left( -\frac{1}{16}\phi_W + \frac{9}{16}\phi_C + \frac{9}{16}\phi_E - \frac{1}{16}\phi_{EE} \right) \right] \\
&= u_e \left[ \frac{h^2}{8} \phi_e'' \right]
\end{aligned}
\quad (39a) $$

and

$$ 
\begin{aligned}
(\mathbf{t}_{e(l)}^D - \mathbf{t}_{e(h)}^D) \phi_{s(l)} &= \Gamma_e \left[ \phi_e^{(p)}(x_e) - \phi_e'(x_e) \right] \\
&= \Gamma_e \left[ \left( \frac{\phi_E - \phi_C}{h} \right) - \phi_e'(x_e) \right] \\
&= \Gamma_e \left[ \left( \frac{\phi_E - \phi_C}{h} \right) - \left( \frac{1}{24h}\phi_W - \frac{9}{8h}\phi_C + \frac{9}{8h}\phi_E - \frac{1}{24h}\phi_{EE} \right) \right] \\
&= \Gamma_e \left[ \frac{h^2}{24} \phi_e''' \right]
\end{aligned}
\quad (39b) $$

are its (weighted) convective and diffusive components, respectively.

Looking back on the estimate derived in Eq. (38), we can see that the expression for $\bar{\tau}_e^{\phi_{(l)}}(x_e)$ is formally equivalent to that of $\bar{\tau}_e^{\phi_{(h)}}|_e(x_e)$ in Eq. (36), with the exception that the facial values $\phi_e''$ and $\phi_e'''$ are extrapolated from the underlying discrete solution, rather than evaluated using a higher-order one (see Figs. 7 and 8). This estimate (for the tau-k component of the face truncation error) is thus less accurate than the one derived in Section 5.3.1, since the higher-order problem does not need to be solved, only formulated.

In addition, even though this error estimator is not able to estimate the tau-0 component of the face truncation error, we can always confidently approximate the cell discretisation and truncation errors (note that this component is nil about the cell, see Eq. (34a)). Nevertheless, as we will see in Section 6.1.1, the tau-k component of this error is, in fact, the source of the (cell) discretisation error, not the sum of both of the aforementioned components.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** 1D distribution of the magnitude of the error for a uniformly spaced grid with 100 CVs (2nd-order accuracy).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Error norm decay for increasingly refined uniformly spaced grids, with $h_0 = 1e^{-02}$ (2nd-order accuracy).</figcaption></figure>
297
6. Adaptive mesh refinement

298
299
300
In addition to estimating the discretisation error, we also desire methods for reducing it when
it is found to be too large. Adaptive mesh refinement strategies are generally classified as one of
the following: *global*, *semi-global*, or *local* refinement.

301
302
303
304
305
306
307
308
In global *p*-refinement (also called uniform *p*-refinement), the degree of every element in the
mesh is increased by the same quantity. Even though in most cases this is the simplest strategy
to implement, it can be very wasteful in the sense that many elements away from the area(s) of
interest are refined. A modification of the global strategy is semi-global *p*-refinement, in which
the degree of every element in one or more selected cross-sections is refined. In certain cases,
this strategy may be implemented as easily as global *p*-refinement and may be less wasteful.
However, complications occur when the desired refinement does not correspond to a natural
cross-section.

309
310
The last refinement strategy is (adaptive) local *p*-refinement, in which only selected elements
in restricted portions of the domain are refined.

311
For (adaptive) local *p*-refinement to be efficient, it is necessary that:

312
• Computation is not expensive<sup>6</sup>.

313
• Neighbouring elements are *not* too different in degree.

314
• The algorithm can be implemented cheaply<sup>7</sup>.

315
316
317
While the first and third requirements are natural for any practical algorithm, the main purpose of the second one is to smooth the changes of the approximate solution from element to element.

318
6.1. *p*-adaptive meshing

319
320
321
322
323
324
An automatic mesh adaptation algorithm should have the capability of detecting important
flow features in the same way that the user does and adapt the mesh accordingly. In *p*-adaptive
meshing, we fix the geometry of the grid and achieve a locally better approximate solution by
modifying (usually increasing) the local formal order of accuracy of the discretisation scheme
(note that since the geometry of the grid remains unchanged upon adaptation, we assume that all
the grids we are working with are geometrically admissible).

325
326
327
Even though the method proposed here is valid for any unstructured grid, extra care must be
taken when dealing with (non-isotropic) irregular grids (see Section 4.3.3). For this reason, in
the numerical examples presented in Section 6.2 we simply use a Cartesian grid<fnref n="8" />.

<footnote n="6"><sup>6</sup>If, upon refinement, the system of equations becomes increasingly dense, then the advantage gained by refining only a small number of elements is partially lost.</footnote>

<footnote n="7">⁷If the (adaptive) process cannot be implemented cheaply, then the actual cost of computing an approximate solution to the (partial) differential equation(s) may become higher with (adaptive) local refinement than with global or semi-global refinement, even in cases where they are extremely wasteful.</footnote>

<footnote n="8">⁸The ideal mesh is a Cartesian distribution, where all the points are equidistant and where all cells are perfect cubes, with $\Delta x = \Delta y = \Delta z$. This grid is associated with the highest possible accuracy of the discretised formulas, where the FVM leads to the same formulas as FDs.</footnote>
### 6.1.1. Adaptation driver

Recall that in Section 5.2 we derived the linear discrete form of the (discretisation) error transport equation (28a). This equation is repeated here for convenience

$$ \bar{\tau}^{\Phi} = -\mathcal{L}_{(p)}(e^{\Phi}), \qquad (40) $$

where $\mathcal{L}_{(p)}(\cdot)$ is the pth-order accurate discrete form of the symbolic differential operator $\mathcal{L}(\cdot)$. According to this equation, the truncation error serves as the local source of the discretisation error [29], and thus reducing the latter should result in a commensurate reduction of the latest.

However, as we have seen in section Section 5.2.1, the tau-0 component of the face truncation error does not directly contribute to the cell discretisation and truncation errors (see Eq. (34a)), and thus we should investigate whether the face truncation (or tau) error in Eq. (33) is the rightful criterion to drive the (adaptive) process, or if it is perhaps its tau-k component alone.

**Example 2.** Consider the 1D steady-state form of Eq. (2) initially discretised with a second-order accurate scheme over the domain $[-1, 1]$, where $u = 1$ [$m/s$] and $\Gamma = 1$ [$m^2/s$]. Dirichlet conditions were prescribed for both boundaries, and the (manufactured) solution to this equation was assumed to be the one used in the previous example (see Eq. (35)).

The 1D distribution of the magnitude of the error before refinement for the uniformly spaced grid with 100 CVs used in Example 1 is shown in Fig. 9.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** 1D distribution of the magnitude of the error before refinement (2nd-order accuracy).</figcaption></figure>

By sorting the magnitude of both criteria (tau and tau-k) in descending order and flagging for refinement the first $\lambda \sim 50\%$ of elements (or faces) from this list, we compare the 1D distribution of the magnitude of the error after refinement with the one prior to it (see error thresholds in Fig. 9a and the respective refinement zones marked in red and yellow, respectively).

The procedure detailed here is different from the one proposed by Syrakos et al. [30] in the sense that we gauge how both criteria fare after the first adaptation cycle and not at the end of the (adaptive) process. As such, flagged elements (or faces) are refined owing to large error estimates and not because they were in violation of the refinement rules proposed in Section 6.1.3 (recall that the sole purpose of the adaptation criterion is to pinpoint regions of the domain deemed inaccurate).
<figure><img src="image_12.png" /><figcaption>**Figure 10.** 1D distribution of the magnitude of the (face) truncation error after refinement.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** 1D distribution of the magnitude of the (cell) error after refinement.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** 1D distribution of the refinement level after refinement.</figcaption></figure>

To decide which of these is the most appropriate criterion to drive the (adaptive) process, we measure the (analytical) error reduction ratio between the error after (e) and before (e₀) the refinement through a global norm over the entire domain Ω.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** $l_1$- and $l_\infty$-error norm reduction ratios.</figcaption></figure>

As shown in Fig. 13, since the sub-problem on the right (10b, 11b and 12b) produces a significantly greater reduction of the (cell) discretisation error, we will use the tau-k component of the face truncation error alone (or perhaps its estimate) to drive the (adaptive) process.

### 6.1.2. Refinement mechanics

Almost all adaptation methods proceed from the fundamental premise of error equidistribution [31]. Equidistribution aims at producing a mesh that contains the same value of (adaptive) refinement parameter (or criterion) over the entire domain, thereby improving the accuracy of the solution. Since this parameter (or criterion) is a stand-in for the face truncation error, this goal spreads the (cell) discretisation error out evenly over the entire domain<fnref n="9" />. In practice, because this principle assumes that all errors are equally important, it is typically somewhat over-conservative, and more often than not, this is certainly not the case. However, without any additional guidance about what is important, equidistribution simply ensures that everything is equally correct.

**Decision criteria.** Bangerth and Rannacher [32] described non-trivial strategies on how to decide the subset of elements (or faces) to be flagged for refinement, based on the adaptation criterion selected in Section 6.1.1.

On the one hand, *fixed-error-reduction* or *fixed-fraction* strategies select subsets of elements whose criteria accumulate to a predefined fraction of their global sum. This strategy is only applicable when the sum of all criteria actually has meaning (e.g., local errors that add up to the global one). Even though this strategy is more sophisticated than the one introduced below, it tends to only refine (or adapt) very few elements whenever smooth rapid solution changes or singularities are encountered.

On the other hand, *fixed-rate* or *fixed-number* strategies select predefined fractions of elements with the highest (or lowest) criterion for adaptation, which not only allows us to control how aggressively each subsequent mesh is refined but also to predict the approximate increase in the computational cost of the (adaptive) process. For these reasons, we will use this strategy in our investigation.

<footnote n="9">⁹The idealised histogram of the refinement parameter (or criterion) looks like a delta function with height equal to the (total) number of faces.</footnote>
**Stopping criteria.** One or more of the following parameters will ultimately define the termination criteria of the (adaptive) process:

* A (remaining) desired error level for the functional of interest.

* An upper limit on the (total) number of DoFs<fnref n="10" />.

* An upper limit on the (total) number of adaptation cycles.

One additional parameter is the maximum level of refinement (note that the conditioning of Vandermonde matrices becomes worse as the degree of the polynomial approximation increases). Hence, if, upon refinement, the local coefficient matrix becomes singular to working precision, then that element (or face) can no longer be refined.

### 6.1.3. Refinement rules

To achieve the second requirement for efficient (adaptive) local p-refinement introduced at the beginning of Section 6, we propose the following rules inspired by those of the same names in Red-Green (h-) refinement [14]. Much like any other adaptation method, the set of rules proposed here is one that we can only hope will drive the (adaptive) process to approximate the optimal mesh, while not compromising error estimation (note that the lack of smoothing between lower- and higher-order regions may lead the algorithm to refine elements (or faces) that would not otherwise be refined had reliable estimates been obtained).

**Neighbour rule (p-version).** For any element (or face) of degree p to be considered eligible for refinement, all the elements (or faces) that belong to the cell(s) that contain(s) the element (or face) flagged in Fig. 14a should be, at least, of degree p (see Fig. 14b).

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Neighbour rule (refinement level: ●-p₁, ●-p₂, with p₂ ≡ p = p₁ + 1).</figcaption></figure>

$$h \sim \text{DoFs}^{-1/d}, \qquad (41)$$

where d is the dimension of the problem.

<footnote n="10">¹⁰The (total) number of DoFs is often used to define a length scale, especially when mesh adaptation is performed by modifying the local formal order of accuracy of the discretisation scheme without changing the geometry of the grid. The grid size is related to the number of DoFs through the following equation [33]:</footnote>
**1-Irregular rule (p-version).** The difference in degree of elements (or faces) in neighbouring cells can be, at most, 1; that is, the degree of all the elements (or faces) that belong to the cell(s) that contain(s) an element (or face) that shares a cell with the flagged element (or face) should be, at least, of degree *p* (see Fig. 15b). Note that this rule is complementary to the previous one, i.e., it provides an extra layer of smoothing.

<figure><img src="image_17.png" /><figcaption>**Figure 15.** 1-Irregular rule (refinement level: ●-p₁, ●-p₂, with p₂ ≡ p = p₁ + 1).</figcaption></figure>

## 6.2. Numerical examples

In this section, we ultimately seek to answer the following set of questions:

* Are the proposed error estimators able to accurately drive the (adaptive) process?

- If so, how does the distribution of the refinement level using the analytical solution compare to that predicted by the aforementioned error estimators?

- If so, how do the magnitude of the analytical discretisation and truncation errors compare to that predicted by these estimators?

* Does the (adaptive) process effectively reduce the spatial discretisation error?

- If so, how efficiently (i.e., how does it compare to global (p-) refinement)?

**Example 3.** Consider the 2D steady-state form of Eq. (2) initially discretised with a second-order accurate scheme over the domain $[-1, 1] \times [-1, 1]$, where $(u, v) = (1, 1)$ [$m/s$] and $(\Gamma_x, \Gamma_y) = (1, 1)$ [$m^2/s$]. Dirichlet conditions were prescribed for all the boundaries, and the (manufactured) solutions to this equation were assumed to be:

$$ \Phi_a(x, y) = \exp(-100x^2y^2) \qquad (42a) $$

$$ \Phi_b(x, y) = \exp(-100x^2y^2) + \sqrt{(x - \pi/8)^2 + y^2} + \sqrt{x^2 + (y - \pi/8)^2}. \qquad (42b) $$

In addition, the (adaptive) process was run under the following conditions:

* Fixed-rate (or fixed-number) parameter: $\lambda = 10\%$.

* Stopping criterion: computational cost ratio increase¹¹ of (approximately) 3.
The 2D distribution of the magnitude of the analytical scalar property and respective volumetric source (or sink) term is shown in Figs. 16 and 17, respectively (note that the white contours lines correspond to regions where the respective functions have a constant value equal to the ticks of the colour bar).

While the smooth rapid change in the middle of the domain makes the solution on the left (16a and 17a) tailor-made for *p*-adaptation, the addition of singularities at $(\pi/8, 0)$ and $(0, \pi/8)$ (16b and 17b) has the opposite effect (recall that the partial derivatives of $\Phi_b(x, y)$ with respect to $x$ and $y$ are respectively unbounded at these points).

<figure><img src="image_18.png" /><figcaption>**Figure 16.** 2D distribution of the magnitude of the analytical scalar property.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 17.** 2D distribution of the magnitude of volumetric source (or sink) term.</figcaption></figure>

The 2D distribution of the magnitude of the analytical (cell) discretisation and truncation errors for these solutions in a Cartesian grid with 10k CVs ($h = 1e-02$) is shown in Figs. 18 and 19, respectively.

As expected, the regions that contain elements (or faces) near the origin of the exponential and two singularities are the ones where the greater error magnitude occurs (see Figs. 18a and 18b), and so it is expected that the algorithm will prioritise their refinement.

<footnote n="11">¹¹The computational cost increase (or “gain”) for the nth adaptation cycle refers to the ratio between the number of DoFs of cycle n and the number of DoFs of the uniform 2nd-order accurate case.</footnote>
<figure><img src="image_20.png" /><figcaption>Figure 18. 2D distribution of the magnitude of the analytical (cell) truncation error (2nd-order accuracy).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 19. 2D distribution of the magnitude of the analytical (cell) discretisation error (2nd-order accuracy).</figcaption></figure>

The 2D distribution of the refinement level (in the last adaptation cycle) for the cases where the value of the adaptation criterion is evaluated using the analytical solution or estimated using the two error estimators derived in Sections 5.3.1 and 5.3.2 is shown in Figs. 20 to 22, respectively.

From a distribution point of view, it is clear that both estimators were able to reliably identify where further resolution was needed, even for the case on the right ($\Phi_b(x, y)$), where one may have initially thought that the two singularities would have prevented the algorithm from doing so.

While for the case on the left (20a, 21a and 22a) only the region near the origin of the exponential is targeted for refinement, the one on the right (20b, 21b and 22b) targets both this region and the one near the two singularities. However, if we were to instead employ a fixed-error-reduction (or fixed-fraction strategy), the algorithm would eventually stop flagging elements (or faces) near the origin of the exponential and strictly flag the ones near the two singularities. Nevertheless, since the proposed method is not able to reduce the magnitude of the error near discontinuities and/or singularities, the only meaningful difference between these strategies is that the latest would become less wasteful.
<figure><img src="image_22.png" /><figcaption>Figure 20. 2D distribution of the refinement level using the analytical solution.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 21. 2D distribution of the refinement level using the estimate provided by the error estimator in Section 5.3.1.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 22. 2D distribution of the refinement level using the estimate provided by the error estimator in Section 5.3.2.</figcaption></figure>

Finally, we present the cycle-by-cycle evolution of both the adaptation criteria and the discretisation error (see Figs. 23 to 26). In these figures, we denote by $\Phi$, $\phi_{(h)}$ and $\phi_{(l)}$ the three cases in Figs. 20 to 22, respectively, and by $\hat{p}_k$ the analytical uniform $k$-degree accurate response.
As for the first case (23 and 24), we can see that the (adaptive) process is able to reduce the spatial discretisation error to similar levels to that of a uniform 14th-order accurate scheme ($\hat{\mathbf{p}}_{13}$) for a computational cost inferior to its 6th-order accurate counterpart ($\hat{\mathbf{p}}_{5}$), i.e., at roughly ~15% of the cost (note that it only makes sense to use methods of this high of an order if the asymptotic range has been reached and provided that the round-off error is small enough to where it does not pollute the solution).

<figure><img src="image_25.png" /><figcaption>Figure 23. $\ell_1$- and $\ell_\infty$-error norm decay (tau-k) norm decay for $\Phi_a(x, y)$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 24. $\ell_1$- and $\ell_\infty$-error norm decay for $\Phi_a(x, y)$.</figcaption></figure>

As for the second case (25 and 26), unlike the previous one, neither estimator was able to provide reliable estimates for the magnitude of the error. In fact, this was to be expected, as the discretisation schemes were not formulated to handle discontinuities and/or singularities. Once the algorithm starts targeting the region of the domain that contains the two singularities and the magnitude of the error in this region surpasses the one that contains the origin of the exponential (i.e., after the first adaptation cycle), both estimators appear to grossly underestimate the magnitude of the *true* error, thereby misleading the user into thinking that the (adaptive) process is accomplishing something that it was not meant to: reduce the magnitude of the error.
<figure><img src="image_27.png" /><figcaption>Figure 25. $l_1$- and $l_\infty$-error (face) truncation error (tau-k) norm decay for $\Phi_b(x, y)$.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 26. $l_1$- and $l_\infty$-discretisation error norm decay for $\Phi_b(x, y)$.</figcaption></figure>

467

# 7. Conclusions

468

The main conclusions drawn from the work developed in this paper include:

469

470

471

472

473

474

475

476

477

* In Section 4, we used a least-squares approach to approximate the local underlying solution about the centroid of the face. Using least-squares methods to compute the coefficients of the polynomial regression model [34] offers more flexibility with regard to the order of accuracy achieved [35] and the discretisation stencil used [36]. Nevertheless, this flexibility comes at a cost, as proper weighting is needed (note that the inverse-distance weighting function in Eq. (15) is simply an “educated guess”; the optimal choice of weights should be made on the basis of the element residuals. However, since we can only compute them once the solution to the (linear) system of equations in Eq. (23) is obtained, this would make for an iterative process).
478

479

480

* In Section 5.2, we used the discrete form of the (discretisation) error transport equation to derive two distinct *a posteriori* error estimators, which were then used to drive the (adaptive) process in Section 6.

481

482

483

- Even though both estimators were able to reliably pinpoint regions of the domain deemed inaccurate, they grossly underestimated the magnitude of both the (face) truncation and (cell) discretisation errors in problems with singularities.

484

485

486

487

- Seeing that both estimators yielded similar results for the numerical examples presented in Section 6.2, we can conclude that the tau-extrapolation estimate derived in Section 5.3.2 is computationally more attractive since the higher-order problem does not need to be solved, only formulated.

488

489

490

491

* As far as the effectiveness of the (adaptive) process goes, we were able to corroborate our initial assumption that, despite having proven accurate for smooth problems, order-based adaptation (or *p*-adaptation) is typically not as successful when applied in the presence of discontinuities and/or singularities [4].

492

493

494

495

- The work of Nguyen [37] in a FE framework suggests that *hp*-adaptive methods provide the best approach to solving this class of problems, where one can use smaller elements of lower degree near the discontinuities and/or singularities and bigger elements of higher degree elsewhere.

496

497

* The present study is only concerned with the steady-state case. Transient problems would come with the following additional complexities:

498

499

500

501

- The truncation error would contain terms that are a function of both the grid spacing (*h*) and the time step ($\Delta t$), and thus the general expressions of the error estimators derived in Sections 5.3.1 and 5.3.2 should be reformulated to account for this additional component (for further reference, see [38]).

502

503

504

- The (adaptive) process should not be re-run for each time step. Instead, it should be run once and coarsened (or unrefined) and refined simultaneously for the subsequent ones, so as to not drastically increase the computational cost of the initial problem.

505

**Declaration of competing interest**

506

507

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

508

**Data availability**

509

No data was used for the research described in this article.

510

**Acknowledgments**

511

512

513

514

The authors would like to acknowledge the financial support received by FCT under the research project: Higher-order Immersed Boundary Method for Moving Body Problems (HIB-forMBP) with the reference PTDC/EME-EME/32315/2017. Institutional acknowledgment: This work was supported by FCT, through IDMEC, under LAETA, project UIDB/50022/2020.
## Appendix A. Least-squares fitting

The term least-squares describes a frequently used approach to solving overdetermined (or inexactely specified) systems of equations in an approximate sense; i.e., rather than solving the equations exactly, we seek to minimise the “energy” (or “sum of squares”) of the residual.

### Appendix A.1. Unconstrained least-squares (ULS) formulation

Consider the following (linear) system of equations

$$Ax = b, \qquad (A.1)$$

where $A \in \mathbb{R}^{m \times n}$ ($m \ge n$) has full (column) rank, $x \in \mathbb{R}^n$, and $b \in \mathbb{R}^m$.

If there is no solution to this system, we can always solve it in the “least-squares” sense. In other words, if a given system with more equations than unknowns cannot be solved because $b$ does not lie in the range of $A$, then the best approximation is one that minimises the “sum of squares” of the residuals ($r = b - Ax$):

$$ \begin{aligned} \min_{x} J(x) &= \min_{x} \|b - Ax\|_{2}^{2} \\ &\equiv \min_{x} \left( \sum_{i=1}^{m} (b_i - a_i^{\top} x)^2 \right). \end{aligned} \qquad (A.2) $$

Expanding the gradient of $J(x)$

$$ \begin{aligned} \nabla_x J(x) &= \nabla_x (b^\top b - 2b^\top Ax + x^\top A^\top Ax) \\ &= \underbrace{\nabla_x (b^\top b)}_{0} - 2\nabla_x (b^\top Ax) + \nabla_x (x^\top A^\top Ax) \\ &= -2A^\top b + [A^\top A + (A^\top A)^\top] x \\ &= -2A^\top b + 2A^\top Ax, \end{aligned} \qquad (A.3) $$

and setting it to zero gives

$$A^\top Ax^* = A^\top b, \qquad (A.4)$$

which is a symmetric (linear) system of equations known as the normal equations.

Since $A$ has full (column) rank (and provided that $b$ does not lie in the null space $N(A^*)$, we can now solve Eq. (A.4) for a unique solution $x^*$ as

$$x^* = (A^\top A)^{-1} A^\top b. \qquad (A.5)$$

In this derivation, we have used the condition $\nabla_x J(x) = 0$, which is a necessary but not sufficient condition for optimality. We found a critical point, but in general, such a point could be a local minimum, a local maximum, or a saddle point. Fortunately, the objective function is convex, which means that any critical point is indeed a local minimum. To show that $J(x)$ is convex, it suffices to compute its Hessian and show that it is positive semi-definite:

$$ \begin{aligned} \forall x, x^\top (2A^\top A)x &= 2(Ax)^\top Ax \\ &= 2\|Ax\|_2^2 \ge 0. \end{aligned} \qquad (A.6) $$
### Appendix A.1.1. Problem conditioning

Even though we have presented the normal equations as the most straightforward approach to solving overdetermined systems of equations, there are several undesirable aspects of this theoretical formulation (note that the normal equations will always be more poorly conditioned than the original system in Eq. (A.1)). In fact, (the $ℓ_2$-norm of) the condition number is squared:

$$ \kappa(\mathbf{A}^\mathrm{T}\mathbf{A}) = \kappa(\mathbf{A})^2. \qquad (A.7) $$

With finite-precision computation, the normal equations can actually become singular and $\mathbf{A}^\mathrm{T}\mathbf{A}$ non-existent. This inconvenience does not make the method of normal equations useless but provides motivation to seek more robust methods to solve least-squares problems.

### Appendix A.1.2. Solving the unconstrained least-squares (ULS) problem via QR decomposition

In this subsection, we consider the QR decomposition of $\mathbf{A} \in \mathbb{R}^{m \times n}$ into a product of an orthogonal matrix $\mathbf{Q} \in \mathbb{R}^{m \times m}$ and an upper triangular matrix $\mathbf{R} \in \mathbb{R}^{m \times n}$ as

$$ \mathbf{A} = \mathbf{QR} = \begin{bmatrix} \tilde{\mathbf{Q}} & \tilde{\mathbf{Q}}_{\perp} \end{bmatrix} \begin{bmatrix} \frac{\tilde{\mathbf{R}}}{\mathbf{0}} \end{bmatrix} = \tilde{\mathbf{Q}}\tilde{\mathbf{R}}, \qquad (A.8) $$

where the columns of $\tilde{\mathbf{Q}}_{\perp}$ span a vector space that is complementary and orthogonal to that spanned by $\tilde{\mathbf{Q}}$.

The ULS problem in Eq. (A.2) may now be written as follows:

$$ \begin{align} \min_{\mathbf{x}} \|\mathbf{b} - \mathbf{A}\mathbf{x}\|_2^2 &= \min_{\mathbf{x}} \|\mathbf{Q}^\mathrm{T}(\mathbf{b} - \mathbf{A}\mathbf{x})\|_2^2 \nonumber \\ &= \min_{\mathbf{x}} \left\| \begin{bmatrix} \tilde{\mathbf{Q}}^\mathrm{T} \\ \tilde{\mathbf{Q}}_{\perp} \end{bmatrix} \mathbf{b} - \begin{bmatrix} \tilde{\mathbf{R}} \\ \mathbf{0} \end{bmatrix} \mathbf{x} \right\|_2^2 \tag{A.9} \\ &= \min_{\mathbf{x}} \|\tilde{\mathbf{Q}}^\mathrm{T}\mathbf{b} - \tilde{\mathbf{R}}\mathbf{x}\|_2^2 + \|\tilde{\mathbf{Q}}_{\perp}^\mathrm{T}\mathbf{b}\|_2^2 \nonumber \\ &\ge \min_{\mathbf{x}} \|\tilde{\mathbf{Q}}^\mathrm{T}\mathbf{b}\|, \nonumber \end{align} $$

where $\tilde{\mathbf{Q}}^\mathrm{T}\mathbf{b}$ is the minimum residual norm of $\|\mathbf{b} - \mathbf{A}\mathbf{x}\|_2$.

Since $\mathbf{A}$ has full (column) rank, then $\tilde{\mathbf{R}}$ is non-singular, and thus we can solve the triangular system $\tilde{\mathbf{R}}\mathbf{x}^* = \tilde{\mathbf{Q}}^\mathrm{T}\mathbf{b}$ for a unique solution $\mathbf{x}^* = \tilde{\mathbf{R}}^{-1}\tilde{\mathbf{Q}}^\mathrm{T}\mathbf{b}$.

### Appendix A.2. Constrained least-squares (CLS) formulation

The (linearly) CLS problem refers to the problem of finding a “least-squares” solution that exactly satisfies $q$ additional constraints of the form

$$ \mathbf{C}\mathbf{x} = \mathbf{d}, \qquad (A.10) $$

where $\mathbf{C} \in \mathbb{R}^{q \times n}$ and $\mathbf{d} \in \mathbb{R}^q$.

Minimisation with constraints can be done using the method of Lagrange multipliers. Consider the Lagrangian function $\mathcal{L}(\mathbf{x}, \boldsymbol{\mu})$ defined as

$$ \mathcal{L}(\mathbf{x}, \boldsymbol{\mu}) = \|\mathbf{b} - \mathbf{A}\mathbf{x}\|_2^2 + \boldsymbol{\mu}^\mathrm{T}(\mathbf{C}\mathbf{x} - \mathbf{d}), \qquad (A.11) $$
where $\mu \in \mathbb{R}^q$ is a set of Lagrange multipliers. This method tells us that if $\mathbf{x}^*$ is a solution to the CLS problem, then there is a set of Lagrange multipliers $\mu^*$ that satisfies

$$ \nabla_{\mathbf{x}} \mathcal{L}(\mathbf{x}, \mu) = \mathbf{0} \qquad (A.12a) $$

$$ \nabla_{\mu} \mathcal{L}(\mathbf{x}, \mu) = \mathbf{0}. \qquad (A.12b) $$

These are the optimality conditions for the CLS problem; any feasible solution must satisfy them. Combining this set of (linear) equations, we can rewrite the optimality conditions as

$$ \left[ \frac{2\mathbf{A}^{\mathrm{T}}\mathbf{A}}{\mathbf{C}} \frac{\mathbf{C}^{\mathrm{T}}}{\mathbf{0}} \right] \left[ \frac{\mathbf{x}^{*}}{\mu^{*}} \right] = \left[ \frac{2\mathbf{A}^{\mathrm{T}}\mathbf{b}}{\mathbf{d}} \right]. \qquad (A.13) $$

The $(n+q) \times (n+q)$ matrix on the LHS of Eq. (A.13) is called the KKT matrix. It is invertible if and only if $\mathbf{C}$ has linearly independent rows and the stacked matrix $\left[ \frac{\mathbf{A}}{\mathbf{C}} \right]$ has linearly independent columns. While the first condition requires $\mathbf{C}$ to be wide (or square), i.e., that there are fewer constraints than variables, the second one depends on both $\mathbf{A}$ and $\mathbf{C}$, and it can be satisfied even when the columns of $\mathbf{A}$ are linearly dependent.

### Appendix A.2.1. Solving the constrained least-squares (CLS) problem via QR decomposition

Assuming that $\mathbf{C}$ has linearly independent rows and the stacked matrix $\left[ \frac{\mathbf{A}}{\mathbf{C}} \right]$ has linearly independent columns, we start by rewriting Eq. (A.13) as

$$ \left[ \frac{\mathbf{A}^{\mathrm{T}}\mathbf{A} + \mathbf{C}^{\mathrm{T}}\mathbf{C}}{\mathbf{C}} \frac{\mathbf{C}^{\mathrm{T}}}{\mathbf{0}} \right] \left[ \frac{\mathbf{x}^{*}}{\nu^{*}} \right] = \left[ \frac{\mathbf{A}^{\mathrm{T}}\mathbf{b}}{\mathbf{d}} \right], \qquad (A.14) $$

where $\nu^* = \frac{1}{2}\mu^* - \mathbf{d}$.

Next, we use QR factorisation

$$ \tilde{\mathbf{Q}}\tilde{\mathbf{R}} = \left[ \frac{\mathbf{A}}{\mathbf{C}} \right] = \left[ \frac{\tilde{\mathbf{Q}}_1}{\tilde{\mathbf{Q}}_2} \right] \tilde{\mathbf{R}} \qquad (A.15) $$

to simplify Eq. (A.14)

$$ \left[ \frac{\tilde{\mathbf{R}}^{\mathrm{T}}\tilde{\mathbf{R}}}{\tilde{\mathbf{Q}}_2\tilde{\mathbf{R}}} \frac{\tilde{\mathbf{R}}^{\mathrm{T}}\tilde{\mathbf{Q}}_2^{\mathrm{T}}}{\mathbf{0}} \right] \left[ \frac{\mathbf{x}^{*}}{\nu^{*}} \right] = \left[ \frac{\tilde{\mathbf{R}}^{\mathrm{T}}\tilde{\mathbf{Q}}_1^{\mathrm{T}}\mathbf{b}}{\mathbf{d}} \right]. \qquad (A.16) $$

Multiplying the first set of $m$ equations by $\tilde{\mathbf{R}}^{-\mathrm{T}}$ and rewriting Eq. (A.16) in terms of $\mathbf{y}^* = \tilde{\mathbf{R}}\mathbf{x}^*$ gives

$$ \left[ \frac{\mathbf{I}}{\tilde{\mathbf{Q}}_2} \frac{\tilde{\mathbf{Q}}_2^{\mathrm{T}}}{\mathbf{0}} \right] \left[ \frac{\mathbf{y}^{*}}{\nu^{*}} \right] = \left[ \frac{\tilde{\mathbf{Q}}_1^{\mathrm{T}}\mathbf{b}}{\mathbf{d}} \right]. \qquad (A.17) $$

Substituting the first set of $m$ equations of Eq. (A.17) into $\tilde{\mathbf{Q}}_2\mathbf{y}^* = \mathbf{d}$ gives the following set of equations in $\nu^*$:

$$ \tilde{\mathbf{Q}}_2 \tilde{\mathbf{Q}}_2^{\mathrm{T}} \nu^* = \tilde{\mathbf{Q}}_2 \tilde{\mathbf{Q}}_1^{\mathrm{T}} \mathbf{b} - \mathbf{d}. \qquad (A.18) $$

Next, we use QR factorisation

$$ \tilde{\mathbf{Q}}_2^{\mathrm{T}} = \tilde{\mathbf{T}}\tilde{\mathbf{U}} \qquad (A.19) $$
to simplify Eq. (A.18) and solve for $v^*$ as

$$v^* = \tilde{U}^{-1} [\tilde{T}^T \tilde{Q}_1^T b - \tilde{U}^{-T} d]. \quad (A.20)$$

Plugging this result into $y^* + \tilde{Q}_2^T v^* = \tilde{Q}_1^T b$ and solving $y^* = \tilde{R}x^*$ for a unique solution $x^*$ we get:

$$x^* = \tilde{R}^{-1} \{[\tilde{Q}_1^T - \tilde{Q}_2^T (\tilde{U}^{-1} \tilde{T}^T) \tilde{Q}_1^T] b + (\tilde{U}^{-1} \tilde{T}^T)^T d\}. \quad (A.21)$$

## Appendix B. Polynomial fit standardisation

When the X- (independent data) values are large and the degree of the polynomial approximation is high enough to produce numbers with drastically different orders of magnitude in the Vandermonde matrix (which is constructed during the fitting procedure), scaling problems may arise. To improve the precision of the computed coefficients (or parameters), the columns of the Vandermonde matrix can be brought to the same order of magnitude by centring the X-data at zero mean and scaling it to unit standard deviation:

$$\hat{X} = \frac{X - \bar{X}}{\sigma}, \quad (B.1)$$

where $\bar{X}$ and $\sigma$ are the mean and standard deviation of the original X-values, respectively.

2D polynomial regression models for curves are given by

$$\phi(x, y) = \sum_{i=0}^{n} \sum_{j=0}^{n-i} c_{ij} x^i y^j, \quad (B.2)$$

where $n$ is the degree of the polynomial fit.

Writing the polynomial function $\hat{\phi}(x, y)$ in terms of $(\hat{x}, \hat{y})$ and of the standardised coefficients (or parameters) $\hat{c}$ gives

$$\begin{aligned}
\hat{\phi}(x, y) &= \sum_{i=0}^{n} \sum_{j=0}^{n-i} \hat{c}_{ij} \hat{x}^i \hat{y}^j \\
&= \sum_{i=0}^{n} \sum_{j=0}^{n-i} \hat{c}_{ij} \left( \frac{x - \bar{x}}{\sigma_x} \right)^i \left( \frac{y - \bar{y}}{\sigma_y} \right)^j \quad (B.3) \\
&= \sum_{i=0}^{n} \sum_{j=0}^{n-i} \frac{\hat{c}_{ij}}{(\sigma_x)^i (\sigma_y)^j} \sum_{k=0}^{i} \sum_{l=0}^{j} \binom{i}{k} \binom{j}{l} x^{i-k} y^{j-l} (-\bar{x})^k (-\bar{y})^l.
\end{aligned}$$

To relate the coefficients (or parameters) of $\hat{\phi}(x, y)$ with those of $\phi(x, y)$, we use the following transformation:

$$c_{ij} = \sum_{k=i}^{n} \sum_{l=j}^{n-k} \binom{k}{k-i} \binom{l}{l-j} \frac{\hat{c}_{kl}}{(\sigma_x)^k (\sigma_y)^l} (-\bar{x})^{k-i} (-\bar{y})^{l-j}. \quad (B.4)$$

In a weighted least-squares formulation, we also need to rewrite the weighting function in terms of $(\hat{x}, \hat{y})$. The inverse-distance weighting function proposed in Section 4.3.1 now becomes

$$\hat{w}_i = \frac{1}{\hat{d}_i^n}, \quad (B.5)$$
where the standardised distance between the centroid of face f and the reference location of the ith element of the discretisation stencil is given by

$$ \begin{aligned} \hat{d}_i &= \sqrt{\left(\frac{x_i - \mu}{\sigma_x} - \frac{x_f - \mu}{\sigma_x}\right)^2 + \left(\frac{y_i - \mu}{\sigma_y} - \frac{y_f - \mu}{\sigma_y}\right)^2} \quad (B.6) \\ &\equiv \sqrt{\left(\frac{x_i - x_f}{\sigma_x}\right)^2 + \left(\frac{y_i - y_f}{\sigma_y}\right)^2}. \end{aligned} $$

As shown in Figs. B.27 and B.28, (the ℓ₂-norm of) the condition number of the Vandermonde D_f and weighted Vandermonde $\sqrt{\mathbf{W}_f}\mathbf{D}_f$ matrices defined in Section 4.3 no longer depends on the grid spacing, only on the degree of the polynomial fit and on the exponent to which the standardised distance is elevated.

<figure><img src="image_29.png" /><figcaption>Figure B.27. (Matrix) D_f conditioning for polynomial fits of different degree.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Figure B.28. (Matrix) $\sqrt{\mathbf{W}_f}\mathbf{D}_f$ conditioning for different inverse-distance weighting functions (linear polynomial fit).</figcaption></figure>

In the context of the problem at hand, the standardisation procedure introduced here is especially useful when dealing with extremely dense grids and/or polynomial fits of very high order.
604 References

605 [1] W. L. Oberkampf, T. G. Trucano, Verification and validation benchmarks, Nuclear Engineering and Design 238
606 (2008) 716-743. Benchmarking of CFD Codes for Application to Nuclear Reactor Safety.

607 [2] J. Abanto, D. Pelletier, A. Garon, J.-Y. Trepanier, M. Reggio, Verification of some Commercial CFD Codes
608 on Atypical CFD Problems, in: 43rd AIAA Aerospace Sciences Meeting and Exhibit, American Institute of
609 Aeronautics and Astronautics, Reno, Nevada, 2005.

610 [3] D. F. Hawken, J. J. Gottlieb, J. S. Hansen, Review of some adaptive node-movement techniques in finite-element
611 and finite-difference solutions of partial differential equations, Journal of Computational Physics 95 (1991) 254-
612 302.

613 [4] Y. Li, S. Premasuthan, A. Jameson, Comparison of Adaptive h and p Refinements for Spectral Difference Methods,
614 in: 40th Fluid Dynamics Conference and Exhibit, American Institute of Aeronautics and Astronautics, Chicago,
615 Illinois, 2010.

616 [5] J. P. P. Magalhães, D. M. Albuquerque, J. M. Pereira, J. C. Pereira, Adaptive mesh finite-volume calculation of 2d
617 lid-cavity corner vortices, Journal of Computational Physics 243 (2013) 365-381.

618 [6] D. M. S. Albuquerque, J. M. C. Pereira, J. C. F. Pereira, Residual Least-Squares Error Estimate for Unstructured
619 h-Adaptive Meshes, Numerical Heat Transfer, Part B: Fundamentals 67 (2015) 187-210.

620 [7] C. Burg, Analytic study of 2D and 3D grid motion using modified Laplacian, International Journal for Numerical
621 Methods in Fluids 52 (2006) 163-197.

622 [8] L. D. Lorenzis, A. Düster, Modeling in Engineering Using Innovative Numerical Methods for Solids and Fluids,
623 CISM International Centre for Mechanical Sciences, 1st ed., Springer, 2021.

624 [9] A. Patra, J. T. Oden, Computational techniques for adaptive hp finite element methods, Finite Elements in Analysis
625 and Design 25 (1997) 27-39. Adaptive Meshing, Part 1.

626 [10] Z. J. Wang, K. Fidkowski, R. Abgrall, F. Bassi, D. Caraeni, A. Cary, H. Deconinck, R. Hartmann, K. Hillewaert,
627 H. T. Huynh, N. Kroll, G. May, P.-O. Persson, B. van Leer, M. Visbal, High-order CFD methods: current status
628 and perspective, International Journal for Numerical Methods in Fluids 72 (2013) 811-845.

629 [11] Z. J. Wang, High-order computational fluid dynamics tools for aircraft design, Philosophical Transactions of the
630 Royal Society A: Mathematical, Physical and Engineering Sciences 372 (2014) 20130318.

631 [12] M. Kronbichler, P.-O. Persson, Efficient High-Order Discretizations for Computational Fluid Dynamics, volume
632 602 of CISM International Centre for Mechanical Sciences, 1st ed., Springer, 2021.

633 [13] J. Loffeld, J. A. F. Hittinger, On the arithmetic intensity of high-order finite-volume discretizations for hyperbolic
634 systems of conservation laws, The International Journal of High Performance Computing Applications 33 (2019)
635 25-52.

636 [14] R. E. Bank, A. H. Sherman, A. Weiser, Some Refinement Algorithms And Data Structures For Regular Local Mesh
637 Refinement, Scientific Computing 1 (1983) 3-17.

638 [15] J. A. White, H. Nishikawa, R. A. Baurle, Weighted Least-squares Cell-Average Gradient Construction Methods For
639 The VULCAN-CFD Second-Order Accurate Unstructured Grid Cell-Centered Finite-Volume Solver, in: AIAA
640 Scitech 2019 Forum, American Institute of Aeronautics and Astronautics, San Diego, California, 2019.

641 [16] P. Tsoutsanis, V. A. Titarev, D. Drikakis, WENO schemes on arbitrary mixed-element unstructured meshes in three
642 space dimensions, Journal of Computational Physics 230 (2011) 1585-1601.

643 [17] M. Dumbser, M. Käser, V. A. Titarev, E. F. Toro, Quadrature-free non-oscillatory finite volume schemes on un-
644 structured meshes for nonlinear hyperbolic systems, Journal of Computational Physics 226 (2007) 204-243.

645 [18] S. Diot, S. Clain, R. Loubère, Improved detection criteria for the Multi-dimensional Optimal Order Detection
646 (MOOD) on unstructured meshes with very high-order polynomials, Computers & Fluids 64 (2012) 43-63.

647 [19] S. Clain, S. Diot, R. Loubère, A high-order finite volume method for systems of conservation laws—Multi-
648 dimensional Optimal Order Detection (MOOD), Journal of Computational Physics 230 (2011) 4028-4050.

649 [20] M. Dumbser, M. Castro, C. Parés, E. F. Toro, ADER schemes on unstructured meshes for nonconservative hyper-
650 bolic systems: Applications to geophysical flows, Computers & Fluids 38 (2009) 1731-1748.

651 [21] A. Jalali, C. F. O. Gooch, Higher-Order Finite Volume Solution Reconstruction on Highly Anisotropic Meshes, in:
652 21st AIAA Computational Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, San
653 Diego, California, 2013.

654 [22] X. Nogueira, L. Cueto-Felgueroso, I. Colominas, F. Navarrina, M. Casteleiro, A new shock-capturing technique
655 based on Moving Least Squares for higher-order numerical schemes on unstructured grids, Computer Methods in
656 Applied Mechanics and Engineering 199 (2010) 2544-2558.

657 [23] M. C. Bessemoulin-Chatard, Development and analysis of finite volume schemes motivated by the preservation of
658 asymptotic behaviors. Application to models from physics and biology., PhD dissertation, Université Blaise Pascal
659 - Clermont-Ferrand II, 2012.

660 [24] C. Roy, Strategies for Driving Mesh Adaptation in CFD (Invited), in: 47th AIAA Aerospace Sciences Meeting in-
661 including The New Horizons Forum and Aerospace Exposition, American Institute of Aeronautics and Astronautics, Orlando, Florida, 2009.

663 [25] A. Hay, M. Visonneau, Error estimation using the error transport equation for finite-volume methods and arbitrary meshes, International Journal of Computational Fluid Dynamics 20 (2006) 463–479.

665 [26] A. Hay, A. Leroyer, M. Visonneau, H-adaptive Navier-Stokes simulations of free-surface flows around moving bodies, Journal of Marine Science and Technology 11 (2006) 1–18.

667 [27] V. Ervin, W. Layton, An Analysis of a Defect-Correction Method for a Model Convection-Diffusion Equation, SIAM Journal on Numerical Analysis 26 (1989) 169–179.

669 [28] N. A. Pierce, M. B. Giles, Adjoint and defect error bounding and correction for functional estimates, Journal of Computational Physics 200 (2004) 769–794.

671 [29] J. H. Ferziger, M. Perić, Computational Methods for Fluid Dynamics, 3rd ed., Springer, 2002.

672 [30] A. Syrakos, G. Efthimiou, J. G. Bartzis, A. Goulas, Numerical experiments on the efficiency of local grid refinement based on truncation error estimates, Journal of Computational Physics 231 (2012) 6725–6753.

674 [31] M. Aftosmis, M. Berger, Multilevel error estimation and adaptive h-refinement for Cartesian meshes with embedded boundaries, in: 40th AIAA Aerospace Sciences Meeting & Exhibit, American Institute of Aeronautics and Astronautics, Reno, Nevada, 2002.

677 [32] W. Bangerth, R. Rannacher, Adaptive Finite Element Methods for Differential Equations, 1st ed., Birkhäuser, 2003.

678 [33] Z. Wang, A perspective on high-order methods in computational fluid dynamics, Science China Physics, Mechanics & Astronomy 59 (2015) 614701.

679 [34] S. Muzaferija, D. Gosman, Finite-Volume CFD Procedure and Adaptive Error Control Strategy for Grids of Arbitrary Topology, Journal of Computational Physics 138 (1997) 766–787.

680 [35] T. J. Barth, D. C. Jespersen, The design and application of upwind schemes on unstructured meshes, in: 27th AIAA Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, Moffett Field, California, 1989.

685 [36] C. Ollivier-Gooch, M. V. Altena, A High-Order-Accurate Unstructured Mesh Finite-Volume Scheme for the Advection–Diffusion Equation, Journal of Computational Physics 181 (2002) 729–752.

687 [37] H. T. Nguyen, p-Adaptive and Automatic h-Adaptive Finite Element Methods for Elliptic Partial Differential Equations, PhD dissertation, University of California, San Diego, 2010.

689 [38] P. M. P. Costa, D. M. S. Albuquerque, A novel approach for temporal simulations with very high-order finite volume schemes on polyhedral unstructured grids, Journal of Computational Physics 453 (2022) 110960.