Source-linked AI summary

Model order reduction of HDG discretizations for harmonic wave propagation

Hélène Barucq, Florian Faucher, Ha Pham, Ishak Tifouti

arXiv:2610.20145v1math.NAcs.NAmath.AP

TL;DR

Repeated frequency sweeps for large-scale time-harmonic wave problems require many costly forward solves. The paper develops and analyzes projection-based model reduction for HDG discretizations, combining adaptive frequency sampling, derivative enrichment, and POD. The reduced model is evaluated on two- and three-dimensional heterogeneous media, with reported high online speedups and error bounds linked to the HDG structure.

  • Problem

    Frequency sweeps require repeated forward solves, and these multi-query evaluations can become prohibitively costly for large-scale simulations, especially in three dimensions.

    Source: §1 · p. 1
  • Method

    The method builds separate volume and trace reduced spaces using residual-driven greedy frequency sampling, derivative enrichment, and POD compression, then projects the non-condensed HDG formulation.

    Source: §1.2 · p. 3, §7 · p. 28
  • Results

    For Marmousi, 65 modes per variable achieve relative errors below 5 × 10^-4 and an online speedup of approximately 2 × 10^3; for a roughly 12-million-volume-unknown 3D problem, 39 modes achieve errors below 10^-4 and speedup of approximately 4 × 10^5.

    Source: §7 · p. 28
  • Takeaways & Limitations

    Numerical results show that lifting decouples trace error from volume-basis approximation error, while a sufficiently rich volume basis makes lifting unnecessary for the problems considered.

    Source: §7 · p. 28

    For the problems considered, lifting is not necessary when the volume basis is sufficiently rich; limited snapshots or distinct volume and trace structures can make the gap terms significant.

    Source: §5.2 · p. 20

Abstract

from arXiv · show

We propose and analyze a projection-based model reduction method for the rapid and reliable solution of time-harmonic acoustic wave propagation problems over a frequency band, with particular emphasis on seismic applications. The high-fidelity model is based on a hybridizable discontinuous Galerkin (HDG) discretization, formulated as a generalized saddle-point system coupling the volume and trace unknowns. The basis construction combines residual-based greedy sampling of the frequency parameter, enrichment with frequency derivatives of the HDG solution, and proper orthogonal decomposition (POD) compression. We establish the well posedness of both the high-fidelity and reduced problems and characterize their stability with inf-sup constants. We investigate the interaction between static condensation, inherent to the HDG discretization, and projection onto the reduced spaces, characterizing its implications for the reduced formulation and approximation errors. Since these operations do not commute in general, we derive a priori error estimates for both the trace and volume variables to quantify the gap between the two operations. We further identify a lifting condition on the volume reduced space that eliminates these gaps and yields a trace error estimate independent of the volume approximation error. Numerical experiments on the two-dimensional Marmousi model and a three-dimensional heterogeneous medium, with more than $10^{7}$ unknowns, illustrate the theoretical estimates and show that the reduced model achieves speedups of several orders of magnitude compared to the high-fidelity solver.

1 Introduction

Repeated frequency-domain solves make large-scale wave simulations costly, motivating model reduction for multi-query evaluation. This work adapts projection-based reduction to HDG’s separate volume and trace variables, combining adaptive sampling and derivative enrichment with analysis of stability and approximation errors.

1 Introduction

  • HDG reduces global coupling by eliminating element-wise volume unknowns, but repeated solves across frequencies remain costly, motivating MOR for large-scale multi-query problems.Prior approaches include POD and greedy sampling of solution fields, as well as derivative-based moment matching for algebraic systems.Source: §1 · p. 1, §1.1 · p. 2

1.2 Contributions

  • The proposed basis combines weak-greedy frequency sampling, frequency-derivative enrichment, and POD compression, with separate spaces for volume and hybrid trace variables.Derivative enrichment adds local information to global adaptive sampling without increasing the number of full-order training solves.Source: §1.2 · p. 3
  • The reduced formulation projects the non-condensed HDG system, retaining flux variables and preserving affine parameter dependence.The formulation uses distinct reduced spaces for volume and trace variables.Source: §1.2 · p. 3
  • The analysis establishes well-posedness for the full and reduced problems and derives a priori estimates quantifying how hybrid approximation error propagates to volume variables.Separate reduced spaces provide flexibility in trial and test functions and allow HDG structure to inform energy identities.Source: §1.2 · p. 3
  • Lifting trace-coupling contributions into the volume space decouples trace error from volume-space approximation error, allowing the trace approximation to inherit the full-order HDG approximation properties.The lifting strategy is an enrichment constructed from the error estimates.Source: §1.2 · p. 3
  • The paper combines established reduction ingredients in an HDG framework and evaluates them on two- and three-dimensional heterogeneous acoustic wave problems.The parametrization is restricted to frequency for a fixed medium; medium variability is left for future work.Source: §1.2 · p. 3

2 Acoustic wave propagation

The paper models damped time-harmonic acoustic waves in two- and three-dimensional media with separate reflecting and absorbing boundary portions. Its first-order formulation uses complex frequency to represent oscillation and damping, with boundary conditions defined using the medium and field traces.

2 Acoustic wave propagation

  • The acoustic domain is Ω⊂R^d for d = 2, 3, with a Lipschitz boundary partitioned into disjoint Neumann and absorbing portions.The density and bulk modulus are piecewise constant, bounded, and positive.Source: §2 · p. 3
  • The complex frequency σ = iω − s uses angular frequency ω and a positive damping parameter s.The medium parameters are assumed bounded and positive, ensuring well-defined bounded boundary traces.Source: §2 · p. 3
  • The first-order acoustic formulation uses reflecting Neumann boundaries and first-order absorbing boundary conditions, with the latter expressed in physical fields independently of frequency.The conditions are well defined using the stated regularity of the boundary traces of pressure and normal velocity.Source: §2 · p. 4

3 HDG discretization

The HDG discretization introduces volume and trace unknowns, which can be condensed to a global trace system or retained in an affine uncondensed system. Under stated damping, stabilization, and coefficient assumptions, the systems are well posed, and the solution map is analytic in frequency.

3 HDG discretization

  • The HDG method introduces a mesh-skeleton trace unknown alongside cell-based pressure and velocity variables.The trace space is defined on mesh faces, while the pressure and velocity spaces are defined on cells.Source: §3 · p. 4

3.2 Condensed formulation

  • Static condensation eliminates local volume variables, leaving a global system whose size is determined by face degrees of freedom.The condensed operator loses affine frequency dependence because it contains the inverse of the local matrix.Source: §3.2 · p. 5

3.3 Uncondensed formulation

  • The reduced formulation retains volume and trace variables independently, preserving affine frequency dependence and using distinct reduced spaces for the two fields.The separate spaces preserve the block structure and support the well-posedness analysis.Source: §3.3 · pp. 5–6

3.4 Well-posedness

  • For s > 0, τ ≥ τmin > 0, and µ ≥ µmin > 0, the HDG system is well posed and local matrices are invertible.The local invertibility makes static condensation well defined.Source: §3.4 · p. 7

3.5 Analyticity of the solution map

  • The solution map is analytic on Dσ because the local and global operators are invertible there.Analyticity permits enrichment of the reduced basis with frequency derivatives.Source: §3.5 · p. 8

4 Reduced basis approximation

Solutions over a frequency range can be approximated by low-dimensional spaces because nearby frequencies typically yield similar solution features. The reduced formulation uses separate low-rank spaces for volume and trace variables and exploits affine frequency dependence.

4 Reduced basis approximation

  • Nearby frequencies typically produce similar solutions, so a small set of snapshots can approximate solutions at other frequencies.The Kolmogorov n-width measures the best possible worst-case error from n-dimensional subspaces.Source: §4 · p. 9
  • The volume and trace variables use separate reduced spaces with ranks nU ≪ NU and nΛ ≪ NΛ.The associated bases are orthonormal in their respective L2 inner products.Source: §4 · p. 9
  • The uncondensed operator is affine in frequency, supporting an offline projection and reuse of the reduced model for online frequency queries.The right-hand side is assumed frequency-independent.Source: §4 · p. 9

4.1 Construction of the reduced basis

The reduced basis is built through adaptive frequency sampling, local derivative enrichment, and separate POD compression of volume and trace snapshots.

4.1 Construction of the reduced basis

  • A weak-greedy loop selects frequencies using a residual indicator, enriches snapshots with solution derivatives, and compresses them by weighted POD.The three ingredients are combined in the basis-construction procedure.Source: §4.1 · p. 9

4.1.1 Enriching the reduced basis with frequency derivatives

  • Frequency derivatives add local slope information to snapshots at each selected training frequency.The derivatives reuse factorizations and each costs one additional back-substitution.Source: §4.1.1 · p. 10

4.1.2 Greedy selection of training points

  • The residual-based indicator selects the least well-approximated training frequency, and termination is based on the true relative L2 error there.The indicator is a computable surrogate for the true error; its threshold is not used directly because it differs by a scaling factor.Source: §4.1.2 · pp. 10–11

4.1.3 POD compression

  • Volume and trace snapshots are compressed separately by POD using their L2 variational norms.Separate compression keeps the reduced dimensions decoupled and preserves the reduced system's block structure.Source: §4.1.3 · p. 11

4.2 Galerkin projection

Galerkin projection of the uncondensed HDG system preserves affine frequency dependence and yields a uniformly stable reduced model. Its online cost depends on reduced dimensions rather than the high-fidelity dimension.

4.2 Galerkin projection

  • Projecting the condensed system destroys explicit affine frequency dependence, while the uncondensed formulation supports an efficient offline/online decomposition.Direct projection of the condensed system would also require elementwise recovery of volume unknowns online.Source: §4.2 · p. 11
  • Online reduced-model cost is independent of N, with inversion costs scaling as O(nU^3) and O(nΛ^3).The reduced operators are projected and stored offline, then assembled cheaply for each query frequency.Source: §4.2 · p. 12
  • The reduced model retains the full-order inf-sup constant, independent of the reduced dimensions, when volume and trace spaces are separate.The separate spaces keep the reduced space invariant under the transformation used in the energy argument.Source: §4.2 · p. 12

4.3 Error analysis

The analysis first establishes a coupled quasi-optimal error bound, then separates trace and volume errors to quantify the effects of noncommuting projection and static condensation. A lifting condition can eliminate the resulting gap terms at selected frequencies, while increasing the volume reduced-space dimension.

4.3.1 Coupled error bound

  • The joint reduced solution is quasi-optimal: its error is bounded by the best approximation of the high-fidelity solution in the joint reduced space.The bound follows from inf-sup stability and continuity of the sesquilinear form.Source: §4.3.1 · p. 12

4.3.2 Decoupled error bounds

  • Projection and static condensation do not generally commute, so the reduced trace solution and the high-fidelity condensed trace solution satisfy different reduced problems.The trace estimate accounts for the gap between the reduced problem and the projection of the condensed high-fidelity problem.Source: §4.3.2 · p. 13
  • The trace error bound combines best approximation in the trace space with the gap between the reduced and projected condensed problems.The decoupled analysis uses coercivity and the Schur-complement structure because classical saddle-point MOR analysis does not apply when C ≠ B^T and L ≠ 0.Source: §4.3.2 · p. 13
  • The volume error is controlled by its best-approximation error and by the trace error propagated through the coupling terms.The coupling constant is Ccoupl = γc/αa.Source: §4.3.2 · p. 14
  • The gap terms are bounded by distances of local volume outputs from the primal reduced space; including the relevant outputs makes the reduced ROM and projected condensed system coincide.The primal lifting condition also makes both gap terms vanish at the selected training frequencies.Source: §4.3.2 · pp. 14–15
  • Lifting can add up to n_lift n_Λ dimensions to the primal space; selecting training frequencies and POD-compressing lifted vectors can control this growth.The added basis vectors are inexpensive to compute, but the increase in reduced dimension may become significant.Source: §4.3.2 · p. 15

4.4 Implementation

The offline construction combines frequency derivatives, residual-driven greedy sampling, weighted POD compression, and Galerkin projection. Its implementation exploits HDG element locality and distributed computation, while online evaluation solves only the reduced system.

4.4 Implementation

  • The offline procedure builds an initial model from sampled frequencies and solution derivatives, then greedily enriches it at the least well-approximated frequency until the selected true relative error meets tolerance.Snapshots and their derivatives are compressed by weighted POD, and the reduced operators are assembled by element-wise projection.Source: §4.4 · pp. 15–16

Implementation.

  • For each selected frequency, static condensation solves globally for trace unknowns; volume unknowns are recovered locally, and derivative solves reuse existing factorizations.No global volume system is assembled or factorized.Source: Implementation · p. 16
  • The implementation constructs reduced operators from element-wise contributions and never assembles the global uncondensed volume operator.MPI distributes elements and reduced contributions, while OpenMP parallelizes element operations; POD keeps volume snapshots distributed rather than gathering them on one process.Source: Implementation · pp. 16–17
  • After reduced Gram matrices are built, the residual indicator is evaluated over the training set using reduced solves and residual norms.Online evaluation assembles and solves the reduced system for each requested frequency.Source: §4.4 · p. 16, Implementation · p. 17

5 Two-dimensional model: Marmousi

The Marmousi experiments evaluate frequency-derivative enrichment, lifting, the residual indicator, greedy convergence, computational cost, accuracy, and frequency-band effects. The final ROM remains accurate on independent test frequencies while providing substantial online speedups.

5 Two-dimensional model: Marmousi

  • The 2D Marmousi study uses a fixed heterogeneous medium and evaluates the ROM construction through basis enrichment, error bounds, indicator behavior, costs, and out-of-sample accuracy.The HDG discretization uses polynomial degree p = 5 over the frequency band [8.5, 10] × 2π.Source: §5 · pp. 17–18

5.1 Enhancement of the ROB via frequency derivatives

  • Frequency-derivative enrichment improves mean projection error by two orders of magnitude and reduces maximum volume ROM error from 5.15 × 10−2 to 9.83 × 10−4.Derivative computation adds 1.97 s, or 8.2% of the 24.08 s training cost; effectivity remains approximately 1.6–2.2 with or without enrichment.Source: §5.1 · p. 19

5.2 Analysis of error bounds

  • With lifting, the trace error stays 1.01 × 10−4 across six volume-basis configurations despite more than an order-of-magnitude variation in volume projection error.The gaps fall to numerical precision; without lifting, sufficiently rich volume snapshot bases capture enough of the trace lifting that explicit enforcement is unnecessary in this test case.Source: §5.2 · p. 20

5.3 Behavior of the error indicator

  • The residual indicator selects one of the two critical frequencies but ranks them oppositely to the true error; its selected frequency still captures 97% of the maximum error.The indicator is intended to identify frequencies with large ROM error, not to estimate relative-error magnitude sharply.Source: §5.3 · p. 21, Figure 2 · p. 21

5.4 Convergence of the greedy algorithm

  • Across three greedy iterations, the relative error at selected points falls from 7.11 × 10−4 to 2.16 × 10−4, while indicator evaluation over 100 candidates costs 3 × 10−3 s per iteration.Both the indicator and true relative error decrease monotonically.Source: §5.4 · p. 21, Figure 3 · p. 22

5.5 Offline and online costs

  • 1952× online speedup offsets offline training after approximately 72 queries; each ROM evaluation averages 3.3 × 10−4 s versus 0.637 s for a full-order solve.The final greedy ROM uses three iterations and reduced bases of size nU = nΛ = 65.Source: Table 3 · p. 22, §5.5 · p. 22

5.6 Accuracy: FOM vs ROM

  • On 50 independent test frequencies, the final ROM keeps maximum relative errors below the 5 × 10−4 greedy tolerance, with mean pressure error 4.41 × 10−5.The relative-to-projection error ratio stays close to 1.8 across fields.Source: Table 4 · p. 22, §5.6 · p. 22

5.7 Influence of the frequency band

  • The frequency-band tests stay below the prescribed maximum-error tolerance, while wider, higher-frequency bands require larger bases and reduce online speedup from 1952× to 879×.The wider bands require 12 and 15 greedy iterations rather than 3, and bases of size 81 and 89 rather than 65.Source: §5.7 · p. 23, Table 5 · p. 23

6 Three-dimensional model

In the 3D heterogeneous test, frequency-derivative enrichment improves the reduced basis, while lifting decouples trace errors from the volume basis. The greedy ROM meets the prescribed error tolerance across tested frequency bands with large online speedups.

6 Three-dimensional model

  • The 3D test uses a degree-3 HDG discretization on 148,422 tetrahedral elements, with 11,873,760 volume unknowns and 3,059,460 trace unknowns.The heterogeneous medium is tested over the frequency band [8, 11] × 2π.Source: §6 · pp. 23–24, Figure 5 · p. 24

6.1 Enhancement of the ROB via frequency derivatives

  • Derivative enrichment lowers the mean projection error for uh from 1.30×10^-3 to 1.80×10^-5 and the maximum ROM error from 9.05 × 10^-3 to 1.01 × 10^-4.The effectivity index remains approximately 1.5–1.8 with or without enrichment, supporting the Galerkin reduced solution's quasi-optimality.Source: §6.1 · p. 24, Table 6 · p. 25

6.2 Analysis of error bounds

  • With lifting, the trace relative error remains 1.8952 × 10^-5 across six volume-basis configurations, while the gap terms fall to numerical precision.Without lifting, trace error varies with the volume basis; as that basis becomes more accurate, gap terms decrease and lifted and unlifted errors become comparable.Source: §6.2 · p. 25, Table 7 · p. 25

6.3 Behavior of the error indicator

  • The residual indicator follows relative-error variations across the band and selects the same maximum-error frequency, ω = 10.92.It is a frequency-selection tool, not a sharp estimate of relative-error magnitude; the relative errors of uh and λh coincide.Source: §6.3 · p. 26, Figure 6 · p. 26

6.4 Convergence of the greedy algorithm

  • After five greedy iterations, the relative error decreases from 8.40 × 10^-4 to 4.85×10^-5, while evaluating the indicator over 100 candidates costs 1.1×10^-3 s per iteration.Both the indicator and relative error decrease monotonically.Source: §6.4 · p. 26, Figure 7 · p. 26

6.5 Offline and online costs

  • The 39-mode ROM costs 1.11 × 10^-4 s per online frequency evaluation versus 48.25 s for a full-order solve, a speedup of approximately 4.3 × 10^5.The offline training cost is recovered after approximately 27 online queries.Source: §6.5 · p. 27, Table 8 · p. 27

6.6 Accuracy: FOM vs ROM

  • On an independent test grid of 50 frequencies over Dω = [8, 11] × 2π, maximum relative errors for all fields remain below greedy tolerance 10^-4, with relative-to-projection error ratios close to 1.6.For pressure, the mean and maximum relative errors are 1.04 × 10^-5 and 4.03 × 10^-5, respectively.Source: Table 9 · p. 27, §6.6 · p. 27

6.7 Influence of the frequency band

  • Across three frequency bands, the greedy ROM keeps maximum test-set relative error below 10^-4 with 31–39 modes and reaches speedups from about 5 × 10^5 to 1.7 × 10^6.The offline cost is dominated by high-fidelity solves (83%–97%), and training is recovered after approximately 20, 27, and 17 online queries, respectively.Source: §6.7 · p. 28, Table 10 · p. 28

7 Conclusion

The paper develops a reduced-order framework tailored to HDG, with stability and error analysis for the interaction between projection and static condensation. Tests in two and three dimensions show accurate, fast reduced solves, while extensions to multiple sources and changing media remain future work.

7 Conclusion

  • The framework combines residual-indicator greedy frequency sampling, frequency-derivative enrichment, and separate POD compression for volume and trace unknowns.It establishes well-posedness for the damped acoustic HDG and reduced problems with the same stability constant, using decoupled reduced bases.Source: §7 · p. 28
  • Projection and static condensation do not commute in general; the analysis gives error estimates and a lifting condition that eliminates the resulting gap.Numerical tests show lifting decouples trace error from volume-basis approximation error, while sufficiently rich volume bases make lifting unnecessary for the tested problems.Source: §7 · p. 28
  • For the 3D problem with approximately 12 million volume unknowns, 39 modes achieve relative errors below 10^-4 with an online speedup of approximately 4 × 10^5.For Marmousi, 65 modes per variable achieve errors below 5 × 10^-4 and an online speedup of approximately 2 × 10^3.Source: §7 · p. 28
  • The current reduced spaces are constructed for a single source and fixed medium; multiple-source configurations and medium variability during inversion remain future-work challenges.The paper also identifies elastic and coupled acoustic–elastic wave propagation as future extensions.Source: §7 · p. 29
Loading 2610.20145v1…