PaperlayerView exact revision on arXiv

Source-linked AI summary

Scalable DPG Multigrid Solver for Helmholtz Problems: A Study on Convergence

Jacob Badger, Stefan Henneking, Socratis Petrides, Leszek Demkowicz

arXiv:2304.01728v3math.NAcs.CEcs.NA

TL;DR

High-frequency wave operators pose stability, cost, and solver challenges, motivating a scalable multigrid preconditioner for DPG discretizations. The solver studies show robust convergence under uniform refinement and linear dependence on frequency, with similar frequency dependence for hp-adaptive refinement.

  • Problem

    High-frequency wave problems need very fine meshes for stability, become prohibitively expensive, and yield highly indefinite systems that can break down standard iterative solvers.

    Source: §1.1 · p. 2

  • Method

    DPG-MG is a multigrid-preconditioned conjugate gradient solver on a refinement hierarchy, using the DPG error indicator to select subsequent meshes.

    Source: §3 · p. 5

  • Results

    The solver demonstrates h- and p-robust convergence with linear dependence on wave frequency under uniform refinement, and a similar linear frequency dependence under hp-adaptive refinement.

    Source: §6 · p. 18

  • Takeaways & Limitations

    The scalable MPI/OpenMP implementation is competitive for high-frequency wave propagation problems despite iteration counts increasing linearly with frequency.

    Source: §6 · p. 18

    The cause of deteriorated convergence on hp-adaptive meshes remains under investigation, with multiple element sizes and Gram-matrix scales present simultaneously.

    Source: §4.3 · p. 12

Abstract

from arXiv · show

This paper presents a scalable multigrid preconditioner targeting large-scale systems arising from discontinuous Petrov-Galerkin (DPG) discretizations of high-frequency wave operators. This work is built on previously developed multigrid preconditioning techniques of Petrides and Demkowicz (Comput. Math. Appl. 87 (2021) pp. 12-26) and extends the convergence results from $\mathcal{O}(10^7)$ degrees of freedom (DOFs) to $\mathcal{O}(10^9)$ DOFs using a new scalable parallel MPI/OpenMP implementation. Novel contributions of this paper include an alternative definition of coarse-grid systems based on restriction of fine-grid operators, yielding superior convergence results. In the uniform refinement setting, a detailed convergence study is provided, demonstrating h and p robust convergence and linear dependence with respect to the wave frequency. The paper concludes with numerical results on hp-adaptive simulations including a large-scale seismic modeling benchmark problem with high material contrast.

1 Introduction

High-frequency wave problems challenge standard numerical methods through stability, cost, and indefiniteness. The paper studies a scalable DPG multigrid solver, building on prior work and examining its convergence under several refinement strategies.

1.1 Background

  • High-frequency wave problems require very fine meshes for stability, become prohibitively expensive, and can defeat standard iterative solvers because their systems are highly indefinite.Existing preconditioners also lose efficiency in heterogeneous media and can be difficult to apply in complex geometries.Source: §1.1 · p. 2

1.2 DPG-MG Solver

  • The ultraweak DPG formulation is presented as less dissipative than other DPG formulations, with dispersion error roughly commensurate to Galerkin discretizations.Prior work used it to define an hp-adaptive multilevel preconditioner for wave problems [50].Source: §1.2 · p. 2
  • The DPG-MG solver uses a condensed global system of trace degrees of freedom, for which constructing stable prolongation between non-nested meshes is a central challenge.Petrides and Demkowicz addressed this challenge for general DPG problems.Source: §1.2 · pp. 2–3

1.3 Direction and Outline

  • The authors developed a scalable hybrid MPI/OpenMP implementation and use it to study convergence under uniform h, uniform p, and hp-adaptive refinement.The paper aims to elucidate solver scaling characteristics and identify aspects of its construction that may be improved.Source: §1.3 · p. 3

2 Ultraweak DPG for Helmholtz

This section formulates time-harmonic acoustics with an impedance boundary condition using a broken ultraweak DPG formulation. It describes the optimal-test-space construction and the trace variables introduced by the broken test space.

DPG method

  • The DPG method constructs stable discretizations using optimal test functions, approximated in practice by inverting the global Riesz map over an enriched discontinuous test space.The test-space polynomial order is uniformly increased by 1 for enrichment, and its element-local computation follows from the broken test space.Source: §2 · p. 3

Broken ultraweak formulation

  • Breaking the test space introduces trace unknowns on the mesh skeleton: normal velocity (flux) and pressure on element boundaries.The normal velocity is discretized as the normal trace of H(div)-conforming elements, and pressure as the continuous trace of H1-conforming elements.Source: §2 · pp. 4–5

Helmholtz problem

  • The model is the first-order mixed form of time-harmonic linear acoustics with an inhomogeneous impedance boundary condition.Here p denotes pressure, u velocity, and ω angular wave frequency; the impedance condition uses Z and the outward normal flux u_n.Source: §2 · p. 4

Broken ultraweak formulation

  • α is a scaling constant, and numerical results throughout the work use α = 1; test functions satisfy Z−1q + vn = 0 on Γ to build in the impedance boundary condition.The impedance boundary condition implicitly implies additional regularity of boundary velocities; the load is well-defined for u0 ∈ H−1/2(Γ) with duality pairing or u0 ∈ L2(Γ) with the L2 pairing.Source: §2 · p. 5

3 DPG-MG Solver

DPG-MG applies multigrid-preconditioned conjugate gradients on a refinement hierarchy, with the solver operating on trace degrees of freedom. Its construction combines two-stage prolongation, macro-grid smoothing, and coarse-grid operators that may be assembled directly or restricted from fine-grid operators.

V-cycle and smoother

  • DPG-MG uses a refinement-generated mesh hierarchy and can use intermediate solutions only to the accuracy needed to produce the next hp-adaptive mesh.The DPG error indicator selects refinements, while mesh-independent stability permits initialization on arbitrarily coarse meshes.Source: §3 · pp. 5–6

Prolongation

  • Prolongation first condenses fine-grid degrees of freedom onto the coarse-grid skeleton, then applies natural inclusion on coarse-grid edges and, in 3D, faces.Natural inclusion is constructed through recursive constrained approximation.Source: §3 · pp. 5–6, Figure 1 · p. 6

V-cycle and smoother

  • The solver applies conjugate-gradient iterations and additive Schwarz smoothing on macro-grid degrees of freedom supported on coarse-grid vertex patches.Vertex patches can create large blocks, especially for high-order discretizations, but avoid added complexity for H(curl) and H(div) variables; the solver applies without modification to well-posed DPG problems with exact-sequence energy spaces.Source: §3 · pp. 6–7

Alternative construction of coarse-grid operators

  • Restriction-based coarse-grid stiffness matrices restore expected convergence at high frequencies, unlike direct assembly, which showed rapidly worsening convergence as frequency increased.The paper reports that restriction typically costs 1%–10% of fine-grid assembly; it defers a rigorous explanation of the difference between the constructions.Source: §3 · p. 7

Scalable MPI/OpenMP implementation in hp3D

  • The solver’s hybrid MPI/OpenMP implementation extends hp3D for distributed computation, using dynamic load balancing across grid levels to support hp-adaptive meshes.The original OpenMP-only implementation was limited by typical node memory; the present paper uses the scalable implementation to study convergence, leaving its scaling characteristics for future publication.Source: §3 · pp. 7–8

4 Convergence Studies

The convergence studies compare coarse-grid systems formed by restriction of fine-grid operators with systems stored from previous meshes, across uniform h-, uniform p-, and hp-adaptive refinements. Restriction gives better uniform-refinement convergence, while hp-adaptive meshes retain linear frequency dependence and require more iterations than uniform h-refinements.

4 Convergence Studies

  • Studies fix a residual-reduction tolerance of 10^7, reset each refined-grid initial solution to zero, and use one pre- and post-smoothing step except for the five-step comparison.The initial single-element mesh has order p = 2; it is solved by MUMPS, and its iterations are not reported.Source: §4 · p. 8

4.1 Problem Setup

  • The model problem propagates a Gaussian beam of waist-radius 0.1 through a homogeneous unit cube with impedance boundary conditions and a manufactured impedance load near the origin.Experiments use Frontera Cascade Lake nodes; timing statistics are deferred to Section 5.Source: §4.1 · p. 8, §4 · p. 8

4.2 Direct Assembly vs. Fine-Grid Restriction for Coarse-Grid Operators

  • Under uniform h-refinement, restriction yields h-robust convergence and fewer iterations than storing coarse-grid systems, which does not demonstrate h-robustness.Iterations increase roughly linearly with frequency in the asymptotic regime.Source: §4.2 · p. 9, Figure 2 · p. 9
  • Five smoothing steps tend to restore h-robustness for stored coarse-grid systems, but the number needed tends to increase with frequency.With restriction, five smoothing steps reduce smoothing iterations by only about a factor of two versus one step; one step is reported qualitatively optimal in the experiments.Source: §4.2 · pp. 9–10, Figure 3 · p. 10
  • Under uniform p-refinement after h-refinement to at least two elements per wavelength, both approaches are p-robust, but restriction requires far fewer iterations.The two-elements-per-wavelength threshold satisfies the Nyquist criterion.Source: §4.2 · p. 10, Figure 4 · p. 11

4.3 hp-Adaptive Refinements

  • Adaptive refinement uses Dörfler marking, h-refining marked elements until their maximum edge-length is below half a wavelength and otherwise p-refining them.The resulting meshes sweep from the beam’s corner into the domain; intermediate solutions are stable and partially resolve the wave.Source: §4.3 · p. 10, Figure 5 · p. 11
  • On hp-adaptive meshes, iterations still increase linearly with frequency, and the maximum is consistently higher than for uniform refinements.The final meshes show considerably fewer iterations and consist exclusively of p-refinements with a consistent, uniform element size.Source: §4.3 · pp. 11–12, Figure 6 · p. 12
  • The proposed explanation for degraded hp-adaptive convergence remains an initial hypothesis involving simultaneous element sizes and their Gram-matrix scales.Final p-refinement meshes have a consistent element size, and their iteration counts decrease considerably.Source: §4.3 · p. 12

5 Application – Seismic Modeling

The DPG-MG solver is applied to a high-contrast seismic subduction-zone benchmark using hp-adaptive meshes. The study reports billion-DOF runs, near-linear per-iteration scaling at fixed frequency, and a comparison with uniform refinement.

hp-adaptive refinements

  • The observed linear increase in iterations with frequency implies O(N^4/3) preasymptotic complexity, yet the solver remains competitive for large-scale high-frequency wave propagation.The reported frequency dependence is worse than the logarithmic increase observed for sweeping-type preconditioners.Source: §5 · p. 12

Benchmark setup

  • The GO_3D_OBS subduction-zone section spans 20 × 102 × 28.3 km, with wavespeeds from 1,500 to 8,500 m/s; runs used 28,672 cores.The model's material data were sampled on a 100 m grid.Source: §5 · p. 13, Figure 7 · p. 13

hp-adaptive refinements

  • Starting from an arbitrarily coarse mesh, seven initial h-refinements resolve the source before hp-refinement continues until the DPG residual falls by a factor of 3 · 10^2 or memory is exhausted.The strategy chooses h-refinement based on wavelength and high-contrast regions, otherwise p-refining up to order 5 before returning to h-refinement.Source: §5 · p. 13
  • The adaptive run ended with over 6.3 million elements and 1.9 billion DOFs, but the upper-right region remained incompletely refined when memory ran out.The implementation stored all 34 adaptive meshes and their coarse-grid operators; the authors suggest reducing the number of grid levels.Source: §5 · p. 14
  • The finest adaptive mesh solved in 210 seconds; the full job across 34 meshes took 3,029 seconds.The full-job time includes assembly and solution.Source: §5 · p. 15
  • At fixed frequency, time per iteration scales roughly linearly with DOFs, while super-linearity appears when DOFs grow sub-geometrically between levels.The authors suggest that reducing the number of grid levels could mitigate super-linear scaling and iteration cost.Source: §5 · p. 16, Figure 10 · p. 16
  • Scalable adaptivity is identified as a key DPG-MG differentiator and, to the authors’ knowledge, novel among Helmholtz solvers.The authors also note that adaptive meshes can attain similar accuracy with a small fraction of the DOFs, though the current implementation adds expense by storing previous grids.Source: §5 · p. 17

Uniform refinements

  • The adaptive mesh reached 1.9 billion DOFs; uniform refinement reached 5.6 billion, using nearly three times as many DOFs for a DPG residual ten times larger.The uniform-refinement comparison reports a time per iteration three times smaller, or nine times smaller when normalized by DOFs.Source: §5 · pp. 14, 18, Figure 11 · p. 17

6 Conclusion

The conclusion reports robust h- and p-convergence and linear dependence on wave frequency with restricted fine-grid coarse operators. The scalable implementation solves a high-contrast seismic benchmark at billion-DOF scale, while accuracy comparisons remain future work.

Conclusion

  • The solver handled the GO_3D_OBS benchmark at 1.9 billion DOFs on hp-adaptive meshes and 5.6 billion DOFs on uniform meshes.The authors describe these as larger-scale high-frequency heterogeneous wave-propagation runs than any they are aware of.Source: §6 · p. 18
  • With coarse-grid operators formed by restricting fine-grid operators, DPG-MG shows h- and p-robust convergence and linear dependence on wave frequency.A similar linear frequency dependence was observed for hp-adaptive refinements.Source: §6 · p. 18
  • The scalable MPI/OpenMP implementation is reported as competitive for high-frequency wave propagation despite iterations increasing linearly with frequency.The solver uses DPG properties including mesh-independent stability, a built-in error indicator, and Hermitian positive-definite discrete systems.Source: §6 · p. 18

Future directions

  • The authors identify refinement-tree refactoring, GPU integration, and fully distributed data structures as near-term ways to improve scaling and performance.They also plan to investigate deteriorated convergence with stored coarse-grid operators and hp-adaptive meshes.Source: §6 · pp. 18–19
Loading 2304.01728v3…