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.
Abstract.
We present a new framework for the fast solution of inhomogeneous elliptic boundary value problems in domains with smooth boundaries. High-order solvers based on adaptive box codes or the fast Fourier transform can efficiently treat the volumetric inhomogeneity but require care to be taken near the boundary to ensure that the volume data is globally smooth. We avoid function extension or cut-cell quadratures near the boundary by dividing the domain into two regions: a bulk region away from the boundary that is efficiently treated with a truncated free-space box code, and a variable-width boundary-conforming strip region that is treated with a spectral collocation method and an accompanying fast direct solver. Particular solutions in each region are then combined with Laplace layer potentials to yield the global solution. The resulting solver has an optimal computational complexity of [Formula: see text] for an adaptive discretization with [Formula: see text] degrees of freedom. With an efficient two-dimensional (2D) implementation we demonstrate adaptive resolution of volumetric data, boundary data, and geometric features across a wide range of length scales, to typically 10-digit accuracy. The cost of all boundary corrections remains small relative to that of the bulk box code. The extension to 3D is expected to be straightforward in many cases because the strip “thickens” an existing boundary quadrature.
Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/danfortunato/fully-adaptive-poisson and in the supplementary materials ( fully-adaptive-poisson-master.zip [170KB]). [Formula: see text]
D. Fortunato, David B. Stein, Alex H. Barnett· SIAM Journal on Scientific C...· 0 citations
In this paper, we propose high-order weighted essentially non-oscillatory (WENO)-based semi-implicit Newton-type Gauss-Seidel Lax-Friedrichs fast sweeping methods for solving the generalized Eikonal equation arising in wave propagation through a moving fluid. Building upon the Newton-type framework of Li and Qian (2020), which updates the solution line-wise using Newton's method with a tridiagonal and strictly diagonally dominant Jacobian, we extend the local solver to fifth-, seventh-, and ninth-order accuracy by incorporating high-order WENO approximations of the spatial derivatives into the numerical Hamiltonian. Three alternating sweeping strategies, namely column-wise, row-wise, and column-row-wise, are considered within the Gauss-Seidel iteration framework. Numerical examples in both two and three spatial dimensions demonstrate the efficiency and accuracy of the proposed schemes.
We present a nodal discontinuous Galerkin method for solving the Bhatnagar-Gross-Krook (BGK) kinetic equation on multi-dimensional, unstructured grids. The method uses implicit, sweep-based solvers and a moment-preserving projection of the Maxwellian source to enable high-order accuracy in time while avoiding restrictive time steps imposed by boundary layers and other geometry-induced features. We verify that the method is correct in the continuum limit by comparing to closed-form and high-order solutions of the Sod shock problem on 2 and 3D unstructured grids. Linear L2 stability is demonstrated for a B-stable diagonally implicit Runge-Kutta method of third order. The solver uses a hybrid parallel scheme based on spatial domain decomposition with local sweeps performed on CPU and GPU hardware. Platform-portability is demonstrated through the development of new GPU-friendly, graph-based sweep algorithms that are implemented using the Kokkos performance portability library and achieve greater than 20 times speedup on NVIDIA H100 GPUs compared to 64-core AMD EPYC 9654 CPUs. Finally, we show results on the Frontier supercomputer at the Oak Ridge Leadership Computing Facility for a boundary value problem with 2.77 trillion phase space degrees of freedom that executed on 1536 nodes utilizing 6144 AMD MI250X GPUs.
T. Evans, Ryan S. Glasby, Cory D. Hauck et al.· 0 citations
Time-fractional diffusion equations describe subdiffusive transport in heterogeneous media, but their numerical treatment is complicated by the nonlocal Caputo derivative and by the weak singularity that the solution develops at the initial time. We study a global-in-time discretization that combines spectral collocation in time—on the fractional power basis {tℓα}ℓ=0N, evaluated at Chebyshev–Gauss–Lobatto nodes, which reproduces the leading terms of the singular expansion of the solution—with a second-order conservative finite-difference stencil in space that uses harmonic averaging of the diffusivity at the cell faces and therefore remains accurate across discontinuous media. The resulting fully discrete problem is a large, nonsymmetric, dense-in-time linear system whose two-norm condition number grows like the inverse square of the spatial mesh size, so that Krylov subspace iteration without preconditioning stalls under refinement. Exploiting the Kronecker sum structure of the discrete operator, we build a preconditioner by fast diagonalization of the spatial factor through the discrete sine transform. For constant diffusivity the preconditioner reproduces the operator exactly and yields a direct solver; for variable diffusivity it is spectrally equivalent to the operator, and we prove that the eigenvalues of the preconditioned system cluster in a disk centered at one whose radius depends only on the coefficient contrast, and not on the mesh, the number of temporal degrees of freedom, or the fractional order. Numerical experiments in one and two space dimensions confirm second-order spatial accuracy and a preconditioned iteration count that stays flat—twelve iterations from M=32 up to M=1024 in one dimension and eleven up to M=256 per direction in two—while the unpreconditioned count grows by more than two orders of magnitude. In time, the accuracy is spectral until round-off in the ill-conditioned Vandermonde matrix of the power basis takes over: the barrier is reached at N=9,10,13 for α=0.3,0.5,0.7, where the attainable error is about 10−6. A benchmark against the L1 scheme on uniform and graded meshes, the Alikhanov L2-1σ scheme and Grünwald–Letnikov convolution quadrature quantifies when the global approach pays: on forced problems and on modes with κλTα≲2 it reaches a prescribed accuracy one to two orders of magnitude faster and with several times less memory, while for strongly damped modes the fractional power basis converges only algebraically and graded time marching is preferable below a relative error of 10−2.
An efficient aggregation-based multigrid approach is developed in this work to solve linear algebraic systems generated by discretizing the Landau-Lifshitz equation in micromagnetics. The discretization couples the first-order backward differentiation formula with first-order extrapolation, and the underlying equation describes magnetization dynamics within micromagnetic simulations. For these linear systems, standard iterative methods and conventional multigrid solvers intended for symmetric problems suffer from deteriorated convergence under mesh refinement and poor computational performance. Existing algorithms struggle to handle the system's intrinsic features, namely its sparse spectral properties and non-symmetric matrix structure. To overcome these limitations, a multigrid framework is constructed, where smoothing and coarse-grid correction operators are customized for the target linear system. Numerical tests confirm the robustness and mesh-independent convergence of the resulting method. Compared with well-established Krylov-subspace solvers and conventional multigrid techniques, the aggregation-based multigrid solver cuts down iteration counts and overall computational cost considerably, while producing faithful representations of magnetization dynamics.
We propose an adaptive multigrid (MG) method for discontinuous Galerkin formulations of elliptic problems using Brandt's full approximation scheme (FAS). Unlike common approaches, this method achieves local $hp$-refinement of hexahedral meshes without the need for hanging nodes. The core component of the FAS-MG method is an overlapping Schwarz smoother, which is optionally accelerated by a Krylov method. This smoother is designed for unstructured curvilinear meshes but maintains a tensor-product structure for fast diagonalization. Numerical experiments demonstrate the exceptional efficiency of the FAS-MG method. Dedicated studies confirm its robustness against high aspect ratios, element deformation, and irregular mesh topology. We also verify its capability for dynamic parallel mesh adaptation using the wave-front benchmark of \v{C}erven\'y, Dobrev, and Kolev (SIAM J. Sci. Comp. 41, 2019). Finally, we present preliminary results of extending the method to incompressible Navier-Stokes problems.
Jörg Stiller· 0 citations
We use cookies to run the site and, with your consent, for analytics and to show ads.
See our Cookie Policy.
