The surrogate matrix methodology: a priori error estimation

TL;DR

This paper rigorously analyzes the surrogate matrix methodology for finite element analysis, demonstrating its efficiency and broad applicability through theoretical error bounds and numerical experiments showing significant computational speedups.

Contribution

It provides the first rigorous a priori error analysis of the surrogate matrix methodology and extends its application to nonlinear and complex geometries.

Findings

Up to twenty-fold reduction in computation time.

Validated convergence for variable coefficient Poisson equation.

Effective in matrix-free environments with multigrid solvers.

Abstract

We give the first mathematically rigorous analysis of an emerging approach to finite element analysis (see, e.g., Bauer et al. [Appl. Numer. Math., 2017]), which we hereby refer to as the surrogate matrix methodology. This methodology is based on the piece-wise smooth approximation of the matrices involved in a standard finite element discretization. In particular, it relies on the projection of smooth so-called stencil functions onto high-order polynomial subspaces. The performance advantage of the surrogate matrix methodology is seen in constructions where each stencil function uniquely determines the values of a significant collection of matrix entries. Such constructions are shown to be widely achievable through the use of locally-structured meshes. Therefore, this methodology can be applied to a wide variety of physically meaningful problems, including nonlinear problems and problems with curvilinear geometries. Rigorous a priori error analysis certifies the convergence of a novel surrogate method for the variable coefficient Poisson equation. The flexibility of the methodology is also demonstrated through the construction of novel methods for linear elasticity and nonlinear diffusion problems. In numerous numerical experiments, we demonstrate the efficacy of these new methods in a matrix-free environment with geometric multigrid solvers. In our experiments, up to a twenty-fold decrease in computation time is witnessed over the classical method with an otherwise identical implementation.

Tables

Table 1. Table 1: Relative H 1 superscript 𝐻 1 H^{1} errors and experimental orders of convergence for fixed h ℎ h and varying q 𝑞 q and H 𝐻 H in the case of problem  65 with the scalar coefficient  66 . Here, the relative H 1 superscript 𝐻 1 H^{1} error with the classical FEM is 1.23 × 10 − 08 1.23E-08 1.23\text{\times}{10}^{-08} .

	$q=1$		$q=2$		$q=3$		$q=4$
$\frac{H}{H_0}$	rel. $H^1$ err.	eoc	rel. $H^1$ err.	eoc	rel. $H^1$ err.	eoc	rel. $H^1$ err.	eoc
\csvreader[late after line=
,late after last line=
\eocFour

Table 2. Table 2: Relative L 2 superscript 𝐿 2 L^{2} errors and experimental orders of convergence for fixed h ℎ h and varying q 𝑞 q and H 𝐻 H in the case of problem  65 with the scalar coefficient  66 . Here, the relative L 2 superscript 𝐿 2 L^{2} error with the classical FEM is 4.10 × 10 − 09 4.10E-09 4.10\text{\times}{10}^{-09} .

	$q=1$		$q=2$		$q=3$		$q=4$
$\frac{H}{H_0}$	rel. $L^2$ err.	eoc	rel. $L^2$ err.	eoc	rel. $L^2$ err.	eoc	rel. $L^2$ err.	eoc
\csvreader[late after line=
,late after last line=
\eocFour

Table 3. Table 3: Relative H 1 superscript 𝐻 1 H^{1} errors and experimental orders of convergence for fixed h ℎ h and varying q 𝑞 q and H 𝐻 H in the case of problem  65 with the tensorial coefficient  67 and curved boundary. The relative H 1 superscript 𝐻 1 H^{1} error with the classical FEM is 1.59 × 10 − 08 1.59E-08 1.59\text{\times}{10}^{-08} .

	$q=1$		$q=2$		$q=3$		$q=4$
$\frac{H}{H_0}$	rel. $H^1$ err.	eoc	rel. $H^1$ err.	eoc	rel. $H^1$ err.	eoc	rel. $H^1$ err.	eoc
\csvreader[late after line=
,late after last line=
\eocFour

Table 4. Table 4: Relative L 2 superscript 𝐿 2 L^{2} errors and experimental orders of convergence for fixed h ℎ h and varying q 𝑞 q and H 𝐻 H in the case of problem  65 with the tensorial coefficient  67 and curved boundary. The relative L 2 superscript 𝐿 2 L^{2} error with the classical FEM is 4.56 × 10 − 09 4.56E-09 4.56\text{\times}{10}^{-09} .

	$q=1$		$q=2$		$q=3$		$q=4$
$\frac{H}{H_0}$	rel. $L^2$ err.	eoc	rel. $L^2$ err.	eoc	rel. $L^2$ err.	eoc	rel. $L^2$ err.	eoc
\csvreader[late after line=
,late after last line=
\eocFour

Table 5. Table 5: Relative L 2 superscript 𝐿 2 L^{2} errors, experimental orders of convergence, and relative time-to-solutions for fixed H 𝐻 H and varying q 𝑞 q and h ℎ h in the case of problem  65 with the tensorial coefficient  67 and curved boundary.

standard

$q=3$

$q=5$

$q=7$

$\frac{h}{H_0}$

$\frac{H_{\mathrm{LS}}}{h}$

DoFs

rel. $L^2$ err.

eoc

rel. $L^2$ err.

eoc

rtts

rel. $L^2$ err.

eoc

rtts

rel. $L^2$ err.

eoc

rtts

\csvreader[late after line=

,late after last line=

\relttsseven

Table 6. Table 6: Relative L 2 superscript 𝐿 2 L^{2} errors, experimental orders of convergence, and relative time-to-solutions for fixed H 𝐻 H , q = 7 𝑞 7 q=7 , and varying h ℎ h in the case of the linearized elasticity problem  69 .

			standard		$q=7$
$\frac{h}{H}$	$\frac{H_{\mathrm{LS}}}{h}$	DoFs	rel. $L^2$ err.	eoc	rel. $L^2$ err.	eoc	rtts
\csvreader[late after line=
,late after last line=
\relttsseven

Equations

$$\text{Find }u\in V\text{ satisfying }$$

$$\text{Find }{u}_h\in V_h\text{ satisfying }$$

$$\text{Find }{\tilde{u}}_h\in V_h\text{ satisfying }$$

$$\mathsf{Au}=\mathsf{f}\text{ and }\tilde{\mathsf{A}}\tilde{\mathsf{u}}=\mathsf{f},$$

$$a(u,v)={\int }_{\Omega }G(x,u(x),v(x))\mathrm{d}x\text{for all }u,v\in V.$$

$$G(x,u(y),v(y))=0\text{whenever }y\notin \mathrm{supp}(u)\cap \mathrm{supp}(v).$$

$$a_1(u,v)={\int }_{\Omega }\nabla u(x{)}^{\top }K(x)\nabla v(x)\mathrm{d}x\text{for all }u,v\in V=H^1(\Omega ).$$

$$a({\varphi }_{x_j},{\varphi }_{x_i})={\int }_{\Omega }G(y,{\varphi }_{x_j}(y),{\varphi }_{x_i}(y))\mathrm{d}y$$

$${\Phi }_i^{\delta }(x)={\int }_{{\Omega }_{\delta }}G(x+y,{\varphi }_{\delta }(y),\varphi (y))\mathrm{d}y\text{for each }\delta \in \mathcal{D}(x_i).$$

$$a({\varphi }_{x_j},{\varphi }_{x_i})=\{\begin{array}{ll}{\Phi }_i^{\delta }(x_i),& \text{if }\delta =x_j-x_i\in \mathcal{D}(x_i),\\ 0,& \text{otherwise.}\end{array}$$

$${\Phi }_i^{\delta }(x_i)={\int }_{{\Omega }_{\delta }}G(x_i+y,\varphi (y),{\varphi }_{\delta }(y))\mathrm{d}y={\int }_{\Omega }G(x_j+y,{\varphi }_{-\delta }(y),\varphi (y))\mathrm{d}y={\Phi }_j^{-\delta }(x_j)$$

$$a({\varphi }_{x_j},{\varphi }_{x_i})=\{\begin{array}{ll}{\Phi }^{\delta }(x_i),& \text{if }\delta =x_j-x_i\in \mathcal{D}(\mathbb{L}),\\ 0,& \text{otherwise.}\end{array}$$

$$a({\varphi }_j,{\varphi }_i)=\{\begin{array}{ll}{\Phi }^{\delta }(x_i),& \text{if }\delta =x_j-x_i\in \{dx,0,dx\},\\ 0,& \text{otherwise,}\end{array}$$

$${\mathbb{X}}_m=\bigcup \{T_m:T\in {\mathcal{T}}_H\}.$$

$$\mathcal{S}(T)=\big{\{}[y_{0},\tfrac{y_{0}+y_{1}}{2},\tfrac{y_{0}+y_{2}}{2}],[\tfrac{y_{0}+y_{1}}{2},y_{1},\tfrac{y_{1}+y_{2}}{2}],[\tfrac{y_{0}+y_{2}}{2},\tfrac{y_{1}+y_{2}}{2},y_{2}],[\tfrac{y_{0}+y_{1}}{2},\tfrac{y_{1}+y_{2}}{2},\tfrac{y_{0}+y_{2}}{2}]\big{\}}.$$

$${\mathcal{S}}^{m+1}(T)=\bigcup \{\mathcal{S}(t):t\in {\mathcal{S}}^m(T)\}\text{ for all }m\ge 1.$$

$${\mathcal{S}}^m({\mathcal{T}}_H)=\bigcup \{{\mathcal{S}}^m(T):T\in {\mathcal{T}}_H\},\text{ for all }m\ge 1.$$

$$a({\varphi }_j,{\varphi }_i)=\{\begin{array}{ll}{\Phi }_T^{\delta }(x_i),& \text{if }\delta =(x_j-x_i)\in \mathcal{D}({\mathring{T}}_m),\\ 0,& \text{otherwise,}\end{array}$$

$$a_{1,H}(u,v)={\int }_{{\Omega }_H}\nabla u(x{)}^{\top }{K}_H(x)\nabla v(x)\mathrm{d}x,K_H=\frac{D{\phi }^{-1}(K\circ \phi )D{\phi }^{-\top }}{|\det \left(D{\phi }^{-1}\right)|}.$$

$$T_{\delta }=\{x\in T:x+y\in T,\text{ for all }y\in {\Omega }_{\delta }=\mathrm{supp}(\varphi )\cap \mathrm{supp}({\varphi }_{\delta })\}.$$

$$G(x+y,{\varphi }_{\delta }(y),\varphi (y))=\sum _{|\alpha |\le l}{c}_{\alpha }(y)(x+y{)}^{\alpha }=\sum _{|\alpha |\le l}\sum _{|\nu |\le \alpha }(\genfrac{}{}{0ex}{}{\alpha }{\nu })c_{\alpha }(y)y^{\nu }{x}^{\alpha -\nu }$$

$${\tilde{\mathsf{A}}}_{ij}=\{\begin{array}{ll}{\int }_{\Omega }G(y,{\varphi }_j(y),{\varphi }_i(y))\mathrm{d}y,& \text{if both }{x}_i\text{ and }{x}_j\in \partial {\mathbb{X}}_m\\ {\tilde{\Phi }}_T^{\delta }(x_i),& \text{if }\delta =(x_j-x_i)\in \mathcal{D}({\mathring{T}}_m)\text{ and }{x}_i\text{ or }{x}_j\in {\mathring{T}}_m,\\ 0,& \text{otherwise.}\end{array}$$

$$\Big{[}\Pi^{\delta}_{T}\Phi^{\delta}_{T}\Big{]}(x)=\Big{[}\Pi^{-\delta}_{T}\Phi^{-\delta}_{T}\Big{]}(x+\delta)\,,$$

$${\tilde{\mathsf{A}}}_{ij}={\Pi }_T^{\delta }{\Phi }_T^{\delta }(x_i)={\Pi }_T^{-\delta }{\Phi }_T^{-\delta }(x_i+\delta )={\Pi }_T^{-\delta }{\Phi }_T^{-\delta }(x_j)={\tilde{\mathsf{A}}}_{ji}.$$

$${\mathsf{A}}_{ij}={\int }_{\Omega }\nabla {\varphi }_i(x{)}^{\top }K(x)\nabla {\varphi }_j(x)\mathrm{d}x$$

$$a_{2}(u,v)=\int_{\Omega}2\mu\epsilon(u):\epsilon(v)+\lambda\hskip 0.75pt\mathrm{div}(u)\,\mathrm{div}(v)\hskip 0.5pt\,\mathrm{d}x\qquad\text{for all }u,v\in\big{[}H^{1}_{0}(\Omega)\big{]}^{n}.$$

$$b(u,v)={\int }_{\Omega }dt|\nabla \tilde{u}{|}^{p-2}\nabla u\cdot \nabla v+uv\mathrm{d}x,\text{for all }u,v\in W^{1,p}(\Omega ).$$

$$m(u,v)={\int }_{\Omega }uv\mathrm{d}x\text{and}{a}_3(u,v)={\int }_{\Omega }|\nabla \tilde{u}{|}^{p-2}\nabla u\cdot \nabla v\mathrm{d}x.$$

$$\mathsf{A}_{ii}-\widetilde{\mathsf{A}}_{ii}=-\sum_{j\neq i}\big{(}\mathsf{A}_{ij}^{\mathrm{ext}}-\widetilde{\mathsf{A}}_{ij}^{\mathrm{ext}}\big{)}=-\sum_{j\neq i}\big{(}\mathsf{A}_{ij}-\widetilde{\mathsf{A}}_{ij}\big{)}\,.$$

$$\underset{v_h\in V_h\cap S}{\sup }a(w_h,v_h)\ge \alpha \Vert w_h{\Vert }_V\text{for all }{w}_h\in V_h.$$

$$\underset{v_h\in V_h\cap S}{\sup }\tilde{a}(w_h,v_h)\ge \tilde{\alpha }\Vert w_h{\Vert }_V\text{for all }{w}_h\in V_h.$$

$$|{\lambda }_k(\mathsf{M})-{\lambda }_k(\mathsf{N})|\le \Vert \mathsf{M}-\mathsf{N}{\Vert }_{\infty }.$$

Full text

{AMS}

65D05, 65M60, 65N30, 65Y05, 65Y20.

1 Introduction

In the field of computational science, major funding initiatives in North America, Europe, and Asia have thrust high performance computing (HPC) to ascendancy. In anticipation of future exascale computers, much work in this discipline involves the deep and careful reconstruction of long-established computing practices. An important characteristic of numerical algorithms aimed for these computers, is the floating-point operation (FLOP) per byte ratio. In order to achieve optimal performance and power efficiency on future machines, the time spent on FLOPs relative to memory transfer needs to be substantial.

Most traditional finite element softwares assemble global stiffness matrices by looping over elements and adding the corresponding local contributions to the global matrix. Storing the resulting sparse matrices requires significantly more memory than just storing the degrees of freedom. However, memory consumption is certainly not the only obstacle at the computational frontier. Indeed, at such scales, the memory traffic and latency involved in loading indices and entries for matrix vector products (MVPs) also presents critical challenges.

Since iterative solvers only require MVPs, it is not necessary to store all of the nonzeros of the global matrix in memory. Instead, it is sufficient to compute the nonzero entries on-the-fly, i.e., matrix-free. Different tactics exist to implement matrix-free methods, but the predominant candidate for low-order finite elements is the element-by-element approach [3, 15, 19, 24, 38], wherein local stiffness matrices are multiplied by local vectors and later added to the global solution vector. These local stiffness matrices may either be stored in memory—which actually requires more memory than storing the global matrix—or computed on-the-fly. When using high-order finite elements, the weak forms can be integrated on-the-fly using standard or reduced quadrature formulas [17, 30, 31, 32, 34]. This is a well-suited tactic for future machines because of its large arithmetic intensity [33].

Significant performance gains are often attributed to a problem, scale, and architecture-specific balance between FLOPs and memory traffic. As a matter of course, exploiting symmetries in a problem or discretization can significantly improve the time to solution. This is often the cause of enormous speed-ups in computations on structured meshes. Likewise, high performance of matrix-free methods can be most easily achieved in homogeneous problems with simple geometries. Nevertheless, most geometries coming from significant real-world problems cannot be adequately approximated by fully-structured meshes. A possible trade-off is to use locally-structured meshes like hierarchical hybrid grids (HHG) where initially unstructured coarse grids are locally refined in a uniform way. This local structure allows the application of stencil-based finite element procedures which operate similar to finite difference methods. By using these grids, efficient stencil-based methods have been successfully applied to a wide range of problems [11, 12, 13, 22, 26]. A related approach, suitable for low-order finite element discretizations of scalar elliptic partial differential equations (PDEs) with variable coefficients, based on scaling of reference stencils is discussed in [7].

In this paper, we revisit the classical lowest-order Bubnov–Galerkin finite element method and analyze a modification of it which is strongly amenable to stencil-based matrix-free computation. In our approach, a macro-mesh, which is not required to have any global structure, is used to triangulate the model geometry. This macro-mesh is then uniformly refined a large number of times, resulting in a fine-scale locally-structured mesh. For each macro-element, a local approximation of the fine-scale global matrix delivers a fine-scale surrogate matrix which maintains the convergence properties of the fine-scale discrete solution, up to the original order of the approximation. A related investigation [10] illustrated the promise of this methodology and provided numerical evidence for the convergence rates which are proven here rigorously. This work considered only Poisson’s equation. Later on, Stokes flow (with variable viscosity) was considered in two follow-up articles [8, 9]. Each of these initial studies focused on the massively parallel high performance computing aspects of their respective methods. These studies used the HHG software framework [11, 12, 13] in their experiments. Here, the finite-element software framework HyTeG [29] is used.

In this paper, we recast the central features of the original work as a methodology complete with a mathematical framework suitable for rigorous analysis. The principal novelty is the mathematical foundation developed here, which can be used to analyze further incarnations of the methodology. In total, we consider three specific mathematical models; namely, the variable coefficient Poisson equation, linear elastostatics, and $p$-Laplacian diffusion. Although our presentation demonstrates that the surrogate matrix methodology applies to each of these models equally well, we only employ a complete a priori analysis of the simplest model, the variable coefficient Poisson equation. In our numerical experiments, we carefully verify the proven a priori convergence rates with the variable coefficient Poisson equation. We also include proof-of-concept demonstrations from numerical experiments with the linear elastostatics and $p$-Laplacian diffusion problems.

2 Notation and outline

Let $V$ be a reflexive Banach space over $\mathbb{R}$, the field of real numbers, and let $V_{h}\subsetneq V$ be a finite-dimensional subspace. Consider a continuous and weakly coercive bilinear form $a:V\times V\to\mathbb{R}$ and a bounded linear functional $F\in V^{\ast}$, the topological dual of $V$.

In this paper, we are concerned with the solutions $u$, $u_{h}$, and $\widetilde{u}_{h}$ of the following three abstract variational problems.

[TABLE]

In 1c, a surrogate bilinear form $\widetilde{a}:V_{h}\times V_{h}\to\mathbb{R}$ has been introduced. In order to properly define $\widetilde{a}(\cdot,\cdot)$, some additional assumptions on $a(\cdot,\cdot)$ are still required; see Section 3.

The discrete variational problems 1b and 1c induce matrix equations for coefficients $\mathsf{u},\widetilde{\mathsf{u}}$ in some $\mathbb{R}^{N}$,

[TABLE]

respectively. In the first case, fix a basis for $V_{h}$, say $\{\phi_{i}\}_{i=1}^{N}$. For this basis, each $(i,j)$-component of the stiffness matrix $\mathsf{A}$ is simply $\mathsf{A}_{ij}=a(\phi_{j},\phi_{i})$. In the following section, we present a methodology to construct a surrogate stiffness matrix $\widetilde{\mathsf{A}}\approx\mathsf{A}$ which can be used in place of the true stiffness matrix $\mathsf{A}$. This methodology stands apart from technical details, such as differences in quadrature formulas. Section 4 provides a short (non-comprehensive) list of examples fitting into our framework. In Section 5, we discuss the incorporation of non-homogeneous boundary conditions and what we hereon refer to as the zero row sum property. In Section 6, sufficient conditions for the discrete stability of surrogate bilinear forms $\widetilde{a}(\cdot,\cdot)$ are briefly discussed. Next, in Section 7, we perform a rigorous a priori error analysis of our approach applied to the variable coefficient Poisson equation. A brief description of our implementation is given in Section 8. Then, in Section 9, we document several numerical experiments. Here, a thorough verification of each error estimate in Section 7 is given. This is complemented by performance measurements for the additional examples.

Throughout this article, we assume that $\Omega\operatorname{\subseteq}\mathbb{R}^{n}$ is a bounded Lipschitz domain. For matrices $\mathsf{M}\in\mathbb{R}^{l\times m}$, define the $\ell^{\infty}$-, and ${\operatorname*{\vphantom{p}max}}$-norms, $\|\mathsf{M}\|_{\infty}=\operatorname*{\vphantom{p}max}_{i}\sum_{j}|\mathsf{M}_{ij}|$ and $\|\mathsf{M}\|_{\operatorname*{\vphantom{p}max}}=\operatorname*{\vphantom{p}max}_{i,j}|\mathsf{M}_{ij}|$. Likewise, for any function $v:\Omega\to\mathbb{R}$, we will use the similar notation, $\|v\|_{0}$, $\|v\|_{1}$, and $\|v\|_{2}$, for the canonical $L^{2}(\Omega)$-, $H^{1}(\Omega)$-, and $H^{2}(\Omega)$-norms, respectively. When dealing with a subset $T\operatorname{\subseteq}\Omega$, denote the related $L^{2}(T)$-, $H^{1}(T)$-, and $H^{2}(T)$-norms by $\|v\|_{0,T}$, $\|v\|_{1,T}$, and $\|v\|_{2,T}$, respectively. For any simplex $T$ and integer $0\leq q<\infty$, we denote the space of polynomials of degree at most $q$ as $\mathcal{P}_{q}(T)$. All remaining notation will be defined as it arises.

3 Surrogate stiffness matrices

In this section, we present the constitutive elements of the surrogate matrix methodology. Our approach here is to gradually introduce the necessary concepts, all the while maintaining a clear sense of generality. In order to arrive at a tractable framework for our problems of interest, we gradually refine the presentation from general $n$-dimensional spaces to only $n=1,2$ or $3$ and from general Banach spaces to only $W^{1,p}(\Omega)$ (or products thereof), where $1<p<\infty$. The intention of proceeding in this way is to indicate that the methodology can be applied to an extremely broad set of problems and, specifically, to most problems where finite element methods are traditionally applied.

3.1 Preliminary assumptions

Given a bounded domain $\Omega\operatorname{\subseteq}\mathbb{R}^{n}$, assume that the true bilinear form can be expressed as

[TABLE]

These assumptions permit us to consider the discretization of most classical differential operators. Indeed, in the assumptions above, the integrand $G(x,u,v)$ may induce distributional derivatives on its second and third arguments. Meanwhile, the first argument can be identified with the spatial argument of any associated variable coefficients. For example, in the weak form of a Poisson-type equation, $-\operatorname{div}(K\nabla u)=f$, with a variable, symmetric positive-definite diffusion tensor $K(x)$ (cf. Section 4.1), we have the bilinear form

[TABLE]

Here, taking any point $x\in\Omega$, the integrand in 3a reduces to $G(x,u,v)={\nabla u^{\top}K(x)\hskip 0.25pt\nabla v}$. Evidently, this $G$ satisfies the sparsity assumption 3b.

3.2 Stencil functions

Let $\phi\in V$ be a test function with compact support in $\Omega$ and, for any fixed $y\in\mathbb{R}^{n}$, define $\phi_{y}(x)=\phi(x-y)$. Now, consider any fixed set of ordered points $\mathbb{X}=\{x_{i}\}$ in $\Omega$ and recall 3a. Assuming that both $\phi_{x_{i}},\phi_{x_{j}}\in V$, observe (via a simple change of variables) that

[TABLE]

where $\delta=x_{j}-x_{i}$ and $\Omega_{\delta}=\operatorname*{supp}(\phi)\cap\operatorname*{supp}(\phi_{\delta})$. In the second equality, passing from an integral over $\Omega$ to an integral over the subset $\Omega_{\delta}\operatorname{\subseteq}\Omega$ follows immediately from the sparsity assumption 3b. For each fixed $x_{i}$, the affine structure of the identity above may be illuminated by collecting each contributing translation $\delta$ into the set $\mathscr{D}(x_{i})=\{x_{j}-x_{i}\,:\,x_{j}\in\mathbb{X},\,a(\phi_{x_{j}},\phi_{x_{i}})\neq 0\}$ and defining a stencil function

[TABLE]

We have just reduced the computation of any $a(\phi_{x_{j}},\phi_{x_{i}})$ to the evaluation of scalar-valued functions enumerated by affine coordinates $(x_{i},x_{j}-x_{i})$. Indeed,

[TABLE]

In the present scenario, there may be a different set of translations $\mathscr{D}(x_{i})$ for every point $x_{i}$. However, if each point is drawn from a point lattice, most of the sets $\mathscr{D}(x_{i})$ are identical. This observation is the subject of the following subsection and a foundational principle in our approach.

Remark 3.1**.**

In the scenario that the bilinear form is symmetric, $a(u,v)=a(v,u)$, it is natural to assume that the stencil functions 6 will inherit a similar symmetry. Indeed, under the equivalent symmetry condition $G(x,u,v)=G(x,v,u)$ almost everywhere, one may easily verify that if $\delta=x_{j}-x_{i}$, then

[TABLE]

*or, equivalently, $\Phi^{\delta}_{i}(x_{i})=\Phi^{-\delta}_{j}(x_{i}+\delta)$. *

3.3 Local stencil functions and locally-structured meshes

An affine point lattice $\mathbb{L}$, from here on referred to only as a lattice, is a regularly spaced array of points in $\mathbb{R}^{n}$ where every point $x_{i}\in\mathbb{L}$ belongs to a neighborhood containing no other points in $\mathbb{L}$. In this paper, each (possibly finite) lattice is determined by a finite linearly independent set of translations in $\mathbb{R}^{n}$; i.e., $\mathbb{L}\operatorname{\subseteq}\{\delta_{0}+a_{1}\delta_{1}+\cdots+a_{l}\delta_{l}\,:\,a_{1},\ldots,a_{l}\in\mathbb{Z}\}$.

Assuming that the test function $\phi\in V$ is sufficiently localized and each point $x_{i}$ is drawn from a lattice $\mathbb{L}\operatorname{\subseteq}\Omega$, then each $\mathscr{D}(x_{i})$ is a subset of a small number of admissible translations $\mathscr{D}(\mathbb{L})=\bigcup\{\mathscr{D}(x_{i}):x_{i}\in\mathbb{L}\}$, determined solely by the lattice structure. In such a scenario, every stencil function is closely related; i.e., $\Phi_{i}^{\delta}(x)=\Phi_{j}^{\delta}(x)$, whenever both are defined. Therefore, it is prudent to drop the subscript and define only one common stencil function $\Phi^{\delta}(x)$ for each $\delta\in\mathscr{D}(\mathbb{L})$. Clearly,

[TABLE]

We are interested in exploiting 9 for solving a wide variety of PDEs with curvilinear geometries. Toward this end, the following examples help motivate our construction further.

Remark 3.2**.**

*From now on, it is useful to let $V$ be a closed subset of $W^{1,p}(\Omega)$, for some $1<p<\infty$. This will allow us to use a basis for $V_{h}\operatorname{\subseteq}V$ consisting of finite element vertex functions [23]. Other problems where V\operatorname{\subseteq}\big{[}W^{1,p}(\Omega)\big{]}^{k}, $k\in\mathbb{N}$, can be handled similarly (see, e.g., Section 4.2), by employing a basis consisting of the same vertex functions in each component. *

3.3.1 The one-dimensional setting

Let $V=H^{1}_{0}(\Omega)$, where $\Omega=(0,1)\operatorname{\subseteq}\mathbb{R}$, and fix a small translation $dx=1/(N+1)$. Consider the scenario where each point, ${x_{i}=x_{i-1}+dx}$, evenly divides $\Omega$ and $\phi$ is the piecewise-linear hat function defined $\phi(x)=\operatorname*{\vphantom{p}max}(1-\frac{|x|}{dx},0)$. Let $V_{h}=\{v\in H^{1}_{0}(\Omega):v|_{t}\in\mathcal{P}_{1}(t),\,t=(x_{i},x_{i+1}),\text{ for each }1\leq i\leq N\}$ and identify each shifted hat function with the standard basis, $\phi_{x_{i}}=\phi_{i}\in V_{h}$. Here, we may define $\mathbb{L}=\{x_{i}\}$. In this case, for each $i\geq 1$, the value $a(\phi_{i},\phi_{j})$ can either be computed directly from (3a), in the standard way, or evaluated using 9, assuming that each $\Phi^{\delta}$ is available at the onset of computation. Note that $\mathscr{D}(x_{1})=\{0,dx\}$ and $\mathscr{D}(x_{N})=\{-dx,0\}$, but $\mathscr{D}(x_{i})=\{-dx,0,dx\}$ for each $2\leq i\leq N-1$. Therefore, $\mathscr{D}(\mathbb{L})=\{-dx,0,dx\}$. Ultimately, defining each structured stencil function $\Phi^{\delta}(x)=\Phi((x_{i},\delta);x)$ from an arbitrary candidate point $x_{i}$, one may verify that

[TABLE]

which is clearly the same format as 9.

Recall Remark 3.1. If $a(\cdot,\cdot)$ is symmetric, then, by 8, $\Phi^{dx}(x_{i})=\Phi^{-dx}(x_{i}+dx)$ and the number of required stencil functions can be reduced to two. In some situations—e.g., when the zero row sum property can be employed (see Section 5)—only a single stencil function is actually required.

3.3.2 Locally-structured meshes with triangles

Let $m\in\mathbb{N}$. Beginning with a scaled Cartesian lattice $\mathbb{L}_{m}=2^{-m}\mathbb{Z}^{n}$, it is useful to define its intersection with the closure of the right-angled reference simplex $\hat{T}=\{\hat{x}\in\mathbb{R}^{n}:\|\hat{x}\|_{1}<1,\,\hat{x}\cdot e_{i}>0,\,\forall i=1,\ldots,n\}$. This simplicial lattice, $\hat{T}_{m}=\mathbb{L}_{m}\cap\overline{\hat{T}}$, can easily be transformed into a similar simplicial lattice $T_{m}$ for any arbitrary simplex $T\operatorname{\subseteq}\Omega$ via an affine transformation (see, e.g., Figure 1). Indeed, fixing the unique $A\in\mathbb{R}^{n\times n}$ and $b\in\mathbb{R}^{n}$ such that $T=\{A\hat{x}+b:\hat{x}\in\hat{T}\}$, the corresponding lattice is clearly $T_{m}=\{A\hat{x}_{i}+b:\hat{x}_{i}\in\hat{T}_{m}\}$. Note that there are no interior points, $\mathring{T}_{m}=T_{m}\cap T=\emptyset$, if $m<2$. We also define the set of boundary points $\partial T_{m}=T_{m}\setminus\mathring{T}_{m}$, which is always non-empty.

Let $\Omega\operatorname{\subseteq}\mathbb{R}^{n}$, where $n=2,3$. The utility of this transformation is evident upon considering a macro-triangulation of $\Omega$, say $\mathcal{T}_{H}$, where each macro-element $T\in\mathcal{T}_{H}$ is endowed with a simplicial lattice $T_{m}$, as defined above. Here, $H=\operatorname*{\vphantom{p}max}_{T\in\mathcal{T}_{H}}H_{T}$, where each $H_{T}=\mathrm{diam}(T)$ denotes the diameter of $T$. Notice that, for any fixed level $m\geq 0$, all interface points $x_{i}\in\partial T_{m}\cap\Omega$ are coincident with an interface lattice point on some neighbouring simplex. We may now define the set of all vertices on level $m$:

[TABLE]

For every $\mathbb{X}_{m}$, there is a corresponding finite element mesh such that each vertex function within a fixed macro-element $T\in\mathcal{T}_{H}$ is self-similar. In the cases $n=2$ or $3$, we are left to define each triangle or tetrahedron whose vertices coincide with points in $\mathbb{X}_{m}$. Let $[y_{0},y_{1},\ldots,y_{k}]\operatorname{\subseteq}\mathbb{R}^{n}$ denote the convex combination of the points $y_{0},y_{1},\ldots,y_{k}\in\Omega$. When $n=2$, the natural construction begins by considering the following uniform subdivision of a triangle $T=[y_{0},y_{1},y_{2}]\in\mathcal{T}_{H}$ into a set of four equal-volume triangles:

[TABLE]

For an illustration of this $n=2$ case, see Fig. 1. For $n=3$, see the construction in [14]. Further subdivisions can then be defined recursively, viz.,

[TABLE]

The set of all vertices in a given subdivision $\mathcal{S}^{m}(T)$ forms an evenly spaced set of points inside $T$. The set of vertices in $\mathcal{S}^{1}(T)$ clearly coincides with $T_{1}$ and one can easily verify from the recursive definition that the set of vertices in $\mathcal{S}^{m}(T)$ also coincides with $T_{m}$. We may finally define the sequence of locally-structured meshes:

[TABLE]

Notice that for each $m$, $\mathbb{X}_{m}$ coincides with the set of all vertices in $\mathcal{S}^{m}(\mathcal{T}_{H})$. Therefore, each point $x_{i}$ in $\mathbb{X}_{m}$ can be identified with a vertex function $\phi_{i}$ supported by only the neighboring elements appearing in $\mathcal{S}^{m}(\mathcal{T}_{H})$; see Fig. 1.

Assume that the fine mesh level, $m$, is chosen large enough that one may find several points in $\mathbb{X}_{m}$ which lie in the interior of some macro-element $T$. If $x_{i}\in\mathring{T}_{m}$ is any such point, then the translation set $\mathscr{D}(x_{i})=\mathscr{D}(\mathring{T}_{m})$ will only contain translations aligned with edges appearing in the original subdivision $\mathcal{S}^{1}(T)$. Identifying $\phi_{x_{i}}=\phi_{i}$ we see that for every $x_{i}\in\mathring{T}_{m}$ and $x_{j}\in T_{m}$,

[TABLE]

where $\Phi^{\delta}_{T}:\mathrm{conv}(\mathring{T}_{m})\to\mathbb{R}$ is a local stencil function for the current level $m$ and macro-element $T$ and $\mathrm{conv}(\mathring{T}_{m})$ is the convex hull of $\mathring{T}_{m}$. In general, notice that $\Phi^{\delta}_{T}\neq\Phi^{2\delta}_{T}$ are two different stencil functions, corresponding to the same direction but different mesh levels.

Consider the bilinear form 4 with a variable diffusion coefficient. A visualization of several corresponding local stencil functions, coming from locally-structured meshes used in our numerical experiments, is given in Figure 2. It is clear from this figure that each $\Phi^{\delta}_{T}$ has the potential to be a smooth function. We now come to the final essential component of our surrogate methodology; the approximation of $\Phi^{\delta}_{T}$.

Remark 3.3**.**

The lattice structure of locally-structured meshes is destroyed under smooth, non-affine transformations $\hat{T}\to\widetilde{T}\operatorname{\subseteq}\Omega$. This offers a possible impediment to our construction in the case of non-polygonal domains $\overline{\Omega}\neq\overline{\Omega}_{H}=\bigcup_{T\in\mathcal{T}_{H}}\overline{T}$. In fact, if a globally continuous transformation $\varphi:\overline{\Omega}_{H}\to\overline{\Omega}$ is available, with $\varphi|_{T}$ a smooth bijection for every $T\in\mathcal{T}_{H}$, then an equivalent method can be found using local pull-backs of $\varphi$. For example, the bilinear form in 4, $a_{1}:H^{1}(\Omega)\times H^{1}(\Omega)\to\mathbb{R}$, simply transforms to $a_{1,H}:H^{1}(\Omega_{H})\times H^{1}(\Omega_{H})\to\mathbb{R}$, where

[TABLE]

*Therefore, from now on, we operate under the assumption that the macro-mesh $\mathcal{T}_{H}$ is geometrically conforming, $\overline{\Omega}=\overline{\Omega}_{H}$. *

3.4 The inherited regularity of stencil functions

Our approach is to locally project each stencil function $\Phi^{\delta}_{T}$ in 14 onto a high-dimensional space of polynomials and later use this projection to compute approximate values of the stiffness matrix $\mathsf{A}$. Let $\mathcal{T}_{H}$ be a shape-regular simplicial mesh and consider a locally-structured mesh of level $m$ subordinate to this mesh, $\mathcal{S}^{m}(\mathcal{T}_{H})$. Before moving on, observe that definition 14 in fact holds for any $x_{i},x_{j}\in T_{m}$ if $x_{i}$ or $x_{j}\in\mathring{T}_{m}$. Therefore, due to the structure of the vertex functions in a locally-structured mesh, the domain of each $\Phi^{\delta}_{T}$ can actually be extended to a set $\overline{T}_{\delta}$ lying between $\mathrm{conv}(\mathring{T}_{m})$ and $\overline{T}$. Indeed, identifying the test function in 6 with the $i$-th vertex function, $\phi(x)=\phi_{i}(x+x_{i})$, for an arbitrary vertex $x_{i}\in\mathring{T}_{m}$, define

[TABLE]

From now on, we assume $\Phi^{\delta}_{T}:\overline{T}_{\delta}\to\mathbb{R}$. See Fig. 3 for a depiction of the sets in a triangular mesh and note that $T_{-\delta}=T_{\delta}+\delta$, for every $\delta\in\mathscr{D}(T_{m})$.

For any $T\in\mathcal{T}_{H}$, let $\mathcal{P}_{q}(T_{\delta})$ denote the space of polynomials of degree at most $q$ on the simplex $\overline{T}_{\delta}$ and let $\Pi^{\delta}_{T}:C^{0}(\overline{T}_{\delta})\to\mathcal{P}_{q}(T_{\delta})$ be an $L^{\infty}$-continuous projection operator, $\Pi^{\delta}_{T}\circ\Pi^{\delta}_{T}=\Pi^{\delta}_{T}$. For each macro-element $T\in\mathcal{T}_{H}$ and level $m\in\mathbb{N}$, define the surrogate stencil function ${\widetilde{\Phi}^{\delta}_{T}\!:{T}_{\delta}\to\mathbb{R}}$ to be the corresponding polynomial projection of $\Phi^{\delta}_{T}$. Namely, {equation+} ~Φ^δ_T

Π^δ_T  Φ^δ_T .

In order to correctly argue that a polynomial approximation of $\Phi^{\delta}_{T}$ is feasible, it is necessary to classify its regularity depending on the problem at hand. In the following proposition, we show, under the modest assumptions above, that if $G(\cdot|_{T},\cdot,\cdot)$ is a polynomial in its first argument, then $\Phi^{\delta}_{T}$ is also a polynomial of the same degree.

Lemma 3.1**.**

*Fix a simplex $T\in\mathcal{T}_{H}$. Assume that the bilinear form $a(\cdot,\cdot)$ in 1a satisfies assumptions 3a and 3b, where the integrand $G(\cdot|_{T},\cdot,\cdot)$ is a polynomial of at most degree $q$ in its first argument. Then, for any locally-structured mesh $\mathcal{S}^{m}(\mathcal{T}_{H})$, as defined above, every local stencil function $\Phi^{\delta}_{T}\in\mathcal{P}_{q}(T_{\delta})$ is a polynomial of the same degree. *

Proof 3.2**.**

Recall definition 6. For every $x_{i},x_{j}\in\mathbb{X}_{m}$, every stencil function $\Phi((x_{i},x_{j}-x_{i}),x)$ is defined with a test function $\phi\in V$. Fixing any arbitrary vertex $x_{i}\in\mathring{T}_{m}$, identify this test function with the $i$-th vertex function, $\phi(x)=\phi_{i}(x+x_{i})$. Let $\alpha=(\alpha_{1},\ldots,\alpha_{n})\in\mathbb{N}_{0}$ denote a standard multi-index, $|\alpha|=\sum_{i=1}^{n}|\alpha_{i}|$ and $x^{\alpha}=x_{\scriptscriptstyle 1}^{\alpha_{1}}\cdots x_{\scriptscriptstyle n}^{\alpha_{n}}$. By assumption, we may express $G(x,\phi_{\delta}(y),\phi(y))=\sum_{|\alpha|\leq l}c_{\alpha}(y)x^{\alpha}$, where each coefficient function $c(y)$ has support only in $\Omega_{\delta}=\operatorname*{supp}(\phi)\cap\operatorname*{supp}(\phi_{\delta})$. Moreover, if $x+y\in T$, then

[TABLE]

*is clearly an equal degree polynomial in the variable $x$. The proof is completed by noting that the integral in 6 is performed only over the variable $y\in\Omega_{\delta}$ and so the stencil function acts like a convolution. Indeed, under the assumption $x\in T$, only the subset of points $x\in T_{\delta}=\{x\in T:x+y\in T\text{ for all }y\in\Omega_{\delta}\}$ guarantee that 17 holds at every point of integration $y$. In this case, by linearity of integration, $\Phi^{\delta}_{T}$ is a member of $\mathcal{P}_{q}(T_{\delta})$, with its coefficients defined by the associated integrals of the $y$-dependent functions in the right-hand side of 17. *

Corollary 3.3**.**

*In the setting of Lemma 3.1, $\widetilde{\Phi}^{\delta}_{T}=\Phi^{\delta}_{T}$. *

Proof 3.4**.**

*Since $\Phi^{\delta}_{T}\in\mathcal{P}_{q}(T_{\delta})$, we immediately see that $\Pi^{\delta}_{T}\,\Phi^{\delta}_{T}=\Phi^{\delta}_{T}$. *

3.5 Surrogate stiffness matrices

The main goal of the entire effort above is to guide us in reducing the vast majority of the finite element assembly process to the evaluation of a small set of functions which can, in fact, be locally approximated by polynomials. As in Section 3.3.2, let $\mathcal{T}_{H}$ be a shape-regular triangulation of a domain $\Omega$ into disjoint macro-simplices $T$.

Our construction of the surrogate matrix $\widetilde{\mathsf{A}}$ is obviously built to exploit the local lattice structure of locally-structured meshes. As argued previously, a surrogate stencil function $\widetilde{\Phi}^{\delta}_{T}$ can be used to approximate any matrix entry $\mathsf{A}_{ij}$ coming from a locally-structured mesh if at least one of the corresponding vertices $x_{i}$ or $x_{j}$ belongs to $\mathring{T}_{m}$. This leaves us to define only the nonzero matrix entries coming from the mutual interaction of vertex functions at the boundaries of the macro-elements. Although these entries could also be approximated by surrogate stencil functions — in this case, these additional functions would be defined on each subsimplex of the macro-mesh $\mathcal{T}_{H}$ — let us assume that they are computed directly. Because the growth of the macro-mesh boundary and interface interactions grow at an order of magnitude less than the interior interactions, computing these matrix entries directly does not affect to the asymptotic performance of the methodology. Finally, letting $\partial\mathbb{X}_{m}=\bigcup_{T\in\mathcal{T}_{H}}\partial T_{m}$ denote the union of all macro-mesh boundary vertices, we define the general surrogate stiffness matrix

[TABLE]

Remark 3.5**.**

Due to the presence of the surrogate stencil functions $\widetilde{\Phi}^{\delta}_{T}$, even if $\mathsf{A}$ is symmetric, $\widetilde{\mathsf{A}}$ will generally not be. However, recalling 8, observe that if $a(\cdot,\cdot)$ is symmetric, then $\Phi^{\delta}_{T}(x)=\Phi^{-\delta}_{T}(x+\delta)$. Therefore, if we use related projection operators

[TABLE]

for each opposing direction $\delta$ and $-\delta$, then $\widetilde{\mathsf{A}}$ will be symmetric. Indeed, if $\delta=x_{j}-x_{i}$, then

[TABLE]

4 Examples

In this section, we present three example problems which easily fit into the framework above.

4.1 The variable coefficient Poisson equation

Consider the Poisson-type equation ${-\mathrm{div}(K\nabla u)=f}$ in $\Omega$, $u=0$ on $\partial\Omega$, with a load $f\in L^{2}(\Omega)$ and a variable, symmetric positive-definite tensor $K$. Furthermore, assume that for each index $a,b$, $K_{ab}\in\mathcal{P}_{q}(\Omega)$. Recall 4 and note that we have already shown that the weak form of this problem can be cast into the framework above.

Assume that, within some set $T\operatorname{\subseteq}\Omega$, each vertex function $\phi_{i}$ is a translation of a fixed test function $\phi(x)=\phi_{i}(x-x_{i})$. Then for each $\phi_{i},\phi_{j}$, the stiffness matrix entry

[TABLE]

can equally well be expressed as the evaluation (at the point $x_{i}$) of a stencil function, which, by Lemma 3.1, is simply a polynomial of the same degree as the diffusion tensor $K$. In the case of locally-structured meshes, there is a locally defined stencil function $\Phi^{\delta}_{T}:\Omega\to\mathbb{R}$ for each macro-element $T$ and level $m$. In this case, each $\Phi^{\delta}_{T}:\overline{T}_{\delta}\to\mathbb{R}$ is a polynomial (of degree at most $q$) on $T_{\delta}$.

4.2 Linearized elasticity

Let $\vec{\nabla}$ and $\mathrm{Div}$ denote the row-wise distributional gradient and divergence, respectively. Now define $\epsilon(u)=\frac{1}{2}[\vec{\nabla}u+(\vec{\nabla}u)^{\top}]$ to be the symmetric gradient operator $\epsilon:H^{1}(\Omega)^{n}\to L^{2}(\Omega)^{n}$, where $n\geq 2$. Consider the following standard PDE model for the displacement $u\in H^{1}_{0}(\Omega)^{n}$ of a linearly elastic isotropic material: $-\mathrm{Div}\,\sigma=\vec{f}$, where the stress $\sigma=2\mu\epsilon(u)+\lambda I\mathrm{div}\,u$ and the load $f\in L^{2}(\Omega)^{n}$.

The weak form of this equation is well known in the literature [20] and the associated bilinear form is simply

[TABLE]

This bilinear form obviously satisfies assumptions 3a and 3b. If we assume that the Lamé parameters $\mu,\lambda:\Omega\to\mathbb{R}$ are piecewise polynomials on a collection of disjoint subdomains $T\in\mathcal{T}_{H}$, then each associated stencil function is also a piecewise polynomial.

4.3 $p$-Laplacian diffusion

For any $1<p<\infty$, let $\Delta_{p}\,u=\operatorname{div}(|\nabla u|^{p-2}\nabla u)$ be the $p$-Laplacian operator. Fix a valid parameter $p$ and consider the nonlinear diffusion equation $\frac{\partial u}{\partial t}-\Delta_{p}u=f$, where $f\in L^{p}(\Omega)$. A simple Euler time-stepping scheme replaces the time derivative $\frac{\partial u}{\partial t}$ by the quotient $\frac{u_{k+1}-u_{k}}{dt}$, where $dt>0$ is a fixed time-step parameter. Choosing backward Euler time-stepping and defining $f_{k}=f(k\cdot dt)$, we arrive at a semi-discrete nonlinear elliptic PDE for the solution variable $u_{k}$, which must be solved at each step $k\in\mathbb{N}$: $u_{k}-dt\hskip 0.5pt\Delta_{p}u_{k}=dt\hskip 0.5ptf_{k}+u_{k-1}$. Upon fixed point linearization of the weak form of this equation, we uncover the following bilinear form:

[TABLE]

Here, the variable coefficient $\widetilde{u}\in W^{1,p}(\Omega)$ is usually identified with the previous solution iteration in the associated fixed point algorithm (cf. Section 9.3).

The bilinear form $b(\cdot,\cdot)$ can easily be placed into the form of 3a and each matrix entry can therefore be superceded by stencil function evaluations, $\Phi^{\delta}_{T}(x)$, as in 6. Alternatively, one may split $b(\cdot,\cdot)$ into a mass term $m(\cdot,\cdot)$ and a $dt$-weighted stiffness term $a_{3}(\cdot,\cdot)$. Specifically, $b(u,v)=m(u,v)+dt\cdot a_{3}(u,v)$, where

[TABLE]

With this observation in hand, we see that $b(u,v)$ may be discretized by a linear combination of independent surrogates; one for $m(\cdot,\cdot)$ and one for $a_{3}(\cdot,\cdot)$ (cf. Section 9.3). In either approach, the variable coefficient $|\nabla\widetilde{u}(x)|^{p-2}$ will generally not remain a polynomial in a subdomain of $\Omega$ and the accuracy of a surrogate stencil function $\widetilde{\Phi}^{\delta}_{T}$ will reflect the local regularity of the solution from the previous iteration, $\widetilde{u}$.

5 Boundary conditions and the zero row sum property

It is generally appropriate to define the surrogate stiffness matrix component-wise by the rule given in 18. Nevertheless, in some problems the operator to be discretized has a kernel which is not guaranteed to be respected by the surrogate. In such scenarios, it is possible that better performance and accuracy can be achieved if elements of this kernel are incorporated into the construction of the surrogate. This occurrence is most easily illustrated with the Poisson example from Section 4.1.

Consider the bilinear form $a_{1}:H^{1}(\Omega)\times H^{1}(\Omega)\to\mathbb{R}$, defined in 4. Define $V_{h}^{\mathrm{ext}}=\{v\in H^{1}(\Omega):v|_{t}\in\mathcal{P}_{1}(t)\text{ for each }t\in\mathcal{S}^{m}(\mathcal{T}_{H})\}$ and $V_{h}=\{v\in H^{1}_{0}(\Omega):v|_{t}\in\mathcal{P}_{1}(t)\text{ for each }t\in\mathcal{S}^{m}(\mathcal{T}_{H})\}\operatorname{\subseteq}V_{h}^{\mathrm{ext}}$. Let the corresponding vertex function bases be $\{\phi_{i}\}\operatorname{\subseteq}\{\phi_{i}^{\mathrm{ext}}\}$, with $\phi_{i}=\phi^{\mathrm{ext}}_{i}$ for $1\leq i\leq N$. In most finite element software, a space like $V_{h}^{\mathrm{ext}}$ is used to impose Dirichlet boundary conditions. Indeed, a “lift” of the Dirichlet data, say $u_{h}^{\mathrm{ext}}=\sum_{i}\mathsf{u}^{\mathrm{ext}}_{i}\phi_{i}^{\mathrm{ext}}$, is generally constructed from a linear combination of the set $\{\phi_{i}^{\mathrm{ext}}\}\setminus\{\phi_{i}\}$. Then, taking $\mathsf{A}^{\mathrm{ext}}_{ij}=a(\phi_{j}^{\mathrm{ext}},\phi_{i})$ for each valid $i,j$, a modified load vector $\mathsf{f}^{\mathrm{ext}}=\mathsf{f}-\mathsf{A}^{\mathrm{ext}}\mathsf{u}^{\mathrm{ext}}$ is used in computation.

It is obvious that $\nabla\,1=0$ and so $a_{1}(1,v)=0$ for any $v\in H^{1}\Omega)$. Therefore, by the partition of unity property $\sum_{i}\phi_{i}^{\mathrm{ext}}=1$, the zero row sum of the matrix $\mathsf{A}^{\mathrm{ext}}$ also vanishes. Namely, $\sum_{j}\mathsf{A}_{ij}^{\mathrm{ext}}=0$. This property may be induced in the corresponding surrogate matrix if we simply define $\widetilde{\mathsf{A}}_{ii}^{\mathrm{ext}}=-\sum_{j\neq i}\widetilde{\mathsf{A}}_{ij}^{\mathrm{ext}}$, for every $x_{i}\in\mathbb{X}_{m}$, where $\widetilde{\mathsf{A}}_{ij}^{\mathrm{ext}}=\mathsf{A}_{ij}^{\mathrm{ext}}$ for every $j$ where $\mathsf{A}_{ij}$ is not defined, and $\widetilde{\mathsf{A}}_{ij}^{\mathrm{ext}}=\widetilde{\mathsf{A}}_{ij}$ otherwise. With this extra condition, the surrogate matrix 18 actually requires one fewer independent stencil function; i.e., $\Phi^{0}_{T}=-\sum_{\delta\in\mathscr{D}(T_{m})\setminus\{0\}}\Phi^{\delta}_{T}$, for every $T\in\mathcal{T}_{H}$. By this definition, although $\widetilde{\mathsf{A}}$ does not satisfy the zero row sum property, the matrix $\mathsf{A}-\widetilde{\mathsf{A}}$ does. Indeed,

[TABLE]

In this way, the stiffness matrix coming from linearized elasticity 22 is similar to the stiffness matrix coming from the Laplacian. Indeed, the zero row sum property can be incorporated into its surrogate via a straight-forward generalization.

6 Analyzing the surrogate discretization

In this section, we define and motivate what we see as some the most essential features in the analysis of our surrogate methods. We begin with a review of discrete stability in the context of 1c. We then touch on the concept of spectral convergence of the surrogate matrix $\widetilde{\mathsf{A}}\to\mathsf{A}$, which helps us motivate the need to control $\|\mathsf{A}-\widetilde{\mathsf{A}}\|_{\operatorname*{\vphantom{p}max}}$. This specific quantity will repeatedly appear in the a priori error analysis in Section 7.

6.1 Discrete stability

Let $S=\{v\in V\,:\,\|v\|_{V}=1\}$ be the surface of the unit ball in $V$. Recall 2 and assume that the discretization $\mathsf{A}\mathsf{u}=\mathsf{f}$ is stable. In the present context, this is equivalent to the existence of a constant $\alpha>0$ such that

[TABLE]

Likewise, in order for the surrogate discretization $\widetilde{\mathsf{A}}\widetilde{\mathsf{u}}=\mathsf{f}$ to be stable, we must show that there exists a constant $\widetilde{\alpha}>0$ such that

[TABLE]

Inequality 27 guarantees that $\widetilde{\mathsf{A}}\widetilde{\mathsf{u}}=\mathsf{f}$ has a unique solution and that $\|\widetilde{u}_{h}\|_{V}\leq\widetilde{\alpha}^{-1}\|F\|_{V^{\ast}}$. Equally important, however, it is a necessary precursor to Strang’s First Lemma, which in some cases can be used, in part, to show that $\widetilde{u}_{h}$ converges to the exact solution $u$ (see, e.g., Section 7.2).

6.2 Spectral convergence

By analyzing the singular values of $\mathsf{A}$ directly, 27 can sometimes be proven by showing that the spectrum of $\widetilde{\mathsf{A}}$ converges to the spectrum of $\mathsf{A}$ at a fast enough rate. The main takeaway from this section is that spectral convergence can be guaranteed by showing that $\widetilde{\mathsf{A}}\to\mathsf{A}$ in the matrix maximum norm, $\|\cdot\|_{\operatorname*{\vphantom{p}max}}$. Before moving on, denote the $k$-smallest eigenvalue of a matrix $\mathsf{M}\in\mathbb{R}^{N\times N}$ as $\lambda_{k}(\mathsf{M})$ and let $\ell(\mathsf{M})=\operatorname*{\vphantom{p}max}_{1\leq i\leq N}\,\#\{\mathsf{M}_{ij}\neq 0\text{ where }1\leq j\leq N\}$ be the maximum number of nonzero components in $\mathsf{M}$, taken across all individual rows.

Proposition 6.1**.**

Let $\mathsf{M},\mathsf{N}\in\mathbb{R}^{N\times N}$ be symmetric matrices. Then, for each $k=1,\ldots,N$, it holds that

[TABLE]

The proof is placed in Appendix A. Because $\mathsf{A}$ and $\widetilde{\mathsf{A}}$ have the same sparsity pattern, the next result follows trivially. Note that when $n=2$, $\ell(\mathsf{A}-\widetilde{\mathsf{A}})\leq 7$, and when $n=3$, $\ell(\mathsf{A}-\widetilde{\mathsf{A}})\leq 15$, whereas $\ell(\mathsf{A})$ can be larger depending on the structure of the macro-mesh.

Corollary 6.2**.**

Let $\mathsf{A}$ and $\widetilde{\mathsf{A}}$ be the true and surrogate stiffness matrices in 2, respectively. If both $\mathsf{A}$ and $\widetilde{\mathsf{A}}$ are real symmetric matrices, then

[TABLE]

Remark 6.3**.**

*As stated above, if $\|\mathsf{A}-\widetilde{\mathsf{A}}\|_{\operatorname*{\vphantom{p}max}}\to 0$ fast enough, then Corollary 6.2 can be used in proving the stability condition 27. However, 29 is generally a very pessimistic bound and, when available, we recommend using more direct means to prove discrete stability (see, e.g., Theorem 7.1). Nevertheless, this result illustrates the importance of controlling $\|\mathsf{A}-\widetilde{\mathsf{A}}\|_{\operatorname*{\vphantom{p}max}}$, which is a central feature in all of the coming analysis. *

6.3 Controlling $\|\mathsf{A}-\widetilde{\mathsf{A}}\|_{\operatorname*{\vphantom{p}max}}$ with the variable coefficient Poisson equation

Before we begin, some new notation is required. For any tensor $K:\Omega\to\mathbb{R}^{n\times n}$, define $\|K\|_{L^{\infty}(\Omega)}=\operatorname*{\vphantom{p}max}_{a,b}\|K_{ab}\|_{L^{\infty}(\Omega)}$, likewise, for any $r\geq 0$, define $|K|_{W^{r+1,\infty}(T)}=\operatorname*{\vphantom{p}max}_{a,b}|K_{ab}|_{W^{r+1,\infty}(T)}$. From now on, the notation $A\lesssim B$ will be used when two mesh-dependent quantities $A$ and $B$ satisfy an inequality $A\leq CB$, where $C$ some positive $H$-independent constant. Likewise, when $A\lesssim B$ and $B\lesssim A$, we write $A\eqsim B$. Recall that for a macro-mesh $\mathcal{T}_{H}$, the diameter of a single element $T\in\mathcal{T}$ is denoted $H_{T}$ and the mesh size is denoted $H=\operatorname*{\vphantom{p}max}_{T\in\mathcal{T}_{H}}H_{T}$. We also denote the fine-scale element diameter $h_{T}=2^{-m}H_{T}$, for each $T\in\mathcal{T}_{H}$ and $h=2^{-m}H$.

Lemma 6.4**.**

Let $\mathsf{A}$ and $\widetilde{\mathsf{A}}$, respectively, be the true and surrogate stiffness matrices corresponding to the bilinear form 4. Namely, let each component of $\mathsf{A}$ be given by 21 and each component of $\widetilde{\mathsf{A}}$ be defined by 18 with $G(x,u(y),v(y)):={\nabla u(y)^{\top}K(x)\nabla v(y)}$. Fix $T\in\mathcal{T}_{H}$ and $0\leq r\leq q$. If $K_{ab}|_{T}\in W^{r+1,\infty}(T)$ for each index $a,b$, then

[TABLE]

Proof 6.5**.**

We prove only 30a, 30b then follows immediately. Recall 18 and fix $T\in\mathcal{T}_{H}$. Let $i$ and $j$ be the indices of the maximal value |\mathsf{A}_{ij}-\widetilde{\mathsf{A}}_{ij}|=\big{\|}\mathsf{A}-\widetilde{\mathsf{A}}\big{\|}_{\operatorname*{\vphantom{p}max},T_{m}}. Next, because the theorem trivially holds in the degenerate case \big{\|}\mathsf{A}-\widetilde{\mathsf{A}}\big{\|}_{\operatorname*{\vphantom{p}max},T_{m}}=0, we proceed under the assumption that $\mathsf{A}_{ij}\neq\widetilde{\mathsf{A}}_{ij}$. Notably, it follows from Corollary 3.3 that if each $K_{ab}|_{T}\in\mathcal{P}_{q}(T)$, then $\widetilde{\mathsf{A}}_{ij}=\mathsf{A}_{ij}$ and, therefore, we find ourselves in the scenario where the diffusion tensor $K|_{T}$ is not a polynomial (of degree at most $q$). Here, we may also freely assume that $i\neq j$ because of 25. Indeed, for each $i$, $|\mathsf{A}_{ii}-\widetilde{\mathsf{A}}_{ii}|\leq\sum_{j\neq i}|\mathsf{A}_{ij}-\widetilde{\mathsf{A}}_{ij}|\leq\ell(\mathsf{A}-\widetilde{\mathsf{A}})\cdot\operatorname*{\vphantom{p}max}_{j\neq i}{|\mathsf{A}_{ij}-\widetilde{\mathsf{A}}_{ij}|}$.

To fix notation in the remainder of the proof, we take $\phi=\phi_{i}$, and express

[TABLE]

for some nonzero $\delta\in\mathscr{D}(T_{m})$. Let $\mathcal{I}_{T}:C^{0}(\overline{T})\to\mathcal{P}_{r}(T)$ be the local Lagrange interpolant and define \big{[}\mathcal{I}_{T}^{n\times n}K\big{]}_{ab}=\mathcal{I}_{T}K_{ab}, for each index $a,b$. Splitting the two matrix entries $\mathsf{A}_{ij}$ and $\widetilde{\mathsf{A}}_{ij}$ into polynomial and non-polynomial parts and rewriting

[TABLE]

we find that

[TABLE]

Upon canceling the first two terms in the expressions above, we arrive at the inequality

[TABLE]

where \beta_{ij}(x)=\int_{\Omega_{\delta}}\nabla\phi_{\delta}(y)^{\top}\big{[}K-\mathcal{I}_{T}^{n\times n}K\big{]}(x+y)\hskip 0.25pt\nabla\phi(y)\hskip 0.5pt\,\mathrm{d}y. Recall that the projection $\Pi^{\delta}_{T}:C^{0}(\overline{T}_{\delta})\to\mathcal{P}(T_{\delta})$ is continuous in the $L^{\infty}(\Omega)$ norm. Therefore,

[TABLE]

A standard scaling argument shows that $\|\nabla\phi_{\delta}\cdot\nabla\phi\|_{L^{1}(\Omega_{\delta})}\lesssim h_{T}^{n-2}$. This, together with the well-known property $\|K_{ab}-\mathcal{I}_{T}K_{ab}\|_{L^{\infty}(T)}\lesssim H^{r+1}_{T}|K_{ab}|_{W^{r+1,\infty}(T)}$, yields the sufficient result

[TABLE]

Remark 6.6**.**

The proof of Lemma 6.4 can be read as a blueprint which extends to the settings of the bilinear forms $a_{2}(\cdot,\cdot)$, $a_{3}(\cdot,\cdot)$, defined in 22 and 24. Indeed, when $a_{2}(\cdot,\cdot)$ is considered, the only significant modification to the proof above is that an interpolation operator $\mathcal{I}_{T}:C^{0}(\overline{T})\to\mathcal{P}_{r}(T)$ must be introduced for each Lamé parameter $\mu$ and $\lambda$. The adaption to the setting $a(\cdot,\cdot)=a_{3}(\cdot,\cdot)$, is obvious. Ultimately,

[TABLE]

*Moreover, the proof may be easily modified to permit surrogates $\widetilde{\mathsf{A}}$ without the zero row sum property. Likewise, scenarios involving fewer derivatives (which generally do not possess the zero row sum property), e.g., $a(\cdot,\cdot)=m(\cdot,\cdot)$, have similar bounds but invoke a different scaling in $h$. *

7 A priori error estimation for the variable coefficient Poisson equation

In this section, we present a thorough analysis of a surrogate discretization of the variable coefficient Poisson equation. Given a load $f\in L^{2}(\Omega)$ and symmetric positive-definite tensor $K:\Omega\to\mathbb{R}^{n\times n}$, the corresponding weak form may be written as follows.

[TABLE]

where $a(u,v)=\int_{\Omega}\nabla u^{\top}K\hskip 0.25pt\nabla v\,\mathrm{d}x$ and $F(v)=\int_{\Omega}f\hskip 0.25ptv\,\mathrm{d}x$. As done in Section 5, define $V_{h}=\{v\in H^{1}_{0}(\Omega):v|_{t}\in\mathcal{P}_{1}(t)\text{ for each }t\in\mathcal{S}^{m}(\mathcal{T}_{H})\}$ and let the corresponding vertex function basis be $\{\phi_{i}\}$.

7.1 Coercivity

In the present setting, observe that each $v_{h}\in V_{h}$ can be expressed as $v_{h}(x)=\sum_{i}v_{h}(x_{i})\phi_{i}(x)$. Therefore, due to the zero row/column sum property 25, we find

[TABLE]

Due to the mutual sparsity of the matrices $\widetilde{\mathsf{A}}$ and $\mathsf{A}$, every nonzero term in the sum above can be rewritten as \big{(}\mathsf{A}_{ij}-\widetilde{\mathsf{A}}_{ij}\big{)}(v_{h}(x_{i})-v_{h}(x_{i}+\delta))\hskip 0.25pt(w_{h}(x_{i})-w_{h}(x_{i}+\delta)), for some nonzero $\delta$. Because $|\delta|\eqsim h$ by construction, one easily arrives at the following upper bound:

[TABLE]

We now arrive at the main result of this subsection.

Theorem 7.1**.**

*Let $0\leq r\leq q$. Assume that $a(\cdot,\cdot)$ is coercive and K\in\big{[}W^{r+1,\infty}(\mathcal{T}_{H})\big{]}^{n\times n}. Then, for any fine enough macro-mesh $\mathcal{T}_{H}$, the surrogate bilinear form $\widetilde{a}:V_{h}\times V_{h}\to\mathbb{R}$ is also coercive. *

Proof 7.2**.**

Let $S=\{v\in H^{1}\,:\,\|v\|_{1}=1\}$ be the surface of the unit ball in $H^{1}$. Recall that since $a(\cdot,\cdot)$ is coercive, there exists a coercivity constant $\alpha>0$ such that $a(v,v)\geq\alpha$ for all $v\in S$. Notice that $\alpha\leq a(v_{h},v_{h})\leq\widetilde{a}(v_{h},v_{h})+|a(v_{h},v_{h})-\widetilde{a}(v_{h},v_{h})|$ for all $v_{h}\in V_{h}\cap S$ and, therefore,

[TABLE]

Here, the second term on the left may be bounded from above using 41 and 6.4 as follows,

[TABLE]

*Thus, for any small enough $H$, we see that $0<\alpha-|a(v_{h},v_{h})-\widetilde{a}(v_{h},v_{h})|\leq\widetilde{a}(v_{h},v_{h})$, as necessary. *

7.2 Convergence of the surrogate solution in the $H^{1}$ norm

The purpose of this subsection is to derive a mesh-dependent upper bound on the error in the surrogate solution $\widetilde{u}_{h}$ of the form $\|u-{\widetilde{u}}_{h}\|_{1}\leq C(K,\Omega,u)\hskip 0.25pth+\widetilde{C}(K,\Omega,u)H^{r+1}$. In doing so, we choose to emphasize the primary difference from the classical $\|u-u_{h}\|_{1}\leq C(K,\Omega,u)\hskip 0.25pth$ error estimate by absorbing the coercivity and continuity constants (which depend on both $K$ and $\Omega$) into the $\lesssim$ symbol. We begin with a particular version of the First Strang Lemma [37].

Lemma 7.3**.**

Let $S=\{v\in H^{1}\,:\,\|v\|_{1}=1\}$ be the surface of the unit ball in $H^{1}$. Assume that $\widetilde{a}:V_{h}\times V_{h}\to\mathbb{R}$ is coercive. The following error estimate holds for the surrogate solution $\widetilde{u}_{h}$ of the variable coefficient model problem 39:

[TABLE]

Theorem 7.4**.**

Let $0\leq r\leq q$ and assume that K\in\big{[}W^{r+1,\infty}(\Omega)\big{]}^{n\times n} is symmetric and positive definite with $\lambda_{1}(K)$ bounded away from zero almost everywhere. Let $u\in H^{1}(\Omega)$ and $\widetilde{u}_{h}\in V_{h}$ be the unique solutions to 1a and 1c, respectively, where $a(u,v)=\int_{\Omega}\nabla u^{\top}K\hskip 0.25pt\nabla v\,\mathrm{d}x$ and $F(v)=\int_{\Omega}f\hskip 0.25ptv\,\mathrm{d}x$. Then, for any sufficiently fine macro-mesh $\mathcal{T}_{H}$, the following upper bound holds:

[TABLE]

Proof 7.5**.**

With the assumptions above, $a(\cdot,\cdot)$ is coercive. Therefore, by Theorem 7.1, if the macro-mesh $\mathcal{T}_{H}$ is taken fine enough, then $\widetilde{a}:V_{h}\times V_{h}\to\mathbb{R}$ is coercive. We now bound the right-hand side of 44. Invoking 41, we find

[TABLE]

for every $w_{h}\in V_{h}$. Setting $w_{h}=\mathcal{S}\mathcal{Z}_{h}u$, the Scott–Zhang interpolant of $u$ [36], we see that

[TABLE]

*In order to finish the proof, recall that $h^{2-n}\|\mathsf{A}-\widetilde{\mathsf{A}}\|_{\operatorname*{\vphantom{p}max}}\lesssim H^{r+1}|K|_{W^{r+1,\infty}(\Omega)}$, by Lemma 6.4.  *

7.3 Convergence of the surrogate solution in the $L^{2}$ norm

In this subsection, we prove an $L^{2}$ error estimate of the form $\|u-{\widetilde{u}}_{h}\|_{0}\leq C(K,\Omega,u)\hskip 0.25pth^{2}+\widetilde{C}(K,\Omega,u)H^{r+1}$. A second result, which elicits accelerated $H$-convergence, is also proved under the additional assumption \sum_{x_{i},\in T^{\delta}_{m}}\big{[}\Phi^{\delta}_{T}-\Pi^{\delta}_{T}\Phi^{\delta}_{T}\big{]}(x_{i})=0, where $T^{\delta}_{m}=T_{m}\cap\overline{T}_{\delta}$. This is a property which naturally arises for the specific class of least-squares projections introduced in Section 8.1. Again, we emphasize the primary differences from the corresponding classical error estimate by absorbing the standard constants into the $\lesssim$ symbol.

Theorem 7.6**.**

Under the conditions of Theorem 7.4, if $\Omega\operatorname{\subseteq}\mathbb{R}^{n}$ is a convex domain, then the following additional upper bound on the error in the surrogate solution holds:

[TABLE]

Proof 7.7**.**

By the triangle inequality, $\|u-\widetilde{u}_{h}\|_{0}\leq\|u-u_{h}\|_{0}+\|u_{h}-\widetilde{u}_{h}\|_{0}$, where $u_{h}\in V_{h}$ is the discrete solution coming from 1b. It can be shown that if $\Omega$ is convex, then $u\in H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$; see, e.g., [27]. It then follows from standard arguments that $\|u-u_{h}\|_{0}\lesssim h^{2}\hskip 0.25pt|u|_{2}$; see, e.g., [16, Theorem 5.7.6]. Therefore, we only need to analyze the term $\|u_{h}-\widetilde{u}_{h}\|_{0}$. Since, $\|u_{h}-\widetilde{u}_{h}\|_{0}\leq\|u_{h}-\widetilde{u}_{h}\|_{1}$, proceeding as in the proof of Theorem 7.4, we quickly arrive at 48a.

In order to prove 48b, first define $w_{h}\in V_{h}$ satisfying $a(w_{h},v_{h})=(u_{h}-\widetilde{u}_{h},v_{h})_{\Omega}$, for all $v_{h}\in V_{h}$. Observe that the exact solution of the problem $a(w,v)=(u_{h}-\widetilde{u}_{h},v)_{\Omega}$, for all $v\in H^{1}_{0}(\Omega)$, belongs to the space $H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$, $\|w\|_{2}\lesssim\|u_{h}-\widetilde{u}_{h}\|_{0}$. Moreover,

[TABLE]

where the final line follows from 40. As remarked previously, each nonzero term in this sum can be written as \big{(}\mathsf{A}_{ij}-\widetilde{\mathsf{A}}_{ij}\big{)}(\widetilde{u}_{h}(x_{i})-\widetilde{u}_{h}(x_{i}+\delta))\hskip 0.25pt(w_{h}(x_{i})-w_{h}(x_{i}+\delta)), for some $T\in\mathcal{T}_{H}$ and nonzero $\delta\in\mathscr{D}(T_{m})$. We can make better use of this expression with the identity $v_{h}(x_{i})-v_{h}(x_{i}+\delta)=\nabla v_{h}(y_{i,\delta})$, wherein each $y_{i,\delta}$ is chosen from the edge connecting $x_{i}$ and $x_{i}+\delta$, and with the relationship \mathsf{A}_{ij}-\widetilde{\mathsf{A}}_{ij}=\big{[}\Phi^{\delta}_{T}-\Pi^{\delta}_{T}\Phi^{\delta}_{T}\big{]}(x_{i}), since $\mathsf{A}_{ij}-\widetilde{\mathsf{A}}_{ij}\neq 0$ and $i\neq j$. With these observations in hand, we have

[TABLE]

If we set the constant above to the following average value, $C:=\frac{1}{\mathrm{vol}(T)}\int_{T}(\nabla u\cdot\delta)(\nabla w\cdot\delta)\,\mathrm{d}x$, then $\|(\nabla u\cdot\delta)(\nabla w\cdot\delta)-C\|_{L^{1}(T)}\lesssim H_{T}\hskip 0.25pt|(\nabla u\cdot\delta)(\nabla w\cdot\delta)|_{W^{1,1}(T)}$. Therefore,

[TABLE]

After summing over all $T\in\mathcal{T}_{H}$ and $\delta\in\mathscr{D}(T_{m})$ and taking into account $|\delta|\eqsim h$, we arrive at the bound

[TABLE]

*The proof is completed by recalling that $\|w\|_{2}\lesssim\|u_{h}-\widetilde{u}_{h}\|_{0}$. *

Remark 7.8**.**

*As stated previously, the error estimates in Theorems 7.4 and 7.6 not do track the constants appearing in the corresponding classical estimates. When $r>0$, none of the constants in either classical estimate depend on the higher order norms $|K|_{W^{r+1,\infty}(\Omega)}$. *

Remark 7.9**.**

*The main ingredients in the error analysis presented in this section are Lemma 6.4, Strang’s First Lemma, and 40, which simply follows from the zero row sum property and symmetry. Simple generalizations of Lemma 6.4, for the bilinear forms $a_{2}(\cdot,\cdot)$ and $a_{3}(\cdot,\cdot)$ have been stated in 38. Meanwhile, 40 simply follows from the zero row sum property of $\widetilde{\mathsf{A}}$ and symmetry. Therefore, we see no reason to doubt that our analysis here can be generalized to the other problems of interest in this paper. Hence, we proceed with numerical verification and demonstration. *

8 Implementation

The performance of a numerical method largely depends on its implementation. Therefore, in this section, we highlight the important features of ours. We use the HyTeG finite-element software framework [29] as the core framework for all the numerical experiments in Section 9. It offers efficient distributed data structures for simplicial meshes in 2D and 3D, which serve as a basis for the implementation of massively parallel fast iterative solvers. Its main concept is based on the idea that a coarse input mesh is split into its geometrical primitives, i.e., vertices, edges, and faces, and each of these primitives is uniformly refined. Because the primitives of the same dimension are decoupled from the others, all primitives of the same dimension may be processed in parallel. This partitioning and the hierarchy of locally structures meshes allows for efficient parallel implementations of geometric multigrid methods. More importantly, these data structures fit perfectly to the concept of macro-elements introduced in Section 3.3. The problems in this paper are mainly solved by employing a geometric multigrid solver using V-cycles with a hybrid Gauss–Seidel smoother. On the coarsest grid, either a preconditioned conjugate gradient method or the direct solver MUMPS [1, 2], as provided by the PETSc interface [4, 5], is used. For improved parallel scalability of the coarse grid solver, agglomeration techniques as provided by PCTELESCOPE [35] are used in runs with many processes.

8.1 Polynomial least squares regression

An important factor in the performance of the surrogate approach is the approximation of the stencil functions $\Phi^{\delta}_{T}$ by polynomials $\widetilde{\Phi}^{\delta}_{T}$. This step in the solver process must be very fast and is usually done in a pre-process step before the actual solve. After various preparatory experiments, we have seen satisfactory performance and accuracy from simply computing $\widetilde{\Phi}^{\delta}_{T}=\Pi^{\delta}_{T}\Phi^{\delta}_{T}$ via solving a simple least-squares problem, which we now describe.

Let $T\in\mathcal{T}_{H}$ be a macro-element and recall that $T_{m}$ is the associated lattice on level $m$. Suppose that $\Phi_{T}^{\delta}$ is the stencil function in direction $\delta\in\mathscr{D}(T_{m})$ which we want to approximate. For the least-squares regression, we fix a level $m_{\mathrm{LS}}$ with $m\geq m_{\mathrm{LS}}\geq 2$ and define the set of least-squares points ${T}{}^{\delta}_{\mathrm{LS}}\coloneqq{T}_{m_{\mathrm{LS}}}\cap\overline{T}_{\delta}$. Furthermore, let $\{p_{k}\}_{k=1}^{M}$ be a basis of $\mathcal{P}_{q}(T)$, the space of polynomials with maximal degree $q$. Assume that $m_{\mathrm{LS}}$ is chosen large enough such that \big{|}{T}{}^{\delta}_{\mathrm{LS}}\big{|}\geq M and introduce the following norm on $\mathcal{P}_{q}(T)$: $\|p\|^{2}_{{T}{}^{\delta}_{\mathrm{LS}}}\coloneqq\sum_{x_{i}\in{T}{}^{\delta}_{\mathrm{LS}}}p(x_{i})^{2}$. The least-squares regression problem, which in turn defines $\Pi^{\delta}_{T}$, is formalized as follows:

[TABLE]

The approximated stencil function is then defined as $\widetilde{\Phi}^{\delta}_{T}\coloneqq\sum_{k=1}^{M}c_{k}p_{k}$. This problem is equivalent to solving the possibly overdetermined linear system of equations $\mathsf{B}\mathsf{c}=\mathsf{f}$ in a least-squares sense, where $\mathsf{B}_{ij}=p_{j}(x_{i})$ and $\mathsf{f}_{i}=\Phi^{\delta}_{T}(x_{i})$ for $1\leq i\leq|{T}{}^{\delta}_{\mathrm{LS}}|$ and $1\leq j\leq N$. The choice of the polynomial basis is arbitrary. However, for an easier implementation, we employ the monomial basis, even knowing that the resulting linear system is ill conditioned. Since it is crucial to get numerically precise results, a stable solver for this problem has to be chosen. For this purpose, we apply the colPivHouseholderQr method from the Eigen 3.3.5 library [28], which offers a good balance between speed and accuracy. Obviously, each of these linear systems is independent of others, therefore they may be solved in parallel.

Remark 8.1**.**

*Taking into account $\widetilde{\Phi}^{\delta}_{T}=\sum_{k=1}^{M}c_{k}p_{k}$, the first order optimality condition for 60 can be stated as \sum_{x_{i}\in{T}{}^{\delta}_{\mathrm{LS}}}\big{[}\Phi^{\delta}_{T}-\Pi^{\delta}_{T}\Phi^{\delta}_{T}\big{]}(x_{i})=0. If $m=m_{\mathrm{LS}}$, then the secondary assumption in Theorem 7.6 is satisfied and we see higher order convergence in $H$, as stated in 48b. Usually, when $m_{\mathrm{LS}}$ is close but not equal to $m$, we see preasymptotic $H$-convergence in between the two estimates given in 48; see Figure 6. *

Remark 8.2**.**

In the case where the bilinear form $a(\cdot,\cdot)$ is symmetric, we need only approximate a single stencil function $\Phi^{\delta}_{T}$ for both directions $\delta$ and $-\delta$. Indeed, as observed in Remark 3.1, the corresponding stencil functions are identical, up to a shift by $\delta$. Furthermore, the symmetry requirement 19, from Remark 3.5, is satisfied with the projection operator defined above. Indeed, one may verify that for every $\delta$, $T_{\delta}=T_{-\delta}-\delta$. Therefore,

[TABLE]

*Thus, on simplicial meshes in 2D, only four instead of seven polynomials per macro-element have to be determined and stored in memory. In some cases, where the zero row sum property holds, the number of required polynomials may be even reduced to three. *

8.2 Fast polynomial evaluation

An even more important factor with respect to the performance of our implementation is the fast evaluation of the surrogate stencil functions $\widetilde{\Phi}^{\delta}_{T}$. Contrary to the computation of each $\widetilde{\Phi}^{\delta}_{T}$, which will happen only once per solve, these evaluations will be made during every matrix-vector multiplication. Therefore, the costs of evaluating the stiffness matrix entries associated to a degree of freedom, may not exceed the costs of evaluating the bilinear forms with the respective ansatz functions. In this case not only the reduction of floating point operations per degree of freedom is of importance, but also the required memory traffic has to be taken into account.

When performing a matrix-vector multiplication in HyTeG, the degrees of freedom in a macro-element are processed in a row-wise fashion as illustrated in the left of Fig. 4. In each row, the stencil function may be interpreted as a 1D function. We assume without loss of generality, that the 1D stencil functions are aligned with the $x$-axis. This property is also inherited by the approximated stencil function.

To further optimize the evaluation of the 1D polynomial, we exploit that the stencil functions have to be evaluated on a line subdivided into uniformly sized intervals of length $h$. Let $(x_{i},y_{j})$ be a vertex node in the lattice and let $p_{y_{j}}(\cdot)\coloneqq\widetilde{\Phi}_{\delta}(\cdot,y_{j})$ be the approximated 1D stencil function associated to row $y_{j}$.

Assuming that we already have evaluated the stencil function $p_{y_{j}}$ at a point $x_{i}$, we want to evaluate it at the next point $x_{i}+h$ as efficiently as possible. Since the grid points are equidistantly distributed, we can use a special case of the divided differences, called forward differences [18, pg. 126].

First, we need $q+1$ helper variables $\left\{\Delta_{x_{0}}^{(k)}\right\}$ for $k\in\{0,1,\ldots,q\}$ which are defined in a preprocessing step as follows:

[TABLE]

The value at position $p_{y_{j}}(x_{0})$ is then given by $\Delta_{x_{0}}^{(0)}$. In order to obtain the value at $p_{y_{j}}(x_{0}+h)$, one has to update all the helper variables in the following way:

[TABLE]

After that, the value of $p_{y_{j}}(x_{0}+h)$ is given by $\Delta_{x_{0}}^{(0)}$. Doing this recursively yields the approximated stencil function values at all mesh points on a single row using only $q+1$ helper variables and $q+1$ floating point additions. When iterating through a row, in general seven polynomial evaluations, one for each direction, are required, cf. right of Fig. 4. In the symmetric case this may be reduced to six polynomial evaluations, since the western stencil weight may be obtained from the previous eastern evaluation. Keep also in mind that in the symmetric case the polynomials of approximated stencil functions in opposite directions are the same but only evaluated at different positions, cf. Remark 8.2. Therefore, $6\cdot(q+1)$ helper variables are required for a single row. For $q=8$, these results in $54\cdot 8$ bytes of memory which fits easily into a modern L1 CPU cache. When moving from one lattice point to another, $6\cdot(q+1)$ vectorizable floating point additions have to be performed, to obtain the updated polynomial evaluations. Furthermore, in our implementation, the polynomial degree is realized as a C++ template parameter, therefore, all loops concerning the evaluation of a polynomial of a certain degree may be optimized at compile time. Since our focus lies in the theoretical analysis of the surrogate approach, thorough performance studies employing performance models should be considered beyond the scope of this paper. Similar performance studies have been carefully completed in [9, 10]. Thus, in the next section, we only report on relative run-times of the surrogate approach compared to the standard method, also implemented on HyTeG, using on-the-fly quadrature of the integrals stemming from the bilinear forms.

9 Numerical experiments

In this section, we numerically verify Theorems 7.4 and 7.6, both related to the variable coefficient Poisson equation. Additionally, we present proof-of-concept results for a linearized elasticity application and a simple $p$-Laplacian diffusion problem. While not covered by the theory, we include these latter examples to demonstrate the breadth of generality of the methodology.

All run-time measurements in this sections were obtained on a machine equipped with two Intel® Xeon® Gold 6136 processors with a nominal base frequency of 3.0 GHz. Each processor has 12 physical cores which results in a total of 24 physical cores. The total available memory of $251\text{\,}\mathrm{GB}$ is split into two NUMA domains, one for each socket. We use the GCC 7.3.0 compiler and specify the following compiler arguments: -O3 -march=native. All the examples in this section were executed in parallel using all available 24 physical cores.

When comparing run times from the standard and the surrogate approaches, many factors are responsible for the relative speed-up of the surrogate approach. Increasing the polynomial order $q$ of the surrogate stencil functions not only increases the run time of a multigrid iteration but also the time spent in the setup phase (i.e., computing each $\Phi_{T}^{\delta}$). The cost of the setup phase, however, is mostly dominated by the sampling level $m_{\mathrm{LS}}$. Therefore, when the ratio of time spent in the iterative solver to the time spent in the setup phase(s) for solving a particular problem is large, the setup cost is almost negligible and we see the best performance. Since the problems in the following subsections differ in complexity and have different ratios of solver to setup time, the observed relative speed-ups are not directly comparable. Nonetheless, the reported speed-ups for all tested examples range between a factor of 14 and 20. Such significant speed-ups are in particular important in case of dynamic or stochastic applications. Most stochastic applications demand an enormous number of deterministic solves resulting quite often in extreme long run times. Having such a surrogate approach at hand can help to make stochastic approaches such as, e.g., multilevel Monte Carlo and its variants, more accessible for complex applications.

9.1 Quantitative benchmark problem

In this subsection, we examine the surrogate method for the variable coefficient Poisson equation which has been described and analyzed above. The strong form of the problem is

[TABLE]

We consider both the bilinear form coming from the scalar coefficient scenario (i.e., $K=k\cdot\mathrm{Id}$), introduced in 4, and the tensorial coefficient scenario, introduced in 15. In the scalar coefficient experiments, we use the unit-square domain $\Omega=(0,1)^{2}$. In the tensorial coefficient experiments, the domain $\Omega$ has a curvilinear boundary.

9.1.1 Scalar coefficient on unit square

In the first benchmark problem ($K=k\cdot\mathrm{Id}$), we take $\Omega=(0,1)^{2}$ and employ the scalar coefficient function

[TABLE]

in problem 65. The manufactured solution $u$ is chosen as $u(x,y)=\sin(x)\sinh(y)$. The restriction of $u$ to the boundary is chosen as Dirichlet datum $g$. The right-hand-side $f$ is directly computed by inserting $u$ into the equation. In this benchmark, we fix the finest mesh size $h$ and report on the errors depending on $H$ and $q$ to show the proven $\mathcal{O}\left(H^{q+1}\right)$ estimate in the $H^{1}$ norm and $\mathcal{O}\left(H^{q+2}\right)$ in the $L^{2}$ norm. For this purpose, $h$ is chosen to be very small in order for the error to be mostly dominated by the surrogate part.

The reference macro-mesh size is given by $H_{0}$, as illustrated in the left of Fig. 5. All finer macro-meshes, with associated mesh sizes $H<H_{0}$, stem from uniformly refining this reference mesh; see middle and right of Fig. 5. The fine mesh, with associated mesh size $h\ll H$, is the 13 times uniformly refined reference macro-mesh, i.e., $h=2^{-13}H_{0}$. This fine mesh, has about $6.71\text{\times}{10}^{7}$ degrees of freedom. The approximation of the stencil functions through least-squares regression, is done on the mesh associated to mesh size $H_{\mathrm{LS}}=2^{-8}H$. Note that this keeps the number of sampling points in each macro-element constant to $32\,639$. Each linear system is solved by applying geometric multigrid V(2,2) iterations until a relative residual of $1\text{\times}{10}^{-13}$ is obtained.

In Tables 1 and 2, the relative $H^{1}$ and $L^{2}$ errors for decreasing mesh sizes $H$ are shown. Both tables show the expected convergence rates. In the case of the $L^{2}$ norm for $q=3$ and $q=4$, the convergence rate deteriorates for small macro-mesh sizes $H$ because the discretization error is dominating the total error.

In order to show the dependence of the least-squares approach on the sampling level, we provide Figure 6. Here, we show two plots of the relative $L^{2}$ errors for fixed $q\in\{1,3\}$, $h=2^{-13}H_{0}$, and varying $H_{\mathrm{LS}}$. From this figure, one can see that it is crucial to tune the sampling level fine enough in order to achieve optimal $\mathcal{O}\left(H^{q+2}\right)$ convergence in the $L^{2}$-norm.

Remark 9.1**.**

*The choice of sampling level $m_{\mathrm{LS}}$ or, equivalently, $H_{\mathrm{LS}}$ is very important, since the cost of the polynomial regression grows exponentially with $m_{\mathrm{LS}}$. However, choosing a too large $H_{\mathrm{LS}}$ may violate the discrete $L^{2}$ projection property required in Theorem 7.6 in order to obtain an increased order of convergence. Therefore, it is crucial to choose a suitable $H_{\mathrm{LS}}$ for an optimal ratio between the accuracy of the solution and the run time of the stencil function approximation. In each of our experiments, setting $H_{\mathrm{LS}}$ two to four times larger than the fine mesh size $h$ yielded satisfactory results with respect to accuracy and run time. *

Furthermore, in Fig. 7 we want to illustrate the dependence of the polynomial degree $q$ and the macro-mesh size $H$ within the surrogate approach. For this purpose, we plot the central true and surrogate stencil functions over the subdomain $\{(x,y)^{\top}\in\Omega:x+y\geq 1\}$ for different pairings of $q$ and $H$. It can be observed, that there is no visible difference of both functions when either the pairing $H=H_{0}$ and $q=8$, or the pairing $H=\nicefrac{{H_{0}}}{{8}}$ and $q=2$ is chosen. Obviously, the quality of $\widetilde{\mathsf{A}}$ can be improved by either increasing $q$ or decreasing $H$. For smooth coefficients $K$, increasing $q$ is the more efficient option, like in the $hp$-FEM context.

9.1.2 Tensor coefficient on domain with curved boundaries

In the second benchmark problem, we study problem 65 with the symmetric and positive definite tensor coefficient

[TABLE]

Moreover, we consider the domain $\Omega$ with the curved boundary illustrated in Fig. 8. In the following scenarios, $a=0.1$ is used as the amplitude of the boundary perturbation. The mapping from the reference unit-square to the perturbed domain is defined by $\varphi$ in 68. To map the coefficient onto the perturbed domain, we replace the coefficient $K$ in 65 by a new coefficient, $K_{0}$, induced by the domain transformation, viz.,

[TABLE]

The manufactured solution $u$ is chosen to be $u(x,y)=\sin(\varphi_{1}(x,y))\sinh(\varphi_{2}(x,y))$. The restriction of $u$ to the boundary is chosen as Dirichlet datum $g$ and the right-hand-side $f$ is directly computed by inserting $u$ into the strong form of the equation 65.

Here, we perform the same verification as in the previous subsection. That is, fixing $h$ and varying $q$ and $H$ with the same meshes and solver settings. In Tables 3 and 4, the relative $H^{1}$ and $L^{2}$ errors for decreasing mesh sizes $H$ are shown. Both tables show the expected convergence rates. In the case of $q=4$, the convergence rate deteriorates for small macro-mesh sizes $H$, because the discretization error is dominating.

Additionally, we present results for fixed $H=2^{-3}H_{0}$ and varying $h$ and $q$. Table 5 shows the relative $L^{2}$ errors and convergence rates of the standard approach and the surrogate approach with $q\in\{3,5,7\}$. Only for $q=7$, the $L^{2}$ error coincides with the errors from the standard approach for all $h$. The relative time-to-solution (rtts) shown for the surrogate approaches is defined as the time-to-solution (tts) including the setup-phase of the surrogate approach divided by the time-to-solution of the standard approach. In the case with the smallest $h$, the surrogate approach took at most only $7\%$ of the time of the standard approach. That is, a speed-up by more than a factor of $14$.

9.2 Linearized elasticity example

In this subsection, we compare a standard method with a surrogate method applied to the linearized elasticity problem presented in Section 4.2. In our surrogate method, we employ the zero row sum property described in Section 5. We choose an annular domain composed of two distinct and concentric materials under uniform pressure loading. The problem is inspired by a similar 3D experiment documented in [25]. Let $\mathcal{B}_{r}\operatorname{\subseteq}\mathbb{R}^{2}$ be the two-dimensional open ball of radius $r$ with the midpoint at the origin. The computational domain is then defined as $\Omega=\mathcal{B}_{R_{\mathrm{out}}}\setminus\mathcal{B}_{R_{\mathrm{in}}}$. We split this domain into two disjoint sets $\Omega_{I}=\mathcal{B}_{R_{\mathrm{mid}}}\setminus\mathcal{B}_{R_{\mathrm{in}}}$ and $\Omega_{O}=\mathcal{B}_{R_{\mathrm{out}}}\setminus\mathcal{B}_{R_{\mathrm{mid}}}$ corresponding to each material. In our experiments, we fix $R_{\mathrm{in}}=$1\text{,}\mathrm{cm}, $R_{\mathrm{mid}}=$1.75\text{\,}\mathrm{cm}, and $R_{\mathrm{out}}=$2\text{,}\mathrm{cm}$$. Refer to the leftmost diagram in Fig. 9 for an illustration of the setup. Here, the macro-elements adjacent to the boundary and the material interface are mapped to the physical geometry by using the transformation described in [39].

Let $r(x)=|x|$. The strong form of the problem is

[TABLE]

Here, the stress tensor $\sigma$ is given by Hooke’s law for isotropic materials as defined in Section 4.2. The unit vector in radial direction is denoted $\hat{e}_{r}$ and the outward pointing unit normal vector is denoted $\hat{n}$. We neglect body forces and therefore set $\vec{f}=\vec{0}$. The displacement is described in polar coordinates where $u_{r}$ is the radial displacement and $u_{\theta}$ is the tangential displacement. In order to make the system uniquely solvable, we enforce the tangential displacement $u_{\theta}$ to be zero at three points; see Fig. 9. The materials are chosen to be cork in the inner domain $\Omega_{I}$ with an A36 steel layer in the outer domain $\Omega_{O}$. Poisson’s ratio and Young’s modulus for this scenario are $E_{I}=$0.02\text{,}\mathrm{GPa}, $E_{O}=$200.0\text{\,}\mathrm{GPa}, $\nu_{I}=0$, and $\nu_{O}=0.26$. The Lamé parameters are obtained by the expressions $\mu=\frac{E}{2(1+\nu)}$ and $\lambda=\frac{E\nu}{(1-2\nu)\cdot(1+\nu)}$. Note that while these expressions induce piecewise-constant Lamé parameters, $\lambda=\lambda(r)$ and $\mu=\mu(r)$, once the smooth domain is mapped to the computational domain depicted on the right of Fig. 9, these parameters will not be piecewise-constant anymore due to the transformation. The pressure on the outer boundary is set to $p_{\mathrm{out}}=$0\text{,}\mathrm{MPa} and on the inner boundary $p_{\mathrm{in}}=$1\text{\,}\mathrm{MPa} is prescribed.

In this particular scenario, there is an analytic solution available for the radial displacement $u_{r}$ which has the form

[TABLE]

The tangential displacement $u_{\theta}$ is zero everywhere due to the symmetry of the problem. In 70, the coefficients $A$, $B$, $C$, and $D$ are uniquely determined by the following system of linear equations:

[TABLE]

where $d_{I}\coloneqq 2\nu_{I}^{2}+\nu_{I}-1$ and $d_{O}\coloneqq 2\nu_{O}^{2}+\nu_{O}-1$. This system is derived after deducing the continuity of the displacement and surface traction at the material interface, then by incorporating the prescribed external forces at the boundaries.

In order to verify the accuracy of the surrogate method, we select the polynomial degree $q=7$ and the macro-mesh $\mathcal{T}_{H}$, illustrated on the right of Fig. 9. Note that the discontinuity in the material parameters lies along the macro-element interfaces and so $\lambda,\mu\in\prod_{T\in\mathcal{T}_{H}}C^{\infty}(T)\subsetneq W^{r+1,\infty}(\mathcal{T}_{H})$, for any $r>0$.

Each linear system is solved by applying geometric multigrid iterations with V(3,3) cycles until the relative residual is reduced by the factor $1\text{\times}{10}^{-7}$. On the coarsest level used in the multigrid hierarchy, we employ MUMPS as a direct solver. In Table 6, we report on the results for varying $h$ and present the relative $L^{2}$ errors for the standard and surrogate approach, respectively. On the finest mesh involving about $1.8\text{\times}{10}^{8}$ degrees of freedom, the surrogate approach required only 5% of the time required by the standard approach while having the same accuracy. That is a speed up by a factor of $20$. Fig. 10 shows the radial displacement $u_{r}$ and the tangential stress $\sigma_{\theta\theta}$ computed with the surrogate approach on the fine mesh $\mathcal{S}^{6}(\mathcal{T}_{H})$, corresponding to $h=2^{-6}H$, with $q=4$. This is illustrated further by the plots in Fig. 11 which allow a visual comparison between $u_{r}$ and $\sigma_{\theta\theta}$ in the surrogate and analytical solutions.

9.3 $p$-Laplacian diffusion example

In this subsection, we consider the time-dependent example introduced in Section 4.3. Here, we solve the non-linear $p$-Laplacian diffusion problem 72, given in strong form as

[TABLE]

The computational domain is set to the unit disk, i.e., $\Omega\coloneqq\mathcal{B}_{1}$ and the right-hand-side is set to a specific constant, $f(x)=2\hat{q}^{\nicefrac{{p}}{{\hat{q}}}}$, where $\hat{q}=\frac{p}{p-1}$. The initial solution is set to $u_{0}(x)=0.1\cdot\left(1-|x|^{2}\right)$. For this particular problem, the stationary limit $u_{\infty}$ has an analytic solution exists with unit magnitude, namely $u_{\infty}(x)=1-|x|^{\hat{q}}$ [6, Example 3.1].

Our discretization follows a standard approach where a mass matrix $\mathsf{M}_{ij}=\int_{\Omega}\phi_{i}\phi_{j}\hskip 0.25pt\,\mathrm{d}x$ and a stiffness matrix $\mathsf{A}_{ij}(\widetilde{\mathsf{u}})=\int_{\Omega}|\nabla\widetilde{u}_{h}|^{p-2}\hskip 0.25pt\nabla\phi_{j}\cdot\nabla\phi_{i}\hskip 0.25pt\,\mathrm{d}x$ are introduced. At this point, $\widetilde{\mathsf{u}}$ is the coefficient vector, in the $\{\phi_{i}\}$ basis, of an arbitrary discrete function $\widetilde{u}_{h}$. The time derivative is discretized by a backward Euler scheme, and the non-linearity in each time step is resolved by Picard fixed-point iterations. Let $\mathsf{u}^{l}_{k}$ be the coefficient vector of the discrete solution at the $k$-th time step and $l$-th fixed point iteration. Employing the bilinear form 23 and fixing a time step size $dt>0$, the discrete problem in each time step $k>0$ and fixed point iteration $l>0$ reads as follows:

[TABLE]

where $\mathsf{u}_{k-1}$ is the final coefficient vector from the previous time step. In each time step, this system is solved multiple times (once for each fixed-point iteration) by the application of five V(2,2) multigrid cycles. The fixed-point iterations continue until the relative increment $\frac{\|\mathsf{u}^{l}_{k}-\mathsf{u}^{l-1}_{k}\|_{2}}{\|\mathsf{u}^{l}_{k}\|_{2}}$ is smaller than the fixed tolerance $1\text{\times}{10}^{-3}$. Then $k$ is incremented and a new $\mathsf{u}_{k-1}=\mathsf{u}_{k-1}^{l}$ is defined.

In our surrogate method, the stencil functions of the stiffness matrix \mathsf{A}\big{(}\mathsf{u}^{k-1}\big{)} are approximated by solving the least-squares problems after every fixed-point iteration, all the while enforcing the zero row sum property (cf. Section 9.1). Meanwhile, the stencil function of the mass matrix $\mathsf{M}$ is only approximated once in a pre-processing step because it does not depend on any free variables in the computation. The time step surrogate polynomials of both operators are then simply summed together to obtain the time step matrix \mathsf{M}+dt\hskip 0.25pt\mathsf{A}\big{(}\mathsf{u}^{k-1}\big{)}. This particular splitting of the surrogate matrices, which reproduces the zero row sum property in the stiffness matrix, allows for faster re-approximation of the time step matrix stencil function and appears to improve the stability of the method. In this example, we did not enforce the symmetry condition featured in Section 3.5. Instead, whenever a vertex $x_{i}\in\mathbb{X}_{m}$ was on the boundary of a macro-element $\partial T_{m}$, we set the surrogate stiffness matrix to the exact value stiffness matrix $\widetilde{\mathsf{A}}_{ij}=\mathsf{A}_{ij}$. This minor asymmetry, is more amenable to computation because there is less data transfer and it did not affect the behavior of our multigrid solver. In fact, this choice improved our results with this problem, which we believe is due to better accuracy in the surrogate near the singularity in the coefficient; i.e., at the origin $x=(0,0)$. The success of this approach suggests that the definition given in 18 may be relaxed in other applications as well.111See [10] for further evidence. In the proximity of this singularity and for $p>2$, the coefficient depending on the solution of the previous fixed-point iteration is getting very close to zero which serves as a challenge for the approximated off-diagonal stencil functions. Depending on the polynomial degree $q$, they might erroneously take on positive values due to overshoots which possibly results in a loss of positive definiteness of the surrogate matrix. However, this drawback could not be observed in the scenario considered in the following example.

The unit-disk is discretized by the macro-mesh $\mathcal{T}_{H}$ featured on the right of Fig. 2. Note that the vertices of the central macro-elements meet at the origin $x=(0,0)$; i.e., exactly where the singularity occurs in the stationary limit $u_{\infty}$. The simulations are conducted on the mesh $\mathcal{S}^{9}(\mathcal{T}_{H})$, which involves about $4.72\text{\times}{10}^{6}$ degrees of freedom. The macro-elements adjacent to the boundary are mapped to the physical geometry by using the mapping described in [39]. Furthermore, the least-squares regressions are carried out on the mesh corresponding to $H_{\mathrm{LS}}=4h$ and the polynomial degree of the approximated stencil functions is fixed to $q=6$. In this scenario, we consider the $p$-Laplacian operator with $p=3$, fix the time step size $\Delta t=$1\text{\times}{10}^{-2}$$, and solve until time $T_{\mathrm{end}}=1$. Figure 12 illustrates the standard and surrogate solutions plotted over the line $[0,1]\times\{0\}$ for different times $t$. In the left of Fig. 13, the surface plot of the surrogate solution is depicted. Since the difference of the solutions is very small, we added in the right of Fig. 13 a surface plot of the absolute difference of the surrogate and standard solution at the final time $t=T_{\mathrm{end}}$. The surrogate approach required only $5.4$% of the time required by the standard approach. That is, a speed-up by more than a factor of $18$.

Appendix A Proofs

Proof A.1** (Proof of Proposition 6.1).**

By the min-max theorem [21], the $k$-th eigenvalue of $\mathsf{M}$ is

[TABLE]

Define $\mathsf{D}=\mathsf{M}-\mathsf{N}$. We first show that $\lambda_{1}(\mathsf{D})\leq\lambda_{k}(\mathsf{M})-\lambda_{k}(\mathsf{N})\leq\lambda_{N}(\mathsf{D})$. Indeed,

[TABLE]

and, likewise,

[TABLE]

*This immediately leads us to the inequality $|\lambda_{k}(\mathsf{M})-\lambda_{k}(\mathsf{N})|\leq\operatorname*{\vphantom{p}max}\{|\lambda_{1}(\mathsf{D})|,|\lambda_{N}(\mathsf{D})|\}$. Now, for at least one $i$, $|\lambda_{N}(\mathsf{D})|-|\mathsf{D}_{ii}|\leq|\lambda_{N}(\mathsf{D})-\mathsf{D}_{ii}|\leq\sum_{j\neq i}|\mathsf{D}_{ij}|$, by the Gershgorin circle theorem. Therefore, $|\lambda_{N}(\mathsf{D})|\leq\sum_{j}|\mathsf{D}_{ij}|\leq\|\mathsf{D}\|_{\infty}$. Similarly, $|\lambda_{1}(\mathsf{D})|\leq\|\mathsf{D}\|_{\infty}$. *

Acknowledgments

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 800898. This work was also partly supported by the German Research Foundation through the Priority Programme 1648 "Software for Exascale Computing" (SPPEXA) and by grant WO671/11-1. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (GCS, www.gauss-centre.eu) for funding this project by providing computing time on the GCS supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de).