The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Discrete Equivalence of Adjoint Neumann–Dirichlet div-grad and grad-div Equations in Curvilinear 3D Domains**

**Yi Zhang, Varun Jain, Artur Palha, and Marc Gerritsma**

### **1 Introduction**

In R*d* , given a bounded domain with Lipschitz boundary *∂* and *σ*ˆ *n* ∈ *H* −1*/*2*(∂)* = tr *H (*div*, )*, *ω* ∈ *H*1*()* solves the Neumann problem,

$$\begin{cases} \begin{aligned} \frac{\partial \boldsymbol{\omega}}{\partial \boldsymbol{\pi}} &= \hat{\boldsymbol{\sigma}}\_{\text{H}} \qquad & \text{on } \partial \Omega \\ -\text{div}\left(\text{grad } \boldsymbol{\omega}\right) + \boldsymbol{\omega} &= 0 \qquad & \text{in } \Omega \end{aligned} \end{cases} , \tag{1}$$

if and only if *σ* ∈ *H (*div*, )* which solves the Dirichlet problem,

$$\begin{cases} \sigma \cdot \mathfrak{n} = \hat{\sigma}\_{\mathfrak{n}} & \text{on } \partial \Omega \\\\ -\text{grad}\,(\text{div } \sigma) + \sigma = 0 & \text{in } \Omega \end{cases}, \tag{2}$$

satisfies *σ* = grad *ω* [3]. This is obvious at the continuous level. The question is whether we can find a set of finite dimensional function spaces such that *σ h* = grad *ωh* holds if *ωh* and *σ h* solve the discrete Neumann and Dirichlet problems respectively. The answer is yes.

Y. Zhang (-) · V. Jain · M. Gerritsma

Delft University of Technology, Delft, Netherlands e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
© The Author(s) 2020

A. Palha Eindhoven University of Technology, Eindhoven, Netherlands e-mail<EMAIL_ADDRESS>
S. J. Sherwin et al. (eds.), *Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018*, Lecture Notes in Computational Science and Engineering 134, https://doi.org/10.1007/978-3-030-39647-3\_15

Throughout this paper, we restrict ourselves to R3. We will first construct the primal polynomial spaces and their algebraic dual representations, and then use them to discretize problems (1) and (2) such that the identity *σ h* = grad *ωh* holds at the discrete level in any curvilinear domain for any polynomial approximation degree. This work extends [7, 9], where similar dual Neumann–Dirichlet problems are considered, to 3-dimensional space. These primal spaces and their algebraic dual representations can be ideal for the so-called mimetic or structure-preserving discretizations [1, 4, 8, 11, 12]. Together with their trace spaces, they can be used for the hybrid finite element methods which first decompose the domains into discontinuous elements then connect them with Lagrange multipliers living in the trace spaces [2, 13, 14].

The outline of this paper is as follows: In Sect. 2, we introduce the construction of polynomial spaces and their algebraic dual representations. The discrete formulations of the Neumann–Dirichlet problems and the proof of their equivalence at the discrete level follow in Sect. 3. A 3-dimensional numerical test case is then presented in Sect. 4. Finally, conclusions are drawn in Sect. 5.

#### **2 Function Spaces**

#### *2.1 Primal Polynomial Spaces*

Let −1 = *ξi* 0 *< ξi* 1 *<* ··· *< ξi I i* = 1, *i* = 1*,* 2*,* 3, being three partitionings of [−1*,* 1]. The associated Lagrange polynomials are

$$h\_j(\boldsymbol{\xi}^i) = \prod\_{m=0, m \neq j}^{I^i} \frac{\boldsymbol{\xi}^i - \boldsymbol{\xi}\_m^i}{\boldsymbol{\xi}\_j^i - \boldsymbol{\xi}\_m^i}, \ j = 0, 1, \dots, I^i.$$

They are polynomials of degree *I i* which satisfy the Kronecker delta property, *hj (ξi k )* = *δjk*. The associated edge functions can be derived as [6],

$$e\_j(\xi^i) = -\sum\_{k=0}^{j-1} \frac{\mathrm{d}h\_k(\xi^i)}{\mathrm{d}\xi^i}, \ j = 1, 2, \cdots, I^i,$$

which are polynomials of degree *I i* − 1. Edge functions also satisfy the Kronecker delta property, but in the integral sense,

$$\int\_{\xi\_{k-1}^{\vec{\ell}}}^{\xi\_k^{\vec{\ell}}} e\_f(\xi^{\vec{\ell}}) \, \mathrm{d}\xi^{\vec{\ell}} = \delta\_{jk}.$$

Consider a reference domain ref|*ξ* 1*,ξ* 2*,ξ* 3 := [−1*,* 1] 3. With the tensor product, we can construct finite dimensional scalar function space P*I* 1*,I* 2*,I* 3 spanned by polynomial basis functions

$$\left\{ h\_i(\xi^1) h\_j(\xi^2) h\_k(\xi^3) \right\},$$

and vector-valued function space L*I* 1*,I* 2*,I* 3 spanned by polynomial basis functions

$$\left\{ e\_l(\xi^1) h\_j(\xi^2) h\_k(\xi^3), \ h\_l(\xi^1) e\_j(\xi^2) h\_k(\xi^3), \ h\_l(\xi^1) h\_j(\xi^2) e\_k(\xi^3) \right\}.$$

Let *ωh* ∈ P*I* 1*,I* 2*,I* 3 be

$$
\omega^h = \sum\_{i=0}^{I^1} \sum\_{j=0}^{I^2} \sum\_{k=0}^{I^3} w\_{i,j,k} h\_i(\xi^1) h\_j(\xi^2) h\_k(\xi^3). \tag{3}
$$

Due to the way of constructing the edge functions, we can easy derive *ρh* = grad *ωh* ∈ L*I* 1*,I* 2*,I* 3 ,

$$
\rho^h = \text{grad } \omega^h = (\rho\_1, \,\rho\_2, \,\rho\_3)^\mathsf{T},
$$

where [6],

$$\rho\_1 = \sum\_{i=1}^{I^1} \sum\_{j=0}^{I^2} \sum\_{k=0}^{I^3} \left( w\_{i,j,k} - w\_{i-1,j,k} \right) e\_i(\xi^1) h\_j(\xi^2) h\_k(\xi^3),$$

$$\rho\_2 = \sum\_{i=0}^{I^1} \sum\_{j=1}^{I^2} \sum\_{k=0}^{I^3} \left( w\_{i,j,k} - w\_{i,j-1,k} \right) h\_i(\xi^1) e\_j(\xi^2) h\_k(\xi^3),$$

$$\rho\_3 = \sum\_{i=0}^{I^1} \sum\_{j=0}^{I^2} \sum\_{k=1}^{I^3} \left( w\_{i,j,k} - w\_{i,j,k-1} \right) h\_i(\xi^1) h\_j(\xi^2) e\_k(\xi^3).$$

Let *ω*, *ρ* be the vectors of expansion coefficients of *ωh*, *ρh*. We can obtain

$$
\underline{\underline{\rho}} = \mathbb{E}\,\underline{\underline{\omega}},\tag{4}
$$

where E is called the incidence matrix. The incidence matrix is very sparse, only consists of ±1 as non-zero entries. If we squeeze, stretch or distort the domain, of course, the polynomial basis functions change, but the incidence matrix will remain the same. It only depends on the topology of the mesh and the numbering of the degrees of freedom. And it is exact. In other words, it introduces no extra error. All these features make it an excellent discrete counterpart of the grad operator. Examples of incidence matrices can be found in [8, 10–12].

For a comprehensive explanation of these polynomial basis functions, we refer to [6]. In isogeometric analysis, tensor-product B-splines with similar properties have been developed, see, for example [5]. For tetrahedral elements, an analogue development can be found in [15].

From (3), we can derive the trace of *ωh*, for example, on the back boundary of ref, b = *ξ* 1 = −1*, ξ* 2*, ξ* 3 ∈ [−1*,* 1] ,

$$\text{tr}\_{\mathbf{b}}\,\omega^{h} = \sum\_{j=0}^{I^{2}}\sum\_{k=0}^{I^{3}}w\_{0,j,k}h\_{0}(-1)h\_{j}(\xi^{2})h\_{k}(\xi^{3})\,.$$

Let *ω*b be the vector of expansion coefficients of trb *ωh*. Clearly, there exists a linear operator Nb such that

$$
\underline{a}\_{\mathsf{b}} = \mathbb{N}\_{\mathsf{b}} \; \underline{a}.
$$

The same processes can be done for other boundaries. If we collect the traces of *ωh* on all boundaries and combine their vectors of expansion coefficients and corresponding linear operators, we can eventually obtain

$$
\underline{a}\_{\underline{v}} = \mathbb{N} \underline{a},
$$

where the matrix N, like E, is sparse and only depends on the topology of the mesh and the numbering of the degrees of freedom. Furthermore, it contains only 1 as non-zero entries. An example of N can be found in [7]. Now, we can conclude that the trace space, *PI* 1*,I* 2*,I* 3 = tr P*I* 1*,I* 2*,I* 3 , is given as

$$P^{I^1, I^2, I^3} := P^{I^2, I^3}\_{-1} \cup P^{I^2, I^3}\_1 \cup P^{I^1, I^3}\_{-1} \cup P^{I^1, I^3}\_1 \cup P^{I^1, I^2}\_{-1} \cup P^{I^1, I^2}\_1,$$

where *PI* 2*,I* 3 −1 is the space spanned by *h*0*(*−1*)hj (ξ* 2*)hk(ξ* 3*)* , *PI* 2*,I* 3 1 is the space spanned by *hI* 1 *(*1*)hj (ξ* 2*)hk(ξ* 3*)* and so on. Notice that the polynomial basis functions in *h*0*(*−1*)hj (ξ* 2*)hk(ξ* 3*)* are exactly the same as those in *hI* 1 *(*1*)hj (ξ* 2*)hk (ξ* 3*)* because *h*0*(*−1*)* = *hI* 1 *(*1*)* = 1. But here we still distinguish them because they represent basis functions at different boundaries.

#### *2.2 Algebraic Dual Polynomial Spaces*

We first consider the space P*I* 1*,I* 2*,I* 3 . Let MP be the symmetric mass matrix, for example,

$$\begin{aligned} \mathbb{M}\mathbb{M}\mathbb{M}\_{l} + j(I^{1}+1) + k(I^{1}+1)(I^{2}+1), \; l+m(I^{1}+1)+n(I^{1}+1)(I^{2}+1) &:=\\ \iiint\_{\Omega\_{\text{ref}}} h\_{l}(\xi^{1})h\_{j}(\xi^{2})h\_{k}(\xi^{3})h\_{l}(\xi^{1})h\_{m}(\xi^{2})h\_{n}(\xi^{3}) \, \mathrm{d}\xi^{1} \, \mathrm{d}\xi^{2} \, \mathrm{d}\xi^{3} \, \mathrm{d}\xi \end{aligned}$$

The associated algebraic dual polynomial representations, or simply dual polynomials, are linear combinations of the polynomial basis functions, or simply primal polynomials, defined in the previous section,

$$\begin{aligned} \left[ \widehat{h\_{0,0,0}}(\xi^1, \xi^2, \xi^3), \dots, \widehat{h\_{I^1,I^2,I^3}}(\xi^1, \xi^2, \xi^3) \right] \\ := \left[ h\_0(\xi^1) h\_0(\xi^2) h\_0(\xi^3), \dots, \, \, h\_{I^1}(\xi^1) h\_{I^2}(\xi^2) h\_{I^3}(\xi^3) \right] \mathbb{M}\_{\mathcal{P}}^{-1} . \end{aligned}$$

These dual polynomials are always well-defined. This is because the primal polynomials are linearly independent. So the mass matrix MP is injective and surjective, therefore invertible. Let the finite dimensional space spanned by *h i,j,k (ξ* 1*, ξ* 2*, ξ* 3*)* be denoted by 3P*I* 1*,I* 2*,I* 3 . We say 3P*I* 1*,I* 2*,I* 3 is the algebraic dual space of the primal space P*I* 1*,I* 2*,I* 3 . Note that P*I* 1*,I* 2*,I* 3 and 3P*I* 1*,I* 2*,I* 3 actually represent the same space. The change of basis functions only leads to a different representation. Therefore, we also call the algebraic dual space a dual representation. Let M3P be the mass matrix of 3P*I* 1*,I* 2*,I* 3 , we can easily see that

$$
\tilde{\mathbb{M}}\_{\mathcal{P}} \mathbb{M}\_{\mathcal{P}} = I,\tag{5}
$$

where I is the identity matrix. Similarly, we can derive the algebraic dual space 3 L*I* 1*,I* 2*,I* 3 of the primal space L*I* 1*,I* 2*,I* 3 . Let M3L and ML be their mass matrices, we have

$$
\widetilde{\mathbb{M}}\_{\mathcal{L}} \mathbb{M}\_{\mathcal{L}} = I. \tag{6}
$$

If *ρh* ∈ L*I* 1*,I* 2*,I* 3 , *σ h*, whose vector of expansion coefficients *σ* satisfies

$$
\underline{\sigma} = \mathbb{M}\_{\mathcal{L}} \underline{\sigma},
\tag{7}
$$

will be the representation of *ρh* in the algebraic dual space 3 L*I* 1*,I* 2*,I* 3 .

To explain how the algebraic dual space of the trace space *PI* 1*,I* 2*,I* 3 is derived, we take *PI* 2*,I* 3 −1 as example. We already know that *PI* 2*,I* 3 −1 is a space spanned by primal polynomials *h*0*(*−1*)hj (ξ* 2*)hk(ξ* 3*)* . With these primal polynomials, we can compute its mass matrix, denoted by Mb. The dual polynomials are then computed by

$$\begin{aligned} \left[ \widehat{h\_{0,0,0}}(-1, \xi^2, \xi^3), \dots, \widehat{h\_{0,I^2,I^3}}(-1, \xi^2, \xi^3) \right] \\ = \left[ \widehat{h\_0(-1)h\_1(\xi^2)h\_1(\xi^3)}, \dots, \widehat{h\_0(-1)h\_{I^2}(\xi^2)h\_{I^3}(\xi^3)} \right] \mathbb{M}\_{\mathsf{b}^{-1}}^{-1} . \end{aligned}$$

The algebraic dual space3*PI* 2*,I* 3 −1 is spanned by dual polynomials *h* -0*,j,k(*−1*, ξ* 2*, ξ* 3*)* . The algebraic dual space of the trace space *PI* 1*,I* 2*,I* 3 eventually can be written as

$$
\widetilde{P}^{I^1,I^2,I^3} = \widetilde{P}^{I^2,I^3}\_{-1} \cup \widetilde{P}^{I^2,I^3}\_{1} \cup \widetilde{P}^{I^1,I^3}\_{-1} \cup \widetilde{P}^{I^1,I^3}\_{1} \cup \widetilde{P}^{I^1,I^2}\_{-1} \cup \widetilde{P}^{I^1,I^2}\_{1}.
$$

The divergence of *σ h* ∈ 3 L*I* 1*,I* 2*,I* 3 can be done with the help of the boundary value *σ*ˆ *h* ∈ 3*PI* 1*,I* 2*,I* 3 . With vector proxies, it can be written as

$$
\underline{\text{div}\,\underline{\sigma}^h} = \mathbb{N}^\mathsf{T}\,\underline{\hat{\underline{\sigma}}^h} - \mathbb{E}^\mathsf{T}\,\underline{\underline{\sigma}}^h.\tag{8}
$$

A detailed introduction of algebraic dual polynomial spaces is given in [9].

#### *2.3 Function Spaces in Curvilinear Domains*

So far, all polynomial spaces are defined only in the reference domain ref|*ξ* 1*,ξ* 2*,ξ* 3 = [−1*,* 1] 3. Consider an arbitrary domain and a C1 diffeomorphism : ref|*ξ* 1*,ξ* 2*,ξ* 3 → |*x*1*,x*2*,x*3 . In , the primal polynomials change. Therefore, the mass matrices will also change. But the process of constructing dual polynomials does not change. And as we mentioned before, the metric-independent incidence matrix E and the matrix N remain the same. The way of converting polynomials in Cartesian domain into those in curvilinear domains follows the general coordinate transformation process, for example, see [16].

From now on, notations mentioned in this section not only refer to the reference domain ref, but also refer to the physical domain .

#### **3 Weak Formulations**

#### *3.1 Discrete Neumann Problem*

With integration by parts, we can derive the weak formulation of the Neumann problem, (1), written as: For given *σ*ˆ ∈ *H* −1*/*2*(∂)*, find *ω* ∈ *H*1*()* such that

$$\left(\text{grad}\,\boldsymbol{\omega},\,\,\text{grad}\,\tilde{\boldsymbol{\omega}}\right)\_{L^{2}} + (\boldsymbol{\omega},\,\,\tilde{\boldsymbol{\omega}})\_{L^{2}} = \left<\text{tr}\,\,\tilde{\boldsymbol{\omega}},\,\,\hat{\boldsymbol{\sigma}}\right>,\,\,\,\forall\tilde{\boldsymbol{\omega}}\in H^{1}(\Omega).\tag{9}$$

Note that on the right hand side, we use ·*,*· to represent the duality pairing between tr *ω*¯ ∈ *H*1*/*2*(∂)* and *σ*ˆ ∈ *H* −1*/*2*(∂)*. We use finite dimensional space P*I* 1*,I* 2*,I* 3 to approximate the space *H*1*()* and use the algebraic dual trace space 3*PI* 1*,I* 2*,I* 3 to approximate the space *H* −1*/*2*(∂)*. Then we obtain

$$\left(\text{grad}\,\boldsymbol{\omega}^{h},\,\text{grad}\,\bar{\boldsymbol{\omega}}^{h}\right)\_{L^{2}}=\bar{\boldsymbol{\omega}}^{h,\,\text{T}}\,\mathbb{E}^{\mathsf{T}}\mathbb{M}\_{\mathcal{L}}\mathbb{E}\,\underline{\boldsymbol{\omega}}^{h},$$

$$\left(\boldsymbol{\omega}^{h},\,\,\bar{\boldsymbol{\omega}}^{h}\right)\_{L^{2}}=\bar{\boldsymbol{\omega}}^{h,\,\mathsf{T}}\,\mathbb{M}\boldsymbol{\rho}\,\,\underline{\boldsymbol{\omega}}^{h},$$

and

$$\int\_{\partial\Omega} \text{tr}\,\bar{\boldsymbol{\phi}}^h \,\hat{\boldsymbol{\sigma}}^h \,\mathrm{d}\Gamma = \underline{\underline{\boldsymbol{\phi}}}^{h,\mathsf{T}} \,\mathbb{N}^{\mathsf{T}} \,\underline{\hat{\underline{\mathbf{o}}}^h}.$$

which eventually leads to the discrete formulation of (9),

$$\mathbb{E}^{\mathsf{T}}\mathbb{M}\_{\mathcal{L}}\mathbb{E}\,\underline{\boldsymbol{\omega}}^{h} + \mathbb{M}\_{\mathcal{P}}\,\underline{\boldsymbol{\alpha}}^{h} = \mathbb{N}^{\mathsf{T}}\,\underline{\hat{\mathsf{o}}}^{h}.\tag{10}$$

#### *3.2 Discrete Dirichlet Problem*

For the Dirichlet problem, (2), the weak formulation is given as: For given *σ*ˆ ∈ *H* −1*/*2*(∂)*, find *σ* ∈ *H (*div*, )*, tr *σ* = ˆ*σ* such that

$$(\operatorname{div}\sigma,\,\operatorname{div}\bar{\sigma})\_{L^2} + (\sigma,\,\bar{\sigma})\_{L^2} = 0,\quad \forall \bar{\sigma} \in H\_0(\operatorname{div},\Omega). \tag{11}$$

We use algebraic dual space 3 L*I* 1*,I* 2*,I* 3 to approximate *H (*div*, )*. With *σ*ˆ *h* ∈ 3*PI* 1*,I* 2*,I* 3 given and (8), we obtain

$$\left(\operatorname{div}\boldsymbol{\sigma}^{h},\operatorname{div}\bar{\boldsymbol{\sigma}}^{h}\right)\_{L^{2}}=-\bar{\underline{\sigma}}^{h,\mathsf{T}}\operatorname{\mathbb{E}\tilde{\mathbb{M}}\_{\mathcal{P}}\left(\mathbb{N}^{\mathsf{T}}\,\underline{\hat{\underline{\sigma}}}^{h}-\mathbb{E}^{\mathsf{T}}\,\underline{\underline{\sigma}}^{h}\right),$$

and

$$\left(\boldsymbol{\sigma}^{h},\,\,\bar{\boldsymbol{\sigma}}^{h}\right)\_{L^{2}}=\underline{\bar{\boldsymbol{\sigma}}^{h,\mathsf{T}}\,\,\widetilde{\mathbb{M}}\_{\mathcal{L}}\,\underline{\underline{\sigma}}^{h}\,\,\cdot$$

Therefore, the discrete formulation of (11) is written as

$$\mathbb{E}\widetilde{\mathbb{M}}\_{\mathcal{P}}\mathbb{E}^{\mathsf{T}}\underline{\mathbf{o}}^{h} + \widetilde{\mathbb{M}}\_{\mathcal{L}}\underline{\mathbf{o}}^{h} = \mathbb{E}\widetilde{\mathbb{M}}\_{\mathcal{P}}\mathbb{N}^{\mathsf{T}}\underline{\hat{\mathbf{o}}}^{h}.\tag{12}$$

#### *3.3 Equivalence Between Discrete Formulations*

Now it is time to check if the equivalence between (1) and (2) holds at the discrete level. In other words, it is time to check if the statement that *ωh* solves (10) if and only if *σ h* = grad *ωh* solves (12) is correct.

From (4) and (7), we know that *σ h*,

$$
\underline{\underline{\sigma}}^h = \mathbb{M} \mathcal{L} \, \underline{\mathbb{E}} \, \underline{\underline{\omega}}^h,\tag{13}
$$

is the vector representation of grad *ωh* in the dual space. If we insert (13) into (12), we obtain

$$\mathbb{E}\widetilde{\mathbb{M}}\mathbb{M}^{\mathbb{E}}\mathbb{M}^{\mathbb{T}}\mathbb{M}^{\mathbb{T}}\mathbb{E}\,\underline{\boldsymbol{\omega}}^{h} + \widetilde{\mathbb{M}}\mathbb{M}^{\mathbb{T}}\mathbb{L}^{\mathbb{E}}\underline{\mathbb{E}}\,\underline{\boldsymbol{\omega}}^{h} = \mathbb{E}\widetilde{\mathbb{M}}\mathbb{M}^{\mathbb{N}}\,\underline{\hat{\mathbf{o}}}^{h}.\tag{14}$$

From (10), we know that

$$\mathbb{E}^{\mathsf{T}}\mathbb{M}\_{\mathbb{C}}\mathbb{E}\,\underline{\boldsymbol{\omega}}^{h} = -\mathbb{M}\boldsymbol{\rho}\,\underline{\boldsymbol{\omega}}^{h} + \mathbb{N}^{\mathsf{T}}\,\underline{\hat{\underline{\sigma}}}^{h}.\tag{15}$$

By inserting (15) into (14), we get

$$\mathbb{E}\widetilde{\mathbb{M}}\mathfrak{\rho}\left(-\mathbb{M}\underline{\rho}\,\,\underline{\omega}^{h}+\mathbb{N}^{\mathsf{T}}\,\underline{\hat{\mathfrak{a}}}^{h}\right)+\widetilde{\mathbb{M}}\_{\mathcal{L}}\mathbb{M}\_{\mathcal{L}}\mathbb{E}\,\underline{\omega}^{h}=\mathbb{E}\widetilde{\mathbb{M}}\mathbb{M}^{\mathsf{T}}\,\underline{\hat{\mathfrak{a}}}^{h}.\tag{16}$$

From (5) and (6), we know that (16) holds, which proves the equivalence.

If the equivalence holds, relation *ωh H*1*()* = *σ h H (*div*,)* should also be satisfied. To prove this, we have

 *σ h* 2 *H (*div*,)* (8) = *σ h,*T M3L *σ h* + NT *σ*ˆ *h* − ET*σ h* T M3P NT *σ*ˆ *h* − ET*σ h* (13) = MLE *ωh* T M3L MLE *ωh* + NT *σ*ˆ *h* − ET MLE *ωh* T M3P NT *σ*ˆ *h* − ET MLE *ωh* (10) = *ωh,*T ETMLE *ωh* + *ωh,*T MPM3PMP *ωh* = *ωh* 2 *H*1*() ,*

where we constantly use (5) and (6) and the fact that mass matrices are symmetric.

#### **4 Numerical Test**

Consider the mapping which maps the Cartesian reference domain ref|*ξ* 1*,ξ* 2*,ξ* 3 := [−1*,* 1] 3 into the physical domain |*x*1*,x*2*,x*3 = [0*,* 1] 3 by

$$\alpha^{i} = \frac{1}{2} + \frac{1}{2} \left( \xi^{i} + c \prod\_{j} \sin(\pi \xi^{j}) \right), \ i = 1, 2, 3.$$

When the deformation coefficient *c* = 0, the domain is Cartesian. Otherwise the domain is curvilinear, meaning that a curvilinear coordinate system parametrizes . Examples of such curvilinear domains in R2 are shown in Fig. 1.

A manufactured solution of the Neumann problem, (1), is

$$a\_{\text{exact}} = e^{\boldsymbol{\chi}^1} + e^{\boldsymbol{\chi^2}} + e^{\boldsymbol{\chi^3}}.$$

Clearly, *σ*exact = grad *ω*exact = *ex*1 *, ex*2 *, ex*3 T solves the Dirichlet problem, (2).

In the domains of different deformation coefficient *c*, with the boundary condition *σ*ˆ = tr *σ*exact imposed, we solve the discrete formulations (10) and (12) using Gauss–Lobatto–Legendre (GLL) polynomial spaces of degree *I* 1 = *I* 2 = *I* 3 = *N*.

The results of the *L*2-error of *σ h* − grad *ωh* are shown in Fig. 2 (Left) where we can see that the relation *σ h* = grad *ωh* is preserved up to the machine precision. With the growth of the polynomial degree, the error increases slowly because of the accumulation of the machine error as the amount of degrees of freedom grows significantly.

In Table 1, the results of the *H*1-norm of *ωh* and *H (*div*)*-norm of *σ h* are presented. It is shown that the relation *ωh H*1*()* = *σ h H (*div*,)* holds for all polynomial degrees irrespective of whether we use the Cartesian domain, *c* = 0, or

**Fig. 1** Curvilinear domains for *c* = 0*.*15 (Left) and *c* = 0*.*3 (Right) in R2. The gray lines illustrate the coordinate lines

**Fig. 2** The *L*2-error of *σ h* − grad *ωh* (Left) and the *p*-convergence of the *H*1-error of *ωh* (Right) for *N* = 2*,* 4*,*··· *,* 20 and *c* = 0*,* 0*.*15*,* 0*.*3

**Table 1** The *H*1-norm of *ωh* and *H (*div*)*-norm of *σ h* for polynomial degree *N* = 2*,* 4*,*··· *,* 20 and deformation coefficient *c* = 0*,* 0*.*15*,* 0*.*3

curvilinear domains, *c* = 0*.*15*,* 0*.*3. It is also seen that the results always converge to the analytical value *ω*exact*H*1 = *σ h H (*div*)* = 6*.*0730653668.The *p*-convergence for the *H*1-error of *ωh*, therefore also for the *H (*div*)*-error of *σ h*, is shown in Fig. 2 (Right), which shows the exponential convergence of the method.

#### **5 Conclusions**

By constructing and using primal polynomial spaces and their algebraic dual representations both in the domain and on the boundary, we successfully preserve the equivalence of the div-grad Neumann problem and the grad-div Dirichlet problem at the discrete level in 3-dimensional curvilinear domains. This suggests the further usage of these spaces to structure-preserving methods and hybrid methods.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **A Conservative Hybrid Method for Darcy Flow**

**Varun Jain, Joël Fisser, Artur Palha, and Marc Gerritsma**

### **1 Introduction**

Hybrid formulations [1, 3, 10] are classical domain decomposition methods which reduce the problem of solving one global system to many small local systems. The local systems can then be efficiently solved independently of each other in parallel.

In this work we present a hybrid mimetic spectral element formulation to solve Darcy flow. We follow [8] which render the constraints on divergence of mass flux, the pressure gradient and the inter-element continuity metric free. The resulting system is extremely sparse and shows a reduced growth in condition number as compared to a non-hybrid system.

This document is structured as follows: In Sect. 2 we define the weak formulation for Darcy flow. The basis functions are introduced in Sect. 3. The evaluation of weighted inner product and duality pairings are discussed in Sect. 4. In Sect. 5 we discuss the formulation of discrete algebraic system. In Sect. 6 we present results for a test case taken from [7].

V. Jain (-) · J. Fisser · A. Palha · M. Gerritsma

Faculty of Aerospace Engineering, TU Delft, Delft, The Netherlands e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
S. J. Sherwin et al. (eds.), *Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018*, Lecture Notes in Computational Science and Engineering 134, https://doi.org/10.1007/978-3-030-39647-3\_16

#### **2 Darcy Flow Formulation**

For ∈ R*d* , where *d* is the dimension of the domain, the governing equations for Darcy flow, are given by,

$$\begin{cases} \mathfrak{u} + \mathbb{A} \, \nabla p = 0 \\ \nabla \cdot \mathfrak{u} \, \qquad = f \end{cases} \quad \text{in } \Omega \quad \text{and} \quad \begin{cases} \partial \Omega &= \Gamma\_D \cup \Gamma\_N \\ p &= \hat{p} \\ \mathfrak{u} \cdot \mathfrak{n} = \hat{\mathfrak{u}}\_n & \text{on } \Gamma\_N \end{cases} ,$$

where, *u* is the velocity, *p* is the pressure, *f* the prescribed RHS term, A is a *d* × *d* symmetric positive definite matrix, *p*ˆ and *u*ˆ *n* are the prescribed pressure and flux boundary conditions, respectively.

#### *2.1 Notations*

For *f, g* ∈ *L*2 *()*, *f, g* denotes the usual *L*2-inner product.

For vector-valued functions in *L*2 we define the weighted inner product by,

$$(\mathfrak{u}, \mathfrak{v})\_{\mathbb{A}^{-1}, \Omega} = \int\_{\Omega} \left( \mathfrak{u}, \mathbb{A}^{-1} \mathfrak{v} \right) d\Omega \,, \tag{1}$$

where *(*· *,*·*)* denotes the pointwise inner product.

Duality pairing, denoted by ·*,*·, is the outcome of a linear functional on *L*2 *()* acting on elements from *L*2 *()*.

Let *K* be a disjoint partitioning of with total number of elements *K*, and *Ki* is any element in *K*, such that, *Ki* ∈ *K*. We define the following broken Sobolev spaces [2], *H* div; *K* = A *i H* div; *Ki* , and *H*1*/*2 *(∂K)* = A *i H*1*/*2 *(∂Ki)*.

#### *2.2 Weak Formulation*

The Lagrange functional for Darcy flow is defined as,

$$\begin{split} \mathcal{L}\left(\mathfrak{u},p,\lambda;f\right) &= \frac{1}{2} \int\_{\partial\Omega\_{K}} \mathfrak{u}^{T} \mathbb{A}^{-1} \mathfrak{u} \,\mathrm{d}\Omega\_{K} - \int\_{\Omega\_{K}} p \left(\nabla \cdot \mathfrak{u} - f\right) \mathrm{d}\Omega\_{K} \\ &+ \int\_{\partial\Omega\_{K} \backslash \Gamma\_{D}} \lambda \left(\mathfrak{u} \cdot \mathfrak{n}\right) \,\mathrm{d}\Gamma + \int\_{\Gamma\_{D}} \hat{p} \left(\mathfrak{u} \cdot \mathfrak{n}\right) \,\mathrm{d}\Gamma - \int\_{\Gamma\_{N}} \lambda \left(\hat{\mathfrak{u}}\_{\mathfrak{n}}\right) \,\mathrm{d}\Gamma \end{split}$$

The variational problem is then given by: For given *f* ∈ *L*2 *(K)*, *p*ˆ ∈ *H*1*/*2*(ΓD)* and *u*ˆ *n* ∈ *H* −1*/*2*(ΓN )* find *u* ∈ *H (*div; *K)*, *p* ∈ *L*2 *(K)*, *λ* ∈ *H* 1 2 *(∂K)*, such that,

$$\begin{cases} \left. \left( \boldsymbol{\mathfrak{v}}, \boldsymbol{\mathfrak{u}} \right) \right|\_{\boldsymbol{\Lambda}^{-1}, \boldsymbol{\Omega}\_{K}} - \left< \nabla \cdot \boldsymbol{\mathfrak{v}}, p \right>\_{\boldsymbol{\Omega}\_{K}} + \left< \left( \boldsymbol{\mathfrak{v}} \cdot \boldsymbol{\mathfrak{n}} \right), \lambda \right>\_{\partial \Omega\_{K} \backslash \Gamma\_{D}} = - \left< \boldsymbol{\mathfrak{v}} \cdot \boldsymbol{\mathfrak{n}}, \hat{p} \right>\_{\Gamma\_{D}} & \forall \, \mathfrak{v} \in H(\operatorname{div}; \Omega\_{K}) \\\\ - \left< q, \nabla \cdot \boldsymbol{\mathfrak{n}} \right>\_{\Omega\_{K}} & = - \left< q, \boldsymbol{f} \right>\_{\Omega\_{K}} & \forall \, \mathfrak{q} \in L^{2}(\Omega\_{K}) \\\\ \left< \mu, \left( \boldsymbol{\mathfrak{n}} \cdot \boldsymbol{\mathfrak{n}} \right) \right>\_{\partial \Omega\_{K} \backslash \Gamma\_{D}} & = \left< \mu, \left< \boldsymbol{\mathfrak{n}}, \boldsymbol{\mathfrak{n}} \right>\_{\Gamma\_{N}} & \forall \, \mu \in H^{\frac{1}{2}}(\partial \Omega\_{K}) \end{cases} \end{cases} \tag{2}$$

#### **3 Basis Functions**

#### *3.1 Primal and Dual Nodal Degrees of Freedom*

Let *ξj* , *j* = 0*,* 1*,...,N*, be the *N* + 1 Gauss–Lobatto–Legendre (GLL) points in *I* ∈ −1*,* 1 . The Lagrange polynomials *hi(ξ )* through *ξj* , of degree *N*, given by,

$$h\_l\left(\xi\right) = \frac{\left(\xi^2 - 1\right)L'\_N\left(\xi\right)}{N\left(N + 1\right)L\_N\left(\xi\_l\right)\left(\xi - \xi\_l\right)},$$

form the 1D primal nodal polynomials which satisfy, *hi(ξj )* = *δij* .

Let *ah* and *bh* be two polynomials expanded in terms of *hi ξ* . The *L*2—inner product is then given by,

$$\left(a^{h}, b^{h}\right)\_{I} = \mathbf{a}^{T} \mathbb{M}^{(0)} \mathbf{b} \,, \quad \text{where} \quad \mathbb{M}^{(0)}\_{i,j} = \int\_{-1}^{1} h\_{l}(\xi) \, h\_{j}(\xi) \, d\xi \,,$$

and, **a** = [a0 a1 *...* a*N* ] and **b** = b0 b1 *...* b*N* are the nodal degrees of freedom. We define the algebraic *dual* degrees of freedom, **a**, such that the duality pairing is simply the vector dot product between primal and dual degrees of freedom,

$$\begin{aligned} \text{for dot product between primal and dual degrees of } \mathbf{f} \\\\ \left\langle a^h, b^h \right\rangle\_I = \widetilde{\mathbf{a}}^T \mathbf{b} := \mathbf{a}^T \mathbb{M}^{(0)} \mathbf{b} \quad \Rightarrow \quad \widetilde{\mathbf{a}} = \mathbb{M}^{(0)} \mathbf{a} \end{aligned}$$

Thus, the dual degrees of freedom are linear functionals of primal degrees of freedom.

#### *3.2 Primal and Dual Edge Degrees of Freedom*

The edge polynomials, for the *N* edges between *N* + 1 GLL points *ξj*−1*, ξj* , of polynomial degree *N* − 1, are defined as [4],

$$e\_j(\xi) = -\sum\_{k=0}^{j-1} \frac{\mathrm{d}h\_k}{\mathrm{d}\xi}(\xi) \ , \quad \text{such that} \quad \int\_{\xi\_{j-1}}^{\xi\_j} e\_i(\xi) = \delta\_{ij} \ .$$

Let *ph* and *qh* be two polynomials expanded in edge basis functions. The inner product in *L*2 space is given by,

$$\left(p^h, q^h\right)\_I = \mathbf{p}^T \mathbb{M}^{(1)} \mathbf{q} \,, \quad \text{where} \quad \mathbb{M}\_{l,j}^{(1)} = \int\_{-1}^1 e\_l(\xi) \, e\_j(\xi) \, d\xi \,, \quad \text{as} \quad \text{for} \quad l \in \mathbb{N}$$

and, **p** = p1 p2 *...* p*N* and **q** = q1 q2 *...* q*N* are the edge degrees of freedom. As before, we define the *dual* degrees of freedom such that,

 *ph, qh I* = **p***T* **q** := **p***T* M*(*1*)* **q** ⇒ **p** = M*(*1*)* **p** *.*

A similar construction can be used for dual degrees of freedom in higher dimensions. For construction of the dual degrees of freedom in 2D see [8] and for 3D see [9].

#### *3.3 Differentiation of Nodal Polynomial Representation*

Let *ah ξ* be expanded in Lagrange polynomials, then

$$\frac{d}{d\xi}a^h\left(\xi\right) = \frac{d}{d\xi}\sum\_{l=0}^{N} \mathbf{a}\_l h\_l\left(\xi\right) = \sum\_{l=1}^{N} \left(\mathbf{a}\_l - \mathbf{a}\_{l-1}\right) e\_l\left(\xi\right) \,. \tag{3}$$

Therefore, taking the derivative of a polynomial involves two steps: First, take the difference of degrees of freedom; second, change of basis from nodal to edge [4].

#### **4 Discrete Inner Product and Duality Pairing**

For 2D domains, the higher dimensional primal basis are constructed using the tensor product of the 1D basis.

For the weak formulation (2) we expand the velocity *uh* in primal edge basis as,

$$\mathfrak{u}^{h}\left(\boldsymbol{\xi},\boldsymbol{\eta}\right) = \sum\_{l=0}^{N} \sum\_{j=1}^{N} \mathfrak{u}\_{\mathfrak{X}l,j} \, h\_{l}\left(\boldsymbol{\xi}\right) \, e\_{j}\left(\boldsymbol{\eta}\right) \, \hat{\mathfrak{z}} + \sum\_{l=1}^{N} \sum\_{j=0}^{N} \mathfrak{u}\_{\mathfrak{Y}l,j} \, e\_{l}\left(\boldsymbol{\xi}\right) \, h\_{j}\left(\boldsymbol{\eta}\right) \, \hat{\mathfrak{z}}\,,\tag{4}$$

where u*x i,j* denotes the flux, ; *u* · *n*, over the vertical edges and u*y i,j* the flux over the horizontal edges, see Fig. 1.

#### *4.1 Weighted Inner Product*

Using (1) and the expansions in (4), the weighted inner product is evaluated as,

$$\left(\boldsymbol{v}^{h},\boldsymbol{u}^{h}\right)\_{\mathbb{A}^{-1},\mathfrak{Q}\_{K}} = \sum\_{K\_{\ell}} \mathbf{v}\_{K\_{\ell}}^{T} \, \mathbb{M}\_{\mathbb{A}^{-1},K\_{\ell}}^{(1)} \, \mathbf{u}\_{K\_{\ell}} \,, $$

**Fig. 1** Discretized domain for *K* = 3 × 3, *N* = 3. The blue dots represent the pressure boundary condition *p*ˆ, and the blue edges represent the velocity boundary condition *u*ˆ *n*

where, **u***Ki* are the degrees of freedom in element *Ki*, and

$$\mathbb{M}^{(1)}\_{\mathbb{A}^{-1}, K\_l} = \int\_{K\_l} \begin{pmatrix} h\_l(\xi) \ e\_j(\eta) \\ e\_l(\xi) \ h\_j(\eta) \end{pmatrix}^T \mathbb{A}^{-1} \begin{pmatrix} \xi, \eta \end{pmatrix} \begin{pmatrix} h\_l(\xi) \ e\_j(\eta) \\ e\_l(\xi) \ h\_j(\eta) \end{pmatrix} \, \mathrm{d}K\_l \ .$$

For mapping of elements please refer to [6].

#### *4.2 Divergence of Velocity*

Divergence of velocity, ∇ · *uh*, is evaluated using (3), but now for 2D,

$$\begin{split} \nabla \cdot \mathbf{u}^{h} &= \frac{\partial}{\partial \mathbf{x}} \sum\_{l=0}^{N} \sum\_{j=1}^{N} \mathbf{u}\_{\mathbf{x}l,j} \, h\_{l}(\xi) e\_{j}(\eta) + \frac{\partial}{\partial \mathbf{y}} \sum\_{l=1}^{N} \sum\_{j=0}^{N} \mathbf{u}\_{\mathbf{y}\_{i,j}} e\_{l}(\xi) h\_{j}(\eta) \\ &= \sum\_{l,j=1}^{N} \left( \mathbf{u}\_{\mathbf{x}l,j} - \mathbf{u}\_{\mathbf{x}l-1,j} + \mathbf{u}\_{\mathbf{y}\_{i,j}} - \mathbf{u}\_{\mathbf{y}\_{i,j-1}} \right) e\_{l}(\xi) \, e\_{j}(\eta) \end{split} \tag{5}$$

For pressure we will use dual degrees of freedom. Therefore the weak constraint on divergence of velocity is a duality pairing evaluated as, **-**

$$\left\langle q^h, \nabla \cdot \mathbf{u}^h \right\rangle\_{\Omega\_K} = \sum\_{K\_l} \widetilde{\mathbf{q}}\_{K\_l}^T \, \mathbb{E}^{2,1} \, \mathbf{u}\_{K\_l} \,, \, $$

where E2*,*1 represents the discrete divergence operator. It is an incidence matrix that is metric-free and topological, and remains the same for each element in *K*. For an extensive discussion on the incidence matrix, see for instance [6]. For an element of degree *N* = 3,

E2*,*1 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ −1 0 0 1 0 0 0 0 0 000 −1 0 0 1 0 0 0 0 0 000 000 −1 0 0 1 0 0 000 0 −1 0 0 1 0 0 0 0 000 000000 −1 0 0 100 0 0 −1 0 0 1 0 0 0 000 0 −1 0 0 1 0 0 0 0 000 0 0 0 −1 0 0 1 0 0 000 0000 −1 0 0 1 0 000 0 0 0 0 −1 0 0 1 0<PHONE_NUMBER> −1 0 010 0 0 0 0 0 −1 0 0 1 000 0 0 −1 0 0 1 0 0 0 000 0 0 0 0 0 0 −1 0 0 100 00000 −1 0 0 1 000 0 0 0 0 0 0 0 −1 0<PHONE_NUMBER>0 −1001 0 0 0 0 0 0 0 0 −1001 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ *.*

#### *4.3 Connectivity Matrix*

The connectivity matrix ensures continuity of the velocity flux across the elements. *λ* is the interface variable defined between the elements, shown as red dots in Fig. 1. *λ* acts as Lagrange multiplier that imposes the continuity constraint given by,

$$\begin{aligned} \text{nge multiplier that imposes the continuity constraint:}\\ \left\langle \boldsymbol{\mu}^{h}, \boldsymbol{\mu}^{h} \cdot \mathbf{n} \right\rangle\_{\partial \Omega\_{K} \backslash \Gamma\_{D}} = \sum\_{K} \widetilde{\boldsymbol{\mu}}\_{K\_{i}}^{T} \boxtimes \mathbf{u}\_{K\_{i}} = \widetilde{\boldsymbol{\mu}}^{T} \boxtimes\_{\mathbb{N}} \mathbf{u}\_{-,h} \end{aligned}$$

where N is the discrete trace operator. It is a sparse matrix that consists of 1, −1 and 0 only. For construction of N please refer to [5]. EN is the assembled N for all elements. For, *K* = 2 × 2, *N* = 2, EN is shown in (6). The matrix size of EN is 8 × 64, but it has only 16 non-zero entities. It is an extremely sparse matrix that is metric free and the location of ± valued entries depend only on the connection between different elements.

EN = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0000 1 0 00000 0 0000 -1 0 000000000 0 000000000 0 0 0 000000000 0 00000 0 00000000 00000 1 00000 0 00000 -1 000000000 0 000000000 0 0 0 000000000 0 00000 0 00000000 00000 0 00000 0 00000 0 000000000 0 00000000 1 0 0 0 00000000 -1 0 00000 0 00000000 00000 0 00000 0 00000 0 000000000 0 000000000 1 0 0 000000000 -1 00000 0 00000000 00000 0 0000 1 0 00000 0 000000000 0 000000000 0 -1 0 000000000 0 00000 0 00000000 00000 0 00000 1 00000 0 000000000 0 000000000 0 0 -1 000000000 0 00000 0 00000000 00000 0 00000 0 00000 0 00000000 1 0 000000000 0 0 0 000000000 0 0000 -1 0 00000000 00000 0 00000 0 00000 0 000000000 1 000000000 0 0 0 000000000 0 00000 -1 00000000 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ *.* (6)

#### **5 Discrete Formulation**

Using the weighted inner product and duality pairings discussed in Sect. 4, we can write the discrete form of weak formulation in (2) as,

$$
\begin{bmatrix}
\mathbb{B} & \mathbb{E}\_{\mathbb{N}}^{\mathcal{T}} \\
\mathbb{E}\_{\mathbb{N}} & 0
\end{bmatrix}
\begin{bmatrix}
\mathbf{x} \\
\mathbf{x}
\end{bmatrix} =
\begin{bmatrix}
\mathbf{F} \\
\mathbf{0}
\end{bmatrix}.
\tag{7}
$$

where, B is an invertible block diagonal matrix given by,

$$
\mathbb{B} = \begin{bmatrix}
\mathbb{M}\_{\mathbb{A}^{-1}, K\_1}^{(1)} & \mathbb{E}^{2, 1} & & & \\
& \mathbb{E}^{2, 1} & \mathbb{0} & & & \\
& & \mathbb{M}\_{\mathbb{A}^{-1}, K\_2}^{(1)} \mathbb{E}^{2, 1} & & & \\
& & \mathbb{E}^{2, 1} & \mathbb{0} & & \\
& & & \ddots & \ddots & \\
& & & & \ddots & \ddots \\
& & & & & \mathbb{M}\_{\mathbb{A}^{-1}, K\_K}^{(1)} \mathbb{E}^{2, 1} \\
& & & & \mathbb{E}^{2, 1} & \mathbb{0}
\end{bmatrix}, \tag{8}
$$

EN is as given in (6), **X** = < *i* ⎡ ⎣ **u p** ⎤ ⎦ *Ki* , and **F** = < *i* ⎡ ⎣ **p**ˆ **f** ⎤ ⎦ *Ki* , where **f** are the expansion coefficients of *f h x,y* = <*N i,j* f*ij ei (x) ej y* .

In (8), the mass matrix M*(*1*)* A−1*,Ki* is the only dense matrix and also the only matrix that changes with each local element, *Ki*. EN is a sparse incidence matrix for the global system and E2*,*1 is a sparse incidence matrix for the local systems that remains the same for each element.

Using the Schur complement method, the global system (7) can be reduced to solve for *λ*, [1],

$$\lambda = \left(\mathbb{E}\_{\mathbb{N}} \mathbb{E}^{-1} \, \mathbb{E}\_{\mathbb{N}} \, ^T\right)^{-1} \cdot \left(\mathbb{E}\_{\mathbb{N}} \mathbb{E}^{-1} \, \mathbf{F}\right) \, . \tag{9}$$

To evaluate *λ* in (9) we need B−1 that can be calculated efficiently by taking inverse of each block of B separately. This part is trivially parallelized. Once the *λ* is determined the solution in each element, *Ki*, can be evaluated independent of each other.

The system (9) solves for interface degrees of freedom between the elements and will always be smaller than the full global system. For a comparison of the size of *λ* system with full system see Table 1 (for 2D), and Table 2 (for 3D). On the left of Tables 1 and 2 we see that, for constant *K*, increasing the order of polynomial basis the growth in size of *λ* system is less than the growth in size of full system. Thus, hybrid formulations are beneficial for high order methods where local degrees of freedom of an element are much higher than interface degrees of freedom.

**Table 1** For 2D

Left: Number of total unknowns as a function of *N*, for *K* = 3 × 3. Right: Number of total unknowns as a function of the number of elements *K*, for *N* = 3

Left: Number of total unknowns as a function of *N*, for *K* = 3 × 3 × 3. Right: Number of total unknowns as a function of the number of elements *K*, for *N* = 3

On the right of Tables 1 and 2 we see that, for constant *N*, the *λ* system is smaller than the full system, although the growth ratio of the size of *λ* and full systems do not change significantly.

#### **6 Results**

In this section we present the results for a test problem from [7] by solving system (7). The domain of the test problem is, ∈ 0*,* 1 2 . The RHS term is defined as,

$$f\_{e\boldsymbol{\alpha}} = \nabla \cdot (-\mathbb{A}\nabla p\_{e\boldsymbol{\alpha}})\ , \quad \text{where } \boldsymbol{\alpha}$$

$$\begin{array}{ll} \mathsf{A} &= \frac{1}{\chi^2 + \chi^2 + a} \begin{pmatrix} 10^{-3} \mathsf{x}^2 + \mathsf{y}^2 + a & \left(10^{-3} - 1\right) \mathsf{x} \mathsf{y} \\\\ \left(10^{-3} - 1\right) \mathsf{x} \mathsf{y} & \mathsf{x}^2 + 10^{-3} \mathsf{y}^2 + a \end{pmatrix}; & a = 0.1 \\\ p\_{e\boldsymbol{x}} &= \sin\left(2\pi\boldsymbol{x}\right) \sin\left(2\pi\boldsymbol{\chi}\right) \end{array}$$

and Dirichlet boundary conditions are imposed along the entire boundary,D = *∂* and N = ∅. We solve this problem on an orthogonal and a curved mesh, see Fig. 2.

The same problem was earlier addressed in [6], but for a method with continuous elements and *primal* basis functions only. For the configuration *K* = 3 × 3, *N* = 6, we compare the sparsity structure of the two approaches in Fig. 3. On the left we see

*,*

**Fig. 2** Mesh configuration: *K* = 3 × 3, *N* = 6. Left: orthogonal. Right: curved

**Fig. 3** Sparsity plots *K* = 3×3, *N* = 6. Left: hybrid elements method. Right: continuous element method

the hybrid formulation, and on the right we see the continuous elements formulation [6]. The number of non zero entries are almost half in the hybrid formulation, 66*,*384, as compared to the continuous element formulation, 117*,*504. Here, the sparsity is due to use of algebraic dual degrees of freedom and is not because of hybridization of the scheme.

In Fig. 4, on the left we compare the growth in condition number, for the *λ* system (9) with full continuous element system, for *N* = 7 on the curved mesh, with increasing number of elements, *K*. We observe similar growth rates for hybrid and continuous formulation, however the condition number for continuous elements formulation is almost O 102 higher. On the right we see the growth in condition number with increasing polynomial degree for *K* = 9 × 9 on the curved mesh. A reduced growth rate in condition number for hybrid formulation is observed. Thus hybrid formulations are beneficial for high order methods.

**Fig. 4** Growth in condition number for hybrid elements in dark line, and continuous elements in dotted line. Left: *h*-refinement; Right: *N*-refinement. 'c' refers to the curved mesh

**Fig. 5** *L*2-error in divergence of velocity: Left: *h*-refinement; Right: *N*-refinement. 'o' refers to the orthogonal mesh and 'c' to the curved mesh

In Fig. 5 we show the *L*2-error for ∇ ·*uh* −*f h*. On the left side as a function of element size, *h* = 1*/* √*K*, and on the right side as a function of polynomial degree of the basis functions. In both cases the maximum error observed is of O 10−12 .

In Fig. 6, on the top two figures we show the error in the *H* div; norm for the velocity; and at the bottom two figures we show the error in *L*2 *()* norm for the pressure. On the left we have *h*-convergence plots, and on the right we have *N*-convergence plots. In all the figures, for the same number of elements, *K*, and polynomial degree, *N*, the error is higher for the curved mesh.

On the left we see that the error decreases with the element size. The slope of error rate of convergence is *N*, which is optimal for both curved and orthogonal meshes. On the right we see exponential convergence of the error with increasing polynomial degree of basis for both orthogonal and curved meshes.

**Fig. 6** Top row: error in *H* div; norm for velocity; Bottom row: *L*2-error in pressure. Left: *h*-refinement; Right: *N*-refinement. 'o' refers to the orthogonal mesh and 'c' to the curved mesh

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **High-Order Mesh Generation Based on Optimal Affine Combinations of Nodal Positions**

**Mike Stees and Suzanne M. Shontz**

#### **1 Introduction**

The advantage of high-order numerical methods for solving partial differential equations is their higher degree of accuracy compared to low-order numerical methods. A major hurdle in the usage of these methods in the presence of complex geometries is the absence of robust high-order mesh generation methods [23]. In other words, these methods need a high-order mesh that accurately captures the features of the geometry to achieve their full potential [1, 10].

The typical approach for generating high-order meshes is to transform a coarse linear mesh [2–6, 9, 11, 12, 14, 16, 17, 19–22, 24]. At a high-level, these transformations usually consist of the following three steps: (1) the low-order mesh is enriched with additional nodes; (2) the new nodes that lie along the boundary of the mesh are moved to the true boundary; (3) the interior nodes are moved based on the boundary deformation. The main challenge of these methods arises from step (2). In particular, the curving of the elements along the boundary can result in invalid mesh elements. With that in mind, these high-order mesh generation methods use different approaches in step (3) in an effort to obtain a valid high-order mesh. Methods for transforming the linear mesh usually fall into two groups. The first

e-mail<EMAIL_ADDRESS>
M. Stees (-)

Department of Electrical Engineering and Computer Science, Information and Telecommunication Technology Center, University of Kansas, Lawrence, KS, USA e-mail<EMAIL_ADDRESS>
S. M. Shontz

Department of Electrical Engineering and Computer Science, Bioengineering Program, Information and Telecommunication Technology Center, University of Kansas, Lawrence, KS, USA

S. J. Sherwin et al. (eds.), *Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018*, Lecture Notes in Computational Science and Engineering 134, https://doi.org/10.1007/978-3-030-39647-3\_17

group of methods transform the mesh based on the solution to a partial differential equation [3, 12, 14, 24]. The second group of methods are based on optimization of an objective function [2, 4, 5, 9, 16–21].

In this paper, we describe an optimization-based approach for generating highorder meshes based on affine combinations of nodal positions. The remainder of this paper is organized as follows. In Sect. 2, we present our new method for highorder mesh generation. In Sect. 3, we demonstrate the performance of our proposed method on several aerospace engineering geometries. Finally, in Sect. 4, we offer some concluding remarks and possible directions for our future work.

#### **2 High-Order Mesh Generation Based on Affine Combinations of Nodal Positions**

In this section, we present our optimization-based method for high-order mesh generation. Our proposed method uses affine combinations of nodal positions to determine the movement of the interior nodes after deforming the boundary. Our method consists of three steps. First, for each interior node in the high-order straightsided mesh, an optimization problem is solved to calculate a set of weights that relates the interior node to its neighbors. Second, the boundary nodes are moved to the true boundary. Third, the new positions of the interior nodes are calculated by solving a linear system of equations using the weights and the new boundary positions. In spirit, this method is similar to the weight-based method that we proposed in [19] with two major differences. The first difference is that we propose an affine combination of nodal positions in this work, as opposed to a convex combination. This change allows us to remove the inequality constraint and logbarrier term, leaving only the equality constraints. We also propose an alternative objective function that when combined with the equality constraints allows us to directly solve the optimization problem via a QR factorization. This change results in simplified computational complexity and faster execution time.

To frame our discussion of the method, we introduce the following notation for the 2D formulation of the problem; the 3D formulation is similar. Let the *x*- and *y*-coordinates of the *i*th interior node be represented as *(xi, yi)*. In addition, define the *x*- and *y*-coordinates of the vertices adjacent to node *i* as {*(xj , yj )* : *j* ∈ *Ni*}, where *Ni* is the set of neighbors of node *i*. For each interior node *i*, this information can be represented as the following linear system, where *wij* are the weights:

$$\sum\_{j \in N\_l} w\_{lj} x\_j = x\_l$$

$$\sum\_{j \in N\_l} w\_{lj} y\_j = y\_l,$$

where

$$N\_l = \{\text{high-order nodes of the patch to which } i \text{ belongs}\}.$$

There are several potential choices for the local neighboring set based on use of the low-order nodes, high-order nodes, or both. We include only the highorder nodes as neighbors, as only the high-order boundary nodes move during the boundary deformations. Using either the low-order nodes or both the low- and high-order nodes would dampen the effect the boundary deformation has on the interior nodes, which might lead to tangling near the boundary. Including additional nodes as neighbors would also result in a less sparse matrix when solving (7). Another important consideration is that while the weight calculation is based on only the local neighbors, the position of an interior node is indirectly affected by the deformation of all the interior nodes through the solution of (7).

Adding the additional constraint that the weights sum to one results in the following linear system *Aw* = *b* for finding an affine combination of the x- and y-coordinates of the vertices adjacent to node *i*:

$$
\begin{bmatrix}
\boldsymbol{x}\_{1} \ \boldsymbol{x}\_{2} \ \dots \ \boldsymbol{x}\_{n} \\
\boldsymbol{y}\_{1} \ \boldsymbol{y}\_{2} \ \dots \ \boldsymbol{y}\_{n} \\
\boldsymbol{1} \ \boldsymbol{1} \ \dots \ \boldsymbol{1} \\
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{w}\_{i1} \\
\boldsymbol{w}\_{i2} \\
\vdots \\
\boldsymbol{w}\_{in}
\end{bmatrix} = \begin{bmatrix}
\boldsymbol{x}\_{i} \\
\boldsymbol{y}\_{i} \\
\boldsymbol{1}
\end{bmatrix},
$$

where *n* = |*Ni*|. Based on the set of neighbors, this linear system will be underdetermined (i.e., *A* = *m* × *n* with *m<n*) in general. If we assume that *A* has full rank, we can find one particular solution to our problem by requiring that *w* has the smallest norm of any solution. This results in the following optimization problem:

$$\min\_{w} ||w||\_2^2 \tag{1}$$

$$\text{subject to } A\,w = b.\tag{2}$$

From the Karush–Kuhn–Tucker (KKT) theory [13], we know that the following conditions must hold for a solution *(w*∗*, λ*∗*)* to our problem to be optimal:

$$\nabla\_w \mathcal{L}(w^\*, \lambda^\*) = 0 \tag{3}$$

$$Aw^\* - b = 0\tag{4}$$

$$
\lambda^\*(Aw^\*-b) = 0.\tag{5}
$$

The Lagrangian of our problem is given by:

$$
\mathcal{L}(w, \lambda) = w^T w - \lambda^T (Aw - b),
$$

where *λ* are the Lagrange multipliers.

Using (3)–(5), we can find the following solution pair *(w*∗*, λ*∗*)* as follows:

$$
\nabla\_w \mathcal{L}(w, \lambda) = 2w - A^T \lambda.
$$

$$
\nabla\_w \mathcal{L}(w^\*, \lambda^\*) = 0 \Rightarrow w^\* = \frac{1}{2} A^T \lambda^\*.
$$

$$
A w^\* - b = 0 \Rightarrow A(\frac{1}{2} A^T \lambda^\*) - b = 0.
$$

$$
\lambda^\* = 2(A A^T)^{-1} b
$$

$$
w^\* = \frac{1}{2} A^T \lambda^\* = \frac{1}{2} A^T 2(A A^T)^{-1} b = A^T (A A^T)^{-1} b.
$$

Although we have verified that *(w*∗*, λ*∗*)* is a stationary point, we cannot yet claim that it is a minimum. To do so, we must investigate ∇2 *w*L*(w*∗*, λ*∗*)*:

$$
\nabla\_w^2 \mathcal{L}(w, \lambda) = 2I\_{|N\_l| \times |N\_l|}
$$

$$
\nabla\_w^2 \mathcal{L}(w^\*, \lambda^\*) = 2I\_{|N\_l| \times |N\_l|}.
$$

From the second-order sufficient conditions, if *w*∗ satisfies (3)–(5) and the following condition is satisfied:

$$z^T \nabla\_w^2 \mathcal{L}(w^\*, \lambda^\*) z > 0, \text{for all } z \in C(w^\*, \lambda^\*), z \neq 0,\tag{6}$$

where *C(w*∗*, λ*∗*)* = {*z* | ∇*wc(w*∗*)T z* = 0} is the critical cone and *c(w)* = *Aw* − *b*, then our solution is a minimum. Since ∇2 *w*L*(w*∗*, λ*∗*)* is symmetric positive definite, the inequality in (6) is satisfied for any choice of *z*. Thus we can conclude that our solution *w*∗ is a minimum of (1)–(2).

Now that we have established that *w*∗ is our solution, we will discuss calculating it via a reduced QR factorization. Suppose that *AT* = *QR*, where *Qn*×*m, Rm*×*m* is upper triangular, and *QT Q* = *Im*×*m*. Substituting in the QR factorization of *AT* into *w*∗, we get the following:

$$\begin{aligned} w^\* = A^+ b = A^T (A A^T)^{-1} b &= \mathcal{Q} R (R^T \mathcal{Q}^T \mathcal{Q} R)^{-1} b &= \mathcal{Q} R (R^T R)^{-1} b \\ &= \mathcal{Q} R R^{-1} R^{-T} b = \mathcal{Q} R^{-T} b. \end{aligned}$$

Rearranging this into linear system form, we have:

$$\mathcal{R}^T \mathcal{Q}^T w^\* = b.$$

If we let *t* = *QT w*∗, then *RT t* = *b* and *w*∗ = *Qt*. Thus calculating *w*∗ involves a QR decomposition of *AT* , solving the lower triangular system *RT t* = *b* by forward substitution, and calculating the matrix-vector product *Qt*.

After calculating the weights, a boundary deformation is applied. The final step is to solve for the new locations of the interior nodes [ ˆ*xI , y*ˆ*I* ] by solving the following global linear system:

$$A\_I[\hat{\mathbf{x}}\_I, \hat{\mathbf{y}}\_I] = -A\_B[\hat{\mathbf{x}}\_B, \hat{\mathbf{y}}\_B],\tag{7}$$

where *x*ˆ*B* and *y*ˆ*B* are the new x- and y-coordinates for the boundary nodes, and *AI* and *AB* contain the weights for the interior nodes and boundary nodes, respectively. In this global linear system, each row of the weight matrix corresponds to an interior node with nonzero entries for the node's neighbors and zero entries for the remainder of the row. The resulting global weight matrix is very sparse with irregular structure. In an effort to shift the nonzero entries closer to the diagonal, we apply the sparse reverse Cuthill–McKee ordering provided in Matlab. In Fig. 1, we show the matrix sparsity plots for the natural node ordering and the updated node ordering for the first two examples in Sect. 3. After applying the matrix reordering, the linear system is solved using a sparse LU factorization.

#### **3 Numerical Experiments**

In this section, we demonstrate the results from applying our method to generate several high-order meshes. We use Gmsh [7, 8, 15, 21] to generate the initial straightsided high-order meshes. Our method then uses this mesh to calculate the weights (step 1). Next, we curve the boundary nodes (step 2) using Gmsh. The positions of the interior nodes in the resulting curved high-order mesh are then updated (step 3) by our method. For each example, we show the mesh which results from our method (with high-order nodes visible), the mesh element distortion for the curved boundary mesh generated using Gmsh, and the distortion for the mesh resulting from our method. When reporting the mesh distortion, we list the minimum distortion, maximum distortion, average distortion computed over all elements (referred to as Avg1 in figures), and average distortion computed over curved elements (referred to as Avg2 in figures). The ideal distortion value is 1, indicating that the element is straight. We also list the execution times needed for steps 1 and 3 of our method (excluding I/O) in Table 1. The code was run using Matlab R2018a, and the wallclock execution times were measured on a machine with 16GB of RAM and a Ryzen 7 1700 CPU. All mesh visualizations and distortion evaluations were done using Gmsh. Our first example is a third-order mesh of a square region around

**Fig. 1** Sparsity plots for the first two examples in Sect. 3: (**a** and **c**) the nonzero entries using the original node ordering; (**b** and **d**) the nonzero entries after applying the sparse reverse Cuthill– McKee ordering

a NACA0012 airfoil. In Fig. 2a, b, we show the mesh resulting from our method and a table of the mesh quality values as measured by the distortion metric. In this example, our method increased the minimum distortion from 0.744 to 0.799, while causing only minor changes in the average distortion.

In our second example, we extrude the NACA0012 airfoil and create a thirdorder mesh of the resulting region. In Fig. 3a, b, we show the mesh resulting from our method, and a table of the mesh quality values. For this example, our method improved the minimum distortion by 0.125, increasing it from 0.317 to 0.442.

**Table 1** The number of elements and the wall clock times for steps 1 and 3 of our method (excluding I/O) for each example

**Fig. 2** NACA0012 airfoil example: (**a**) the mesh resulting from our method and (**b**) the mesh quality as measured by the element distortion metric

(b)

Our third and final example is a second-order mesh of an Airbus A319 aircraft. Unlike our previous examples, this geometry resulted in tangled elements after curving the boundary. Although our method still increased the minimum quality, it was not able to untangle the mesh. To address this, we applied the high-order

**Fig. 3** Extruded NACA0012 airfoil example: (**a**) the mesh resulting from our method and (**b**) the mesh quality as measured by the element distortion metric

(b)

(c)

regularization scheme available in Gmsh as a post-processing step. Aside from changing the target Jacobian range to 0.3–2 and fixing the boundary nodes, all other parameters were left at their default values. The untangling for the original mesh took 14.14 s, while untangling the mesh resulting from our method required only 1.64 s. In Fig. 4a–c, we show the surfaces of the mesh resulting from our method, a view of a cut through the interior volume mesh, and a table of the mesh quality values. In Fig. 4c, we list the distortion for the mesh after curving the boundary, the distortion for the mesh resulting from our method, and the distortions of both meshes after applying the regularization scheme in Gmsh.

Aside from the third test case, all of these examples were relatively straightforward. In each case, our method increased the minimum distortion when compared to only curving the nodes along the boundary. While additional testing is necessary to confirm this, our results for the third example seem to indicate that our method could be used to reduce the severity of the mesh tangling and thus simplify the work for an untangling method during post-processing.

**Fig. 4** Airbus A319 example: (**a**) surfaces of the mesh resulting from our method; (**b**) a view of a cut through the interior volume mesh, and (**c**) the quality of the mesh with only boundary curving, the quality of the mesh resulting from our method, and the quality of both meshes after applying the regularization scheme available in Gmsh

#### **4 Concluding Remarks and Future Work**

We have presented a new optimization-based method for generating high-order meshes. Our examples have shown that the proposed method based on affine combinations of nodal positions tends to improve the quality of the most distorted elements, while causing minor changes to the least distorted elements. While our approach is optimization-based, we have demonstrated that the optimization problem can be solved directly using a QR factorization as opposed to the typical iterative optimization approach. This change results in lessened computational complexity and reduced execution time.

As part of our future work, we will consider other definitions for the set of neighbors of an interior node. We will also investigate other aspects of the linear system including other node reordering schemes and solvers. Finally, we will apply the untangling method that we proposed in [20] after extending it to 3D.

**Acknowledgements** The work of the first author was funded in part by the Madison and Lila Self Graduate Fellowship and NSF CCF grant 1717894. The work of the second author was supported in part by NSF grants CCF 1717894 and OAC 1808553. The authors wish to thank the anonymous referee for his or her careful reading of the paper and for the helpful suggestions which strengthened it.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Sparse Spectral-Element Methods for the Helically Reduced Einstein Equations**

**Stephen R. Lau**

#### **1 Introduction**

To model the inspiral and merger of binary objects (blackholes or neutron stars), many researchers have been solving the Einstein equations numerically. Such simulation involves both the construction of gravitational initial data at time *t*0 and its subsequent evolution to a final time *tF* ' *t*0. Interpretation of experimental detections of gravitational waves relies on numerical simulation. Moreover, detection of weak signals is facilitated by statistical techniques alongside "template banks" of numerically generated signals. We consider a nonstandard problem, solution of the Einstein equations reduced by helical symmetry, as described by Beetle, Bromley, Hernández, and Price (BBHP) [1, 2]. Heuristically, helical reduction is a data+evolution synthesis. Although solutions to the BBHP equations are ultimately unphysical, they may approximate the early phase of inspiral and serve as reduced order models. Moreover, they may provide excellent "trial data" (the starting point for the construction of *t* = *t*0 initial data). Finally, they would address bewitching mathematical issues concerning exact helical symmetry in general relativity.

We consider a spectral element approach [3–8] for solving the BBHP equations. Although the equations involve a mixed-typed operator *L*, we solve them via relaxation using a Broyden-Krylov approach. The computational domain which surrounds the compact objects is split into 11+ subdomains (blocks, spherical shells, and cylindrical shells with classical spectral expansions thereon). To rapidly solve the linear systems arising in our scheme, we have developed sparse modal methods based on the application of spectral integration matrices, extending ideas originally described in the 1990s by Coutsias, Hagstrom, Hesthaven, and Torres.

S. R. Lau (-)

Mathematics and Statistics, University of New Mexico, Albuquerque, NM, USA e-mail<EMAIL_ADDRESS>
© The Author(s) 2020

S. J. Sherwin et al. (eds.), *Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018*, Lecture Notes in Computational Science and Engineering 134, https://doi.org/10.1007/978-3-030-39647-3\_18

We use preconditioned GMRES to solve these systems, with standard domain decomposition methods. In addition, we have developed fast methods for inversion of subdomain approximations of *L*, either via modal-based preconditioning or direct schemes.

#### **2 Background**

This section first describes the helically reduced wave equation (HRWE), a model for the helically reduced Einstein equations. Our HRWE description fixes ideas.

#### *2.1 Helically Reduced Wave Equation*

The wave equation is *ψ* = 0, where = −*∂*2 *t* + *Δ***x**. Assume that *ψ* rotates rigidly with rate *Ω*. With the *z*-axis as the rotation axis, *ψ* then depends on time *t* only through *ϕ* = *φ* − *Ωt*, where *φ* is the azimuthal angle. Via the *∂t* → −*Ω∂ϕ* replacement, *ψ* = 0 becomes the HRWE *Lψ* = 0, where *L* = *Δ***x**˜ − *Ω*2*(x∂*˜ *y*˜ − ˜*y∂x*˜*)*2 in terms of co-rotating coordinates *(t ,*˜ *x,*˜ *y,*˜ *z)*˜ = *(t, x* cos *Ωt* + *y* sin *Ωt,y* cos *Ωt* − *x* sin *Ωt, z)*.

 Its boundary *∂D* = *∂S***−** *I* ∪ *∂S***−** *II* ∪ *∂S***+** out is the union of two inner spheres (the **−**'s) and one outer sphere (the **+**). We consider the mixed-type problem

$$L\psi = 0,\qquad \psi = f^- \text{ on } \partial S\_I \cup \partial S\_{II}, \qquad \left(\partial\_r - \Omega \,\partial\_\psi + r^{-1}\right)\psi = f^+ \text{ on } \partial S\_{\text{out}}^\dagger,\tag{1}$$

**Fig. 1** 2-center domain decomposition. The left panel depicts the *z*-cross section of the inner subdomains with *S*out suppressed. The right panel depicts all subdomains, although for visualization the outer radius for *S*out has been chosen rather small. In this work the blocks *B*2 and *B*4 are absent. (**a**) Inner domain decomposition. (**b**) Double cross section

with inner Dirichlet conditions and an outer radiation boundary condition on *ψ*. For simplicity here, we have put on *∂S***+** out a simple Sommerfeld condition; in practice, *f* **+** is a nonlocal function of *ψ* enforcing an exact Dirichlet-to-Neumann map [4].

A class of solutions to (1) stems from Liénard-Wiechert potentials. Indeed, consider *Lψ* = −16*πγMδ(***x**˜ − **x**˜*p)*. The "particle" location **x**˜*p* = *(*±*a,* 0*,* 0*)* is the center of either *∂S***−** *II* or *∂S***−** *I* ; whence the equation is homogeneous on *D*. The mass *M* and relativistic factor *γ* = *(*1 − *v*2*)*−1*/*2 are constants, with *v* = *aΩ <* 1 so the particle moves subluminally in the *(t , x, y, z)* frame. The retarded solution to this problem is

$$\psi(\tilde{x}, \tilde{y}, z) = \psi(\rho \cos \varphi, \rho \sin \varphi, z) = \frac{4\gamma M}{\lambda \mp v\rho \sin(\varphi + \mathcal{Q}\lambda)} =: \gamma^2 \frac{4M}{\mathcal{R}}.\tag{2}$$

Evaluation of this expression involves a numerical component: solution of the fixed point equation *λ* = *z*2 + *ρ*2 + *a*2 ∓ 2*aρ* cos*(ϕ* + *Ωλ)*1*/*2 for (the retarded time) *λ*.

#### *2.2 Helically Reduced Einstein Equations*

We consider the vacuum Einstein equations in Landau-Lifshitz form. Write the (densitized contravariant) metric tensor as g*μν* = *ημν* − *hμν* , where *ημν* = diag*(*−1*,* 1*,* 1*,* 1*)* is the flat metric and *hμν* is the *metric perturbation*. Assume the *harmonic gauge condition ∂νhμν* = 0. Then the vacuum Einstein equations can be expressed as

$$
\Box h^{\mu\upsilon} = S^{\mu\upsilon\kappa\beta}\_{\tau\phi\gamma\alpha}(\mathfrak{g}) \frac{\partial h^{\tau\phi}}{\partial x^{\kappa}} \frac{\partial h^{\gamma\alpha}}{\partial x^{\beta}} + h^{\alpha\beta} \frac{\partial^2 h^{\mu\upsilon}}{\partial x^{\alpha} \partial x^{\beta}},\tag{3}
$$

where *Sμνκβ τφγ α(*g*)* depends on g*μν* and its inverse g*μν* (but not on derivatives of either). Einsteins equations are then a constrained system of 10 nonlinear wave equations.

The BBHP reduction of the Einstein equations is similar to the one outlined for the wave equation. Technically, it assumes the existence of a *Killing vector field*, but we give a brief and heuristic description of their approach. The challenge is that the perturbations *hμν* themselves are not "helical scalars"; helical reduction is therefore not tantamount to the replacement *hμν* → *Lhμν* . However, BBHP have introduced helical scalars *ψA* through which the reduction can be carried out. These are

$$\begin{cases} \psi^{(nn)} = h^{\text{II}}, & \psi^{(n0)} = \sqrt{2}h^{\text{Iz}}, \quad \psi^{(00)} = \sqrt{\frac{1}{2}}(h^{\text{xx}} + h^{\text{yy}} + h^{\text{zz}}) \\\\ \psi^{(20)} = -\sqrt{\frac{1}{6}}(h^{\text{xx}} + h^{\text{yy}} - 2h^{\text{zz}}), & \psi^{(n1)} = e^{\text{i}\Omega t} \left( -h^{\text{Ix}} + \text{i}h^{\text{Iy}} \right) \\\\ \psi^{(21)} = e^{\text{i}\Omega t} \left( -h^{\text{x}\varepsilon} + \text{i}h^{\text{y}\varepsilon} \right), & \psi^{(22)} = e^{\text{i}\Omega t} \left[ \frac{1}{2}(h^{\text{x}\varepsilon} - h^{\text{yy}}) - \text{i}h^{\text{xy}} \right]. \end{cases} (4)$$

These 4 real and 3 complex quantities contain the 10 degrees of freedom in the *hμν* . We express this transformation as *ψA* = *e*i*μ(A)ΩtMA μνhμν* , where *A* runs over the (tensor-spherical-harmonic) labels, and *μ(A)* is 0,1, or 2. Contraction of *MA μν* on (3), with subsequent helical reduction based on the action of on *ψA*, yields

$$\begin{aligned} &L\psi^A - 2\text{i}\mu(A)\mathfrak{L}^2\partial\_{\psi}\psi^A + \mu^2(A)\mathfrak{L}^2\psi^A = \\ &h^{a\beta}\psi^A\_{,\alpha\beta} - 2\text{i}\mu(A)\mathfrak{L}h^{\dagger\beta}\psi^A\_{,\beta} - \mu^2(A)\mathfrak{L}^2h^{\dagger\iota}\psi^A + M^A\_{\ \mu\upsilon}S^{\mu\upsilon\kappa\beta}\_{\text{\tiny tr\beta\forall\alpha}}(\mathfrak{g})h^{\dagger\phi}{}\_{,\kappa}h^{\prime\nu\alpha}{}\_{,\beta}.\end{aligned} \tag{5}$$

Here *L* is the operator appearing in the HRWE. Similar to the boundary conditions appearing in (1), the boundary conditions we adopt for (5) are

$$\psi^A = (f^A)^\top \text{ on } \partial S\_I \cup \partial S\_{II}, \qquad \left(\partial\_r - \mathcal{Q}\,\partial\_\theta + r^{-1}\right) e^{\mathbf{i}\mu(A)\varphi} \psi^A = (f^A)^\top \text{ on } \partial S\_{\text{out}}^\bot. \tag{6}$$

Again for simplicity, here we have a Sommerfeld condition on *∂S***+** out, but in practice use a nonlocal outgoing condition based on an exact Dirichlet-to-Neumann map.

Price has written down the analog of (2) for linearized gravity, i.e. (5) when the right-hand side of the equation is set to zero. This solution may be viewed as an exact solution to *hμν* = −16*πT μν* , where the stress energy tensor *T μν* corresponds to a massive point particle in an eternal circular orbit (as discussed, such a point source is excised from our domain *D*). Price's solution is analogous to the electromagnetic solution given by G. A. Schott, and it is given by the leading 1*/*R terms in the appendix expressions (12). Unfortunately, *∂νhμν* = 0 for this solution.

#### **3 Sparse Spectral Element Methods**

This section summarizes our numerical methods. It is necessarily impressionistic, as even an incomplete presentation of details would take too much space.

#### *3.1 Overview*

We split *D* into subdomains, here with the minimal configuration of 11 subdomains: blocks *B*1*,*3*,*5; cylinders *C*1*,*2*,*3*,*4*,*5; an inner shell *SI* around "particle" *I* ; an inner shell *SII* around "particle" *II* ; and an outer shell *S*out. This corresponds to a "binary blackhole" (BBH) domain with two excised inner balls. For a "binary neutron star" (BNS) domain, we further split both *SI* and *SII* into two overlapping concentric spheres (a stellar surface then resides in each overlap [8]), and fill in the excised regions with two extra blocks *B*2*,*4. Figure 1 depicts a BNS domain. Pioneered by Pfeiffer et al. [9], such decompositions are used in the EllipticSolver of SpEC [10].

The unknowns in our approach are the modal expansion coefficients associated with subdomain expansions in terms of classical (Chebyshev, Fourier, and spherical harmonic) basis functions. As described below, we make extensive use of integration matrices to achieve sparse representations of the relevant operators. Before presenting details, we first address how we handle the nonlinearities in (5). Let *ψ* 3*A* (the vector of unknowns) be a concatenation of the modal coefficients from all 11 subdomains. Then, as sketched below, upon approximation (5) becomes

$$
\mathcal{A}\mathcal{B}\widetilde{\mathcal{Q}}\widetilde{\Psi}^A = \mathcal{A}\widetilde{\mathfrak{g}}^A,\tag{7}
$$

where *BL* approximates the operator on the left-hand side of (5), with *B* representing the action of "integration preconditioning" (see below) on each subdomain. For linearized gravity, with the right-hand side of (5) set to zero, the vector 3*gA* is zero, save for select entries related to the inner Dirichlet values *(f A)***−** of *ψA* on *∂S***−** *I* ∪ *∂S***−** *II* . For the full Einstein equations 3*gA* depends on *ψ* 3*A*, and its evaluation relies on spectral analysis/synthesis (forward/backward transform) and numerical differentiation on each subdomain. We then view

$$
\ddot{\boldsymbol{\Psi}}\_k^A = (\boldsymbol{\cir}\boldsymbol{\mathcal{X}})^{-1}\boldsymbol{\mathcal{X}}\widetilde{\mathbf{g}}^A(\ddot{\boldsymbol{\Psi}}\_{k-1}^B) \tag{8}
$$

as a fixed-point equation, accelerating its convergence with the Broyden algorithm.

This approach relies on approximation and inversion of the operator appearing on the left-hand side of (5). Reference [4] is a detailed account of the case *μ(A)* = 0, i.e. the HRWE. For *μ(A)* = 1 or 2, the operator mixes the *UA* and *V A* in *ψA* = *UA* + i*V A*. We have not yet described our treatment of this scenario, but note that it relies on Schur-complement techniques. Here we describe only the *μ(A)* = 0 case.

#### *3.2 Integration Preconditioning and Other Key Aspects*

We use "integration preconditioning" [11] in order to achieve sparse linear systems. Especially for the cylinders and inner shells, the details are formidable. We convey the basic ideas with the Laplacian *Δ*, rather than *L*, on a cube (suppressing tildes on the co-rotating coordinates). Let *u(x, y, z)* ≈ <*Nx i*=0 <*Ny j*=0 <*Nz k*=0 3*uijkTi(x)Tj (y)Tk(z)* obey *Δu* = *g* on C = [−1*,* 1] 3. An approximation of the Poisson equation is

$$(D\_\chi^2 \otimes I\_\chi \otimes I\_\mathfrak{z} + I\_\chi \otimes D\_\chi^2 \otimes I\_\mathfrak{z} + I\_\chi \otimes I\_\chi \otimes D\_\mathfrak{z}^2)\widetilde{\mathbf{u}} = \widetilde{\mathbf{g}},\tag{9}$$

with 3**u** the vector of 3*uijk* with appropriate ordering, and *D*2 representing double differentiation in the modal Chebyshev basis. Let *B*2 [2] represent double integration in the modal basis, where the [2] indicates that the first two rows have all zero entries. To (9) we apply the "preconditioner" B = *B*2 *x*[2] ⊗ *B*2 *y*[2] ⊗ *B*2 *z*[2] , thereby reaching

$$\mathbb{P}\left(I\_{\mathbf{x}\{2\}}\otimes B\_{\mathbf{y}\{2\}}^2\otimes B\_{\mathbf{z}\{2\}}^2 + B\_{\mathbf{x}\{2\}}^2\otimes I\_{\mathbf{y}\{2\}}\otimes B\_{\mathbf{z}\{2\}}^2 + B\_{\mathbf{x}\{2\}}^2\otimes B\_{\mathbf{y}\{2\}}^2\otimes I\_{\mathbf{z}\{2\}}\tilde{\mathbf{u}} = \mathfrak{B}\widetilde{\mathbf{g}}.\qquad(10)$$

The system (10) is sparse, and the number of *empty rows* (all 0's) equals the number of *tau-conditions* to be enforced, auxiliary equations enforcing, say, *u*|*∂*C = *f* .

"Integration preconditioning" of (9) results in the sparse system (10), but the issue of condition number is subtle. Indeed, passage from (9) to (10) arguably *worsens conditioning*. For this reason, integration preconditioning has been viewed as bad for PDE; we call the technique *integration sparsification*. Conditioning issues are then surmounted either by further *genuine* preconditioning on top of the sparsification or fast direct solves. Our use of 2-center domains with sparse modaltau methods features: (a) sparse representation of *L* on subdomains; (b) "gluing" of conforming and overlapping subdomains; (c) modal-based preconditioning of subdomain solves; (d) "fast" direct-solves on blocks and *S*out; (e) standard global preconditioning; (f) low-rank treatment of stellar surfaces; (g) nonlocal domain reduction.

#### *3.3 Current Complexity Estimates*

We aim to solve the linear systems (say approximating the HRWE or the Helmholtz equation posed on the 2-center domain *D*) arising in our problems at *demonstrably* sub-quadratic complexity; indeed, we believe that an order- 5 3 complexity is achievable. This is the complexity associated with matrix-vector multiplication; despite our sparse representations, the "gluing" of overlapping subdomains in our decomposition of *D* prevents realization of a linear-complexity matrix-vector product.

To document progress towards our goal, we summarize our current complexity estimates associated with solution of the HRWE on each subdomain type. These solves serve as part of our preconditioner for the global GMRES solution of the HRWE on *D*. For this discussion, let N represent the total number of modes on a given subdomain; e.g. N = *(Nx* + 1*)(Ny* + 1*)(Nz* + 1*)* for the block considered above.

For *S*out let M be the matrix which represents*r*2*L* (the *r*2 factor here is explained in [4]) and includes inserted tau-vectors to enforce boundary conditions. Ignoring tau-vectors, M is block diagonal in the spectral space of spherical harmonics indexed by *(, m)*. View its elements as M*mk, m k* , where *mk* is a "clumped index' and *k, k* are Chebyshev indices. Then, apart from tau-vectors, M*mk, mk* = 0, unless = and *m* = *m* . Moreover, each *(Nr* + 1*)* ×*(Nr* + 1*)* block M*mk,mk* is itself sparse and banded. While these desirable structures are somewhat spoiled by the tau-conditions, through the use of the Woodbury formula and band solvers, for *S*out we are able to directly invert M (i.e. solve the HRWE on *S*out) at *O(*N*)* cost.

For the inner shells, *SI* and *SII* , our representation M of *r*2*L* is only block banded. Indeed, the centers of these shells lie off the rotation axis, and *r*2*L* mixes spherical harmonic modes. Its spectral representation [4] is remarkably complicated, and relies on identities for spherical harmonics found in the treatise [12]. We solve the HRWE on inner shells via preconditioned GMRES, with a modal block-Jacobi preconditioner defined by inversion of the diagonal blocks M*mk,mk*. Apart from tau-conditions, these blocks are again sparse and banded, and therefore amenable to the fast methods alluded to in the last bullet. Construction and reuse of this block-Jacobi preconditioner therefore has *O(*N*)* cost. Moreover, we have empirically observed (see [4]) that such preconditioning yields low and essentially *resolution independent* iteration counts. While more analysis is needed, from a practical standpoint solution of the HRWE on an inner shell has an *O(*N*)* cost.

The situation on blocks is worse. Part of our global preconditioner for solving the HRWE on *D* involves inversion of the Poisson problem on blocks as an approximation to the HRWE (and inversion of the Helmholtz equation when the spin index *μ(A)* = 0). This works extremely well, likely due to the fact that *Ω* 1 and the blocks are close to the rotation center. In any case, we solve the Poisson/Helmholtz problem on a block via a direct approach [13] based on a rank-augmenting generalization of the Woodbury formula. This direct method is empirically well-conditioned and low-memory. Moreover, it has an *O(*N2*)* set-up cost with a small constant, followed by an *O(*N4*/*3*)* cost for subsequent solves. If possible, we hope to reduce the set-up cost to an *O(*N5*/*3*)* complexity.

The situation on cylinders is worst of all, although to date we have not focused much attention on these subdomains. We solve the HRWE on cylinders (or the collection of cylinders) via GMRES, with modal-based preconditioners that empirically yield resolution-independent iteration counts. Application of the preconditioner currently involves an *O(*N7*/*3*)* set-up cost, followed by an *O(*N5*/*3*)* cost for subsequent solves. Here, we believe improvement is possible.

#### **4 Numerical Tests**

Our decomposition of *D* is from Table IV of [4], except here *∂S*out has *r*out = 15. For that table *∂S***−** *I* has radius *rI,*min = 0*.*4 and center *(x*˜*I* = −0*.*9*,* 0*,* 0*)*, and *∂S***−** *II* has radius *rII ,*min = 0*.*3 and center *(x*˜*II* = 1*.*0*,* 0*,* 0*)*. The coordinates *X,* 3 *Y ,* 3 *Z* 3 in [4] are *y,*˜ *z,*˜ *x*˜ here. The subdomain truncations (number of modes) adopted here are nearly the same as those in [4]; however, here we only record the total number of modes over all subdomains. Unless otherwise stated, *Ω* = 0*.*075, *MI* = 0*.*05, and *MII* = 0*.*1.

#### *4.1 Comparison with Exact Solutions*

Our first test uses PriceLG, the aforementioned solution for linearized gravity due to Price. Here the solution is a superposition of two point-sources with the above masses. Each source point's contribution to the helical scalars is defined by the leading 1*/*R term in (12), and these expressions seed inner Dirichlet conditions on *∂S***−** *I* ∪ *∂S***−** *II* . The outer boundary conditions are nonlocal conditions which are exact for this solution. Table 1 lists errors for *ψ(nn)*; errors for the other scalars are similar. These are relative *L*2-errors (against the exact solution) computed on both *B*3 and *S*out via interpolation onto uniform reference grids. Since the problem is linear, each line of table corresponds to a single GMRES solve performed in parallel on 10 nodes. For the table middle about 0.4% of the matrix entries are nonzero. The last table line has *∂νhμν*rms (3.9420e-12,1.7310e-03,8.9752e-05,9.0076e-16), with the rms calculation taken over the (relatively coarse) dual-nodal subdomain grids. The first and last components of *∂νhμν* vanish for the exact PriceLG solution.

Our next test is SchwarzH, the Schwarzschild metric in harmonic coordinates:

$$h^{ll} = -1 + r^{-2}(r+M)^3/(r-M), \quad h^{jk} = r^{-2}M^2\nu^j\nu^k. \tag{11}$$

Here the radius *r* and the direction cosines *ν* = *(*sin *θ* cos *φ,*sin *θ* sin *φ,* cos *θ )* are chosen relative to a point *(x*0*, y*0*, z*0*)* = (-0.9+1.37e-3,-1.6854e-4,2.9985 e-3) which is off-center but close to the center of *∂S***−** *I* . The mass is *M* = 0*.*05, and for this choice the horizon of the blackhole lies inside of (but is not concentric with) *∂S***−** *I* . For this test *Ω* = 0, and the exact solution (11) seeds inner Dirichlet boundary conditions on both *∂S***−** *I* and *∂S***−** *II* . On *∂S***+** out rather than radiation conditions, we adopt an inhomogeneous Neumann condition based on the exact solution. Table 2 lists

**Table 1** PriceLG solution

Here relative *L*2 errors are listed only for *ψ(nn)*. The lowest resolution run has zero initial iterate; afterwards the initial iterate stems from the previous solution. Respectively, tGMRES and iGMRES are the tolerance and iteration number for the GMRES solve

**Table 2** SchwarzH solution

As for the PriceLG test, only errors for *ψ(nn)* are listed

**Table 3** Full GR with the same boundary conditions as PriceLG

errors with the same meanings as before. Now *∂νhμν*rms for each *μ* is comparable to the corresponding table errors; i.e. the gauge constraint converges to zero. Table 2 also lists the number iBROY of iterations performed by the Broyden solver to achieve the tolerance tBROY. Each Broyden iteration itself involves a linear solve via GMRES. Each GMRES tolerance is tBROY*/*10, the same as the corresponding line in Table 1.

#### *4.2 Gauge Constraint Tests*

Our next two tests explore to what extent the gauge constraint *∂νhμν* = 0 is satisfied for the Einstein problem (5) and (6) in a binary scenario. For the first test we redo the PriceLG test, except now with the Einstein equations. The inner and outer boundary conditions are exactly as before. Table 3 again lists errors for *ψ(nn)* with the same meanings as before. Errors are computed against the finest-resolution numerical solution. The table also lists the *L*2-norm of the nonlinear residual. The last table line has *∂νhμν*rms (3.2031e-03 5.1309e-03 4.5881e-03 4.6711e-03). That is, the gauge constraint does not converge to zero. Since the harmonic gauge is not satisfied, we cannot view these as solutions to the Einstein equations!

Presumably, the violation of the harmonic gauge in the preceding example stems from the fact that *∂νhμν* = 0 for the Price solution used to fix the inner boundary conditions. Beetle, Bromley, and, Hernández, and Price have given a refined set of inner boundary conditions, based on the near-field asymptotics of a moving Schwarzschild blackhole and meant to improve on the point-particle boundary conditions. The appendix lists (our understanding of) these conditions. With the hope that these refined boundary will result in lower gauge errors, we have performed the previous test with them. Convergence of the numerical solution is similar, but *∂νhμν*rms has comparable size and still does not converge to zero.

#### **5 Conclusions and Acknowledgments**

Our tentative conclusion for helically symmetric BBH models is that violation of the gauge constraint stems from imperfect inner boundary conditions. We have also found a persistent gauge error in our BNS models, despite being several orders of magnitude smaller in that context. The BNS model, with stars in place of excised regions, involves no inner boundary conditions. We believe that in this context it is the outer boundary conditions which give rise to the constraint violation. Likely, at both the inner and outer boundaries some of the helical scalars need to be fixed using the gauge constraint itself (similar to "constraint preserving boundary conditions" used in evolution codes). The author is grateful for correspondence with Richard H. Price, and for assistance from UNM's Center for Advanced Research Computing.

#### **Appendix: Beetle, Bromley, Hernández, and Price Inner Boundary Conditions**

This appendix lists expansions for the helical scalars which somewhat generalize the ones in [1]. We need two expansions, one for a "particle" at *(*−*aI ,* 0*,* 0*)*, and one at *(aII ,* 0*,* 0*)*. Each has its own mass *MI,II* . The top choice of ± or ∓ refers to *II* and the bottom to *I* . For both *aI* and *aII* , we define *v* = *aΩ*, *γ* = *(*1 − *v*2*)*−1*/*2, R = *γ λ* ∓ *vγρ* sin*(ϕ* + *Ωλ)*, *KR* = −*ρ* cos *ϕ* ± *a* cos*(Ωλ)* ± *vγ*R sin*(Ωλ)*, and *KI* = *ρ* sin *ϕ*±*a* sin*(Ωλ)*∓*vγ*R cos*(Ωλ)*. For *λ* see after (2). Assuming *M/*R 1, we have

$$
\psi^{(nn)} \sim \chi^2 \left(\frac{4M}{\mathcal{R}} + \frac{7M^2}{\mathcal{R}^2}\right) + \frac{M^2(\lambda - \chi\mathcal{R})^2}{\mathcal{R}^4} \tag{12a}
$$

$$
\psi^{(n0)} \sim 0 \cdot \left(\frac{4M}{\mathcal{R}} + \frac{7M^2}{\mathcal{R}^2}\right) + \sqrt{2} \frac{M^2 (\lambda - \mathcal{Y}\mathcal{R}) z}{\mathcal{R}^4} \tag{12b}
$$

$$\psi^{(00)} \sim \frac{v^2 \chi^2}{\sqrt{3}} \left( \frac{4M}{\mathcal{R}} + \frac{7M^2}{\mathcal{R}^2} \right) + \frac{M^2 \left[ \lambda^2 - 2\gamma \lambda \mathcal{R} + (2 + v^2 \chi^2) \mathcal{R}^2 \right]}{\sqrt{3} \mathcal{R}^4} \tag{12c}$$

$$\psi^{(20)} \sim -\frac{v^2 \nu^2}{\sqrt{6}} \left(\frac{4M}{\mathcal{R}} + \frac{7M^2}{\mathcal{R}^2}\right) - \frac{M^2 \left[\lambda^2 - 3z^2 - 2\gamma\lambda\mathcal{R} + (2 + v^2\nu^2)\mathcal{R}^2\right]}{\sqrt{6}\mathcal{R}^4} \tag{12d}$$

$$\psi^{(nl)} \sim \pm i v \gamma^2 e^{i\Omega\lambda} \left(\frac{4M}{\mathcal{R}} + \frac{7M^2}{\mathcal{R}^2}\right) + \frac{M^2(\lambda - \gamma\mathcal{R})(K^R + iK^I)}{\mathcal{R}^4} \tag{12e}$$

$$
\psi^{(21)} \sim 0 \cdot \left(\frac{4M}{\mathcal{R}} + \frac{7M^2}{\mathcal{R}^2}\right) + \frac{M^2 z (K^R + iK^I)}{\mathcal{R}^4} \tag{12f}
$$

$$\Psi^{(22)} \sim -\frac{1}{2}v^2\gamma^2 e^{2i\Omega\lambda} \left(\frac{4M}{\mathcal{R}} + \frac{7M^2}{\mathcal{R}^2}\right) + \frac{M^2}{2\mathcal{R}^4} [(K^R)^2 - (K^I)^2 + 2iK^R K^I].\tag{12g}$$

Worried about a possible sign discrepancy with the results in [1], we have also considered (12a–g) with all correction terms (those with R4 in the denominator) flipped by a sign.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Spectral Analysis of Isogeometric Discretizations of 2D Curl-Div Problems with General Geometry**

**Mariarosa Mazza, Carla Manni, and Hendrik Speleers**

### **1 Introduction**

In this paper we focus on isogeometric Galerkin discretizations of the weighted curl-div operator

$$\mathcal{L}\_{\alpha,\beta}\mathfrak{u} := \alpha \nabla \times \nabla \times \mathfrak{u} - \beta \nabla \nabla \cdot \mathfrak{u}, \quad 0 < \alpha, \beta. \tag{1}$$

This parameter-dependent operator appears in several problems, including the Stokes equation and Maxwell equations [2]. Moreover, containing a weighting of the curl and div operators, it captures the essential features of the so-called Alfvénlike operator [14], which is of interest in magnetohydrodynamics [15]. We note that L*α,β* can be seen as a weighted Laplacian for vector fields (equivalently, Hodge Laplace for 1-forms). Indeed, when *α* = *β* = 1, it is equal to the standard (negative) vector Laplace operator, i.e.,

$$\nabla \times \nabla \times \mathfrak{u} - \nabla \nabla \cdot \mathfrak{u} = -\nabla^2 \mathfrak{u} \dots$$

We assume that (1) is defined on a sufficiently smooth domain ∈ R2 that can be described through a geometry map *G* : [0*,* 1] 2 → , and we consider homogeneous

M. Mazza (-)

C. Manni · H. Speleers

© The Author(s) 2020

Max Planck Institute for Plasma Physics, Garching, Germany e-mail<EMAIL_ADDRESS>
Department of Mathematics, University of Rome Tor Vergata, Rome, Italy e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
S. J. Sherwin et al. (eds.), *Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018*, Lecture Notes in Computational Science and Engineering 134, https://doi.org/10.1007/978-3-030-39647-3\_19

Dirichlet (no-slip) boundary conditions, i.e., *u* = 0 on *∂*. This leads us to the following variational formulation

$$\alpha(\mathcal{L}\_{\mathfrak{u},\beta}\mathfrak{u},\mathfrak{v}) = \mathfrak{a}(\nabla \times \mathfrak{u}, \nabla \times \mathfrak{v}) + \beta(\nabla \cdot \mathfrak{u}, \nabla \cdot \mathfrak{v}), \quad \mathfrak{u}, \mathfrak{v} \in \left(H\_0^1(\Omega)\right)^2. \tag{2}$$

We refer the reader to [3, 15] for a discussion about well-posedness.

To find an approximate solution of the problem L*α,βu* = *f* , with the stated boundary conditions, we consider the variational formulation (2) in a finite dimensional vector space V*h* ⊂ *H*1 0 *()*2 , i.e.,

$$\alpha(\mathcal{L}\_{a,\beta}\mathfrak{u}\_h, \mathfrak{v}\_h) = \alpha(\nabla \times \mathfrak{u}\_h, \nabla \times \mathfrak{v}\_h) + \beta(\nabla \cdot \mathfrak{u}\_h, \nabla \cdot \mathfrak{v}\_h), \quad \mathfrak{u}\_h, \mathfrak{v}\_h \in \mathbb{V}\_h. \tag{3}$$

We focus on isogeometric analysis (IgA) as discretization technique, where the approximation space V*h* is chosen to be composed of vector fields whose components are linear combinations of tensor-product B-splines mapped according to *G*.

The discretization (3) leads to solving linear systems, which turn out to be severely ill-conditioned and require ad hoc fast solvers for a proper treatment [4, 6, 15]. This requires a deep understanding of the spectral properties of the related matrices. They depend on many factors: the problem parameters *α*, *β*, the basic curl and div operators, the mesh-size, the degree of the B-spline approximation, and the map *G* used to describe the geometry of the computational domain.

In this paper we provide a spectral study of these matrices using the theory of (multilevel block) Toeplitz [13, 17, 19] and generalized locally Toeplitz [10–12] sequences. More precisely, we show that such matrices admit a spectral distribution which can be described in terms of a so-called *spectral symbol*. We determine this spectral symbol and we reveal its dependence on the characteristic parameters of the problem listed above. The spectral analysis presented in this paper extends the results of [15] to the case of non-trivial geometry and relies on the spectral theory developed for isogeometric discretizations of elliptic problems in [7, 8]. We also refer the reader to [16] for a spectral analysis of the curl-curl operator.

The remainder of the paper is organized as follows. In Sect. 2 we introduce notations and definitions relevant for our spectral analysis, and we recall the basics of B-splines. In Sect. 3 we detail the IgA discretization matrices and we perform a spectral analysis of them. We numerically illustrate those results in Sect. 4. Finally, we conclude the paper in Sect. 5.

#### **2 Preliminaries**

In this section we collect some preliminary tools on spectral analysis and IgA discretizations. In particular, we recall the formal definition of spectral distribution for a general matrix-sequence and the definition of (cardinal) B-splines.

#### *2.1 Spectral Distribution*

Throughout the paper, we follow the standard convention for operations with multiindices (see e.g. [9, 18]). Given a multi-index *n* := *(n*1*,...,nd )* ∈ N*d* , we say *n* → ∞ if *ni* → ∞, *i* = 1*,...,d*. Let C0*(*C*)* be the set of continuous functions *F* : C → C with compact support.

**Definition 1** Let *f* : *D* → C*s*×*s* be a measurable matrix-valued function, defined on a measurable set *D* ⊂ R*q* with *q* ≥ 1, 0 *< μq (D) <* ∞, where *μq* is the Lebesgue measure. Let {*An*}*n* be a matrix-sequence with dim*(An)* =: *dn* and *dn* → ∞ as *n* → ∞. Then, {*An*}*n* is distributed like the pair *(f, D)* in the sense of the eigenvalues, denoted by {*An*}*n* ∼*λ (f, D)*, if the following limit relation holds for all *F* ∈ C0*(*C*)*:

$$\lim\_{n \to \infty} \frac{1}{d\_{\mathfrak{n}}} \sum\_{j=1}^{d\_{\mathfrak{n}}} F(\lambda\_j(A\_{\mathfrak{n}})) = \frac{1}{\mu\_q(D)} \int\_D \frac{\sum\_{l=1}^s F(\lambda\_l(f(t)))}{s} d\mathfrak{t},\tag{4}$$

where *λj (An)*, *j* = 1*,...,dn* are the eigenvalues of *An* and *λi(f )*, *i* = 1*,...,s* are the eigenvalues of *f* . We say that *f* is the (spectral) symbol of the matrix-sequence {*An*}*n*.

If *f* is smooth enough and the matrix-size of *An* is sufficiently large, then the limit relation (4) has the following informal meaning: a first set of *dn/s* eigenvalues of *An* is approximated by a sampling of *λ*1*(f )* on a uniform equispaced grid of the domain *D*, a second set of *dn/s* eigenvalues of *An* is approximated by a sampling of *λ*2*(f )* on a uniform equispaced grid of the domain *D*, and so on, up to few outliers.

In general, understanding whether a matrix-sequence admits a symbol and how to compute it is not an easy task. On the other hand, any "reasonable" approximation of partial differential equations by local methods leads to matrix-sequences that are in the so-called generalized locally Toeplitz (GLT) algebra, and so admit a symbol [10–12]. The IgA discretization of our curl-div problem (3) fits in this frame.

#### *2.2 B-Splines*

For *p* ≥ 0 and *n* ≥ 1, consider the uniform knot sequence

$$\xi\_1 = \dots = \xi\_{p+1} := 0 < \xi\_{p+2} < \dots < \xi\_{p+n} < 1 =: \xi\_{p+n+1} = \dots = \xi\_{2^p + n + 1},$$

where *ξi*+*p*+1 := *i n* , *i* = 0*,...,n*. This knot sequence allows us to define *n* + *p* B-splines of degree *p*. Let *χI* denote the characteristic function on the interval *I* .

**Definition 2** The B-splines of degree *p* over a uniform mesh of [0*,* 1], consisting of *n* intervals, are denoted by *Np i* : [0*,* 1] → R, *i* = 1*,...,n* + *p*, and defined recursively as follows: for 1 ≤ *i* ≤ *n* + 2*p*,

$$N\_i^0(\mathbf{x}) := \chi\_{[\xi^i, \xi\_{i+1})}(\mathbf{x});$$

for 1 ≤ *k* ≤ *p* and 1 ≤ *i* ≤ *n* + 2*p* − *k*,

$$N\_l^k(\mathbf{x}) := \frac{\mathbf{x} - \xi\_l}{\xi\_{l+k} - \xi\_l} N\_l^{k-1}(\mathbf{x}) + \frac{\xi\_{l+k+1} - \mathbf{x}}{\xi\_{l+k+1} - \xi\_{l+1}} N\_{l+1}^{k-1}(\mathbf{x}),$$

where a fraction with zero denominator is assumed to be zero.

It is well known (see e.g. [1]) that the B-splines *Np i* , *i* = 1*,...,n* + *p*, form a basis, and

$$N\_l^p(0) = N\_l^p(1) = 0, \quad i = 2, \ldots, n+p-1. \tag{5}$$

The central B-splines *Np i* , *i* = *p* + 1*,...,n*, are uniformly shifted and scaled versions of a single shape function, the so-called cardinal B-spline *φp* : R → R,

$$\phi\_0(t) := \chi\_{[0,1)}(t), \quad \phi\_p(t) := \frac{t}{p}\phi\_{p-1}(t) + \frac{p+1-t}{p}\phi\_{p-1}(t-1), \quad p \ge 1.$$

More precisely, we have

$$N\_i^{p}(\mathbf{x}) = \phi\_p(n\mathbf{x} - i + p + 1), \quad i = p + 1, \ldots, n.$$

The cardinal B-spline *φp* is a *Cp*−1 function which is locally supported on the interval [0*, p* + 1].

Finally, we recall the definition of tensor-product B-splines.

**Definition 3** The tensor-product B-splines of bi-degree *p* := *(p*1*, p*2*)* over a uniform mesh of [0*,* 1] 2, consisting of *n* := *(n*1*, n*2*)* intervals in each direction, are denoted by *Np i* : [0*,* 1] 2 → R, *i* = **1***,..., n* + *p*, and defined as

$$N\_{\boldsymbol{\ell}}^{\boldsymbol{p}} := N\_{\boldsymbol{i}\_1}^{\boldsymbol{p}\_1} \otimes N\_{\boldsymbol{i}\_2}^{\boldsymbol{p}\_2},$$

where **1** := *(*1*,* 1*)* and *i* := *(i*1*, i*2*)* ∈ N2.

We define the tensor-product spline space S*p n* as

$$\mathbb{S}\_{\mathfrak{n}}^{p} := \text{span}\left\{ N\_{l}^{p} : l = \mathbf{2}, \dots, n + p - 1 \right\}.\tag{6}$$

Note that all the elements of this space vanish at the boundary of [0*,* 1] 2; see (5). Hence, the space incorporates homogeneous Dirichlet boundary conditions.

#### **3 Spectral Analysis of Isogeometric Discretizations in 2D**

Suppose that the physical domain can be described by a global geometry map, *G* := [*G*1*, G*2] *T* , *G* : B → , which is invertible in the parametric domain B := [0*,* 1] 2 and satisfies *G(∂)* B = *∂*. Let

$$\mathbb{V}\_{h} = \text{span}\left\{ \phi\_{l\_{1},l\_{2}}^{p,1}, \phi\_{j\_{1},j\_{2}}^{p,2} : i\_{l}, j\_{l} = 2, \dots, n+p-1; \ l = 1,2 \right\},\tag{7}$$

where

$$
\phi\_{i\_1,i\_2}^{p,1} := \begin{bmatrix} \varphi\_{i\_1,i\_2} \\ 0 \end{bmatrix}, \quad \phi\_{j\_1,j\_2}^{p,2} := \begin{bmatrix} 0 \\ \varphi\_{j\_1,j\_2} \end{bmatrix},
$$

and for *kl* ∈ {*il, jl*}, *l* = 1*,* 2,

$$
\varphi\_{k\_1,k\_2}(\mathbf{x}\_1,\mathbf{x}\_2) := \hat{\varphi}\_{k\_1,k\_2}(\mathbf{G}^{-1}(\mathbf{x}\_1,\mathbf{x}\_2)) = \hat{\varphi}\_{k\_1,k\_2}(\hat{\mathbf{x}}\_1,\hat{\mathbf{x}}\_2), \quad (\mathbf{x}\_1,\mathbf{x}\_2) = \mathbf{G}(\hat{\mathbf{x}}\_1,\hat{\mathbf{x}}\_2).
$$

Then, we set *ϕ*ˆ*k*1*,k*2 = *Np k*1 ⊗ *Np k*2 , i.e., the tensor-product B-splines in (6). For simplicity of notation, we have taken *n*1 = *n*2 = *n* and *p*1 = *p*2 = *p*. Also note that

$$\begin{split} \nabla \varphi\_{k\_{1},k\_{2}} &= (J\_{G})^{-T} \nabla (N\_{k\_{1}}^{p} \otimes N\_{k\_{2}}^{p}) \\ &= \frac{1}{\det(J\_{G})} \begin{bmatrix} \frac{\partial G\_{2}}{\partial \hat{\boldsymbol{x}}\_{2}} (N\_{k\_{1}}^{p})' \otimes N\_{k\_{2}}^{p} - \frac{\partial G\_{2}}{\partial \hat{\boldsymbol{x}}\_{1}} N\_{k\_{1}}^{p} \otimes (N\_{k\_{2}}^{p})' \\ -\frac{\partial G\_{1}}{\partial \hat{\boldsymbol{x}}\_{2}} (N\_{k\_{1}}^{p})' \otimes N\_{k\_{2}}^{p} + \frac{\partial G\_{1}}{\partial \hat{\boldsymbol{x}}\_{1}} N\_{k\_{1}}^{p} \otimes (N\_{k\_{2}}^{p})' \end{bmatrix}. \end{split}$$

where

$$J\_G := \begin{bmatrix} \frac{\partial G\_1}{\partial \hat{\boldsymbol{x}}\_1} \frac{\partial G\_1}{\partial \hat{\boldsymbol{x}}\_2} \\ \frac{\partial G\_2}{\partial \hat{\boldsymbol{x}}\_1} \frac{\partial G\_2}{\partial \hat{\boldsymbol{x}}\_2} \end{bmatrix}.$$

In the following, we start by discussing the coefficient matrices arising from the IgA discretization of a generalized Poisson problem. Then, we construct the coefficient matrices related to the IgA discretization of our curl-div problem (3) using (7), and we perform a spectral analysis.

#### *3.1 Matrices Related to a Generalized Poisson Problem*

Let us focus on the following bivariate generalized Poisson operator:

$$
\mathcal{L}\!Ku := -\nabla \cdot K\nabla u,\tag{8}
$$

where *K* : → R2×2, and consider homogeneous Dirichlet boundary conditions, i.e., *u* = 0 on *∂*. From [8] we know that the Galerkin discretization of (8) using one component of the space (7) leads to the coefficient matrix A*p,K n,G* defined by

$$\left[\mathcal{R}\_{\mathfrak{n},G}^{p,K}\right]\_{l,j} := \int\_{\widehat{\Omega}} \left[\nabla (N\_{j\_1+1}^p \otimes N\_{j\_2+1}^p)^T \, K\_{\mathcal{G}} \, \nabla (N\_{l\_1+1}^p \otimes N\_{l\_2+1}^p)\right] |\det(J\_{\mathcal{G}})|,$$

where

$$K\_{\mathbf{G}} := (J\_{\mathbf{G}})^{-1} K(\mathbf{G}) (J\_{\mathbf{G}})^{-T}.$$

It has been proved in [8] that such matrices admit a spectral distribution according to Definition 1. To this end, let us define

$$H\_p := \begin{bmatrix} \mathfrak{s}\_p \otimes \mathfrak{m}\_p & \mathfrak{a}\_p \otimes \mathfrak{a}\_p \\ \mathfrak{a}\_p \otimes \mathfrak{a}\_p & \mathfrak{m}\_p \otimes \mathfrak{s}\_p \end{bmatrix},$$

with

$$\begin{aligned} \mathfrak{m}\_p(\theta) &:= \phi\_{2^{p+1}}(p+1) + 2 \sum\_{k=1}^p \phi\_{2^{p+1}}(p+1-k) \cos(k\theta), \\ \mathfrak{a}\_p(\theta) &:= -2 \sum\_{k=1}^p \phi'\_{2^{p+1}}(p+1-k) \sin(k\theta), \\ \mathfrak{s}\_p(\theta) &:= -\phi\_{2^{p+1}}'(p+1) - 2 \sum\_{k=1}^p \phi\_{2^{p+1}}''(p+1-k) \cos(k\theta). \end{aligned}$$

**Theorem 1** *Let G be a regular geometry map, i.e., G* ∈ *C*1*(*[0*,* 1] 2*) and* det*(JG)* = 0 *in* [0*,* 1] 2*, and let K be a symmetric matrix. Then, the matrix-sequence* {A*p,K n,G* }*n with n* = *(n, n) is distributed, in the sense of the eigenvalues, like the function*

$$f\_{\mathbf{G}}^{p,K}(\hat{\mathbf{x}},\boldsymbol{\theta}) := \left[\mathbb{I}\ \mathbf{l}\right] \left(\left|\det(J\_{\mathbf{G}}(\hat{\mathbf{x}}))\right| K\_{\mathbf{G}}(\hat{\mathbf{x}}) \circ H\_{p}(\boldsymbol{\theta})\right) \left[\mathbb{I}\ \mathbf{l}\right]^{T},\tag{9}$$

*where x*ˆ ∈ [0*,* 1] 2*, θ* ∈ [−*π, π*] 2*, and* ◦ *is the Hadamard matrix product.*

We refer the reader to [8, 9] for a detailed discussion about the symbol (9).

#### *3.2 Matrices Related to Our Curl-Div Problem*

We can reformulate (1) in 2D as

$$\mathcal{L}\_{\alpha,\beta}u = \alpha \begin{bmatrix} \frac{\partial^2 u\_2}{\partial x\_1 \ge \varepsilon\_2} - \frac{\partial^2 u\_1}{\partial x\_2^2} \\ \frac{\partial^2 u\_1}{\partial x\_1 \ge \varepsilon\_2} - \frac{\partial^2 u\_2}{\partial x\_1^2} \end{bmatrix} - \beta \begin{bmatrix} \frac{\partial^2 u\_2}{\partial x\_1 \ge \varepsilon\_2} + \frac{\partial^2 u\_1}{\partial x\_1^2} \\ \frac{\partial^2 u\_1}{\partial x\_1 \ge \varepsilon\_2} + \frac{\partial^2 u\_2}{\partial x\_2^2} \end{bmatrix},\tag{10}$$

where *u(x*1*, x*2*)* := [*u*1*(x*1*, x*2*), u*2*(x*1*, x*2*)*] *T* . When discretizing the weak form (3) using the space (7) we arrive at the 2 × 2 block matrix

$$
\mathcal{H}\_{\mathfrak{n},G}^{p,\alpha,\beta} := \alpha \begin{bmatrix}
\mathcal{R}\_{\mathfrak{n},11}^{p,\text{curl}} \, \mathcal{R}\_{\mathfrak{n},12}^{p,\text{curl}} \\
\mathcal{R}\_{\mathfrak{n},21}^{p,\text{curl}} \, \mathcal{R}\_{\mathfrak{n},22}^{p,\text{curl}}
\end{bmatrix} + \beta \begin{bmatrix}
\mathcal{R}\_{\mathfrak{n},11}^{p,\text{div}} \, \mathcal{R}\_{\mathfrak{n},12}^{p,\text{div}} \\
\mathcal{R}\_{\mathfrak{n},21}^{p,\text{div}} \, \mathcal{R}\_{\mathfrak{n},22}^{p,\text{div}}
\end{bmatrix}.
$$

The blocks related to the curl-curl operator *(*∇×·*,* ∇×·*)* are given by

! A*p,*curl *n,*11 " *i,j* = B ! −*∂G*1 *∂x*ˆ2 *(Np j*1+1*)* ⊗ *Np j*2+1 + *∂G*1 *∂x*ˆ1 *Np j*1+1 ⊗ *(Np j*2+1*)* " ! −*∂G*1 *∂x*ˆ2 *(Np i*1+1*)* ⊗ *Np i*2+1 + *∂G*1 *∂x*ˆ1 *Np i*1+1 ⊗ *(Np i*2+1*)* " 1 | det*(JG)*| *,* ! A*p,*curl *n,*12 " *i,j* = − B ! *∂G*2 *∂x*ˆ2 *(Np j*1+1*)* ⊗ *Np j*2+1 − *∂G*2 *∂x*ˆ1 *Np j*1+1 ⊗ *(Np j*2+1*)* " ! −*∂G*1 *∂x*ˆ2 *(Np i*1+1*)* ⊗ *Np i*2+1 + *∂G*1 *∂x*ˆ1 *Np i*1+1 ⊗ *(Np i*2+1*)* " 1 | det*(JG)*| *,* ! A*p,*curl *n,*22 " *i,j* = B ! *∂G*2 *∂x*ˆ2 *(Np j*1+1*)* ⊗ *Np j*2+1 − *∂G*2 *∂x*ˆ1 *Np j*1+1 ⊗ *(Np j*2+1*)* " ! *∂G*2 *∂x*ˆ2 *(Np i*1+1*)* ⊗ *Np i*2+1 − *∂G*2 *∂x*ˆ1 *Np i*1+1 ⊗ *(Np i*2+1*)* " 1 | det*(JG)*| *,*

and A*p,*curl *n,*21 = A*p,*curl *n,*12 . Note that all those blocks are symmetric matrices. Similarly, the blocks related to the div-div operator *(*∇·*,* ∇·*)* are given by (see also (10))

$$\mathcal{R}^{p, \text{div}}\_{\mathfrak{n}, 11} = \mathcal{R}^{p, \text{curl}}\_{\mathfrak{n}, 22}, \quad \mathcal{R}^{p, \text{div}}\_{\mathfrak{n}, 12} = \mathcal{R}^{p, \text{div}}\_{\mathfrak{n}, 21} = -\mathcal{R}^{p, \text{curl}}\_{\mathfrak{n}, 12}, \quad \mathcal{R}^{p, \text{div}}\_{\mathfrak{n}, 22} = \mathcal{R}^{p, \text{curl}}\_{\mathfrak{n}, 11}.$$

In the next subsection we compute the symbol of the matrix-sequence {A*p,α,β n,G* }*n*.

#### *3.3 Spectral Symbol of Curl-Div Matrices* A*p,α,β n,G*

We are now ready for the main contribution of the paper: we show that the matrixsequence {A*p,α,β n,G* }*n* admits a spectral distribution according to Definition 1. This extends the symbol computation in [15] to the case of non-trivial geometry.

**Theorem 2** *Let G be a regular geometry map, i.e., G* ∈ *C*1*(*[0*,* 1] 2*) and* det*(JG)* = 0 *in* [0*,* 1] 2*. Then, the matrix-sequence* {A*p,α,β n,G* }*n with n* = *(n, n) is distributed, in the sense of the eigenvalues, like the* 2 × 2 *matrix-valued function*

$$f\_{G}^{p,\alpha,\beta}(\hat{\mathfrak{x}},\theta) := \alpha f\_{G}^{p,\text{curl}}(\hat{\mathfrak{x}},\theta) + \beta f\_{G}^{p,\text{div}}(\hat{\mathfrak{x}},\theta),\tag{11}$$

*where x*ˆ ∈ [0*,* 1] 2*, θ* ∈ [−*π, π*] 2*, and*

$$\begin{split} f\_{\boldsymbol{G}}^{p,\operatorname{curl}}(\hat{\mathbf{x}},\boldsymbol{\theta}) &:= \frac{1}{|\det(J\_{\boldsymbol{G}}(\hat{\mathbf{x}}))|} \boldsymbol{J}\_{\boldsymbol{G}}(\hat{\mathbf{x}}) \, \boldsymbol{P} \, H\_{p}(\boldsymbol{\theta}) \, \boldsymbol{P}^{T} \left(J\_{\boldsymbol{G}}(\hat{\mathbf{x}})\right)^{T}, \quad \boldsymbol{P} := \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix}, \\\ f\_{\boldsymbol{G}}^{p,\operatorname{div}}(\hat{\mathbf{x}},\boldsymbol{\theta}) &:= |\det(J\_{\boldsymbol{G}}(\hat{\mathbf{x}}))| \, (J\_{\boldsymbol{G}}(\hat{\mathbf{x}}))^{-T} \, H\_{p}(\boldsymbol{\theta}) \, (J\_{\boldsymbol{G}}(\hat{\mathbf{x}}))^{-1}. \end{split}$$

*Proof* From (10) it follows that the block A*p,*curl *n,*11 corresponds to the isogeometric discretization of −*∂*2*u*1 *∂x*2 2 . By means of a direct computation we can verify that Theorem 1, with

$$K = \begin{bmatrix} 0 \ 0 \\ 0 \ 1 \end{bmatrix}.$$

ensures that the matrix-sequence {A*p,*curl *n,*11 }*n* is distributed in the sense of the eigenvalues like the entry *(*1*,* 1*)* of the matrix *f p,*curl *G* . The same argument (using a suitable matrix *K*) can also be applied to the remaining blocks. Then, it can be checked that all the considered blocks satisfy the hypotheses of [10, Theorem 5], which implies that A*p,α,β n,G* is similar, via a proper permutation matrix, to a matrix T*p,α,β n,G* such that the matrix-sequence {T*p,α,β n,G* }*n* has its symbol given by (11). %&

In the context of IgA, the geometry map *G* is expressed in terms of the same B-spline basis as used for the discretization space. However, as can be seen from the proof, the spectral result in the above theorem holds for any (smooth enough) geometry map.

Finally, we remark that the *p*-dependence of the symbol in (11) is completely captured by the matrix *Hp(θ)*. As described in Sect. 3.1 this matrix also appears in the symbol expression of a generalized Poisson problem; its properties have been discussed in [5, 8].

#### **4 Numerical Example**

In this section we numerically illustrate the spectral results obtained in Sect. 3.3, using the same test problem as in [15, Sect. 5]. More precisely, we consider (3) defined on a quarter of an annulus,

$$\mathfrak{U}\mathfrak{D} = \{ (\mathbf{x}\_1, \mathbf{x}\_2) \in \mathbb{R}^2 \colon r^2 < \mathbf{x}\_1^2 + \mathbf{x}\_2^2 < \mathbf{R}^2, \ \mathbf{x}\_1 > \mathbf{0}, \ \mathbf{x}\_2 > \mathbf{0} \}, \quad r = 1, \quad \mathbf{R} = 4, \dots$$

with

$$\mathbf{G}(\hat{\boldsymbol{\alpha}}\_{1},\hat{\boldsymbol{\alpha}}\_{2}) = \begin{cases} \boldsymbol{x}\_{1} = \left[r + \hat{\boldsymbol{\alpha}}\_{1}(\boldsymbol{R} - r)\right] \cos\left(\frac{\pi}{2}\hat{\boldsymbol{\alpha}}\_{2}\right) \\ \boldsymbol{x}\_{2} = \left[r + \hat{\boldsymbol{\alpha}}\_{1}(\boldsymbol{R} - r)\right] \sin\left(\frac{\pi}{2}\hat{\boldsymbol{\alpha}}\_{2}\right) \end{cases} \quad (\hat{\boldsymbol{\alpha}}\_{1},\hat{\boldsymbol{\alpha}}\_{2}) \in \left[0, 1\right]^{2} \dots$$

Let us fix *n* := *(n, n)* ∈ N2, *p* := *(p, p)* ∈ N2 and *m* ∈ N2 such that *m*2 = *n* + *p* − **2**. We start by defining two equispaced grids on [0*,* 1] 2 and [0*, π*] 2:

$$\text{If } x\_j := \frac{j}{m - 1}, \quad \theta\_k := \frac{k\pi}{m - 1}, \quad j, k = \mathbf{0}, \dots, m - 1.$$

Then, we denote by *\*i* the set of all evaluations of *λi(f p,α,β G )* on := {*(xj , θ k)*, *j, k* = **0***,..., m* − **1**} for a fixed *i* ∈ {1*,* 2}. Note that it suffices to consider only [0*, π*] 2 because the symbol (11) is symmetric on [−*π, π*] 2, and hence also its eigenvalue functions.

In Fig. 1 we numerically check relation (11) by comparing the eigenvalues of A*p,α,β n,G* with the values collected in *\** = {*\**1*, \**2}, ordered in ascending way, for *α* = 1 and *β* = 0*.*1. We observe that, in a complete agreement with the theory, the considered sampling of *λi(f p,α,β G )*, *i* = 1*,* 2, describes quite accurately the behavior of the eigenvalues of A*p,α,β n,G* , also for relatively small matrix-sizes, up to few outliers.

#### **5 Conclusions**

We have analyzed the spectral properties of matrix-sequences arising from isogeometric Galerkin methods for weighted curl-div operators on general planar domains, considering a non-trivial geometry map. More precisely, we have shown that an (asymptotic) spectral distribution exists and it is compactly described by a 2 × 2 spectral symbol. In other words, the eigenvalues of the matrices we are dealing with can be approximated accurately by a uniform sampling of the two eigenvalue functions of the 2 × 2 symbol matrix. The symbol depends on the characteristic parameters of the problem and on the geometry of the physical domain. Its formal

**Fig. 1** Comparison of the eigenvalues of A*p,α,β n,G* (open circle) with *\** = {*\**1*, \**2} collecting uniform samples of *λi(f p,α,β G )*, *i* = 1*,* 2 (asterisk), ordered in ascending way, varying both *n* and *p*, and fixing *α* = 1 and *β* = 0*.*1. (**a**) *p* = 3, *n* = 15. (**b**) *p* = 3, *n* = 35. (**c**) *p* = 4, *n* = 14. (**d**) *p* = 4, *n* = 34. (**e**) *p* = 5, *n* = 13. (**f**) *p* = 5, *n* = 33

structure nicely mimics the structure of the differential problem. The numerical results show a very good matching between the true eigenvalues and the estimates provided by the symbol, already for relatively small matrix-sizes.

The convergence of iterative solvers for linear systems strongly depends on the spectral behavior of the corresponding coefficient matrices. Since the symbol gives a precise description of the spectrum of the curl-div matrix A*p,α,β n,G* , it could be helpful in the design of good preconditioners that lead to better performance than current solution strategies, like the one in [15, Sect. 5].

**Acknowledgements** This work was partially supported by the INdAM research group GNCS, by the MIUR-DAAD Joint Mobility 2017 Programme through the project "ATOMA", and by the MIUR Excellence Department Project awarded to the Department of Mathematics, University of Rome Tor Vergata (CUP E83C18000100006).

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Performance of Preconditioners for Large-Scale Simulations Using Nek5000**

**N. Offermans, A. Peplinski, O. Marin, E. Merzari, and P. Schlatter**

### **1 Introduction**

The preconditioning of elliptic problems characterized by the propagation of information at infinite speed over the domain is a numerically challenging task. We study the case of the Poisson equation arising from the numerical resolution of the incompressible Navier–Stokes equations by operator splitting. We consider Nek5000, a code based on the spectral element method, as our framework. The current preconditioning strategy is based on an additive Schwarz method, which combines a domain decomposition method [5] and a so-called coarse grid problem [10]. The first step consists in solving directly local overlapping Poisson problems and is easily parallelizable. The second step corresponds to a Poisson-like problem over the whole domain and is hard to scale because of its relatively low number of degrees of freedom and the bottleneck induced by global communication.

A scalable solver for the coarse grid problem is critical to ensure strong scaling of the code. Existing strategies include a direct solution method similar to a Cholesky decomposition, called *XXT* [14], and an algebraic multigrid (AMG) solver [6]. While the first choice works well for relatively small problems (typically *<*100*,*000 spectral elements on *<*10*,*000 cores), the second option is preferred for large scale simulations. The current AMG solver, which we will denote as the *in-house AMG*, is fast and scales well [11]. It has been shown that the use of AMG can speed up

e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
N. Offermans (-) · A. Peplinski · P. Schlatter

Linné FLOW Centre and Swedish e-Science Research Centre (SeRC), KTH Mechanics, Stockholm, Sweden

O. Marin · E. Merzari Argonne National Laboratory, Lemont, IL, USA e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
S. J. Sherwin et al. (eds.), *Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018*, Lecture Notes in Computational Science and Engineering 134, https://doi.org/10.1007/978-3-030-39647-3\_20

large-scale simulations by up to 10%. In addition, the *XXT* has been designed for optimal performance on a number of cores which is a power of 2, whereas the AMG is insensitive to this parameter.

However, the AMG solver requires a setup phase, performed once for each mesh, by an external and serial code. Besides inducing an unwanted overhead, it also limits the use of the *in-house AMG* solver in the framework of mesh refinement, which is the main motivation for this work. Therefore, we propose to replace the *in-house AMG* by *BoomerAMG*, a parallel AMG solver for arbitrary unstructured grids from the *hypre* library for linear algebra [1, 4, 8]. *BoomerAMG* offers a number of parallel algorithms for the coarsening, interpolation and smoothing steps of the AMG setup, to accommodate various types of problems, meshes and architectures. The *BoomerAMG* solver will be tested in terms of scalability and time to solution.