It is well known that numerical solutions of Helmholtz boundary-value problems suffer from pollution (dispersion) errors that require much finer grids than expected for high frequencies. In this article a dispersion error analysis is presented that leads to an explicit rule-of-thumb formula to estimate the points-per-wavelength required for a p-th order accurate finite difference approximation for Helmholtz problems that accounts for pollution errors. The rule of thumb only depends on the wave-number (or frequency) of the periodic forcing, the size of the domain, and the desired relative error tolerance. Several numerical examples on simple and complex domains and with closed or open boundaries show that the rule of thumb is a useful predictor of the error. Having such a rule is important in practice to either avoid under-resolving the problem and obtaining the wrong answer, or over-resolving the problem resulting in a much more expensive computation.

We present a numerical study of eigenvector deflation as a means of accelerating the WaveHoltz method for solving the Helmholtz equation. For energy-conserving (Dirichlet or Neumann) boundary conditions the WaveHoltz fixed-point iteration converges slowly at high frequency, requiring approximately 𝒪(ω^2d) iterations in d dimensions. We show that deflating the eigenvectors whose eigenvalues lie nearest the driving frequency substantially reduces iteration counts, and we examine two ways of incorporating the eigenvectors: direct eigenvector deflation (DEVD), in which the forcing and iterate are projected against the deflation set, and augmented-Krylov eigenvector deflation (AUKED) using deflated conjugate gradient (DCG), augmented GMRES (AGMRES), and augmented (recycled) BICGSTAB (ABICGSTAB). The required eigenpairs can be computed efficiently with the EigenWave approach, and we demonstrate, in two dimensions, that when the number of deflation vectors grows quadratically with ω the asymptotic convergence rate remains essentially constant. Because the eigenvectors on structured grids are naturally represented as matrices, we further apply SVD-based compression to reduce their storage. Numerical experiments on single curvilinear grids discretized with summation-by-parts operators, and on overset grids illustrate the robustness and efficiency of the approach, with the deflated solver breaking even against the undeflated solver after as few as two right-hand sides, when accounting for the cost of precomputing the eigenvectors.

An algorithm, named EigenWave, is described to compute eigenvalues and eigenvectors of elliptic boundary value problems. EigenWave is based on the recently developed WaveHoltz scheme and solves a related time-dependent wave equation as part of an iteration involving an initial condition of the time-dependent problem. At each iteration, the solution to the wave equation is filtered in time resulting in updated initial data. The time-filter is designed to promote the contributions of eigenmodes in the solution near a chosen target frequency (target eigenvalue). The ability to choose an arbitrary target frequency enables the computation of eigenvalues anywhere in the spectrum, without the need to invert an indefinite matrix, as is common with other approaches. Furthermore, the iteration can be embedded within a matrix-free Arnoldi algorithm, which enables the efficient computation of multiple eigenpairs near the target frequency. For efficiency, the time-dependent wave equation can be solved with implicit time-stepping and only about 10 time-steps per-period are needed, independent of the mesh spacing. When the (definite) implicit time-stepping equations are solved with a multigrid algorithm, the cost of the resulting EigenWave scheme scales linearly with the number of spatial grid points N as the mesh is refined, giving an optimal O(N) algorithm. The approach is demonstrated by finding eigenpairs of the Laplacian in complex geometry using finite-difference approximations on overset grids. Results in two and three space dimensions are presented using second-order and fourth-order accurate approximations.

We develop efficient and high-order accurate solvers for the Helmholtz equation on complex geometry. The solvers are based on the WaveHoltz algorithm, which computes solutions of the Helmholtz equation by time-filtering solutions of the wave equation. The approach avoids the need to invert an indefinite matrix which can cause convergence difficulties for many iterative solvers for indefinite Helmholtz problems. Complex geometry is treated with overset grids which use Cartesian grids throughout most of the domain together with curvilinear grids near boundaries. The basic WaveHoltz fixed-point iteration is accelerated using GMRES and also by a deflation technique using a set of precomputed eigenmo des. The solution of the wave equation is solved efficiently with implicit time-stepping using as few as five time-steps per period, independent of the mesh size. When multigrid is used to solve the implicit time-stepping equations, the cost of the resulting WaveHoltz scheme scales linearly with the total number of grid points N (at fixed frequency) and is thus optimal in CPU time and memory usage as the mesh is refined. Numerical results are given for problems in two and three dimensions, to second- and fourth-order accuracy, and they show the potential of the approach to solve a wide range of large-scale problems.

We develop and analyze a new approach for simultaneously computing multiple solutions to the Helmholtz equation for different frequencies and different forcing functions. The new Multi-Frequency WaveHoltz (MFWH) algorithm is an extension of the original WaveHoltz method and both are based on time-filtering solutions to an associated wave equation. With MFWH, the different Helmholtz solutions are computed simultaneously by solving a single wave equation combined with multiple time filters. The MFWH algorithm defines a fixed-point iteration which can be accelerated with Krylov methods such as GMRES. The solution of the wave equation can be efficiently solved with either explicit time-stepping or implicit time-stepping using as few as five time-steps per period. When combined with an O(N) solver for the implicit equations, such a multigrid, the scheme has an O(N) solution cost when the frequencies are fixed and the number of grid points N increases. High-order accurate approximations in space are used together with second-order accurate approximations in time. We show how to remove time discretization errors so that the MFWH solutions converge to the corresponding solutions to the discretized Helmholtz problems. Numerical results are given using second-order accurate and fourth-accurate discretizations to confirm the convergence theory.

New implicit and implicit-explicit time-stepping methods for the wave equation in second-order form are described with application to two and three-dimensional problems discretized on overset grids. The implicit schemes are single step, three levels in time, and based on the modified equation approach. Second and fourth-order accurate schemes are developed and they incorporate upwind dissipation for stability on overset grids. The fully implicit schemes are useful for certain applications such as the WaveHoltz algorithm for solving Helmholtz problems where very large time-steps are desired. Some wave propagation problems are geometrically stiff due to localized regions of small grid cells, such as grids needed to resolve fine geometric features, and for these situations the implicit time-stepping scheme is combined with an explicit scheme: the implicit scheme is used for component grids containing small cells while the explicit scheme is used on the other grids such as background Cartesian grids. The resulting partitioned implicit-explicit scheme can be many times faster than using an explicit scheme everywhere. The accuracy and stability of the schemes are studied through analysis and numerical computations.

The formulation of finite difference approximations is a classical problem in numerical analysis. In this article, we consider difference approximations that are based on a series expansion in powers of the second undivided difference. Each additional term in the series increases the order of accuracy by two. These expansions are useful in a variety of contexts such as in the development of modified equation schemes, the design of high-order accurate energy stable discretizations, and error analysis of certain finite element or finite difference schemes. Here, we provide closed form expressions for the coefficients in the series expansions for derivatives of all orders. We also provide some short recursions defining the series coefficients, and formulae for the stencil coefficients in standard difference approximations. The series expansions are used to show some useful properties of the Fourier symbols of difference approximations and to derive rules of thumb for the number of points-per-wavelength needed to achieve a given error tolerance when solving wave propagation problems involving higher spatial derivatives.

High order accurate Hermite methods for the wave equation on curvilinear domains are presented. Boundaries are treated using centered compatibility conditions rather than more standard one-sided approximations. Both first-order-in-time (FOT) and second-order-in-time (SOT) Hermite schemes are developed. Hermite methods use the solution and multiple derivatives as unknowns and achieve space-time orders of accuracy 2m-1 (FOT) and 2m (SOT) for methods using (m+1)^d degree of freedom per node in d dimensions. The compatibility boundary conditions (CBCs) are based on taking time derivatives of the boundary conditions and using the governing equations to replace the time derivatives with spatial derivatives. These resulting constraint equations augment the Hermite scheme on the boundary. The solvability of the equations resulting from the compatibility conditions are analyzed. Numerical examples demonstrate the accuracy and stability of the new schemes in two dimensions.

We develop efficient and high-order accurate solvers for the Helmholtz equation on complex geometry. The schemes are based on the WaveHoltz algorithm which computes solutions of the Helmholtz equation by time-filtering solutions of the wave equation. The approach avoids the need to invert an indefinite matrix which can cause convergence difficulties for many iterative solvers for indefinite Helmholtz problems. Complex geometry is treated with overset grids which use Cartesian grids throughout most of the domain together with curvilinear grids near boundaries. The basic WaveHoltz fixed-point iteration is accelerated using GMRES and also by a deflation technique using a set of precomputed eigenmodes. The solution of the wave equation is solved efficiently with implicit time-stepping using as few as five time-steps per period, independent of the mesh size. The time-domain solver is adjusted to remove dispersion errors in time and this enables the use of such large time-steps without degrading the accuracy. When multigrid is used to solve the implicit time-stepping equations, the cost of the resulting WaveHoltz scheme scales linearly with the total number of grid points N (at fixed frequency) and is thus optimal in CPU-time and memory usage as the mesh is refined. A simple rule-of-thumb formula is provided to estimate the number of points-per-wavelength required for a p-th order accurate scheme which accounts for pollution (dispersion) errors. Numerical results are given for problems in two and three space dimensions, to second and fourth-order accuracy, and they show the potential of the approach to solve a wide range of large-scale problems.

In this study, we employ a hydroelastic analysis to investigate the motion response of large ship hulls, treating them as either Euler–Bernoulli or Timoshenko beams to consider the influence of shear effects. To enhance clarity, we provide a detailed derivation of the equation of motion within the framework of Timoshenko beams. This work solves forward-speed radiation and diffraction problems for flexible bodies, utilizing linearized potential flow theory including generalized modes. Two common base-flow models, the Neumann-Kelvin and double-body base flows, are included in the solver. The solution is numerically implemented in the high-order finite difference and open-source seakeeping solver Oceanwave3D-seakeeping. The numerical implementation involves the discretization of the geometry using overlapping, boundary-fitted grids, which has been validated by three examples involving a barge and two Wigley hulls. The influence of the Doppler shift due to forward speed on the hydroelastic motion response is also discussed. Through the integration of hydroelastic analysis using potential flow theory and advanced numerical techniques, this work contributes to a deeper understanding of the complex interaction between large ship hulls and waves, offering valuable insights for the maritime industry.

We consider the numerical solution of Poisson’s equation on structured grids using geometric multigrid with nonstandard coarse grids and coarse-level operators. We are motivated by the problem of developing high-order accurate numerical solvers for elliptic boundary value problems on complex geometry using overset grids. For flexibility in grid generation, we would like to consider lower-order accurate coarse-level approximations, and coarsening factors other than two. We show that second-order accurate coarse-level approximations are very effective for fourth- or sixth-order accurate fine-level finite-difference discretizations. We study the use of different Galerkin and non-Galerkin coarse-level operators. Using local Fourier analysis (LFA) we choose the smoothing parameter ω and the coarse-level operators to optimize the overall multigrid convergence rate. We show that the results based on LFA for periodic problems also hold for more general boundary conditions provided these are discretized using compatibility conditions . Numerical results for Poisson’s equation on a sample overset grid show that our multigrid solver is many times faster, and uses less memory, than selected Krylov solvers and an algebraic multigrid solver. We also study grid coarsening by a general factor and show that good convergence rates are retained for a range of coarsening factors around two. We ask the question of which coarsening factor leads to the most efficient multigrid algorithm.

Overlapping grids are a powerful tool for representing the complex geometry inherent to many linear and nonlinear electromagnetic scattering problems. However, inter-grid interpolation is known to be potentially destabilizing in the context of wave phenomena. Upwind dissipation has become one of the most powerful tools to cure this instability, but unfortunately some of the original formulations were difficult to implement and expensive, particularly for complex material models. Here we present an optimized upwind scheme whose formulation enables a simple modular implementation so that it can be easily be applied even to complex models, and is computationally efficient in comparison to previous schemes.

Efficient finite-difference schemes for the numerical solution of the time-dependent equations of incompressible linear elasticity on complex geometry are presented. The schemes are second-order accurate in space and time, and solve the equations in displacement-pressure form. The algorithms use a fractional-step approach in which the time-step of the displacements is performed separately from the solution of a Poisson problem to update the pressure. Complex geometry is treated with curvilinear overset grids. Compatibility boundary conditions and an upwind dissipation for wave equations in second-order form are included to ensure stability for overset grids, and for the difficult case of traction (free-surface) boundary conditions. A divergence damping term is added to keep the dilatation small. The stability of the schemes is studied with GKS mode analysis. Exact eigen-mode solutions are obtained for several problems involving rectangular, cylindrical and spherical configurations, and these are valuable as benchmark problems. The accuracy and stability of the approach in two and three space dimensions is illustrated by comparing the numerical solutions with the exact solutions of the benchmark problems.

We describe an algorithm to easily and efficiently incorporate upwinding into finite-difference schemes for solving wave equations in second-order form and apply this scheme to solve problems on complex geometry using overset grids. Upwinding can be added to an existing discretization, such as a centered and dissipation-free scheme, as a modular corrector stage, and takes the form of a special artificial dissipation. This new upwind predictor-corrector scheme significantly improves the run-time performance compared to our original formulation, with typical speedups of factors of ten or more. As with the original upwind formulation, theory and numerical results show that the new algorithm remains robust and stable even for the difficult cases of overset grids with “thin” boundary fitted grids, where nondissipative schemes are generally unstable. Numerical results simulating Maxwell’s equations in second-order form to second- and fourth-order accuracy are used to assess the run-time performance of the new scheme.

We describe a new approach to derive numerical approximations of boundary conditions for high-order accurate finite-difference approximations. The approach, called the Local Compatibility Boundary Condition (LCBC) method, uses boundary conditions and compatibility boundary conditions derived from the governing equations, as well as interior and boundary grid values, to construct a local polynomial, whose degree matches the order of accuracy of the interior scheme, centered at each boundary point. The local polynomial is then used to derive a discrete formula for each ghost point in terms of the data. This approach leads to centered approximations that are generally more accurate and stable than one-sided approximations. Moreover, the stencil approximations are local since they do not couple to neighboring ghost-point values which can occur with traditional compatibility conditions. The local polynomial is derived using continuous operators and derivatives which enables the automatic construction of stencil approximations at different orders of accuracy. The LCBC method is developed here for problems governed by second-order partial differential equations, and it is verified for a wide range of sample problems, both time-dependent and time-independent, in two space dimensions and for schemes up to sixth-order accuracy.

We describe a fourth-order accurate finite-difference time-domain scheme for solving dispersive Maxwell’s equations with nonlinear multi-level carrier kinetics models. The scheme is based on an efficient single-step three time-level modified equation approach for Maxwell’s equations in secondorder form for the electric field coupled to ODEs for the polarization vectors and population densities of the atomic levels. The resulting scheme has a large CFL-one time-step. Curved interfaces between different materials are accurately treated with curvilinear grids and compatibility conditions. A novel hierarchical modified equation approach leads to an explicit scheme that does not require any nonlinear iterations. The hierarchical approach at interfaces leads to local updates at the interface with no coupling in the tangential directions. Complex geometry is treated with overset grids. Numerical stability is maintained using high-order upwind dissipation designed for Maxwell’s equations in second-order form. The scheme is carefully verified for a number of two and three-dimensional problems. The resulting numerical model with generalized dispersion and arbitrary nonlinear multi-level system can be used for many plasmonic applications such as for ab initio time domain modeling of nonlinear engineered materials for nanolasing applications, where nano-patterned plasmonic dispersive arrays are used to enhance otherwise weak nonlinearity in the active media.

Efficient and high-order accurate FDTD schemes are developed for nonlinear dispersive models in active multi-level media with material interfaces in 2D and 3D arbitrary geometry. The considered nonlinear systems are multi-level generalizations of 2-level optical Maxwell-Bloch equations. Arbitrary numbers of atomic levels and macroscopic polarization vectors can be accounted for in the multi-level atomic (MLA) system. The MLA models are suitable for modeling active media with various properties, such as lasing, saturable absorption, reverse saturable absorption, etc. Composite overlapping grids are employed to handle complex geometry with material interfaces, using locally conforming curvilinear grids for curved boundaries and interfaces and Cartesian grids for the rest. The developed high-order schemes allow compact stencils in time integration, and efficient point-wise update of the numerical solutions.

Bianisotropic homogenization for optically thick metasurfaces is demonstrated. The proposed method can homogenize arbitrarily thick metasurface unit cells by dividing the unit cell in propagation direction and treat individual layers. The accuracy of the retrieved parameters is tested with a BA-compatible Maxwell solver. The results indicate a high accuracy for the retrieved BA parameters.

A high-order accurate scheme for solving the time-domain dispersive Maxwell's equations and material interfaces is described. Maxwell's equations are solved in second-order form for the electric field. A generalized dispersive material (GDM) model is used to represent a general class of linear dispersive materials and this model is implemented in the time-domain with the auxiliary differential equation (ADE) approach. The interior updates use our recently developed second-order and fourth-order accurate single-stage three-level space-time finite-difference schemes, and this paper extends these schemes to treat interfaces between different dispersive materials. Composite overlapping grids are used to treat complex geometry, with Cartesian grids generally covering most of the domain and local conforming grids representing curved boundaries and interfaces. Compatibility conditions derived from the interface jump conditions and governing equations are used to derive accurate numerical interface conditions that define values at ghost points. Although some compatibility conditions couple the equations for the ghost points in tangential directions due to mixed-derivatives, it is shown how to decouple the equations to avoid solving a larger system of equations for all ghost points on the interface. The stability of the interface approximations is studied with mode analysis and it is shown that the schemes retain close to a CFL-one time-step restriction. Numerical results are presented in two and three space dimensions to confirm the accuracy and stability of the schemes. The schemes are verified using exact solutions for a planar interface, a disk in two dimensions, and a solid sphere in three dimensions.

We consider the numerical solution of Poisson's equation on structured grids using geometric multigrid with nonstandard coarse grids and coarse level operators. We are motivated by the problem of developing high-order accurate numerical solvers for elliptic boundary value problems on complex geometry using overset grids. Overset grids are typically dominated by large Cartesian background grids and thus fast solvers for Cartesian grids are highly desired. For flexibility in grid generation we would like to consider coarsening factors other than two, and lower-order accurate coarse-level approximations. We show that second-order accurate coarse-level approximations are very effective for fourth- or sixth-order accurate fine-level finite difference discretizations. We study the use of different Galerkin and non-Galerkin coarse-level operators. We use red-black smoothers with a relaxation parameter $\omega$. Using local Fourier analysis we choose $\omega$ and the coarse-level operators to optimize the overall multigrid convergence rate. Motivated by the use of red-black smoothers in one dimension that can result in a direct solver for the standard second-order accurate discretization to Poisson's equation, we show that this direct-solver property can be extended to two dimensions using a rotated grid that results from red-black coarsening. We evaluate the use of red-black coarsening in more general settings. We also study grid coarsening by a general factor and show that good convergence rates are retained for a range of coarsening factors near two. We ask the question of which coarsening factor leads to the most efficient algorithm.