High-order Finite Element Discretizations Research Articles

Statically Condensed Iterated Penalty Method for High Order Finite Element Discretizations of Incompressible Flow

Statically Condensed Iterated Penalty Method for High Order Finite Element Discretizations of Incompressible Flow

A dynamical core based on a discontinuous Galerkin method for higher-order finite-element sea ice modeling

Abstract. The ability of numerical sea ice models to reproduce localized deformation features associated with fracture processes is key for an accurate representation of the ice dynamics and of dynamically coupled physical processes in the Arctic and Antarctic. Equally key is the capacity of these models to minimize the numerical diffusion stemming from the advection of these features to ensure that the associated strong gradients persist in time, without the need to unphysically re-inject energy for re-localization. To control diffusion and improve the approximation quality, we present a new numerical core for the dynamics of sea ice that is based on higher-order finite-element discretizations for the momentum equation and higher-order discontinuous Galerkin methods for the advection. The mathematical properties of this core are discussed, and a detailed description of an efficient shared-memory parallel implementation is given. In addition, we present different numerical tests and apply the new framework to a benchmark problem to quantify the advantages of the higher-order discretization. These tests are based on Hibler's viscous–plastic sea ice model, but the implementation of the developed framework in the context of other physical models reproducing a strong localization of the deformation is possible.

End-to-end GPU acceleration of low-order-refined preconditioning for high-order finite element discretizations

In this article, we present algorithms and implementations for the end-to-end GPU acceleration of matrix-free low-order-refined preconditioning of high-order finite element problems. The methods described here allow for the construction of effective preconditioners for high-order problems with optimal memory usage and computational complexity. The preconditioners are based on the construction of a spectrally equivalent low-order discretization on a refined mesh, which is then amenable to, for example, algebraic multigrid preconditioning. The constants of equivalence are independent of mesh size and polynomial degree. For vector finite element problems in H(curl) and H(div) (e.g., for electromagnetic or radiation diffusion problems), a specially constructed interpolation–histopolation basis is used to ensure fast convergence. Detailed performance studies are carried out to analyze the efficiency of the GPU algorithms. The kernel throughput of each of the main algorithmic components is measured, and the strong and weak parallel scalability of the methods is demonstrated. The different relative weighting and significance of the algorithmic components on GPUs and CPUs is discussed. Results on problems involving adaptively refined nonconforming meshes are shown, and the use of the preconditioners on a large-scale magnetic diffusion problem using all spaces of the finite element de Rham complex is illustrated.

Textbook Efficiency: Massively Parallel Matrix-Free Multigrid for the Stokes System

We employ textbook multigrid efficiency (TME), as introduced by Achi Brandt, to construct an asymptotically optimal monolithic multigrid solver for the Stokes system. The geometric multigrid solver builds upon the concept of hierarchical hybrid grids (HHG), which is extended to higher-order finite-element discretizations, and a corresponding matrix-free implementation. The computational cost of the full multigrid (FMG) iteration is quantified, and the solver is applied to multiple benchmark problems. Through a parameter study, we suggest configurations that achieve TME for both, stabilized equal-order, and Taylor-Hood discretizations. The excellent node-level performance of the relevant compute kernels is presented via a roofline analysis. Finally, we demonstrate the weak and strong scalability to up to $147,456$ parallel processes and solve Stokes systems with more than $3.6 \times 10^{12}$ (trillion) unknowns.

Model and computational advancements to full vectorial Maxwell model for studying fiber amplifiers

We present both modeling and computational advancements to a unique three-dimensional discontinuous Petrov–Galerkin finite element model for the simulation of laser amplification in a fiber amplifier. Our model is based on the time-harmonic Maxwell equations, and it incorporates both amplification via an active dopant and thermal effects via coupling with the heat equation. As a full vectorial finite element simulation, this model distinguishes itself from other fiber amplifier models that are typically posed as an initial value problem and make significantly more approximations. Our model supports co-, counter-, and bi-directional pumping configurations, as well as inhomogeneous and anisotropic material properties. The long-term goal of this modeling effort is to study nonlinear phenomena that prohibit achieving unprecedented power levels in fiber amplifiers, along with validating typical approximations used in lower-fidelity models. The high-fidelity simulation comes at the cost of a high-order finite element discretization with many degrees of freedom per wavelength. This is necessary to counter the effect of numerical pollution due to the high-frequency nature of the wave simulation. To make the computation more feasible, we have developed a novel longitudinal model rescaling, using artificial material parameters with the goal of preserving certain quantities of interest. Our numerical tests demonstrate the applicability and utility of this scaled model in the simulation of an ytterbium-doped, step-index fiber amplifier that experiences laser amplification and heating. We present numerical results for the nonlinear coupled Maxwell/heat model with up to 240 wavelengths.

Numerical Analysis for an Electro‐Magneto‐Mechanical Phenomenon with High‐Order Accurate Methods

AbstractThis paper establishes an axisymmetric model for a levitation device. Therein, the Maxwell equations are combined with the balance of linear momentum. Different possible formulations to describe the MAXWELL equations are presented and compared and discussed in the example. A high order finite element discretization using GALERKIN's method in space and the generalized NEWMARK‐α method in time are developed for the electro‐magneto‐mechanical approach. Several studies on spatial and temporal discretization with respect to convergence will be investigated. In addition, the boundary influences and the domain size with respect to the levitation device are also examined.

Preconditioning the Mass Matrix for High Order Finite Element Approximation on Tetrahedra

A preconditioner for the mass matrix for high order finite element discretization on tetrahedra is presented and shown to give a condition number that is independent of both the mesh size and the polynomial order of the elements. The preconditioner is described in terms of a new, high order basis which has the usual property whereby individual functions are associated with distinct geometric entities of the tetrahedron. It is shown that the basis enjoys the novel property that the resulting mass matrix is spectrally equivalent to its own diagonal with constants independent of h and p. Although the exposition is based on an explicit basis, the preconditioner can be applied to any choice of basis. In particular, the basis can be used to specify a basis-independent additive Schwarz method, meaning that, in order to apply the preconditioner to an alternative basis, one only need to implement an appropriate change of basis.

High-Order Accurate Methods for the Numerical Analysis of a Levitation Device

This paper establishes different axisymmetric and two-dimensional models for a levitation device. Therein, the Maxwell equations are combined with the balance of linear momentum. Different possible formulations to describe the Maxwell equations are presented and compared and discussed in the example. A high order finite element discretization using Galerkin’s method in space and the generalized Newmark $$-\alpha$$ method in time are developed for the electro-magneto-mechanical approach. Several studies on spatial and temporal discretization with respect to convergence will be investigated. In addition, the boundary influences and the domain size with respect to the levitation device are also examined.

High-order matrix-free incompressible flow solvers with GPU acceleration and low-order refined preconditioners

We present a matrix-free flow solver for high-order finite element discretizations of the incompressible Navier-Stokes and Stokes equations with GPU acceleration. For high polynomial degrees, assembling the matrix for the linear systems resulting from the finite element discretization can be prohibitively expensive, both in terms of computational complexity and memory. For this reason, it is necessary to develop matrix-free operators and preconditioners, which can be used to efficiently solve these linear systems without access to the matrix entries themselves. The matrix-free operator evaluations utilize GPU-accelerated sum-factorization techniques to minimize memory movement and maximize throughput. The preconditioners developed in this work are based on a low-order refined methodology with parallel subspace corrections, as described for diffusion problems in [1]. The saddle-point Stokes system is solved using block-preconditioning techniques, which are robust in mesh size, polynomial degree, time step, and viscosity. For the incompressible Navier-Stokes equations, we make use of projection (fractional step) methods, which require Helmholtz and Poisson solves at each time step. The performance of our flow solvers is assessed on several benchmark problems in two and three spatial dimensions.

A global digital volume correlation algorithm based on higher-order finite elements: Implementation and evaluation

We propose a DVC technique that is based on higher-order finite-element discretization of the displacement field and a global optimization procedure. We use curvature penalization to suppress non-physical fluctuations of the displacement field and resulting erroneous strain concentrations. The performance of the proposed method is compared to the commercial code Avizo using trabecular bone images and found to perform slightly better in most cases.In addition, we stress that the performance of a DVC method needs to be evaluated using double scans (zero strain), virtual deformation (imposed deformation) and real deformation. Double scans give insight into the presence of noise and artifacts whereas virtual deformation benchmarks allows evaluation of the performance without noise and artifacts. Investigation of the performance for actual deformed heterogeneous materials is needed for evaluation with noise, artifacts and non-zero strains.We show that both decreasing the resolution of the displacement field (increasing subvolume size) as well as (increasing) curvature penalization (regularization) have a similar effect on the performance of evaluated DVC methods: Decreasing the detrimental effect of noise, artifacts and interpolation errors, but also decreasing the sensitivity of a DVC method to displacement peaks, discontinuities and strain concentrations. The needed amount of regularization is a trade-off between accuracy and precision of the estimated strain fields and their resolution.The obtainable accuracy and precision of the estimated displacement fields are influenced by interpolation errors in the DVC procedure and the relative amount of detail, noise and artifacts in the images. Errors in the displacement field are typically magnified during the strain calculation. Based on the tests and subvolume sizes (16–50 voxels) in this study, the expected order of magnitude of the accuracy and precision is 0.1 micro-voxels and 1 milli-voxels for the displacements and 0.1 and 1 milli-strains of the strain fields.

Residual stresses in metal deposition modeling: Discretizations of higher order

This article addresses the research question if and how the finite cell method, an embedded domain finite element method of high order, may be used in the simulation of metal deposition to harvest its computational efficiency. This application demands for the solution of a coupled thermo-elasto-plastic problem on transient meshes within which history variables need to be managed dynamically on non-boundary conforming discretizations. To this end, we propose to combine the multi-level hp-method and the finite cell method. The former was specifically designed to treat high-order finite element discretizations on transient meshes, while the latter offers a remedy to retain high-order convergence rates also in cases where the physical boundary does not coincide with the boundary of the discretization. We investigate the performance of the method at two analytical and one experimental benchmark.

Coupling a NURBS contact interface with a higher order finite element discretization for contact problems using the mortar method

In this paper, a contact problem between two bodies, discretized by finite elements, is solved by adding an auxiliary NURBS layer between the bodies. The advantages of a smooth contact formulation in a NURBS approach are combined with simple mesh generation procedures for the bodies discretized with finite elements. Mesh tying conditions are used to couple the NURBS layer with the finite element discretization. The NURBS layer is the master side for contact and mesh tying. Mesh tying is enforced either using pointwise or mortar type approaches. Frictionless 2D and 3D contact problems are considered using small deformations. The contact problem is discretized with the mortar method and a penalty approach is used to enforce the contact constraints. A robust element-based quadrature is applied for mortar tying and contact discretizations, thus avoiding computationally expensive segmentation.

Higher order FEM for the obstacle problem of the [formula omitted]-Laplacian—A variational inequality approach

We consider higher order finite element discretizations of a nonlinear variational inequality formulation arising from an obstacle problem with the p-Laplacian differential operator for p∈(1,∞). We prove an a priori error estimate and convergence rates with respect to the mesh size h and in the polynomial degree q under assumed regularity. Moreover, we derive a general a posteriori error estimate which is valid for any uniformly bounded sequence of finite element functions. All our results contain the known results for the linear case of p=2. We present numerical results on the improved convergence rates of adaptive schemes (mesh size adaptivity with and without polynomial degree adaptation) for the singular case of p=1.5 and for the degenerated case of p=3.

Configurational forces in electronic structure calculations using Kohn-Sham density functional theory

We derive the expressions for configurational forces in Kohn-Sham density functional theory, which correspond to the generalized variational force computed as the derivative of the Kohn-Sham energy functional with respect to the position of a material point $\textbf{x}$. These configurational forces that result from the inner variations of the Kohn-Sham energy functional provide a unified framework to compute atomic forces as well as stress tensor for geometry optimization. Importantly, owing to the variational nature of the formulation, these configurational forces inherently account for the Pulay corrections. The formulation presented in this work treats both pseudopotential and all-electron calculations in single framework, and employs a local variational real-space formulation of Kohn-Sham DFT expressed in terms of the non-orthogonal wavefunctions that is amenable to reduced-order scaling techniques. We demonstrate the accuracy and performance of the proposed configurational force approach on benchmark all-electron and pseudopotential calculations conducted using higher-order finite-element discretization. To this end, we examine the rates of convergence of the finite-element discretization in the computed forces and stresses for various materials systems, and, further, verify the accuracy from finite-differencing the energy. Wherever applicable, we also compare the forces and stresses with those obtained from Kohn-Sham DFT calculations employing plane-wave basis (pseudopotential calculations) and Gaussian basis (all-electron calculations). Finally, we verify the accuracy of the forces on large materials systems involving a metallic aluminum nanocluster containing 666 atoms and an alkane chain containing 902 atoms, where the Kohn-Sham electronic ground state is computed using a reduced-order scaling subspace projection technique (P. Motamarri and V. Gavini, Phys. Rev. B 90, 115127).

Electromagnetic Analysis Using High-Order Numerical Schemes in Space and Time

The current paper establishes different axisymmetric and two-dimensional models for an electrostatic, magnetostatic and electromagnetic induction process. Therein, the Maxwell equations are combined in a monolithic solution strategy. A higher order finite element discretization using Galerkin’s method in space as well as in time is developed for the electromagnetic approach. In addition, time integration procedures of the Runge–Kutta family are evolved. Furthermore, the residual error is introduced to open an alternative way for a numerically efficient estimation of the time integration accuracy of the Galerkin time integration method. Runge–Kutta methods are enriched by the embedded error estimate. A family of electrostatic, magnetostatic and electromagneto dynamic boundary and initial boundary value problems with existing analytical solutions are introduced, which will serve as benchmark examples for numerical solution procedures.

High‐Order Numerical Methods for Electromagnetic Induction

AbstractThe present paper establishes an axisymmetric benchmark model of a conducting loop, which implies an electromagnetic induction. Therein, the fully coupled MAXWELL equations are demonstrated in a monolithic solution strategy. This dynamic problem is solved with a high order finite element discretization using GALERKIN's method in space as well as in time. Furthermore, high order RUNGE‐KUTTA time integration methods are analyzed. Studies regarding an h error estimation and the order of convergence are examined. (© 2017 Wiley‐VCH Verlag GmbH & Co. KGaA, Weinheim)

Efficient resonance computations for Helmholtz problems based on a Dirichlet-to-Neumann map

We present an efficient procedure for computing resonances and resonant modes of Helmholtz problems posed in exterior domains. The problem is formulated as a nonlinear eigenvalue problem (NEP), where the nonlinearity arises from the use of a Dirichlet-to-Neumann map, which accounts for modeling unbounded domains. We consider a variational formulation and show that the spectrum consists of isolated eigenvalues of finite multiplicity that only can accumulate at infinity. The proposed method is based on a high order finite element discretization combined with a specialization of the Tensor Infinite Arnoldi method (TIAR). Using Toeplitz matrices, we show how to specialize this method to our specific structure. In particular we introduce a pole cancellation technique in order to increase the radius of convergence for computation of eigenvalues that lie close to the poles of the matrix-valued function. The solution scheme can be applied to multiple resonators with a varying refractive index that is not necessarily piecewise constant. We present two test cases to show stability, performance and numerical accuracy of the method. In particular the use of a high order finite element discretization together with TIAR results in an efficient and reliable method to compute resonances.

Numerical Methods for Power-Law Diffusion Problems

In this paper, we consider numerical methods for nonlinear diffusion problems where the diffusion term follows a power law, e.g., $p$-Laplace-type problems. In the first part, we present continuous higher order finite element discretizations for the model problem and we derive error estimates. In the second part, we discuss Newton iterative methods based on residual-based line-search and error-oriented globalization, which are employed for the numerical solution of the produced nonlinear algebraic system. Third, we formulate the original problem as a saddle point problem in the frame of augmented Lagrangian techniques and present two iterative methods for its solution. We conduct a systematic investigation of all solution algorithms. These algorithms are compared with respect to computational cost and their efficiency. Numerical results demonstrating the theoretical error estimates are also presented in five examples.

Weighted BFBT Preconditioner for Stokes Flow Problems with Highly Heterogeneous Viscosity

We present a weighted BFBT approximation (w-BFBT) to the inverse Schur complement of a Stokes system with highly heterogeneous viscosity. When used as part of a Schur complement-based Stokes preconditioner, we observe robust fast convergence for Stokes problems with smooth but highly varying (up to 10 orders of magnitude) viscosities, optimal algorithmic scalability with respect to mesh refinement, and only a mild dependence on the polynomial order of high-order finite element discretizations ($Q_k \times P_{k-1}^{disc}$, order $k \ge 2$). For certain difficult problems, we demonstrate numerically that w-BFBT significantly improves Stokes solver convergence over the widely used inverse viscosity-weighted pressure mass matrix approximation of the Schur complement. In addition, we derive theoretical eigenvalue bounds to prove spectral equivalence of w-BFBT. Using detailed numerical experiments, we discuss modifications to w-BFBT at Dirichlet boundaries that decrease the number of iterations. The overall algorithmic performance of the Stokes solver is governed by the efficacy of w-BFBT as a Schur complement approximation and, in addition, by our parallel hybrid spectral-geometric-algebraic multigrid (HMG) method, which we use to approximate the inverses of the viscous block and variable-coefficient pressure Poisson operators within w-BFBT. Building on the scalability of HMG, our Stokes solver achieves a parallel efficiency of 90% while weak scaling over a more than 600-fold increase from 48 to all 30,000 cores of TACC's Lonestar 5 supercomputer.

Schwarz type preconditioners for the neutron diffusion equation

Domain decomposition is a mature methodology that has been used to accelerate the convergence of partial differential equations. Even if it was devised as a solver by itself, it is usually employed together with Krylov iterative methods improving its rate of convergence, and providing scalability with respect to the size of the problem.In this work, a high order finite element discretization of the neutron diffusion equation is considered. In this problem the preconditioning of large and sparse linear systems arising from a source driven formulation becomes necessary due to the complexity of the problem. On the other hand, preconditioners based on an incomplete factorization are very expensive from the point of view of memory requirements. The acceleration of the neutron diffusion equation is thus studied here by using alternative preconditioners based on domain decomposition techniques inside Schur complement methodology. The study considers substructuring preconditioners, which do not involve overlapping, and additive Schwarz preconditioners, where some overlapping between the subdomains is taken into account.The performance of the different approaches is studied numerically using two-dimensional and three-dimensional problems. It is shown that some of the proposed methodologies outperform incomplete LU factorization for preconditioning as long as the linear system to be solved is large enough, as it occurs for three-dimensional problems. They also outperform classical diagonal Jacobi preconditioners, as long as the number of systems to be solved is large enough in such a way that the overhead of building the preconditioner is less than the improvement in the convergence rate.