Discretization of PDEs with variable coefficients using locally adaptive sparse grids

Elliptic partial differential equations with variable coefficients can be discretized on sparse grids. With prewavelets being L2‐orthogonal, one can apply the Ritz‐Galerkin discretization to obtain a linear equation system with O(N(logN)d−1)$O(N(\log N)^{d-1})$ unknowns. However, for several applications like partial differential equations with corner singularities or the high‐dimensional Schrödinger equation, locally adaptive grids are needed to obtain optimal convergence. Therefore, we introduce a new kind of locally adaptive sparse grid and a corresponding algorithm that allows solving the resulting finite element discretization equation with optimal complexity. These grids are constructed by local tensor product grids to generate adaptivity but still maintain a local unidirectional approach. First simulation results of a two‐dimensional Helmholtz problem are presented.


INTRODUCTION
Discretizing PDEs by finite elements on full grids leads to a computational amount of (  ), where  is the number of grid points in one direction .This so called curse of dimensionality can be weakened by sparse grids which construct a subspace of the classical finite element space on full grids, see ref. [1].With prewavelets being  2 -orthogonal, one can apply the Ritz-Galerkin discretization to obtain a linear equation system with ((log ) −1 ) unknowns [2].The interpolation error increases for sparse grids compared to full grids by only a small logarithmic factor, that is, ‖ −    ()‖ ∞ = (ℎ 2 log(ℎ −1 ) −1 ) with  = 2  , ℎ = 2 − and  being the level of refinement, see ref. [1].
For several applications like PDEs with corner singularities or the high-dimensional Schrödinger equation, locally adaptive grids are needed to obtain optimal convergence.There have been several approaches on solving PDEs on adaptive sparse grids.However, they were either restricted to constant coefficients, see refs.[3][4][5], or dimension  ≤ 3, see refs.[6][7][8].Therefore, we introduce a new kind of locally adaptive sparse grid that allows solving the resulting finite element discretization equation with optimal complexity by using an efficient algorithm based on ref. [2].

SPARSE GRID BASIC NOTATION AND DISCRETIZATION
The aim of this section is to give an overview about finite elements on sparse grids, prewavelets and a discretization of Problem 2.1.It is mainly based on refs.[2] and [9].All of these concepts shall be used later to discretize Problem 2.1 on locally adaptive sparse grids.
Let  ≥ 2 be the dimension of space and Ω = (0, 1)  .Furthermore, let  ∈  2 (Ω),  ∈ ( ∞ (Ω)) × and  ∈  ∞ (Ω),  ≥ 0 be given.Let us define By using tensor products of one-dimensional finite elements, one can construct sparse grids.To this end, define the one-dimensional grid where   = { |  = 1, … , 2  − 1} for  ∈ ℕ.The complementary index set and corresponding grid points are defined by To describe the position and the corresponding basis function of a grid point  ∈ , we define the depth function , the index function  and the mesh width ℎ  ∶= () by In Section 3 we will need these functions to describe adaptive sparse grids.We call   the space of piecewise linear functions of mesh size 2 − with  lin , being the corresponding nodal basis function at point 2 − .By using these functions, we define prewavelets  , ,  ∈ Ξ  by (see Figure 1) F I G U R E 1 Prewavelet basis functions of levels 1, 2 and 3.
To construct a tensor product basis, let us use a multi-index  = ( 1 , … ,   ) ∈ ℕ  0 and  = ( 1 , … ,   ) ∈   with Applying these abbreviations, we define the set of tensor product indices where  = ( 1 , … ,   ).Furthermore, we can define a uniformly refined sparse grid   ⊂  as A regular two-dimensional sparse grid is depicted in Figure 2.With these constructions the tensor product vector spaces  lin and  prew   can be defined as prew The spaces  lin .
The corresponding prewavelet-based algorithm for the matrix-vector multiplication is presented in ref. [2].Numerical results, such as convergence of Discretization 2.4, for a three-dimensional problem on a curvilinear bounded domain and for a six-dimensional problem with variable coefficients can be found in refs.[2,10] and [9].

ADAPTIVE SPARSE GRIDS
Our aim is to extend the concepts and ideas of refs.[2] and [9] as presented in Section 2 to locally adaptive sparse grids.
The algorithm for the matrix-vector multiplication presented in ref. [2] is based on a unidirectional principle, global, matrix-free and parallelizable.All of these properties lead to an efficient algorithm with ((log ) −1 ) computations and an optimal convergence with respect to the approximation properties of sparse grids.To maintain this optimality of the algorithm, we need to construct locally adaptive sparse grids for prewavelets as described in the following section.
In order to describe grid points corresponding to prewavelet basis functions in an adaptive manner, we need to define multidimensional three-and five-point stencils.
The grid points defined by Three() are also called parents of .They are necessary for the basis transformation from nodal basis to prewavelet basis.This motivates our first criterion for the adaptive grid construction.Since the algorithm for the matrix-vector multiplication is supposed to use algorithmic paths globally, for example, stencil operators on all grid points of a fixed depth , we need to make sure that it is either locally full or locally sparse (Figure 3).We call  full (,) the locally full prewavelet grid w.r.t. the grid points  and .The next definition combines Definition 3.3 and the the five-point stencil.
Definition 3.4 (locally adaptive sparse grid for prewavelets).Let  be an adaptive sparse grid.Let  ⊃  be the smallest grid such that Five() ⊂  ∀ ∈ .
Furthermore, for all ,  ∈  the following statement must be fulfilled: The pair ( , ) is called locally adaptive sparse grid for prewavelets.
Notation 3.5.Let ( , ) be a locally adaptive sparse grid for prewavelets.From now on we will use the following notation:  ad ∶= ,  active ∶=  and  hanging ∶=  ad ⧵  active .
Note that only active nodes,  active , will describe degrees of freedom.
Definition 3.4 states that for all grid points, the grid points of the five-point stencil have to be part of the grid, but not necessarily as active nodes, here degrees of freedom.Furthermore, it requires the grid to be either locally full or locally sparse.
Remark 3.6.Another property is the so called local tensor product structure.This grid structure is only needed for the basic transformation from prewavelets to nodal basis functions.Hence it isn't necessary for the matrix-vector multiplication and can be omitted when solving eigenvalue problems.Further details about this structure shall be published in the future.

DISCRETIZING PDEs ON LOCALLY ADAPTIVE SPARSE GRIDS FOR PREWAVELETS
Using the definition of locally adaptive sparse grids for prewavelets, we can define a Galerkin discretization of Problem 2.1.

NUMERICAL RESULTS
As a first proof of concept, one numerical result is presented.To this end, consider the two-dimensional Helmholtz problem with variable coefficient The right-hand side  is chosen such that (, ) ∶= sin(5) sin() is a solution of (2).We solve this problem, by using Discretization 4.2 and applying the conjugate gradient method with a simple diagonal Jacobi preconditioner.
The condition number of the stiffness matrix, including the preconditioner, is (1), see the multilevel theory in ref. [11].Hence, even with higher DOF only a few CG-iterations are needed to obtain a small algebraic error, see Table 1.Furthermore, the number of hanging nodes stays relatively small in comparison to the number of DOF. Figure 7 shows that the convergence of Problem (2) on locally adaptive sparse grids is similar to the convergence of Problem (2) with (, ) = 0 on regular sparse grids.Since the solution  is smooth, we wouldn't expect faster convergence on adaptive sparse grids, as regular sparse grids are constructed as a good fit for smooth solutions.For solutions with corner-singularities, we would expect faster convergence on adaptive sparse grids than on regular sparse grids.Nevertheless, this shows that Discretization 4.2 with variable coefficients does not introduce any remarkable additional errors.In Figure 6, we can see a corresponding locally adaptive sparse grid with local tensor product structure, to a solution of Problem 2. The algorithm detects the frequency of the solution , see Figure 5, and applies higher refinement in -direction.

CONCLUSION
We introduced a new kind of locally adaptive sparse grid for prewavelets, along with a corresponding discretization that can be carried out by an optimal algorithm, as presented in ref. [2].As a first proof of concept, we solved a simple twodimensional Helmholtz problem.The simulation results showed the expected behavior, as for example, a bounded number of CG-iterations and no additional errors when solving a PDE with variable coefficients.In the future, we want to extend the theory to overcome inhomogeneous boundary conditions.This should widen the range of examples, such as problems with corner singularities.Furthermore, performance optimization of algorithms on sparse grids and application to highdimensional eigenvalue problems shall be implemented.

A C K N O W L E D G M E N T S
Open access funding enabled and organized by Projekt DEAL.

F
I G U R E 3 A locally adaptive sparse grid for prewavelets with tensor product structure.FI G U R E 4Example of a locally full grid (A) and a locally sparse grid (B) and an application of Discretization 4.2.

Figure
Figure 4A depicts an example of  full (,) .The five-point stencil Five() is necessary for the correct evaluation of the bilinear form,  ( ,   ) = ∑ ∈Five()

F I G U R E 5
(, ) = sin(5) sin() on an adaptive sparse grid.F I G U R E 6A locally adaptive sparse grid for the solution of Problem(2).F I G U R E 7 Discretization error  ∞ ∶= ‖ −  prew  active ‖ ∞ .
Let  ∈ ℝ be constant and  ∈ ℝ × a constant diagonal-matrix.Then, for all indices ,  ′ ,  ∈   and  ′ ∈   ′ (1)are equal.However, this equality does not hold for the adaptive versions of these spaces.Furthermore, the  2 -orthogonality of prewavelets(1)is not fulfilled for nodal basis functions  lin , .Hence, two different names for  lin   are used. I this paper we will extend the definition of  prew   to locally adaptive sparse grids  ad ⊂ .Now, the Galerkin discretization of Problem 2.1 using prewavelets on a regular sparse grid   is: Discretization 2.2.Find  prew   .Applying the  2 -orthogonality of prewavelets, that is, property (1), one can easily show the following lemma.Lemma 2.3.Ω (∇ , )  ∇  ′ , ′ +  ,   ′ , ′  = 0. Using Lemma 2.3 as a motivation leads to the following discretization.
Furthermore, we can extend Discretization 2.4 by using Definition 3.3 as a criterion for the  2 -orthogonality of prewavelets.Figure4shows how this discretization can be employed in a local approach for two grid points  and .Degrees of freedom (DOF), number of hanging nodes, number of conjugate gradient-iterations (CG and the discretization error. Lt   be a regular sparse grid and ,  ∈   .If  full (,) ⊂   , then | max((), ())| ≤  + 1.From Lemma 4.3 it follows that applying Discretization 4.2 on a regular sparse grid   leads to Discretization 2.4.