Accurate simulation of three-dimensional time-harmonic wave propagation over many wavelengths requires control of phase error and efficient solution of large indefinite systems. We develop a three-grid solver based on the compact 27-point interpolated optimized finite-difference (IOFD) discretization. Its wavenumber-dependent stencil supports a fine-grid resolution of six points per shortest wavelength and an unshifted physical correction on the \(2h\) grid at only three points per shortest wavelength. The method retains unshifted IOFD operators on the \(h\) and \(2h\) grids, while a complex-shifted \(2h\)--\(4h\) auxiliary cycle preconditions a factorization-free iterative approximation of the coarse inverse. Restricting the shift to this auxiliary cycle preserves the propagative character of the coarse correction. Comparison with the outgoing Green function confirms phase and relative-amplitude accuracy on a sequence of meshes up to \(6144^3\), with the largest problem spanning approximately 1024 wavelengths per coordinate. The same fixed solver configuration retains robust convergence across smooth, discontinuous, high-contrast, and geophysical velocity models and exhibits scalable parallel performance. In particular, a problem spanning approximately 340 wavelengths in each coordinate direction is solved in 18.1 seconds on just 64 NVIDIA A100 GPUs.
Three-dimensional time-harmonic Maxwell simulations generate massive complex indefinite systems whose mesh coarsening is strictly limited by phase accuracy. Although matrix-free finite element kernels utilize GPU throughput efficiently, standard multilevel solvers are ultimately bottlenecked by the memory and communication costs of exact coarse-grid factorizations. We present a fully matrix-free, factorization-free three-grid preconditioner for curl-conforming N{\'e}delec discretizations with perfectly matched layers (PML) and optimally blended quadrature. The method employs an outer FGMRES to solve the unshifted fine-grid equation, while an intermediate-grid correction is computed by a fixed-work FGMRES preconditioned with a complex-shifted $2h$--$4h$ cycle. This strategically confines the complex shift to an auxiliary preconditioner, preserving the physical Maxwell operator. A local Fourier analysis derives the blended Maxwell branches and compatible edge transfers, identifying robust shift and Jacobi damping parameters. Validated against the analytical Maxwell Green tensor, our approach demonstrates extreme scalability: using a single solver configuration, both homogeneous and highly heterogeneous systems with approximately 10.89 billion complex edge unknowns are solved in 42.0--72.0 seconds on just 64 NVIDIA A100 GPUs.
This article presents a fast, high-order Nystr\"om solver for the two-dimensional Lippmann--Schwinger equation arising from time-harmonic scattering by penetrable, piecewise-smooth heterogeneous media. Relying on high-order evaluation of the Newtonian potential on unstructured grids adapted to interfaces of discontinuity, the methodology achieves high-order accuracy using existing fast algorithms such as the fast multipole method. As an iterative method the solver exhibits rapid convergence when coupled to a preconditioning strategy that exploits a class of structured-grid solvers---fast solution methods offering quasi-linear time and memory complexity and nearly-constant iteration counts, but long limited in accuracy. The preconditioning strategy couples the Nystr\"om discretization---given by high-order quadrature nodes over a (curved) unstructured mesh conforming to the support of the spatially varying contrast---to a uniform Cartesian grid underlying the fast preconditioner via a pair of transfer operators. The resulting preconditioner inherits the frequency-robust behavior of its Cartesian counterpart without sacrificing the geometric flexibility and high-order accuracy of the unstructured discretization. We prove that invertibility of the proposed preconditioner holds under explicit conditions on the mesh sizes and on the Cartesian preconditioner. Numerical experiments demonstrate that the preconditioned system requires significantly fewer GMRES iterations than its unpreconditioned counterpart, with iteration counts almost independent of mesh size and wavenumber, and illustrate the method's robustness for inhomogeneities with piecewise-smooth refractive indices and jump discontinuities across interfaces.
T. Anderson, Juan Burbano-Gallegos, Luiz M. Faria et al.· 0 citations
High-precision magnetohydrodynamics (MHD) simulations are essential for understanding the mechanism of solar eruptions and coronal mass ejections (CMEs), where the slow accumulation and impulsive release of magnetic energy are inherent to the corona's extremely low resistivity. However, achieving such simulations with conventional second-order schemes requires extremely high grid resolutions, making routine calculations prohibitively expensive. Here, we present a new code AMR–CESE4–MHD by implementing, for the first time, a fourth-order conservation-element and solution-element (CESE) scheme on general curvilinear geometries with adaptive mesh refinement (AMR) grid for MHD simulations. We propose a smoothness indicator based on second-order partial derivatives to guide the mesh refinement. For solution prolongation and restriction that are required during the grid refinement and coarsening, we use Taylor series expansions based on spatial derivatives available in the CESE solution to preserve solution accuracy during grid changes. We also propose a magnetic divergence control method by correcting the spatial derivatives of the magnetic field so that the magnetic-field divergence and its spatial derivatives evaluated from these auxiliary variables vanish. A series of benchmark tests validate the implementation. The scheme achieves its designed fourth-order accuracy, and the AMR criteria automatically refine near discontinuities. The derivative-control method keeps the ∇·B errors comparable to the combined approach of Powell source term and Marder diffusion term, while showing better conservation properties. The code produces consistent results in both Cartesian and curvilinear coordinates, and is ready for applications in solar eruption and CME modeling.
Chaowei Jiang, Ling Zhang, Xueshang Feng et al.· The Physics of Fluids· 0 citations
We consider plane-wave discontinuous Galerkin (PWDG) approximations of the Helmholtz equation with adaptive local propagation directions. The directions are chosen by minimizing a weighted residual on the mesh skeleton. We study two formulations: in Part A the PWDG system is solved for the coefficients at fixed directions, while in Part B the same residual is minimized jointly over coefficients and directions, with the coefficients eliminated by variable projection. Complex angles unify propagating and evanescent Trefftz waves. The discrete system is normalized in the Trefftz-DG norm, and a local Cauchy-trace Gramian is used to remove numerically dependent directions. On straight edges the residual integrals are exact. For all-Dirichlet problems, the residual equals the squared DG error and provides local adaptive indicators. We prove local quadratic growth of the reduced direction functional near an identifiable zero-residual solution of fixed rank. Numerically, an exact circular DtN test shows that the dominant phase direction of a Hankel wave can be recovered accurately even when a small plane-wave fan gives a larger field error. We also separate the effects of the coefficient basis, trace cutoff, and trace-spectrum arithmetic on the high-$p$ DtN floor. Finally, for finite sums of plane waves, residual-based ENRICH--MOVE continuation recovers all directions to roundoff for $M=1,\ldots,19$. At $M=20$ the automatic birth step enters a false basin, whereas a nearby birth again reaches roundoff. The results indicate that direction adaptation is most effective for solutions of low directional complexity.
Accurate finite-difference modeling of seismic wave propagation with irregular free surfaces remains challenging. The standard staggered-grid (SSG) finite-difference method exhibits significant numerical dispersion when simulating surface waves, necessitating dense spatial sampling to maintain accuracy. In contrast, the Fully Staggered Grid (FSG) or Lebedev scheme, enhanced with Mimetic Finite Difference (MFD) operators, enables high-precision modeling with relatively few grid points (approximately eight points per minimum wavelength). However, FSG requires simultaneous computations on multiple staggered layouts, resulting in a substantial increase in computational cost. In many seismic imaging applications, using an FSG everywhere in the model domain is unnecessary and often inefficient. Motivated by this observation and inspired by overset grid techniques, we propose a hybrid-grid strategy to balance accuracy and efficiency. Near irregular free-surface boundaries, FSG is applied in general curvilinear coordinates to ensure accurate surface wave simulation. In deeper regions where topographic effects diminish, a more efficient SSG scheme is used. An overset grid strategy is introduced to ensure smooth wavefield transitions between the two regions. Numerical experiments demonstrate that the proposed method preserves computational accuracy while significantly reducing runtime, achieving a nearly 50% reduction in computation time. Furthermore, the use of a staggered grid eliminates source injection errors in the transition zone that typically arise when applying the MacCormack scheme in collocated grids, thereby simplifying the source implementation strategy.
Zixiao Zhang, P. Yong, Jianping Huang· Geophysics· 0 citations
We analyze coefficient design in a one-parameter family of explicit palindromic electric--magnetic splittings for the time-domain Maxwell equations. After fourth-order staggered spatial discretization, the Fourier amplification matrix depends on the single scalar $g_2=a(1-2a)/2$. We prove that $a=1/4$ is the unique real coefficient maximizing the spectral CFL interval, with threshold $s_*=12/(7\sqrt d)$. We then identify a real-coefficient obstruction to higher phase accuracy: cancellation of the leading temporal phase defect requires $g_2=1/12$, whereas every real member satisfies $g_2\le 1/16$. The resulting complex-conjugate coefficients give fourth-order temporal phase accuracy for each fixed semidiscrete Fourier mode and have threshold $6\sqrt3/(7\sqrt d)$, while the complete field update remains globally second order in time. For real Maxwell data, the physical output is the real projection of the complex trajectory; this projection is branch independent and preserves the second-order error bound. We further give an exactly equivalent doubled real-arithmetic realization, which clarifies the role of the auxiliary imaginary component without changing the numerical method. A semidiscrete convergence result and numerical experiments confirm the distinction between stability optimization and phase optimization.
Hui Duan, Hongliang Li, Lunzhong Guo· 0 citations
We use cookies to run the site and, with your consent, for analytics and to show ads.
See our Cookie Policy.
