Parallel implementation of the projector augmented plane wave method for charged systems

A parallel implementation of the projector augmented plane wave (PAW) method with the applications to several transition metal complexes is presented. A unique aspect of our PAW code is that it can treat both charged and neutral cluster systems. We discuss how this is achieved via accurate numerical treatment of the Coulomb Green's function with free space boundary conditions. The strategy for parallelizing the PAW code is based on distributing the plane wave basis across processors. This is a versatile approach and is implemented using a parallel three-dimensional Fast Fourier Transformation (FFT). We report parallel performance analysis of our program and of the three-dimensional FFT's and discuss large-scale parallelization issues of the PAW code. Using a series of transition metal monoxides and dioxides, as well as two iron aqueous complexes, it is shown that a free space PAW code can give structural parameters and energies in good accord with Gaussian based methods.     

Flow-induced waterway in a heterogeneous granular material

In a heterogeneous granular material, viscous flow concentrates to regions with lower particle number density, or higher permeability region, denoted here by “macroscopic cavity”. This in turn enhances the normal stress toward the fluid region on the upstream boundary, which destroys the boundary if the local stress exceeds a certain magnitude. The latter may further enhance the concentration of flow into the cavity region, which is repeated to form a large scale fluidized region toward upstream direction. These processes have been elucidated in our previous experiment using glass beads layer confined between two parallel plane walls. In this paper, a numerical simulation taking into account of the global flow field on the basis of the generalized Darcy?s equation as well as the local stick-slip equation on the basis of the Newton?s equation of motion is performed. Our numerical simulation can successfully reproduce our experimental findings by imposing a suitable pressure gradient and frictional coefficient. The present numerical method can be applied to more general distribution of cavities, including three-dimensional ones, which may predict the formation of long underground waterways or creation of the passage of blood flows (angiogenesis).     

Perturbation and initial Reynolds number effects on transition attainment of supercritical, binary, temporal mixing layers

Two- and three-dimensional (2D and 3D) numerical simulations are performed for a heptane/nitrogen (thermodynamically) supercritical mixing layer initially perturbed at different wavelengths, including the most unstable incompressible wavelength. Simulations are performed with spanwise (and streamwise, for 3D) perturbations available in the literature (for direct numerical simulations of turbulent flow) superimposed on the mean flow, and the domain length is four times the perturbation wavelength. The 2D simulations are undertaken to ascertain that perturbations having the shortest unstable wavelength obtained from a linear inviscid stability analysis are unstable for the viscous non-linear flow. For 3D layers, the purpose of the perturbations is to accelerate the growth of the layer in order to attain transitional Reynolds numbers, as well as to generate structures similar to those that have been observed in spatial mixing layers. The goal of the 3D simulations is to ascertain whether perturbing the mixing layer at different wavelengths, in contrast to the most unstable incompressible wavelength as had previously been done, will affect the transition to turbulence. In particular, we inquire whether perturbing the layer at smaller wavelengths, which requires a smaller domain, will reduce the computational time. It is found that transition can be obtained at different perturbation wavelengths, provided that they are longer than the shortest unstable wavelength as determined by the 2D linear inviscid stability analysis, and provided that the initial Reynolds number is proportionally increased as the wavelength is decreased. The transitional states thus obtained display different dynamic and mixture characteristics, and show strong departures from perfect gas, ideal mixtures. The smaller wavelength perturbations were found to have similar computational requirements for transition attainment.     

Boundary element stress analysis for copper-based dies in pressure die casting

It has recently been demonstrated that, for pressure die casting, high rates of heat extraction are possible with the use copper–alloyed dies suitably protected with a thermally sprayed steel layer. The structural integrity of dies of this type is paramount as the thermally sprayed layer has the potential to de-bond necessitating repair or retirement of the dies.In this paper, an efficient three-dimensional elastostatic stress model for the pressure die casting process is described. The collocation based boundary element method is used for the prediction of transient stress fields over a thermally stabilised casting cycle. A peculiar feature of the pressure die casting process is that transient thermal penetration into the die is limited to regions close to the surface of the die cavity and nozzle. The presence of a transient thermal field necessitates the evaluation of domain integrals in the boundary element stress formulation. Two methods for evaluating these integrals are presented in the paper: (i) a simplex method on a mesh local to the cavity surface; (ii) a modified reciprocity method utilising Gaussian radial-basis functions, interpolating on a perturbed thermal field. The simplex method is shown to be superior, providing high accuracy, stability and computational efficiency. The dies present a multi-domain environment for stress predictions making it necessary to utilise a suitably constructed coarse preconditioner to enhance numerical stability and provide for efficient computation. A multiplicative Schwarz method is presented that enables parameter matrix accelerated GMRES to be applied on each domain. Numerical experiments are performed to demonstrate the computational effectiveness of the approach. Predicted strain fields are compared with strain gauge measurements obtained on a purpose built rig designed to be representative of the casting process.     

Parallel TREE code for two-component ultracold plasma analysis

The TREE method has been widely used for long-range interaction N-body problems. We have developed a parallel TREE code for two-component classical plasmas with open boundary conditions and highly non-uniform charge distributions. The program efficiently handles millions of particles evolved over long relaxation times requiring millions of time steps. Appropriate domain decomposition and dynamic data management were employed, and large-scale parallel processing was achieved using an intermediate level of granularity of domain decomposition and ghost TREE communication. Even though the computational load is not fully distributed in fine grains, high parallel efficiency was achieved for ultracold plasma systems of charged particles. As an application, we performed simulations of an ultracold neutral plasma with a half million particles and a half million time steps. For the long temporal trajectories of relaxation between heavy ions and light electrons, large configurations of ultracold plasmas can now be investigated, which was not possible in past studies.     

An efficient parallel algorithm with application to computational fluid dynamics

When solving time-dependent partial differential equations on parallel computers using the nonoverlapping domain decomposition method, one often needs numerical boundary conditions on the boundaries between subdomains. These numerical boundary conditions can significantly affect the stability and accuracy of the final algorithm.In this paper, a stability and accuracy analysis of the existing methods for generating numerical boundary conditions will be presented, and a new approach based on explicit predictors and implicit correctors will be used to solve convection-diffusion equations on parallel computers, with application to aerospace engineering for the solution of Euler equations in computational fluid dynamics simulations. Both theoretical analyses and numerical results demonstrate significant improvement in stability and accuracy by using the new approach.     

基于二维结构化网格的可压缩流体并行算法研究

基于二维/轴对称高精度可压缩多相流计算流体力学方法 MuSiC-CCASSIM的结构化网格部分,设计了区域并行分解方法;针对各处理器边界数据的通信,设计了阻塞式通信与非阻塞式通信并行算法;为了减少通信开销,设计了MPI/OpenMP混合并行优化算法。在天河二号超级计算机上进行了测试,每个核固定网格规模为625*250,最多调用8 192核。测试数据表明,采用MPI/OpenMP混合并行算法、纯MPI非阻塞式通信并行算法和纯MPI阻塞式通信并行算法的程序的平均并行效率分别达到86%、83%和77%,三种算法都具有良好的可扩展性。     

基于GPU的GRAPES数值预报系统中RRTM模块的并行化研究

GRAPES(Global and Regional Assimilation and Prediction System)是由中国气象科学研究院自主研究开发的中国新一代数值天气预报系统,由于其处理的数据量非常庞大以及对实时性的要求较高,因此一直是并行计算领域研究的热点。首次运用GPU(图形处理器)通用计算及CUDA技术对CRAPES_Meso。模式中物理过程的RRTM(快速辐射传输模式)长波辐射模块进行并行化处理。在性能分析的基础上,针对GPU体系结构的特点,从代码优化、存储器优化、编译选项等方面对程序性能进行优化,并取得了14X倍的加速比。经过测试表明,长波辐射RRTM模块在GPU上并行计算过程正确、稳定而且有效,并为GRAPES系统未来在GPU平台上的并行化发展奠定了一定的基础。     

A finite element procedure for radio-frequency sheath–plasma interactions in the ion cyclotron range of frequencies

A numerical code that solves self-consistent radio-frequency (RF) sheath–plasma interactions in the ion cyclotron range of frequencies is developed based on a nonlinear finite element technique. The present code solves for plasma waves based on the cold plasma model subject to a sheath boundary condition. The finite element procedure has been implemented for one- and two-dimensional analyses using simplified models for the poloidal plane of a tokamak. The results show good accuracy and generalize previous analytical calculations. The present algorithmic approach shows promise for developing a code to predict RF sheath potentials in the scrape-off layer of RF-heated fusion experiments.     

Multi-scale simulations of plasma with iPIC3D

The implicit Particle-in-Cell method for the computer simulation of plasma, and its implementation in a three-dimensional parallel code, called iPIC3D, are presented. The implicit integration in time of the Vlasov–Maxwell system, removes the numerical stability constraints and it enables kinetic plasma simulations at magnetohydrodynamics time scales. Simulations of magnetic reconnection in plasma are presented to show the effectiveness of the algorithm.     

An easily implemented task-based parallel scheme for the Fourier pseudospectral solver applied to 2D Navier-Stokes turbulence

An efficient parallel scheme is proposed for performing direct numerical simulation (DNS) of two-dimensional Navier-Stokes turbulence at high Reynolds numbers. We illustrate the resulting numerical code by displaying relaxation to states close to those that have been predicted by statistical-mechanical methods which start from ideal (Euler) fluid mechanics. The validation of these predictions by DNS requires unusually long computation times on single-cpu workstations, and suggests the use of parallel computation. The performance of our MPI Fortran 90 code on the SGI Origin 3800 is reported, together with its comparison with another parallel method. A few computational results that illustrate tests of the statistical-mechanical predictions are presented.     

Scalability analysis of parallel Particle-In-Cell codes on computational grids

We have performed benchmarks of two three-dimensional parallel Particle-In-Cell (PIC) codes that are similar but have quite different communication patterns on different computational Grids. An electrostatic code with only electrons based on the three-dimensional skeleton PIC code employs the FFT Poisson solver that uses collective communication patterns. Another is the TRISTAN (TRI-dimensional STATNford) code parallelized with MPI, an electromagnetic full particle code, which uses a field solver that only requires point-to-point neighbor communication patterns. We present the mpptest benchmarks on cluster-based computational Grids, where both the basic point-to-point communication patterns and the basic collective communication patterns used in these PIC codes are tested. The results of these benchmarks clearly allow us to quantify and understand the scalability of both communication patterns on the Grids. The present results show that the parallelized TRISTAN code (without all-to-all collective communication) is more scalable than the parallelized skeleton PIC code (with all-to-all collective communication), in cluster-based computational Grid systems where communication performances is poor.     

Gas-kinetic numerical study of complex flow problems covering various flow regimes

The Boltzmann simplified velocity distribution function equation, as adapted to various flow regimes, is described on the basis of the Boltzmann–Shakhov model from the kinetic theory of gases in this study. The discrete velocity ordinate method of gas-kinetic theory is studied and applied to simulate complex multi-scale flows. On the basis of using the uncoupling technique on molecular movements and collisions in the DSMC method, the gas-kinetic finite difference scheme is constructed by extending and applying the unsteady time-splitting method from computational fluid dynamics, which directly solves the discrete velocity distribution functions. The Gauss-type discrete velocity numerical quadrature technique for flows with different Mach numbers is developed to evaluate the macroscopic flow parameters in the physical space. As a result, the gas-kinetic numerical algorithm is established for studying the three-dimensional complex flows with high Mach numbers from rarefied transition to continuum regimes. On the basis of the parallel characteristics of the respective independent discrete velocity points in the discretized velocity space, a parallel strategy suitable for the gas-kinetic numerical method is investigated and, then, the HPF (High Performance Fortran) parallel programming software is developed for simulating gas dynamical problems covering the full spectrum of flow regimes. To illustrate the feasibility of the present gas-kinetic numerical method and simulate gas transport phenomena covering various flow regimes, the gas flows around three-dimensional spheres and spacecraft-like shapes with different Knudsen numbers and Mach numbers are investigated to validate the accuracy of the numerical methods through HPF parallel computing. The computational results determine the flow fields in high resolution and agree well with the theoretical and experimental data. This computing, in practice, has confirmed that the present gas-kinetic algorithm probably provides a promising approach for resolving hypersonic aerothermodynamic problems with the complete spectrum of flow regimes from the gas-kinetic point of view for solving the mesoscopic Boltzmann model equation.     

Concurrent finite element analysis of periodic boundary value problems

A parallel finite element approach for analyzing micromechanical problems with periodic unit cells is discussed. The method uses a direct solution strategy so that general periodic boundary conditions can be treated using a two-step domain decomposition strategy. The speedup results show a good performance of the method on coarse-grained problems, i.e. for cases where the computational work done on the substructures that are treated in parallel is relatively large compared to the total amount of computational work. Application examples using crystal-plasticity on an array of planar crystals and a metal matrix composite are used to show that the overall response of these materials is rather strongly dependent on the constraint imposed on the unit cell so that a correct treatment of the periodic boundary conditions is required to accurately predict the macroscopic response of a periodic material even though a unit cell with a large number of grains or fibers is used.     

Comparing methods for the modelling of boundary-driven streaming in acoustofluidic devices

Numerical simulations of acoustic streaming flows can be used not only to explain the complex phenomena observed in acoustofluidic manipulation devices, but also to predict and optimise their performances. In this paper, two numerical methods based on perturbation theory are compared in order to demonstrate their viability and applicability for modelling boundary-driven streaming flows in acoustofluidic systems. It was found that the Reynolds stress method, which predicts the streaming fields from their driving terms, can effectively resolve both the inner and outer streaming fields and can be used to demonstrate the driving mechanisms of a broad range of boundary-driven streaming flows. However, computational efficiency typically limits its useful application to two-dimensional models. We highlight the close relationship between the classical boundary-driven streaming vortices and the rotationality of the Reynolds stress force field. The limiting velocity method, which ignores the acoustic boundary layer and solves the outer streaming fields by applying the ‘limiting velocities’ as boundary conditions, is more computationally efficient and can be used for predicting three-dimensional outer streaming fields and provide insight into their origins, provided that the radius of curvature of the channel surfaces is much greater than the acoustic boundary layer thickness (\(\delta_{v}\)). We also show that for the limiting velocity method to be valid the channel scales must exceed a value of approximately 100 \(\delta_{v}\) (for an error of ~5% on the streaming velocity magnitudes) for the case presented in this paper. Comparisons of these two numerical methods can provide effective guidance for researchers in the field of acoustofluidics on choosing appropriate methods to predict boundary-driven streaming fields in the design of acoustofluidic particle manipulation devices.     

Analysis and shape optimization in incompressible flows with unstructured grids

The aim of this study is to develop and validate numerical methods that perform shape optimization in incompressible flows using unstructured meshes. The three-dimensional Euler equations for compressible flow are modified using the idea of artificial compressibility and discretized on unstructured tetrahedral grids to provide estimates of pressure distributions for aerodynamic configurations. Convergence acceleration techniques like multigrid and residual averaging are used along with parallel computing platforms to enable these simulations to be performed in a few minutes. This computational frame-work is used to analyze sail geometries. The adjoint equations corresponding to the “incompressible” field equations are derived along with the functional form of gradients. The evaluation of the gradients is reduced to an integral around the boundary to circumvent hurdles posed by adjoint-based gradient evaluations on unstructured meshes. The reduced gradient evaluations provide major computational savings for unstructured grids and its accuracy and use for canonical and industrial problems is a major contribution of this study. The design process is driven by a steepest-descent algorithm with a fixed step-size. The feasibility of the design process is demonstrated for three inverse design problems, two canonical problems and one industrial problem.     

Parallel three-dimensional direct simulation Monte Carlo method and its applications

The development of a parallel three-dimensional direct simulation Monte Carlo (DSMC) method using unstructured cells is reported. Variable hard sphere molecular model and no time counter method are used for the molecular collision kinetics, while the cell-by-cell ray-tracing technique is implemented for particle movement. Developed serial code has been verified by comparing the results of a supersonic corner flow with those of Bird’s three-dimensional structured DSMC code. In addition, a benchmark test is performed for an orifice expanding flow to verify the parallel implementation of DSMC method by comparing with available experimental data. Static physical domain decomposition is used to distribute the workload among multiple processors by considering the estimated particle weighting distribution. Two-step multi-level graph partitioning technique is used to perform the required domain decomposition. Completed code is then applied to compute a hypersonic flow over a sphere (external flow) and the flow field of a spiral drag pump (internal flow), respectively. Results of the former are in good agreement with previous numerical results using axisymmetric DSMC method and experimental data. Results of the latter also agree well with previous numerical results.     

A new 3D parallel SPH scheme for free surface flows

We propose a new robust and accurate SPH scheme, able to track correctly complex three-dimensional non-hydrostatic free surface flows and, even more important, also able to compute an accurate and little oscillatory pressure field. It uses the explicit third order TVD Runge-Kutta scheme in time, following Shu and Osher [Shu C-W, Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J Comput Phys 1988;89:439-71], together with the new key idea of introducing a monotone upwind flux for the density equation, thus removing any artificial viscosity term. For the discretization of the velocity equation, the non-diffusive central flux has been used. A new flexible approach to impose the boundary conditions at solid walls is also proposed. It can handle any moving rigid body with arbitrarily irregular geometry. It does neither produce oscillations in the fluid pressure in proximity of the interfaces, nor does it have a restrictive impact on the stability condition of the explicit time stepping method, unlike the repellent boundary forces of Monaghan [Monaghan JJ. Simulating free surface flows with SPH. J Comput Phys 1994;110:399-406]. To asses the accuracy of the new SPH scheme, a 3D mesh-convergence study is performed for the strongly deforming free surface in a 3D dam-break and impact-wave test problem providing very good results.Moreover, the parallelization of the new 3D SPH scheme has been carried out using the message passing interface (MPI) standard, together with a dynamic load balancing strategy to improve the computational efficiency of the scheme. Thus, simulations involving millions of particles can be run on modern massively parallel supercomputers, obtaining a very good performance, as confirmed by a speed-up analysis. The 3D applications consist of environmental flow problems, such as dam-break flows and impact flows against a wall. The numerical solutions obtained with our new 3D SPH code have been compared with either experimental results or with other numerical reference solutions, obtaining in all cases a very satisfactory agreement.     

基于标记语言的跨平台并行编程框架设计

大量遗留的串行代码需要进行并行化改造,而并行程序复杂性及并行计算平台多样性导致改造成本较高.为此,设计了一种基于标记语言的三层并行编程框架,完成了从串行程序层到并行中间代码层、并行中间代码层到目标并行编程语言程序层的二个转换阶段.采用对串行代码进行语言标记的方法来实现并行中间代码层,该代码层实际是共享存储、分布式存储并行平台编程语言的一种抽象.该框架还实现了一种性能标记方法,可用于并行参数自动寻优.用于雷达数据处理的实验结果表明,实现了对应并行代码的生成,且并行加速比与人工实现的并行代码相当.     

Numerical experiments with the lid driven cavity flow problem

The centerline velocity profiles obtained from the solution of the two- and three-dimensional representations of the lid driven cavity flow problem are compared for different Reynolds numbers. Two configurations were used in this study: a unit cavity and a cavity with an aspect ratio of 2. The Reynolds numbers ranged from 100 to 5000 for all of the configurations studied. A new method of extending the Jacobi collocation technique called spectral difference is developed in this paper together with a unique computational grid. In addition, an iterative method for solving the pressure problem is also developed. This new numerical method allowed the calculation of three-dimensional Navier-Stokes equations to be performed in computers with very modest computational capabilities such as workstations.