On Integral Equation and Least Squares Methods for Scattering by Diffraction Gratings

In this paper we consider the scattering of a plane acoustic or electromagnetic wave by a one-dimensional, periodic rough surface. We restrict the discussion to the case when the boundary is sound soft in the acoustic case, perfectly reflecting with TE polarization in the EM case, so that the total field vanishes on the boundary. We propose a uniquely solvable first kind integral equation formulation of the problem, which amounts to a requirement that the normal derivative of the Green's representation formula for the total field vanish on a horizontal line below the scattering surface. We then discuss the numerical solution by Galerkin's method of this (ill-posed) integral equation. We point out that, with two particular choices of the trial and test spaces, we recover the so-called SC (spectral-coordinate) and SS (spectral-spectral) numerical schemes of DeSanto et al., Waves Random Media, 8, 315-414, 1998. We next propose a new Galerkin scheme, a modification of the SS method that we term the SS^∗ method, which is an instance of the well-known dual least squares Galerkin method. We show that the SS^∗ method is always well-defined and is optimally convergent as the size of the approximation space increases. Moreover, we make a connection with the classical least squares method, in which the coefficients in the Rayleigh expansion of the solution are determined by enforcing the boundary condition in a least squares sense, pointing out that the linear system to be solved in the SS^∗ method is identical to that in the least squares method. Using this connection we show that (reflecting the ill-posed nature of the integral equation solved) the condition number of the linear system in the SS^∗ and least squares methods approaches infinity as the approximation space increases in size. We also provide theoretical error bounds on the condition number and on the errors induced in the numerical solution computed as a result of ill-conditioning. Numerical results confirm the convergence of the SS^∗ method and illustrate the ill-conditioning that arises.     

Integral Equation Method for a Non-Selfadjoint Steklov Eigenvalue Problem

We propose a numerical method for a non-selfadjoint Steklov eigenvalueproblem of the Helmholtz equation. The problem is formulated using boundary integrals. The Nyström method is employed to discretize the integral operators, whichleads to a non-Hermitian generalized matrix eigenvalue problems. The spectral indicator method (SIM) is then applied to calculate the (complex) eigenvalues. The convergence is proved using the spectral approximation theory for (non-selfadjoint) compactoperators. Numerical examples are presented for validation.     

A Fourth-Order Kernel-Free Boundary Integral Method for Interface Problems

This paper presents a fourth-order Cartesian grid based boundary integral method (BIM) for heterogeneous interface problems in two and three dimensional space, where the problem interfaces are irregular and can be explicitly given by parametric curves or implicitly defined by level set functions. The method reformulates the governing equation with interface conditions into boundary integral equations (BIEs) and reinterprets the involved integrals as solutions to some simple interface problems in an extended regular region. Solution of the simple equivalent interface problems for integral evaluation relies on a fourth-order finite difference method with an FFT-based fast elliptic solver. The structure of the coefficient matrix is preserved even with the existence of the interface. In the whole calculation process, analytical expressions of Green’s functions are never determined, formulated or computed. This is the novelty of the proposed kernel-free boundary integral (KFBI) method. Numerical experiments in both two and three dimensions are shown to demonstrate the algorithm efficiency and solution accuracy even for problems with a large diffusion coefficient ratio.     

A Highly Scalable Boundary Integral Equation and Walk-on-Spheres (BIE-WOS) Method for the Laplace Equation with Dirichlet Data

In this paper, we study a highly scalable communication-free parallel domain boundary decomposition algorithm for the Laplace equation based on a hybrid method combining boundary integral equations and walk-on-spheres (BIE-WOS)method, which provides a numerical approximation of the Dirichlet-to-Neumann(DtN) mapping for the Laplace equation. The BIE-WOS is a local method on theboundary of the domain and does not require a structured mesh, and only needs acovering of the domain boundary by patches and a local mesh for each patch for a local BIE. A new version of the BIE-WOS method with second kind integral equations isintroduced for better error controls. The effect of errors from the Feynman-Kac formulabased path integral WOS method on the overall accuracy of the BIE-WOS method isanalyzed for the BIEs, especially in the calculation of the right hand sides of the BIEs.For the special case of flat patches, it is shown that the second kind integral equationof BIE-WOS method can be simplified where the local BIE solutions can be given inclosed forms. A key advantage of the parallel BIE-WOS method is the absence of communications during the computation of the DtN mapping on individual patches ofthe boundary, resulting in a complete independent computation using a large numberof cluster nodes. In addition, the BIE-WOS has an intrinsic capability of fault tolerance for exascale computations. The nearly linear scalability of the parallel BIE-WOSmethod on a large-scale cluster with 6400 CPU cores is verified for computing the DtNmapping of exterior Laplace problems with Dirichlet data for several domains.     

A Continuous Finite Element Method with Homotopy Vanishing Viscosity for Solving the Static Eikonal Equation

We develop a second-order continuous finite element method for solving the static Eikonal equation. It is based on the vanishing viscosity approach with a homotopy method for solving the discretized nonlinear system. More specifically, the homotopy method is utilized to decrease the viscosity coefficient gradually, while Newton’s method is applied to compute the solution for each viscosity coefficient. Newton’s method alone converges for just big enough viscosity coefficients on very coarse grids and for simple 1D examples, but the proposed method is much more robust and guarantees the convergence of the nonlinear solver for all viscosity coefficients and for all examples over all grids. Numerical experiments from 1D to 3D are presented to confirm the second-order convergence and the effectiveness of the proposed method on both structured or unstructured meshes.     

A Sparse Grid Discrete Ordinate Discontinuous Galerkin Method for the Radiative Transfer Equation

The radiative transfer equation is a fundamental equation in transport theory and applications, which is a 5-dimensional PDE in the stationary one-velocity case, leading to great difficulties in numerical simulation. To tackle this bottleneck, we first use the discrete ordinate technique to discretize the scattering term, an integral with respect to the angular variables, resulting in a semi-discrete hyperbolic system. Then, we make the spatial discretization by means of the discontinuous Galerkin (DG) method combined with the sparse grid method. The final linear system is solved by the block Gauss-Seidal iteration method. The computational complexity and error analysis are developed in detail, which show the new method is more efficient than the original discrete ordinate DG method. A series of numerical results are performed to validate the convergence behavior and effectiveness of the proposed method.     

Numerical Method for the Deterministic Kardar-Parisi-Zhang Equation in Unbounded Domains

We propose an artificial boundary method for solving the deterministic Kardar-Parisi-Zhang equation in one-, two- and three- dimensional unbounded domains. The exact artificial boundary conditions are obtained on the artificial boundaries. Then the original problems are reduced to equivalent problems in bounded domains. A finite difference method is applied to solve the reduced problems, and some numerical examples are provided to show the effectiveness of the method.     

An Adaptive Finite Element Method with a Modified Perfectly Matched Layer Formulation for Diffraction Gratings

For numerical simulation of one-dimensional diffraction gratings both in TE and TM polarization, an enhanced adaptive finite element method is proposed in this paper. A modified perfectly matched layer (PML) formulation is proposed for the truncation of the unbounded domain, which results in a homogeneous Dirichlet boundary condition and the corresponding error estimate is greatly simplified. The a posteriori error estimates for the adaptive finite element method are provided. Moreover, a lower bound is obtained to demonstrate that the error estimates obtained are sharp.     

Fast Multipole Accelerated Boundary Integral Equation Method for Evaluating the Stress Field Associated with Dislocations in a Finite Medium

In this paper, we develop an efficient numerical method based on the boundary integral equation formulation and new version of fast multipole method to solve the boundary value problem for the stress field associated with dislocations in a finite medium. Numerical examples are presented to examine the influence from material boundaries on dislocations.     

On the Five-Moment Maximum Entropy System of One-Dimensional Boltzmann Equation

The maximum entropy moment system extends the Euler equation to nonequilibrium gas flows by considering higher order moments such as the heat flux. This paper presents a systematic study of the maximum entropy moment system of Boltzmann equation. We consider a hypothetical one-dimensional gas and study a five-moment model. A numerical algorithm for solving the optimization problem is developed to produce the maximum entropy distribution function from known moments, and the asymptotic behaviour of the system around the singular region known as the Junk’s line, as well as that near the boundary of the realizability domain is analyzed. Furthermore, we study the properties of the system numerically, including the behaviour of the system around the Maxwellian and within the interior of the realizability domain, and properties of its characteristic fields. Our studies show the higher order entropy-based moment models to differ significantly from the Euler equations. Much of this difference comes from the singularity near the Junk’s line, which would be removed if a truncation of the velocity domain is employed.     

A High-Accurate Fast Poisson Solver Based on Harmonic Surface Mapping Algorithm

Poisson's equations in a cuboid are frequently solved in many scientific andengineering applications such as electric structure calculations, molecular dynamicssimulations and computational astrophysics. In this paper, a fast and highly accuratealgorithm is presented for the solution of the Poisson's equation in a cuboidal domainwith boundary conditions of mixed type. This so-called harmonic surface mappingalgorithm is a meshless algorithm which can achieve a desired order of accuracy byevaluating a body convolution of the source and the free-space Green's function withina sphere containing the cuboid, and another surface integration over the spherical surface. Numerical quadratures are introduced to approximate the integrals, resultingin the solution represented by a summation of point sources in free space, which canbe accelerated by means of the fast multipole algorithm. The complexity of the algorithm is linear to the number of quadrature points, and the convergence rate can bearbitrarily high even when the source term is a piecewise continuous function.     

High Order Finite Difference Hermite WENO Fixed-Point Fast Sweeping Method for Static Hamilton-Jacobi Equations

In this paper, we combine the nonlinear HWENO reconstruction in [43] andthe fixed-point iteration with Gauss-Seidel fast sweeping strategy, to solve the staticHamilton-Jacobi equations in a novel HWENO framework recently developed in [22].The proposed HWENO frameworks enjoys several advantages. First, compared withthe traditional HWENO framework, the proposed methods do not need to introduceadditional auxiliary equations to update the derivatives of the unknown function $phi$.They are now computed from the current value of $phi$ and the previous spatial derivatives of $phi$. This approach saves the computational storage and CPU time, which greatlyimproves the computational efficiency of the traditional HWENO scheme. In addition,compared with the traditional WENO method, reconstruction stencil of the HWENOmethods becomes more compact, their boundary treatment is simpler, and the numerical errors are smaller on the same mesh. Second, the fixed-point fast sweeping methodis used to update the numerical approximation. It is an explicit method and doesnot involve the inverse operation of nonlinear Hamiltonian, therefore any Hamilton-Jacobi equations with complex Hamiltonian can be solved easily. It also resolves someknown issues, including that the iterative number is very sensitive to the parameter $ε$ used in the nonlinear weights, as observed in previous studies. Finally, to furtherreduce the computational cost, a hybrid strategy is also presented. Extensive numerical experiments are performed on two-dimensional problems, which demonstrate thegood performance of the proposed fixed-point fast sweeping HWENO methods.     

An Adaptive Moving Mesh Discontinuous Galerkin Method for the Radiative Transfer Equation

The radiative transfer equation models the interaction of radiation with scattering and absorbing media and has important applications in various fields in scienceand engineering. It is an integro-differential equation involving time, frequency, spaceand angular variables and contains an integral term in angular directions while being hyperbolic in space. The challenges for its numerical solution include the needsto handle with its high dimensionality, the presence of the integral term, and the development of discontinuities and sharp layers in its solution along spatial directions.Its numerical solution is studied in this paper using an adaptive moving mesh discontinuous Galerkin method for spatial discretization together with the discrete ordinatemethod for angular discretization. The former employs a dynamic mesh adaptationstrategy based on moving mesh partial differential equations to improve computational accuracy and efficiency. Its mesh adaptation ability, accuracy, and efficiency aredemonstrated in a selection of one- and two-dimensional numerical examples.     

A Diagonal Sweeping Domain Decomposition Method with Source Transfer for the Helmholtz Equation

In this paper, we propose and test a novel diagonal sweeping domain decomposition method (DDM) with source transfer for solving the high-frequency Helmholtz equation in$\mathbb{R}^n$. In the method the computational domain is partitioned into overlapping checkerboard subdomains for source transfer with the perfectly matched layer (PML) technique, then a set of diagonal sweeps over the subdomains are specially designed to solve the system efficiently. The method improves the additive overlapping DDM [43] and the L-sweeps method [50] by employing a more efficient subdomain solving order. We show that the method achieves the exact solution of the global PML problem with $2^n$ sweeps in the constant medium case. Although the sweeping usually implies sequential subdomain solves, the number of sequential steps required for each sweep in the method is only proportional to the $n$-th root of the number of subdomains when the domain decomposition is quasi-uniform with respect to all directions, thus it is very suitable for parallel computing of the Helmholtz problem with multiple right-hand sides through the pipeline processing. Extensive numerical experiments in two and three dimensions are presented to demonstrate the effectiveness and efficiency of the proposed method.     

Random Batch Particle Methods for the Homogeneous Landau Equation

We consider in this paper random batch particle methods for efficiently solving the homogeneous Landau equation in plasma physics. The methods are stochastic variations of the particle methods proposed by Carrillo et al. [J. Comput. Phys.: X 7: 100066, 2020] using the random batch strategy. The collisions only take place inside the small but randomly selected batches so that the computational cost is reduced to $\mathcal{O}(N)$ per time step. Meanwhile, our methods can preserve the conservation of mass, momentum, energy and the decay of entropy. Several numerical examples are performed to validate our methods.     

The New Operator Marching Method on Calculating the Electromagnetic Scattered Fields from the Periodic Structures

In this work, the operator marching methods based on the modifiedNeumann-to-Dirichlet (NtD) map are developed to calculate the scattered fields fromthe periodic structures, especially diffraction gratings with metallic material. For thegrating structures coated with the metallic material, the modified NtD map operatormarching scheme is numerically stable for calculating the electromagnetic scatteredfields. The modified NtD map operators are constructed by the integral equation(IE) method in each homogeneous medium of the layered medium structures, andavoid the blow up of the condition numbers from the challenging metallic materialwith the complex refractive index. For the calculations of the scattered fields fromthe diffraction grating structures, the modified NtD map operator marching methodtakes advantage of the periodic structure features, and avoids the calculations of thecomplicated quasi-periodic Green's function. Numerical results demonstrate that theproposed method achieves high accuracy and low computational workload for thescattered fields from both the dielectric and metallic diffraction grating structures.     

Central Discontinuous Galerkin Methods for the Generalized Korteweg-de Vries Equation

In this paper, we develop central discontinuous Galerkin (CDG) finite element methods for solving the generalized Korteweg-de Vries (KdV) equations in onedimension. Unlike traditional discontinuous Galerkin (DG) method, the CDG methods evolve two approximate solutions defined on overlapping cells and thus do notneed numerical fluxes on the cell interfaces. Several CDG schemes are constructed, including the dissipative and non-dissipative versions. L²error estimates are establishedfor the linear and nonlinear equation using several projections for different parameterchoices. Although we can not provide optimal a priori error estimate, numerical examples show that our scheme attains the optimal (k+1)-th order of accuracy when usingpiecewise k-th degree polynomials for many cases.     

Fast Spectral Collocation Method for Surface Integral Equations of Potential Problems in a Spheroid

This paper proposes a new technique to speed up the computation of the matrix of spectral collocation discretizations of surface single and double layer operators over a spheroid. The layer densities are approximated by a spectral expansion of spherical harmonics and the spectral collocation method is then used to solve surface integral equations of potential problems in a spheroid. With the proposed technique, the computation cost of collocation matrix entries is reduced from O(M²N⁴) to O(MN⁴), where N²is the number of spherical harmonics (i.e., size of the matrix) and M is the number of one-dimensional integration quadrature points. Numerical results demonstrate the spectral accuracy of the method.     

A Fast Accurate Boundary Integral Method for Potentials on Closely Packed Cells

Boundary integral methods are naturally suited for the computation of harmonic functions on a region having inclusions or cells with different material properties. However, accuracy deteriorates when the cell boundaries are close to each other. We present a boundary integral method in two dimensions which is specially designed to maintain second order accuracy even if boundaries are arbitrarily close. The method uses a regularization of the integral kernel which admits analytically determined corrections to maintain accuracy. For boundaries with many components we use the fast multipole method for efficient summation. We compute electric potentials on a domain with cells whose conductivity differs from that of the surrounding medium. We first solve an integral equation for a source term on the cell interfaces and then find values of the potential near the interfaces via integrals. Finally we use a Poisson solver to extend the potential to a regular grid covering the entire region. A number of examples are presented. We demonstrate that increased refinement is not needed to maintain accuracy as interfaces become very close.     

Numerical Integrators for Dispersion-Managed KdV Equation

In this paper, we consider the numerics of the dispersion-managed Korteweg-de Vries (DM-KdV) equation for describing wave propagations in inhomogeneous media. The DM-KdV equation contains a variable dispersion map with discontinuity,which makes the solution non-smooth in time. We formally analyze the convergenceorder reduction problems of some popular numerical methods including finite difference and time-splitting for solving the DM-KdV equation, where a necessary constraint on the time step has been identified. Then, two exponential-type dispersion-map integrators up to second order accuracy are derived, which are efficiently incorporated with the Fourier pseudospectral discretization in space, and they can convergeregardless of the discontinuity and the step size. Numerical comparisons show the advantage of the proposed methods with the application to solitary wave dynamics andextension to the fast & strong dispersion-management regime.