High-order numerical method for the nonlinear Helmholtz equation with material discontinuities in one space dimension

The nonlinear Helmholtz equation (NLH) models the propagation of electromagnetic waves in Kerr media, and describes a range of important phenomena in nonlinear optics and in other areas. In our previous work, we developed a fourth order method for its numerical solution that involved an iterative solver based on freezing the nonlinearity. The method enabled a direct simulation of nonlinear self-focusing in the nonparaxial regime, and a quantitative prediction of backscattering. However, our simulations showed that there is a threshold value for the magnitude of the nonlinearity, above which the iterations diverge.In this study, we numerically solve the one-dimensional NLH using a Newton-type nonlinear solver. Because the Kerr nonlinearity contains absolute values of the field, the NLH has to be recast as a system of two real equations in order to apply Newton’s method. Our numerical simulations show that Newton’s method converges rapidly and, in contradistinction with the iterations based on freezing the nonlinearity, enables computations for very high levels of nonlinearity.In addition, we introduce a novel compact finite-volume fourth order discretization for the NLH with material discontinuities. Our computations corroborate the design fourth order convergence of the method.The one-dimensional results of the current paper create a foundation for the analysis of multidimensional problems in the future.     

On the use of discrete wavelet transform for solving integral equations of acoustic scattering

In this paper, the discrete wavelet transform (DWT) is used to solve the dense system of equations which arises from integral equation of acoustic scattering. The DWT using appropriate wavelet family for acquiring larger sparsification of the system matrix is used to obtain a sparse approximation to the transformed matrix that is used in place of the original matrix in an iterative solver. Alternatively DWT is also used to design sparse preconditioners for an iterative method. Also, DWT-based preconditioners are constructed to accelerate iterative Krylov subspace methods. Convergence rates and number of operations are discussed for each case.

     

A domain decomposition method for solving the three-dimensional time-harmonic Maxwell equations discretized by discontinuous Galerkin methods

We present here a domain decomposition method for solving the three-dimensional time-harmonic Maxwell equations discretized by a discontinuous Galerkin method. In order to allow the treatment of irregularly shaped geometries, the discontinuous Galerkin method is formulated on unstructured tetrahedral meshes. The domain decomposition strategy takes the form of a Schwarz-type algorithm where a continuity condition on the incoming characteristic variables is imposed at the interfaces between neighboring subdomains. A multifrontal sparse direct solver is used at the subdomain level. The resulting domain decomposition strategy can be viewed as a hybrid iterative/direct solution method for the large, sparse and complex coefficients algebraic system resulting from the discretization of the time-harmonic Maxwell equations by a discontinuous Galerkin method.     

Three-dimensional electromagnetic nonlinear inversion in layered media by a hybrid diagonal tensor approximation: Stabilized biconjugate gradient fast Fourier transform method

This paper presents an efficient three-dimensional nonlinear electromagnetic inversion method in a multilayered medium for radar applications where the object size is comparable to the wavelength. In the first step of this two-step inversion algorithm, the diagonal tensor approximation is used in the Born iterative method. The solution of this approximate inversion is used as an initial guess for the second step in which further inversion is carried out using a distorted Born iterative method. Since the aim of the second step is to improve the accuracy of the inversion, a full-wave solver, the stabilized biconjugate-gradient fast Fourier transform algorithm, is used for forward modelling. The conjugate-gradient method is applied at each inversion iteration to minimize the functional cost. The usage of an iterative solver based on the FFT algorithm and the developed recursive matrix method combined with an interpolation technique to evaluate the layered medium Green's functions rapidly, makes this method highly efficient. An inversion problem with 32 768 complex unknowns can be solved with 1% relative error by using a simple personal computer. Several numerical experiments for arbitrarily located source and receiver arrays are presented to show the high efficiency and accuracy of the proposed method.     

Three-dimensional electromagnetic nonlinear inversion in layered media by a hybrid diagonal tensor approximation: Stabilized biconjugate gradient fast Fourier transform method

This paper presents an efficient three-dimensional nonlinear electromagnetic inversion method in a multilayered medium for radar applications where the object size is comparable to the wavelength. In the first step of this two-step inversion algorithm, the diagonal tensor approximation is used in the Born iterative method. The solution of this approximate inversion is used as an initial guess for the second step in which further inversion is carried out using a distorted Born iterative method. Since the aim of the second step is to improve the accuracy of the inversion, a full-wave solver, the stabilized biconjugate-gradient fast Fourier transform algorithm, is used for forward modelling. The conjugate-gradient method is applied at each inversion iteration to minimize the functional cost. The usage of an iterative solver based on the FFT algorithm and the developed recursive matrix method combined with an interpolation technique to evaluate the layered medium Green's functions rapidly, makes this method highly efficient. An inversion problem with 32 768 complex unknowns can be solved with 1% relative error by using a simple personal computer. Several numerical experiments for arbitrarily located source and receiver arrays are presented to show the high efficiency and accuracy of the proposed method.     

一种求解三维热辐射输运方程的整体预处理迭代方法及并行计算

为三维灰体热辐射输运方程的隐式离散纵标方法发展一个整体预处理迭代方法并研制并行程序.该方法采用组装线性代数方程组策略,同时求出所有离散方向上的辐射强度.借助预处理的Krylov子空间迭代法,避免复杂网格上扫描算法可能遇到的死锁问题,能够提高健壮性和计算效率.空间离散上采用一阶迎风有限体积格式.数值实验测试变形六面体网格...     

Preconditioning based on Calderon’s formulae for periodic fast multipole methods for Helmholtz’ equation

Solution of periodic boundary value problems is of interest in various branches of science and engineering such as optics, electromagnetics and mechanics. In our previous studies we have developed a periodic fast multipole method (FMM) as a fast solver of wave problems in periodic domains. It has been found, however, that the convergence of the iterative solvers for linear equations slows down when the solutions show anomalies related to the periodicity of the problems. In this paper, we propose preconditioning schemes based on Calderon’s formulae to accelerate convergence of iterative solvers in the periodic FMM for Helmholtz’ equations. The proposed preconditioners can be implemented more easily than conventional ones. We present several numerical examples to test the performance of the proposed preconditioners. We show that the effectiveness of these preconditioners is definite even near anomalies.     

Efficient algebraic multigrid for migration–diffusion–convection–reaction systems arising in electrochemical simulations

The article discusses components and performance of an algebraic multigrid (AMG) preconditioner for the fully coupled multi-ion transport and reaction model (MITReM) with nonlinear boundary conditions, important for electrochemical modeling. The governing partial differential equations (PDEs) are discretized in space by a combined finite element and residual distribution method. Solution of the discrete system is obtained by means of a Newton-based nonlinear solver, and an AMG-preconditioned BICGSTAB Krylov linear solver. The presented AMG preconditioner is based on so-called point-based classical AMG. The linear solver is compared to a standard direct and several one-level iterative solvers for a range of geometries and chemical systems with scientific and industrial relevance. The results indicate that point-based AMG methods, carefully designed, are an attractive alternative to more commonly employed numerical methods for the simulation of complex electrochemical processes.     

A Solver for Helmholtz System Generated by the Discretization of Wave Shape Functions

An interesting discretization method for Helmholtz equations was introduced in B. Després [1]. This method is based on the ultra weak variational formulation (UWVF) and the wave shape functions, which are exact solutions of the governing Helmholtz equation. In this paper we are concerned with fast solver for the system generated by the method in [1]. We propose a new preconditioner for such system, which can be viewed as a combination between a coarse solver and the block diagonal preconditioner introduced in [13]. In our numerical experiments, this preconditioner is applied to solve both two-dimensional and three-dimensional Helmholtz equations, and the numerical results illustrate that the new preconditioner is much more efficient than the original block diagonal preconditioner.     

Discontinuous Galerkin method coupled with an integral representation for solving the three-dimensional time-harmonic Maxwell equations

In this paper, we present a mathematical and numerical studies of the three-dimensional time-harmonic Maxwell equations. The problem is solved by a discontinuous Galerkin DG method coupled with an integral representation. This study was completed by some numerical tests to justify the effectiveness of the proposed approach. The numerical simulation was done by an iterative solver implemented in FORTRAN.     

Three-dimensional elliptic solvers for interface problems and applications

Second-order accurate elliptic solvers using Cartesian grids are presented for three-dimensional interface problems in which the coefficients, the source term, the solution and its normal flux may be discontinuous across an interface. One of our methods is designed for general interface problems with variable but discontinuous coefficient. The scheme preserves the discrete maximum principle using constrained optimization techniques. An algebraic multigrid solver is applied to solve the discrete system. The second method is designed for interface problems with piecewise constant coefficient. The method is based on the fast immersed interface method and a fast 3D Poisson solver. The second method has been modified to solve Helmholtz/Poisson equations on irregular domains. An application of our method to an inverse interface problem of shape identification is also presented. In this application, the level set method is applied to find the unknown surface iteratively.     

On the HLLC Riemann solver for interface interaction in compressible multi-fluid flow

In this work, the HLLC Riemann solver, which is much more robust, simpler and faster than iterative Riemann solvers, is extended to obtain interface conditions in sharp-interface methods for compressible multi-fluid flows. For interactions with general equations of state and material interfaces, a new generalized Roe average is proposed. For single-phase interactions, this new Roe average does not introduce artificial states and satisfies the U-property exactly. For interactions at material interfaces, the U-property is satisfied by introducing ghost states for the internal energy. A number of numerical tests suggest that the proposed Riemann solver is suitable for general equations of state and has an accuracy comparable to iterative Riemann solvers, while being significantly more robust and efficient.     

Event-based parareal: A data-flow based implementation of parareal

Parareal is an iterative algorithm that, in effect, achieves temporal decomposition for a time-dependent system of differential or partial differential equations. A solution is obtained in a shorter wall-clock time, but at the expense of increased compute cycles. The algorithm combines a fine solver that solves the system to acceptable accuracy with an approximate coarse solver. The critical task for the successful implementation of parareal on any system is the development of a coarse solver that leads to convergence in a small number of iterations compared to the number of time slices in the full time interval, and is, at the same time, much faster than the fine solver. Very fast coarse solvers may not lead to sufficiently rapid convergence, and slow coarse solvers may not lead to significant gains even if the number of iterations to convergence is satisfactory. We find that the difficulty of meeting these conflicting demands can be substantially eased by using a data-driven, event-based implementation of parareal. As a result, tasks for one iteration do not wait for the previous iteration to complete, but are started when the needed data are available. For given convergence properties, the event-based approach relaxes the speed requirements on the coarse solver by a factor of ～K, where K is the number of iterations required for a converged solution. This may, for many problems, lead to an efficient parareal implementation that would otherwise not be possible or would require substantial coarse solver development. In addition, the framework used for this implementation executes a task when the data dependencies are satisfied and computational resources are available. This leads to improved computational efficiency over previous approaches that pipeline or schedule groups of tasks to a particular processor or group of processors.     

二维三温热传导方程求解中的非线性迭代初值选取

由于二维三温热传导方程具有很强的非线性特性,因此采用全隐格式对该方程离散后,所得非线性代数方程组的求解将变得非常困难.针对二维三温热传导方程离散所得非线性代数方程组的迭代求解,提出了一种有效的选取初值的方法.对两种不同性质的介质进行数值实验,结果表明,所设计的初值选取方法不仅大大提高了计算效率,而且能够降低非线性解法器对时间步长的影响.     

A well-conditioned augmented system for solving Navier–Stokes equations in irregular domains

An augmented method based on a Cartesian grid is proposed for the incompressible Navier–Stokes equations in irregular domains. The irregular domain is embedded into a rectangular one so that a fast Poisson solver can be utilized in the projection method. Unlike several methods suggested in the literature that set the force strengths as unknowns, which often results in an ill-conditioned linear system, we set the jump in the normal derivative of the velocity as the augmented variable. The new approach improves the condition number of the system for the augmented variable significantly. Using the immersed interface method, we are able to achieve second order accuracy for the velocity. Numerical results and comparisons to benchmark tests are given to validate the new method. A lid-driven cavity flow with multiple obstacles and different geometries are also presented.     

Fast multi-orbital equation of motion impurity solver for dynamical mean field theory

We propose a fast multi-orbital impurity solver for dynamical mean field theory (DMFT). Our DMFT solver is based on the equations of motion (EOMs) for local Green's functions and is constructed by generalizing from the single-orbital case to the multi-orbital case with the inclusion of the inter-orbital hybridizations and applying a mean field approximation to the inter-orbital Coulomb interactions. The two-orbital Hubbard model is studied using this impurity solver within a large range of parameters. The Mott metal-insulator transition and the quasiparticle peak are well described. A comparison of the EOM method with the quantum Monte Carlo method is made for the two-orbital Hubbard model and good agreement is obtained. The developed method hence holds promise as a fast DMFT impurity solver in studies of strongly correlated systems.     

A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources

In this paper, we propose a general time-discrete framework to design asymptotic-preserving schemes for initial value problem of the Boltzmann kinetic and related equations. Numerically solving these equations are challenging due to the nonlinear stiff collision (source) terms induced by small mean free or relaxation time. We propose to penalize the nonlinear collision term by a BGK-type relaxation term, which can be solved explicitly even if discretized implicitly in time. Moreover, the BGK-type relaxation operator helps to drive the density distribution toward the local Maxwellian, thus naturally imposes an asymptotic-preserving scheme in the Euler limit. The scheme so designed does not need any nonlinear iterative solver or the use of Wild Sum. It is uniformly stable in terms of the (possibly small) Knudsen number, and can capture the macroscopic fluid dynamic (Euler) limit even if the small scale determined by the Knudsen number is not numerically resolved. It is also consistent to the compressible Navier–Stokes equations if the viscosity and heat conductivity are numerically resolved. The method is applicable to many other related problems, such as hyperbolic systems with stiff relaxation, and high order parabolic equations.     

The Riemann Problem for the one-dimensional,free-surface Shallow Water Equations with a bed step: Theoretical analysis and numerical simulations

In this paper, the solution of the Riemann Problem for the one-dimensional, free-surface Shallow Water Equations over a bed step is analyzed both from a theoretical and a numerical point of view. Particular attention has been paid to the wave that is generated at the location of the bed discontinuity. Starting from the classical Shallow Water Equations, considering the bed level as an additional variable, and adding to the system an equation imposing its time invariance, we show that this wave is a contact wave, across which one of the Riemann invariants, namely the energy, is not constant. This is due to the fact that the relevant problem is nonconservative. We demonstrate that, in this type of system, Riemann Invariants do not generally hold in contact waves. Furthermore, we show that in this case the equations that link the flow variables across the contact wave are the Generalized Rankine–Hugoniot relations and we obtain these for the specific problem. From the numerical point of view, we present an accurate and efficient solver for the step Riemann Problem to be used in a finite-volume Godunov-type framework. Through a two-step predictor–corrector procedure, the solver is able to provide solutions with any desired accuracy. The predictor step uses a well-balanced Generalized Roe solver while the corrector step solves the exact nonlinear system of equations that consitutes the problem by means of an iterative procedure that starts from the predictor solution. In order to show the effectiveness and the accuracy of the proposed approach, we consider several step Riemann Problems and compare the exact solutions with the numerical results obtained by using a standard Roe approach far from the step and the novel two-step algorithm for the fluxes over the step, achieving good results.     

A fast parallel Poisson solver on irregular domains applied to beam dynamics simulations

We discuss the scalable parallel solution of the Poisson equation within a Particle-In-Cell (PIC) code for the simulation of electron beams in particle accelerators of irregular shape. The problem is discretized by Finite Differences. Depending on the treatment of the Dirichlet boundary the resulting system of equations is symmetric or ‘mildly’ nonsymmetric positive definite. In all cases, the system is solved by the preconditioned conjugate gradient algorithm with smoothed aggregation (SA) based algebraic multigrid (AMG) preconditioning. We investigate variants of the implementation of SA-AMG that lead to considerable improvements in the execution times. We demonstrate good scalability of the solver on distributed memory parallel processor with up to 2048 processors. We also compare our iterative solver with an FFT-based solver that is more commonly used for applications in beam dynamics.     

High-order quadratures for the solution of scattering problems in two dimensions

We construct an iterative algorithm for the solution of forward scattering problems in two dimensions. The scheme is based on the combination of high-order quadrature formulae, fast application of integral operators in Lippmann–Schwinger equations, and the stabilized bi-conjugate gradient method (BI-CGSTAB). While the FFT-based fast application of integral operators and the BI-CGSTAB for the solution of linear systems are fairly standard, a large part of this paper is devoted to constructing a class of high-order quadrature formulae applicable to a wide range of singular functions in two and three dimensions; these are used to obtain rapidly convergent discretizations of Lippmann–Schwinger equations. The performance of the algorithm is illustrated with several numerical examples.