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

TL;DR

This work uses the proximal Galerkin algorithm, a recently introduced mesh-independent algorithm, to obtain a high-order finite element solver for variational problems with pointwise inequality constraints, and applies the resulting algorithm to solve obstacle problems withhp-adaptivity, a gradient-type constrained problem, and the thermoforming problem, an example of an obstacle-type quasi-variational inequality.

Abstract

We leverage the proximal Galerkin algorithm (Keith and Surowiec, Foundations of Computational Mathematics, 2024, DOI: 10.1007/s10208-024-09681-8), a recently introduced mesh-independent algorithm, to obtain a high-order finite element solver for variational problems 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 \leq \varphi$, and gradient-type, $|\nabla u| \leq \varphi$, constraints. We apply the resulting algorithm to solve obstacle problems with $hp$-adaptivity, 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.

Figures (8)

• Figure 1: The spy plots of the scaled stiffness matrix $A_\alpha$, its lower Cholesky factor $L^{\text{chol}}_\alpha$ or reverse lower Cholesky factor $L^{\text{rchol}}_\alpha$, the Gram matrix $B$, the exponential block $D_\psi$, and $D_\psi$ with permuted rows and columns to reveal the hidden block diagonal structure for the Newton linearization of the pG obstacle subproblem \ref{['eq:pG']} when using: (top row) a one-dimensional $p$-FEM discretization on a uniform mesh with five cells and degree $p=10$ on each cell and (bottom row) a two-dimensional tensor-product $p$-FEM discretization with 25 cells and partial degree $p=5$ on each cell.
• Figure 2: (Obstacle problem). We consider $d\in\{1,2\}$ and fix $\alpha =1$, $\beta=0$, $D_\psi \equiv 0$ and consider the dense Schur complement $S$ as defined \ref{['eq:schur2']}, in the context of the obstacle subproblem \ref{['eq:pG']}, and its sparse preconditioner $\hat{S}$ as defined in \ref{['eq:schur-preconditioner']}. We measure the growth in the number of $\hat{S}$-left-preconditioned GMRES iterations to solve $S \boldsymbol{\mathbf{x}} = \boldsymbol{1}$ to a relative error of $10^{-6}$ (i.e. $\|\boldsymbol{\mathbf{r}}_k\|_{\ell^2}/\|\boldsymbol{\mathbf{r}}_0\|_{\ell^2} \leq 10^{-6}$) with respect to partial degree $p$ in (a) and (e) and $1/h$ in (c) and (g). We observe no growth in 1D with respect to $p$ and polylogarithmic growth in the other cases. We also measure the sparse LU factorization time for $\hat{S}$ with respect to partial degree $p$ in (b) and (f) and $1/h$ in (d) and (h).
• Figure 3: (Gradient-type constraint). We consider $d=2$ and fix $\alpha =1$, $\beta=10^{-5}$, $D_\psi \equiv 0$ and consider the dense Schur complement $S$ as defined \ref{['eq:schur2']}, in the context of the pG gradient-type subproblem \ref{['eq:pG:gradient']}, and its sparse preconditioner $\hat{S}$ as defined in \ref{['eq:schur-preconditioner:gradient']}. We measure the growth in the number of $\hat{S}$-left-preconditioned GMRES iterations to solve $S \boldsymbol{\mathbf{x}} = \boldsymbol{1}$ to a relative error of $10^{-6}$ with respect to partial degree $p$ in (a) and $1/h$ in (c). We observe slower than logarithmic growth in both cases. We also measure the sparse LU factorization time for $\hat{S}$ with respect to partial degree $p$ in (b) and $1/h$ in (d).
• Figure 4: (Left) The solution to the 1D obstacle problem with the setup \ref{['eq:osc-data-setup']}. (Right) Convergence of 7 refinement strategies. $hp$-adaptive utilizes \ref{['alg:adaptivity']} ($\delta =0.7$, $\sigma = 0.8$), $h$-adaptive is \ref{['alg:adaptivity']} with $\sigma = 0$ ($\delta = 0.3$), $h$-uniform implies the mesh size is uniformly halved with each refinement and $p$-uniform is where the polynomial degree is incremented by one with each refinement. Any number associated with a data point is the average time taken per Newton iteration, measured in milliseconds, via a direct sparse factorization. In the order listed in the legend, the first two strategies converge at a rate of 1, the next 3 at a rate slightly faster than 3/2 and the final two at a rate of 5 and 10, respectively (in the regime of $h$ and $p$ considered).
• Figure 5: (Left) The solution to the obstacle problem with the setup \ref{['eq:oscillatory-obstacle-setup']}. (Right) A 1D slice through the solution at $y=1/2$.

Theorems & Definitions (8)

• Definition 3.1: $L^2$-conforming hierarchical $p$-FEM space $\Psi_{h,\boldsymbol{p}}$
• Lemma 3.2: Diagonal mass matrix
• Definition 3.3: $H^1$-conforming hierarchical $p$-FEM space $U_{h,\boldsymbol{p}}$
• Remark 3.4
• Remark 4.1: Other Krylov methods
• Remark 5.1: pG error indicator
• Remark 7.1: Comparing solver times
• Remark 7.2: Thermoforming QVI