Data structures and algorithms for high-dimensional structured adaptive mesh refinement

https://doi.org/10.1016/j.advengsoft.2014.12.001Get rights and content

Highlights

  • Structured anisotropic refinement of high-dimensional domains.

  • Recursive bisection of hyperrectangular mesh elements.

  • Linearized kd-tree for storage of the hierarchical mesh decomposition.

  • Examples and scalability studies of meshes in up to 6 dimensions.

  • Complete scheme scales better than nlogn, although worst case is O(n2).

Abstract

Spatial discretization of high-dimensional partial differential equations requires data representations that are of low overhead in terms of memory and complexity. Uniform discretization of computational domains quickly grows out of reach due to an exponential increase in problem size with dimensionality. Even with spatial adaptivity, the number of mesh data points can be unnecessarily large if care is not taken as to where refinement is done. This paper proposes an adaptive scheme that generates the mesh by recursive bisection, allowing mesh blocks to be arbitrarily anisotropic to allow for fine structures in some directions without over-refining in those directions that suffice with less refinement. Within this framework, the mesh blocks are organized in a linear kd-tree with an explicit node index map corresponding to the hierarchical splitting of internal nodes. Algorithms for refinement, coarsening and 2:1 balancing of a mesh hierarchy are derived. To demonstrate the capabilities of the framework, examples of generated meshes are presented and the algorithmic scalability is evaluated on a suite of test problems. In conclusion, although the worst-case complexity of sorting the nodes and building the node map index is n2, the average runtime scaling in the studied examples is no worse than nlogn.

Introduction

Structured adaptive mesh refinement (SAMR) is an active area of research within the scientific computing community [1]. By adjusting the resolution of the computational mesh dynamically to features in the solution or the computational domain, widely varying scales of resolution can be represented simultaneously. In effect, the computational efficiency of a simulation is improved, possibly by orders of magnitude, allowing for larger computations and/or shorter execution times due to a reduction in the number of gridpoints [2]. For problems in two and three dimensions, there are efficient algorithms and data structures available, relying on quad/octrees for structuring the mesh blocks [2], [3], [4], [5]. However, extending quadtrees and octrees to higher dimensional trees is problematic since they yield a fan-out of 2D nodes at every branch, which leads to an exponential increase in the potential number of tree nodes to handle. The contribution of this paper is a framework capable of generating and propagating meshes of arbitrary dimensionality. Our approach is based on recursive bisection and generates far fewer mesh nodes compared to 2D-trees of corresponding refinement.

In order to construct a practical refinement scheme that works well even in higher dimensions, the framework presented in this paper is built on a structured block-based refinement strategy [1] allowing blocks to be refined anisotropically. The mesh nodes and their mutual relationships are maintained in a kd-tree [6]. With anisotropic refinement, a block is not restricted to be refined equally in all dimensions, potentially leading to a more efficient discretization in terms of the number of created blocks since mesh blocks are refined only in the dimensions in which they would benefit from finer resolution. The grid is refined successively by dividing blocks in half, dimension by dimension. If a block needs refinement in more than one dimension, this is done by subsequent division in several steps. This paper does not consider details regarding error estimation and how to determine when refinement/coarsening is required. For the anisotropic refinement strategy to be useful though, the error estimator must be able to detect the discretization error per dimension. An example of such an error estimator is given in [7].

A kd-tree is a binary tree representation of a hierarchical subdivision of a D-dimensional hyperrectangle by recursive bisection [6], [8], [9], [10], [11]. The interior nodes of a kd-tree represent hyperplane cuts, aligned with the Cartesian coordinate axes, and the leaf nodes contain the actual data. In a general kd-tree, a cut can be placed anywhere along the split dimension of a block. However, by restricting the cuts to always be placed in the middle of the block (which we do by imposing a halving of the blocks on each refinement) the scheme is simplified significantly. Furthermore, implementations of kd-trees usually assign split dimensions to nodes cyclically, such that a node at tree level l is split in dimension (lmodD). In the implementation presented here, this restriction is relaxed and nodes are allowed to be split arbitrarily without any intermediate refinement in the other dimensions. Thus, blocks can become as elongated as is needed and in principle there is no restriction on the aspect ratio of the blocks.

It is often motivated for reasons of efficiency and memory requirements that a pointer representation of a tree structure should be avoided. By storing locality information in each mesh node and structuring the nodes in a linear order according to this information, the internal structure of a tree is available implicitly and no pointers are needed for searching and navigating it [2], [3], [4], [12]. With the leaf nodes stored in linear order (e.g. the Morton order space-filling curve [13]), tree search is replaced by binary search, which is further advantageous in terms of search complexity; a tree representing an adaptively refined mesh is potentially very unbalanced with a search complexity approaching O(n) in the number of leaf nodes, whereas binary search in the linear representation is always O(log2n) [4]. This paper extends previous work by other authors [2], [4] and builds upon linear Morton order trees. The data structures and algorithms are generalized to enable an extension to higher dimensionalities. However, the anisotropic node refinement does not map directly to an efficient consecutive order of elements solely from the locality information in the leaf nodes. In order to alleviate this, the linear tree representation is extended with a lightweight representation of the internal node structure (a node map index), corresponding to a hierarchical Morton order of nodes. The node map index is constructed at the same time as the nodes are sorted and stored in a compact array representation for efficient traversal.

The remainder of this paper is organized as follows. In Section 2, some related work is discussed. A key concept in our implementation is the notion of location codes, which is described in Section 3. In Section 4 the construction of a hierarchical Morton order index is derived together with some implementation details. The algorithms for tree search are discussed in Section 5, followed by the adaptive mesh algorithms in Section 6. Finally, Section 7 gives some results of our implementation and Section 8 concludes the paper.

Section snippets

Related work

To our knowledge, there are no dynamically adaptive frameworks available that support anisotropic SAMR on hyperrectangular domains in D dimensions. Klöfkorn and Nolte described the implementation of the SPGrid interface of the DUNE framework [14], [15], which allows anisotropic structured grids of higher dimensionality but only handles static meshes. An alternative approach to tackling higher dimensions is to discretize the domain into structured simplex meshes [16], [17], [18], [19]. In [11],

Location codes

A mesh block represents a hyperrectangle (orthotope) subdomain of the discretized D-dimensional space. The location code [2], [4] of an orthotope is a unique identifier that encodes the location and refinement of the corresponding mesh element. The location code has two components; (1) a coordinate in space with respect to the orthotope anchor, and (2) refinement level information (per dimension). The anchor of an orthotope is chosen by convention to be the lowermost corner in every direction,

Hierarchical Morton order

In order to understand why it is problematic to omit the internal node hierarchy in anisotropic domains, consider the example in Fig. 4. If there is no information available regarding the hierarchical structure of node subdivisions in the kd-tree, all we can get is an ordering based directly on the location code information, corresponding to the Morton ordering in Fig. 4(a). Here, the ordering of atoms clearly violates the order of the structure of the mesh blocks, which is an undesirable

Search

The basic algorithms binarySearchMatching and binarySearchOverlapping perform binary search of the linear kd-tree with a location code as search key (the search routines are not listed here since they are straightforward to implement). Given a search key, it is translated to a linear index which is used for comparisons in the binary search, whereas the search key itself is used to indicate a search hit. In binarySearchMatching an identical match to the search key is looked for, whereas

Adaptive mesh algorithms

This section describes the key algorithms for mesh adaptation; refinement, coarsening and 2:1 balancing. All of these algorithms are of worst-case complexity O(n2), where n is the number of leaf nodes in the tree (thus equal to the number of mesh elements). The actual run-time complexity will depend on how unbalanced the refinement tree is and the number of mesh nodes that need to be updated. In practical computations on localized solutions, we expect only a small portion of the mesh nodes to

Mesh propagation

In order to study how the mesh adaptation performs, examples of different dimensionalities (2D and 4D) are constructed and studied. In each experiment, a function is generated and propagated in time, and the mesh is evolved with it. Here, we do not actually solve a PDE; instead, an analytic function to which we know the values at time t is used, and the mesh points are updated with the exact function values in each time instance. This approach eliminates the impact of all exterior factors, such

Conclusion

In this paper, an efficient scheme for managing adaptively refined meshes with anisotropic features was presented. The framework is built around a linear representation of kd-trees, with algorithms for refinement, coarsening and 2:1 balancing. A linear kd-tree data structure is proposed that can be seen as a generalization of the quadtree/octree data structures for higher dimensionalities, with added flexibility in terms of layout and shape of grid elements. Compared to a binary tree

References (26)

  • J.P. Jessee et al.

    An adaptive mesh refinement algorithm for the radiative transport equation

    J Comput Phys

    (1998)
  • J. Rantakokko et al.

    Parallel structured adaptive mesh refinement

  • C. Burstedde et al.

    p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees

    SIAM J Sci Comput

    (2011)
  • I. Gargantini

    An effective way to represent quadtrees

    Commun ACM

    (1982)
  • H. Sundar et al.

    Bottom-up construction and 2:1 balance refinement of linear octrees in parallel

    SIAM J Sci Comput

    (2008)
  • S.F. Frisken et al.

    Simple and efficient traversal methods for quadtrees and octrees

    J Graph Tools

    (2002)
  • J.L. Bentley

    Multidimensional binary search trees used for associative searching

    Commun ACM

    (1975)
  • Nissen A, Kormann K, Grandin M, Virta K. Stable difference methods for block-oriented adaptive grids. J Sci Comput...
  • R.A. Finkel et al.

    Quad trees: a data structure for retrieval on composite keys

    Acta Inform

    (1974)
  • J.H. Freidman et al.

    An algorithm for finding best matches in logarithmic expected time

    ACM Trans Math Softw

    (1977)
  • G. Gan et al.

    Data clustering: theory, algorithms and applications

    (2007)
  • K. Weiss et al.

    Bisection-based triangulations of nested hypercubic meshes

  • C. Chuanbo et al.

    Linear binary tree

  • Cited by (5)

    View full text