Hierarchical proximal Galerkin: a fast $hp$-FEM solver for variational problems with pointwise inequality constraints

Abstract

We leverage the proximal Galerkin algorithm (Keith and Surowiec, Foundations of Computational Mathematics, 2024), a recently introduced mesh-independent algorithm, to obtain a high-order finite element solver for variational problems, posed on tensor-product domains, with pointwise inequality constraints. This is achieved by discretizing the saddle point systems, arising from the latent variable proximal point method, with the hierarchical $p$-finite element basis. This results in discretized sparse Newton systems that admit a simple and effective block preconditioner. The solver can handle both obstacle-type, $u\le \phi$, and gradient-type, $|\nabla u|\le \phi$, constraints. We apply the resulting algorithm to solve obstacle problems with $hp$-adaptivity, a three-dimensional obstacle problem, a gradient-type constrained problem, and the thermoforming problem, an example of an obstacle-type quasi-variational inequality. We observe $hp$-robustness in the number of Newton iterations and only mild growth in the number of inner Krylov iterations to solve the Newton systems. Crucially we also provide wall-clock timings that are faster than low-order discretization counterparts.