Harmonic Surface Mapping Algorithm for Electrostatic Potentials in an Atomistic/Continuum Hybrid Model for Electrolyte Solutions

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.     

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 Simple Solver for the Two-Fluid Plasma Model Based on PseudoSpectral Time-Domain Algorithm

We present a solver of 3D two-fluid plasma model for the simulation of short-pulse laser interactions with plasma. This solver resolves the equations of the two-fluid plasma model with ideal gas closure. We also include the Bhatnagar-Gross-Krook collision model. Our solver is based on PseudoSpectral Time-Domain (PSTD) method to solve Maxwell's curl equations. We use a Strang splitting to integrate Euler equations with source term: while Euler equations are solved with a composite scheme mixing Lax-Friedrichs and Lax-Wendroff schemes, the source term is integrated with a fourth-order Runge-Kutta scheme. This two-fluid plasma model solver is simple to implement because it only relies on finite difference schemes and Fast Fourier Transforms. It does not require spatially staggered grids. The solver was tested against several well-known problems of plasma physics. Numerical simulations gave results in excellent agreement with analytical solutions or with previous results from the literature.     

A Fast Time Splitting Finite Difference Approach to Gross–Pitaevskii Equations

We propose an idea to solve the Gross–Pitaevskii equation for dark structures inside an infinite constant background density $ρ_∞$=${|ψ_∞|}^2$, without the introduction of artificial boundary conditions. We map the unbounded physical domain $\mathbb{R}^3$ into the bounded domain ${(−1,1)}^3$ and discretize the rescaled equation by equispaced 4th-order finite differences. This results in a free boundary approach, which can be solved in time by the Strang splitting method. The linear part is solved by a new, fast approximation of the action of the matrix exponential at machine precision accuracy, while the nonlinear part can be solved exactly. Numerical results confirm existing ones based on the Fourier pseudospectral method and point out some weaknesses of the latter such as the need of a quite large computational domain, and thus a consequent critical computational effort, in order to provide reliable time evolution of the vortical structures, of their reconnections, and of integral quantities like mass, energy, and momentum. The free boundary approach reproduces them correctly, also in finite subdomains, at low computational cost. We show the versatility of this method by carrying out one- and three-dimensional simulations and by using it also in the case of Bose–Einstein condensates, for which $ψ$→0 as the spatial variables tend to infinity.     

Explicit Computation of Robin Parameters in Optimized Schwarz Waveform Relaxation Methods for Schrödinger Equations Based on Pseudodifferential Operators

The Optimized Schwarz Waveform Relaxation algorithm, a domain decomposition method based on Robin transmission condition, is becoming a popular computational method for solving evolution partial differential equations in parallel. Along with well-posedness, it offers a good balance between convergence rate, efficient computational complexity and simplicity of the implementation. The fundamental question is the selection of the Robin parameter to optimize the convergence of the algorithm. In this paper, we propose an approach to explicitly estimate the Robin parameter which is based on the approximation of the transmission operators at the subdomain interfaces, for the linear/nonlinear Schrödinger equation. Some illustrating numerical experiments are proposed for the one- and two-dimensional problems.     

Kinetic Slip Boundary Condition for Isothermal Rarefied Gas Flows Through Static Non-Planar Geometries Based on the Regularized Lattice-Boltzmann Method

The simulation of rarefied gas flows through complex porous media is challenging due to the tortuous flow pathways inherent to such structures. The Lattice Boltzmann method (LBM) has been identified as a promising avenue to solve flows through complex geometries due to the simplicity of its scheme and its high parallel computational efficiency. It has been proposed to model the stress-strain relationship with the extended Navier-Stokes equations rather than attempting to directly solve the Boltzmann equation. However, a regularization technique is required to filter out non-resolved higher-order components with a low-order velocity scheme. Although slip boundary conditions (BCs) have been proposed for the non-regularized multiple relaxation time LBM (MRT-LBM) for planar geometries, previous slip BCs have never been verified extensively with the regularization technique. In this work, following an extensive literature review on the imposition of slip BCs for rarefied flows with the LBM, it is proven that earlier values for kinetic parameters developed to impose slip BCs are inaccurate for the regularized MRT-LBM and differ between the D2Q9 and D3Q15 schemes. The error was eliminated for planar flows and good agreement between analytical solutions for arrays of cylinders and spheres was found with a wide range of Knudsen numbers.     

Directional $\mathcal{H}^2$ Compression Algorithm: Optimisations and Application to a Discontinuous Galerkin BEM for the Helmholtz Equation

This study aimed to specialise a directional $\mathcal{H}^2 (\mathcal{D}\mathcal{H}^2)$ compression to matrices arising from the discontinuous Galerkin (DG) discretisation of the hypersingular equation in acoustics. The significant finding is an algorithm that takes a DG stiffness matrix and finds a near-optimal $\mathcal{D}\mathcal{H}^2$ approximation for low and high-frequency problems. We introduced the necessary special optimisations to make this algorithm more efficient in the case of a DG stiffness matrix. Moreover, an automatic parameter tuning strategy makes it easy to use and versatile. Numerical comparisons with a classical Boundary Element Method (BEM) show that a DG scheme combined with a $\mathcal{D}\mathcal{H}^2$ gives better computational efficiency than a classical BEM in the case of high-order finite elements and $hp$ heterogeneous meshes. The results indicate that DG is suitable for an auto-adaptive context in integral equations.     

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 the boundary of the domain and does not require a structured mesh, and only needs a covering 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 is introduced for better error controls. The effect of errors from the Feynman-Kac formula based path integral WOS method on the overall accuracy of the BIE-WOS method is analyzed 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 equation of BIE-WOS method can be simplified where the local BIE solutions can be given in closed 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 of the boundary, resulting in a complete independent computation using a large number of 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-WOS method on a large-scale cluster with 6400 CPU cores is verified for computing the DtN mapping of exterior Laplace problems with Dirichlet data for several domains.     

A Non-Ersatz Material Approach for the Topology Optimization of Elastic Structures Based on Piecewise Constant Level Set Method

The topology optimization of a linearized elasticity system with the area (volume) constraint is investigated. A non-ersatz material approach is proposed. By introducing a fixed background domain, the linearized elasticity system is extended into the background domain by a characteristic function. The piecewise constant level set (PCLS) method is applied to represent the original material region and the void region. A quadratic function of PCLS function is proposed to replace the characteristic function. The functional derivative of the objective functional with respect to PCLS function is derived, which is zero in the void region and nonzero in the original material region. A penalty gradient algorithm is proposed. Four numerical experiments of 2D and 3D elastic structures with different boundary conditions are presented, illustrating the validity of the proposed algorithm.