General Field Evaluation in High-Order Meshes on GPUs

TL;DR

This work tackles the challenge of evaluating PDE solutions at arbitrary points within large-scale, high-order curved meshes distributed across MPI. It presents a robust framework that combines global and processor-local maps, axis-aligned and oriented bounding boxes, and Newton's method with a trust region to identify the overlapping element and compute reference-space coordinates, with specialized GPU kernels to accelerate Find and Interpolate on tensor-product elements. Key contributions include the intra- and inter-processor data structures for candidate element/rank selection, AABB/OBB bounding strategies, a trust-region Newton inversion, and extensions to surface meshes along with a GPU-optimized implementation. The approach enables efficient general field evaluation for applications such as Lagrangian particle tracking, r-adaptivity remapping, and interdomain data exchange, delivering substantial speedups over CPU-only workflows and enabling scalable, high-fidelity simulations on modern HPC architectures.

Abstract

Robust and scalable function evaluation at any arbitrary point in the finite/spectral element mesh is required for querying the partial differential equation solution at points of interest, comparison of solution between different meshes, and Lagrangian particle tracking. This is a challenging problem, particularly for high-order unstructured meshes partitioned in parallel with MPI, as it requires identifying the element that overlaps a given point and computing the corresponding reference space coordinates. We present a robust and efficient technique for general field evaluation in large-scale high-order meshes with quadrilaterals and hexahedra. In the proposed method, a combination of globally partitioned and processor-local maps are used to first determine a list of candidate MPI ranks, and then locally candidate elements that could contain a given point. Next, element-wise bounding boxes further reduce the list of candidate elements. Finally, Newton's method with trust region is used to determine the overlapping element and corresponding reference space coordinates. Since GPU-based architectures have become popular for accelerating computational analyses using meshes with tensor-product elements, specialized kernels have been developed to utilize the proposed methodology on GPUs. The method is also extended to enable general field evaluation on surface meshes. The paper concludes by demonstrating the use of proposed method in various applications ranging from mesh-to-mesh transfer during r-adaptivity to Lagrangian particle tracking.

Figures (17)

• Figure 1: Schematic showing a finite element ($p=3$) mapped from the reference space to physical space.
• Figure 2: Key ingredients used to determine candidate elements that overlap a given point ${\bf x}^*$. (a) $p=3$ mesh with $N_E=64$ partitioned on to 4 MPI ranks. Different colors represent different ranks. (b) Global Cartesian mesh $\mathcal{M}_G$ spanning the domain of $\mathcal{M}$, (c) partitioned on to 4 MPI ranks. $\mathcal{M}_G$ is used to construct a map from its elements to MPI ranks corresponding to intersecting elements of $\mathcal{M}$. (d) The process local Cartesian mesh $\mathcal{M}_L$ on one of the ranks, which is used to map elements of $\mathcal{M}_L$ to corresponding intersecting elements of $\mathcal{M}$ on that rank. (e) Axis-aligned bounding boxes (AABBs) and (f) oriented bounding boxes (OBBs) for the elements of $\mathcal{M}_L$.
• Figure 3: (a) $\phi_i(r)$, $i=1 \dots 4$, Lagrange interpolants. (b) Cubic 1D function $u(r)$ with nodal coefficients $u_i=\{-2.4,0.8,1.4,-1\}$ in \ref{['eq_1D_function']}.
• Figure 4: Lagrange interpolants $\phi_i(r)$, $i=1\dots 4$ at GLL points (•) along with the piecewise linear bounds constructed with $M=8$ Chebyshev interval points (x). (a) $\phi_1(r)$, (b) $\phi_2(r)$, (c) $\phi_3(r)$, and (d) $\phi_4(r)$.
• Figure 5: Piecewise linear bounds on $u(r)$ constructed using (a) $M=6$, (b) $M=8$, and (c) $M=12$ Chebyshev interval points. Higher number of interval points lead to tighter bounds around $u(r)$.