Source-linked AI summary
An adaptive multigrid solver for DPG methods with applications in linear acoustics and electromagnetics
Socratis Petrides, Leszek Demkowicz
TL;DR
Wave-problem discretizations face stability and pollution challenges, and their preconditioning remains difficult. This paper develops an adaptive multigrid preconditioner for DPG systems and reports stable, reliable performance on adaptive hp-meshes, with uniform convergence under certain conditions.
Problem
High-frequency wave discretizations can lose stability or suffer pollution, while effective preconditioning of the resulting highly indefinite systems remains challenging.
Source: §1.1 · pp. 1–2Method
The preconditioner operates on DPG trace spaces on the mesh skeleton and is integrated with the DPG adaptive mesh-refinement process.
Source: §4.1 · pp. 11–12Results
The method is reported stable and reliable for adaptive hp-meshes, with uniform convergence in frequency, polynomial order, and discretization size when the coarse grid is sufficiently good.
Source: §6 · p. 27Takeaways & Limitations
Integrating the iterative solver with adaptive DPG refinement is presented as a promising approach, including in the pre-asymptotic regime.
Source: §3.3 · p. 11The theoretical analysis is for the one-level additive preconditioner in acoustics, and the interpolation-norm bound is not expected to be sharp.
Source: §3.3 · p. 11, §6 · p. 27
Abstract
from arXiv · showhide
We propose an adaptive multigrid preconditioning technology for solving linear systems arising from Discontinuous Petrov-Galerkin (DPG) discretizations. Unlike standard multigrid techniques, this preconditioner involves only trace spaces defined on the mesh skeleton, and it is suitable for adaptive hp-meshes. The key point of the construction is the integration of the iterative solver with a fully automatic and reliable mesh refinement process provided by the DPG technology. The efficacy of the solution technique is showcased with numerous examples of linear acoustics and electromagnetic simulations, including simulations in the high-frequency regime, problems which otherwise would be intractable. Finally, we analyze the one-level preconditioner (smoother) for uniform meshes and we demonstrate that theoretical estimates of the condition number of the preconditioned linear system can be derived based on well established theory for self-adjoint positive definite operators.
1. INTRODUCTION
Wave simulations face stability and pollution challenges at high frequencies, while scalable preconditioning for wave problems remains difficult. This paper develops a DPG preconditioning approach informed by prior DPG work and examines its behavior for wave problems.
1.1. Background.
- In two and three dimensions, wave discretizations can lose stability unless the mesh resolves the frequency, and higher-order discretizations reduce but do not eliminate pollution.The Nyquist criterion requires a minimum number of elements or points per wavelength; the authors note no known pollution-free method in these dimensions [4].Source: §1.1 · p. 1
- Wave-problem preconditioners remain challenging because elliptic multigrid methods fail on highly indefinite operators, and coarse-grid solves may need expensive fine meshes.Prior work on DPG wave preconditioners reported convergence behavior sensitive to Schwarz overlap or omitted mesh-size dependence studies.Source: §1.1 · p. 2, §1.2 · p. 2
1.3. Outline.
- Earlier DPG preconditioning studies addressed elliptic problems, while wave-problem studies lacked theoretical results in the authors' account.The paper proceeds from a DPG and acoustics overview to one-level analysis, multigrid construction, and acoustics and electromagnetics experiments.Source: §1.2 · p. 2, §1.3 · p. 2
1.2. Preconditioning the DPG system.
- For DPG wave systems, this work constructs a preconditioner from the wave bilinear form and numerically examines dependence on polynomial order, frequency, patch size, and mesh size.The analysis covers a one-level additive smoother and is intended to inform multilevel construction.Source: §1.2 · p. 2
2. THE DPG METHOD AND LINEAR ACOUSTICS
The paper introduces the ideal Petrov–Galerkin and ultraweak DPG formulations for linear acoustics, then explains how broken test spaces and practical finite-dimensional test spaces support computation and preconditioner analysis.
2.1. The ideal Petrov–Galerkin method: a brief outline.
- The ideal Petrov–Galerkin method uses a trial-to-test operator to construct optimal test functions and minimize the residual in the test-space dual norm.The resulting system is symmetric positive definite and guarantees uniform, unconditional discrete stability.Source: §2.1 · p. 3
- The norm of the error representation function provides a built-in a posteriori indicator of the Galerkin error in the energy norm.The method solves simultaneously for the original unknown and the error representation function.Source: §2.1 · p. 4
2.2. The ultraweak formulation – linear acoustics.
- The ultraweak acoustics formulation transfers regularity to the test space, while introducing trace unknowns to eliminate test-function boundary conditions.The group variable comprises pressure and velocity; the trace variables are pressure and normal-velocity traces.Source: §2.2 · pp. 4–5
2.3. The Discontinuous Petrov–Galerkin method: “breaking” the test space.
- Breaking the test space defines fluxes and traces on the mesh skeleton, and the energy norm's elementwise structure permits single-element norm-equivalence analysis.The global energy norm squared is the sum of the element contributions.Source: §2.3 · p. 5
2.4. The practical DPG method.
- Practical DPG approximates the Riesz-operator inversion with an enriched finite-dimensional test space because exact inversion requires solving an infinite-dimensional boundary-value problem.The test-space dimension is chosen much larger than the trial-space dimension; cited studies report stability is not significantly affected.Source: §2.4 · p. 6
3. ANALYSIS OF THE ONE-LEVEL PRECONDITIONER
The section analyzes a one-level additive Schwarz smoother for the ultraweak DPG acoustics formulation using subspace correction theory for self-adjoint positive definite operators. It verifies the theorem’s assumptions and studies how the interpolation norm and solver convergence vary with discretization and frequency.
3.1. Additive Schwarz preconditioner and the subspace correction theory.
- Subspace correction theory gives a condition-number estimate for the preconditioned operator when strengthened Cauchy–Schwarz and stable-decomposition assumptions hold.The preconditioned matrix is P = AB, with the theorem relating its condition number to the constants β and α.Source: §3.1 · pp. 6–7
3.2.1. Setup.
- The DPG analysis uses the energy norm of the field-and-trace variable and partitions the domain into vertex patches with finite overlap.Local energy subspaces restrict trace variables to the mesh skeleton within each patch.Source: §3.2.1 · p. 7
3.2.2. Strengthened Cauchy–Schwarz inequality.
- Broken test spaces make the trial-to-test operator local, so disjoint patches have zero coupling and finite overlap bounds the strengthened Cauchy–Schwarz constant.The bound is expressed through the spectral radius of the patch-coupling matrix.Source: §3.2.2 · pp. 7–8
3.2.3. Stable Decomposition.
- The stable-decomposition argument constructs patch contributions for both the L2 field space and the trace space, with the proof presented in two dimensions and the three-dimensional case described as analogous.The construction relies on a partition of unity subordinate to the vertex-patch cover.Source: §3.2.3 · p. 8
- The stable-decomposition proof uses an operator-dependent skeleton interpolant; its patch constant is frequency-independent, while the interpolation constant is expected to grow mildly with acoustics frequency.The resulting bound also uses the global boundedness-below constant, which is frequency-independent for the acoustics operator with impedance boundary conditions.Source: §3.2.3 · pp. 8–10
3.3. Numerical study of ∥ˆΠ∥ˆ
- CG convergence is reported independent of polynomial order and mesh size, while iteration counts grow mildly with frequency; unresolved high-frequency cases use fewer iterations.The experiments use plane waves with impedance boundary conditions and uniform h-refinement; the authors identify the pre-asymptotic behavior as potentially useful for adaptive refinement.Source: Table 1 · p. 10, Table 2 · p. 10, Table 3 · p. 10
- The computed interpolation norm tracks solver convergence, but the authors do not claim its bound is sharp.They present it as useful insight into preconditioner effectiveness and investigate its dependence on frequency, mesh size, and polynomial order.Source: §3.2.3 · p. 10, §3.3 · p. 11
- Static condensation eliminates element-interior field variables, leaving a smaller trace-only system whose condition-number bound matches that of the original operator.The stated equivalence follows from a Schur complement result.Source: §3.3 · p. 11
4. EXTENSION TO MULTIGRID
The extension couples an additive Schwarz smoother multiplicatively with a coarse-grid solve, extending the approach to three-dimensional and multilevel simulations. Its implementation combines inter-grid transfer operators, smoothing, and a coarse-grid solver.
4. EXTENSION TO MULTIGRID
- The multilevel extension couples an additive Schwarz smoother multiplicatively with a coarse-grid solve and extends earlier two-dimensional acoustics work to three-dimensional and multilevel simulations.The one-level preconditioner's dependence on subdomain count and overlap creates a trade-off between iteration counts and local-solve work as h decreases.Source: §4 · p. 11
4.1. Discussion on implementation.
- The h-refinement transfer composes natural inclusion from the coarse grid to a macro-grid with Schur-complement extension to the fine grid.The macro-grid is formed by eliminating trace variables off the coarse-grid skeleton; this handles new fine-grid interfaces that have no coarse-grid ancestors.Source: §4.1 · pp. 11–12
- For p-refinement with hierarchical shape functions, new degrees of freedom are set to zero, and restriction is the transpose of prolongation.Source: §4.1 · p. 12
- The additive Schwarz smoother uses patches supported on coarse-grid vertex basis functions one level down, with smoothing solves performed on macro-grids.Figure 4.2 illustrates the patch construction.Source: §4.1 · p. 12, Figure 4.2 · p. 12
- The coarse-grid correction uses the MUMPS multifrontal sparse direct solver on a grid whose field variables are eliminated and whose size is usually significantly smaller than the fine-grid system.Source: §4.1 · p. 13
5. NUMERICAL RESULTS
The numerical-results chapter reports large acoustics and electromagnetic simulations using CG preconditioned by multigrid, including uniform and adaptive refinement settings.
5. NUMERICAL RESULTS
- The chapter presents numerical experiments for acoustics and Maxwell equations with both uniform and adaptive refinement.Source: §5 · p. 13
- The simulations use the Conjugate Gradient solver preconditioned with multigrid technology.Source: §5 · p. 13
- All simulations ran on a single node with 24 cores and 256G of RAM using the hp3D finite element library.Source: §5 · p. 13
5.1. Time harmonic Maxwell equations.
The section presents time-harmonic Maxwell simulations using DPG multigrid, comparing it with standard multigrid before testing adaptive refinements on singular and high-frequency scattering problems.
5.1.1. Comparison with standard multigrid methods.
- At ω = 10, DPG iteration counts grow mildly and converge to the true solution, while GMRES can be unreliable on pre-asymptotic coarse grids.Uniform convergence is recovered for both methods when the coarse grid is fine enough.Source: Table 5 · p. 14, §5.1.1 · p. 14
5.1.2. Fichera “oven” problem.
- In the Fichera “oven” problem, the residual tends to zero and preconditioned CG converges uniformly with discretization size throughout adaptive refinement.The solution has no known exact form, and the problem was too large for a direct solve.Source: §5.1.2 · p. 18
5.1.3. Scattering of a Gaussian beam from a cube.
- The strategy h-refines marked elements unless their size is less than half a wavelength, when it p-refines them; elements adjacent to scatterer corners and edges are forced to h-refine as needed.The simulation uses ω = 50π, corresponding to approximately 40 wavelengths inside the computational domain.Source: §5.1.3 · p. 18
- In the Gaussian-beam case, the preconditioned CG iteration count remains under control throughout adaptive refinement.The DPG residual is driven to zero during the process.Source: §5.1.3 · p. 19
5.2. Linear acoustics equations.
The section demonstrates DPG multigrid on three-dimensional linear-acoustics scattering problems, including a plane wave interacting with a sphere and a cube.
5.2.1. Scattering of a plane wave from a sphere.
- The sphere case uses a plane wave with ω = 35π, approximately 30 wavelengths in the domain, and hard boundary conditions on the sphere.A homogeneous impedance condition truncates the outer cube boundary.Source: §5.2.1 · pp. 21–22
- With a zero initial guess and residual tolerance below 10⁻⁵, preconditioned CG converged in 12 iterations for a uniform-refinement system of approximately 14 million degrees of freedom; sparse multifrontal solution was infeasible at ≈350 GB.The fine grid consists of 40,960 quartic hexahedra, and the coarse grid consists of 640 quartic hexahedra.Source: §5.2.1 · p. 22
5.2.2. Scattering of a plane wave from a cube.
- For plane-wave scattering from a cube, the DPG residual converges and preconditioned CG shows robust convergence throughout adaptive refinement.Refinements occur near scatterer singularities and in the PML region.Source: §5.2.2 · p. 24
5.3. Scattering of a Gaussian beam from a cube.
A high-frequency Gaussian beam scatters from a cube in a PML-truncated domain. Adaptive hp-refinement resolves the scatterer and PML features, while multigrid-preconditioned CG iterations remain controlled.
5.3. Scattering of a Gaussian beam from a cube.
- At ω = 140π, approximately 85 wavelengths occupy the domain; elements are h-refined until smaller than half a wavelength, then p-refined, except near singularities and in the PML.The exception regions are restricted to h-refinement, and adaptation stops when the DPG residual drops by an order of magnitude.Source: §5.3 · p. 24
- The solver uses prolongation of the previous mesh solution as its initial iterate, ten smoothing steps per multigrid level, and skips smoothing where the local patch residual is near zero.The run uses at most eight levels and terminates CG when the residual norm falls below 10^-3.Source: §5.3 · p. 25
- The DPG residual tends to zero only after the decay in the PML region is resolved, while preconditioned CG iterations remain under control throughout adaptation.The meshes and solutions show that scatterer singularities must be resolved before the wave can propagate.Source: §5.3 · p. 25
6. CONCLUSION
The work presents an efficient, robust preconditioner for DPG linear systems, suitable for adaptive hp-meshes and reliable in acoustics and electromagnetic simulations. Integrating it with adaptive refinement supports efficient high-frequency problems that would be intractable on uniform meshes, while the paper also develops one-level smoother theory.
6. CONCLUSION
- The preconditioner is stable and reliable for adaptive hp-meshes, and under a sufficiently good coarse grid can converge uniformly with frequency, polynomial order, and discretization size.The reported simulations cover acoustics and electromagnetics.Source: §6 · p. 27
- Integrating the iterative solver with adaptive refinement allowed efficient solution of high-frequency problems that would be intractable with uniform meshes.Source: §6 · p. 27
- The analysis presents theoretical results for the one-level additive preconditioner and examines its frequency dependence for the two-dimensional acoustics problem; multilevel analysis remains future work.The paper also identifies distributed-memory implementation as ongoing work.Source: §6 · p. 27
APPENDIX A. COMPUTING THE INTERPOLATION NORM
The appendix describes a projection-based interpolation operator on the mesh skeleton and computes its continuity norm through a generalized eigenvalue problem. It builds the interpolant from vertex and edge contributions, with practical calculations using polynomial approximations to minimum-energy extensions.
APPENDIX A. COMPUTING THE INTERPOLATION NORM
- The continuity constant is the square root of the maximum generalized eigenvalue, computed from the interpolation matrix and the Gram matrix for the trace-space inner product.Source: Appendix A · pp. 27–28
- The Gram matrix uses inner products of minimum-energy extensions, which are approximated in practice by polynomial extensions of sufficiently large order.The numerical results use polynomial order 2 for these extensions.Source: Appendix A · p. 28
- The interpolant matches the input at element vertices, then projects the vertex-vanishing remainder onto edge bubble spaces; the final interpolant sums the vertex and edge parts.Source: Appendix A · p. 28
- Increasing the polynomial approximation order further produced no notable change in the numerical results.Source: Appendix A · p. 28
- The matrix representation is assembled by applying the operator to polynomial basis functions, using vertex values and minimization for edge bubbles.Both sides of the edge-bubble minimization can be computed with the Gram matrix.Source: Appendix A · pp. 28–29