This article deals with the efficient and accurate computation of the electrostatic forces between charged, spherical dielectric particles undergoing mutual polarisation. We use the spectral Galerkin boundary integral equation framework developed by Lindgren et al. (J. Comput. Phys. 371 (2018): 712-731) and subsequently analysed in two earlier contributions of the authors to propose a linear scaling in cost algorithm for the computation of the approximate forces. We establish exponential convergence of the method and derive error estimates for the approximate forces that do not explicitly depend on the number of dielectric particles $N$. Consequently, the proposed method requires only $\mathcal{O}(N) $ operations to compute the electrostatic forces acting on $N$ dielectric particles up to any given and fixed relative error.

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.

In this paper, we propose accurate numerical boundary conditions for atomic simulations of twin boundary. The heterogeneity of the lattice structure induces physical reflection across the twin boundary. When numerical boundary and the twin boundary coincide, the goal is to reproduce the correct amount of physical reflection. In particular, we consider waves periodic in the direction parallel to the twin boundary and reduce the problem into a complex-valued chain motion. Using Laplace transform, we design time history kernel (THK) treatment. We further design matching boundary conditions (MBC) by reproducing physical reflection at long wave limit and a specific wave number. Reflection analysis and numerical tests demonstrate the effectiveness of the proposed THK and MBC treatments.

In this paper, we propose a class of numerical methods based on discrete-velocity vector-BGK models for the incompressible Navier-Stokes equations. By analyzing a splitting method with Maxwell iteration, we show that the usual lattice Boltzmann discretization of the vector-BGK models provides a good numerical scheme. Moreover, we establish the stability of the numerical scheme. The stability and second-order accuracy of the scheme are validated through numerical simulations of the two-dimensional Taylor-Green vortex flows. Further numerical tests are conducted to exhibit some potential advantages of the vector-BGK models, which can be regarded as competitive alternatives of the scalar-BGK models.

Long time simulations are needed in the numerical study of the Zeldovich-Neumann-Döring model, in which the quality resolving the dynamics of the detonation front is crucial. The numerical error introduced from the inappropriate outflow boundary condition and the mesh resolution are two main factors qualitatively affecting the dynamics of the detonation front. In this paper we improve the numerical framework in [15] by introducing the Strang splitting method and a new $h$-adaptive method with a feature based $a$ $posteriori$ error estimator. Then a cheap numerical approach is proposed to sharply estimate a time period, in which the unphysical influence on the detonation front can be avoided effectively. The sufficiently dense mesh resolution can be guaranteed around the detonation front and in the reaction zone by the proposed $h$-adaptive method. The numerical results show that the proposed method is sufficiently robust even for long time calculations, and the quality dynamics of the detonation can be obtained by the proposed numerical approach.

We study the cosmological Burgers model, as we call it, which is a nonlinear hyperbolic balance law (in one and two spatial variables) posed on an expanding or contracting background. We design a finite volume scheme that is fourth-order in time and second-order in space, and allows us to compute weak solutions containing shock waves. Our main contribution is the study of the asymptotic structure of the solutions as the time variable approaches infinity (in the expanding case) or zero (in the contracting case). We discover that a saddle competition is taking place which involves, on one hand, the geometrical effects of expanding or contracting nature and, on the other hand, the nonlinear interactions between shock waves.

In the immersed interface method, an object in a flow is formulated as a singular force, and jump conditions caused by the singular force are incorporated into numerical schemes to compute the flow. Previous development of the method considered only smooth objects. We here extend the method to handle non-smooth rigid objects with sharp corners in 2D incompressible viscous flows. We represent the boundary of an object as a polygonal curve moving through a fixed Cartesian grid. We compute necessary jump conditions to achieve boundary condition capturing on the object. We incorporate the jump conditions into finite difference schemes to solve the flow on the Cartesian grid. The accuracy, efficiency and robustness of our method are tested using canonical flow problems. The results demonstrate that the method has second-order accuracy for the velocity and first-order accuracy for the pressure in the infinity norm, and is extremely efficient and robust to simulate flows around non-smooth complex objects.

In the finite difference WENO (weighted essentially non-oscillatory) method, the final scheme on the whole stencil was constructed by linear combinations of highest order accurate schemes on sub-stencils, all of which share the same total count of grid points. The linear combination method which the original WENO applied was generalized to arbitrary positive-integer-order derivative on an arbitrary (uniform or non-uniform) mesh, still applying finite difference method. The possibility of expressing the final scheme on the whole stencil as a linear combination of highest order accurate schemes on WENO-like sub-stencils was investigated. The main results include: (a) the highest order of accuracy a finite difference scheme can achieve and (b) a sufficient and necessary condition that the linear combination exists. This is a sufficient and necessary condition for all finite difference schemes in a set (rather than a specific finite difference scheme) to have WENO-like linear combinations. After the proofs of the results, some remarks on the WENO schemes and TENO (targeted essentially non-oscillatory) schemes were given.

Simulating charged many-body systems has been a computational demanding task due to the long-range nature of electrostatic interaction. For the multi-scale model of electrolytes which combines the strengths of atomistic/continuum electrolyte representations, a harmonic surface mapping algorithm is developed for fast and accurate evaluation of the electrostatic reaction potentials. Our method reformulates the reaction potential into a sum of image charges for the near-field, and a charge density on an auxiliary spherical surface for the far-field, which can be further discretized into point charges. Fast multipole method is used to accelerate the pairwise Coulomb summation. The accuracy and efficiency of our algorithm, as well as the choice of relevant numerical parameters are demonstrated in detail. As a concrete example, for charges close to the dielectric interface, our method can improve the accuracy by two orders of magnitudes compared to the Kirkwood series expansion method.

Via adopting the matrix exponential (ME) approach, efficient and tight complex-frequency-shifted perfectly matched layer (CFS-PML) formulations using the matched Z-transform (MZT) method are well organized and proposed to truncate the monolayer black phosphorous (BP) metasurface problems in the infrared regime. Due to introducing ME approach, the proposed MZT-ME-PML develops into the compact first-order differential matrix form used to the circumvent convolution manipulations, the variable replacements, or the formulation rearrangements during the mathematical derivations. Furthermore, the MZT-ME-PML can not only flexibly truncate arbitrary materials due to using the DB constitutive relation, but efficiently attenuate evanescent waves and reduce late-time reflections with the CFS scheme. Besides, in view of the discretization effect on Terahertz (THz) BP metasurface implementation, the MZT technique is added deliberately into this proposal in that it is more suitable for simulating fine structures. Finally, the ME-based proposal can maintain high absorption when it comes closer to terminate finite-difference time-domain (FDTD) domains so that smaller physical regions can be applied, leading to palpable improvements in memory and CPU time. The three-dimensional (3D) numerical simulations on monolayer BP THz problems have been carried out to illustrate the validity and accuracy of the proposal.

Based on Transition Path Theory (TPT) for Markov jump processes [1, 2], we develop a general approach for identifying and calculating Transition States (TS) of stochastic chemical reacting networks. We first extend the concept of probability current, originally defined on edges connecting different nodes in the configuration space [2], to each sub-network. To locate sub-networks with maximal probability current on the separatrix between reactive and non-reactive events, which will give the Transition States of the reaction, constraint optimization is conducted. We further introduce an alternative scheme to compute the transition pathways by topological sorting, which is shown to be highly efficient through analysis.

Within the mode-coupling theory (MCT) of the glass transition, we reconsider the numerical schemes to evaluate the MCT functional. Here we propose nonuniform discretizations of the wave number, in contrast to the standard equidistant grid, in order to decrease the number of grid points without losing accuracy. We discuss in detail how the integration scheme on the new grids has to be modified from standard Riemann integration. We benchmark our approach by solving the MCT equations numerically for mono-disperse hard disks and hard spheres and by computing the critical packing fraction and the nonergodicity parameters. Our results show that significant improvements in performance can be obtained employing a nonuniform grid.