Feynman’s path integral reformulates the quantum Schrödinger differential equation to be an integral equation. It has been being widely used to compute internuclear quantum-statistical effects on many-body molecular systems. In this Review, the molecular Schrödinger equation will first be introduced, together with the Born-Oppenheimer approximation that decouples electronic and internuclear motions. Some effective semiclassical potentials, e.g., centroid potential, which are all formulated in terms of Feynman’s path integral, will be discussed and compared. These semiclassical potentials can be used to directly calculate the quantum canonical partition function without individual Schrödinger’s energy eigenvalues. As a result, path integrations are conventionally performed with Monte Carlo and molecular dynamics sampling techniques. To complement these techniques, we will examine how Kleinert’s variational perturbation (KP) theory can provide a complete theoretical foundation for developing non-sampling/non-stochastic methods to systematically calculate centroid potential. To enable the powerful KP theory to be practical for many-body molecular systems, we have proposed a new path-integral method: automated integration-free path-integral (AIF-PI) method. Due to the integration-free and computationally inexpensive characteristics of our AIF-PI method, we have used it to perform ab initio path-integral calculations of kinetic isotope effects on proton-transfer and RNA-related phosphoryl-transfer chemical reactions. The computational procedure of using our AIF-PI method, along with the features of our new centroid path-integral theory at the minimum of the absolute-zero energy (AMAZE), are also highlighted in this review.

A globally hyperbolic moment system upto arbitrary order for the Wigner equation was derived in [6]. For numerically solving the high order hyperbolic moment system therein, we in this paper develop a preliminary numerical method for this system following the NRxx method recently proposed in [8], to validate the moment system of the Wigner equation. The method developed can keep both mass and momentum conserved, and the variation of the total energy under control though it is not strictly conservative. We systematically study the numerical convergence of the solution to the moment system both in the size of spatial mesh and in the order of the moment expansion, and the convergence of the numerical solution of the moment system to the numerical solution of the Wigner equation using the discrete velocity method. The numerical results indicate that the high order moment system in [6] is a valid model for the Wigner equation, and the proposed numerical method for the moment system is quite promising to carry out the simulation of the Wigner equation.

The numerical simulation of non conservative system is a difficult challenge for two reasons at least. The first one is that it is not possible to derive jump relations directly from conservation principles, so that in general, if the model description is non ambiguous for smooth solutions, this is no longer the case for discontinuous solutions. From the numerical view point, this leads to the following situation: if a scheme is stable, its limit for mesh convergence will depend on its dissipative structure. This is well known since at least [1]. In this paper we are interested in the “dual” problem: given a system in non conservative form and consistent jump relations, how can we construct a numerical scheme that will, for mesh convergence, provide limit solutions that are the exact solution of the problem. In order to investigate this problem, we consider a multiphase flow model for which jump relations are known. Our scheme is an hybridation of Glimm scheme and Roe scheme.

Parallel computing has become an important subject in the field of computer science and has proven to be critical when researching high performance solutions. The evolution of computer architectures (multi-core and many-core) towards a higher number of cores can only confirm that parallelism is the method of choice for speeding up an algorithm. In the last decade, the graphics processing unit, or GPU, has gained an important place in the field of high performance computing (HPC) because of its low cost and massive parallel processing power. Super-computing has become, for the first time, available to anyone at the price of a desktop computer. In this paper, we survey the concept of parallel computing and especially GPU computing. Achieving efficient parallel algorithms for the GPU is not a trivial task, there are several technical restrictions that must be satisfied in order to achieve the expected performance. Some of these limitations are consequences of the underlying architecture of the GPU and the theoretical models behind it. Our goal is to present a set of theoretical and technical concepts that are often required to understand the GPU and its massive parallelism model. In particular, we show how this new technology can help the field of computational physics, especially when the problem is data-parallel. We present four examples of computational physics problems; n-body, collision detection, Potts model and cellular automata simulations. These examples well represent the kind of problems that are suitable for GPU computing. By understanding the GPU architecture and its massive parallelism programming model, one can overcome many of the technical limitations found along the way, design better GPU-based algorithms for computational physics problems and achieve speedups that can reach up to two orders of magnitude when compared to sequential implementations.

We provide an overview of current techniques and typical applications of numerical bifurcation analysis in fluid dynamical problems. Many of these problems are characterized by high-dimensional dynamical systems which undergo transitions as parameters are changed. The computation of the critical conditions associated with these transitions, popularly referred to as ‘tipping points’, is important for understanding the transition mechanisms. We describe the two basic classes of methods of numerical bifurcation analysis, which differ in the explicit or implicit use of the Jacobian matrix of the dynamical system. The numerical challenges involved in both methods arementioned and possible solutions to current bottlenecks are given. To demonstrate that numerical bifurcation techniques are not restricted to relatively low-dimensional dynamical systems, we provide several examples of the application of the modern techniques to a diverse set of fluid mechanical problems.

In this paper we present recent developments concerning a Cell-Centered Arbitrary Lagrangian Eulerian (CCALE) strategy using the Moment Of Fluid (MOF) interface reconstruction for the numerical simulation of multi-material compressible fluid flows on unstructured grids in cylindrical geometries. Especially, our attention is focused here on the following points. First, we propose a new formulation of the scheme used during the Lagrangian phase in the particular case of axisymmetric geometries. Then, the MOF method is considered for multi-interface reconstruction in cylindrical geometry. Subsequently, a method devoted to the rezoning of polar meshes is detailed. Finally, a generalization of the hybrid remapping to cylindrical geometries is presented. These explorations are validated by mean of several test cases using unstructured grid that clearly illustrate the robustness and accuracy of the new method.

This paper investigates the chaotic lid-driven square cavity flows at extreme Reynolds numbers. Several observations have been made from this study. Firstly, at extreme Reynolds numbers two principles add at the genesis of tiny, loose counterclockwise- or clockwise-rotating eddies. One concerns the arousing of them owing to the influence of the clockwise- or counterclockwise currents nearby; the other, the arousing of counterclockwise-rotating eddies near attached to the moving (lid) top wall which moves from left to right. Secondly, unexpectedly, the kinetic energy soon reaches the qualitative temporal limit’s pace, fluctuating briskly, randomly inside the total kinetic energy range, fluctuations which concentrate on two distinct fragments: one on its upper side, the upper fragment, the other on its lower side, the lower fragment, switching briskly, randomly from each other; and further on many small fragments arousing randomly within both, switching briskly, randomly from one another. As the Reynolds number Re → ∞, both distance and then close, and the kinetic energy fluctuates shorter and shorter at the upper fragment and longer and longer at the lower fragment, displaying tall high spikes which enlarge and then disappear. As the time t → ∞ (at the Reynolds number Re fixed) they recur from time to time with roughly the same amplitude. For the most part, at the upper fragment the leading eddy rotates clockwise, and at the lower fragment, in stark contrast, it rotates counterclockwise. At Re=109 the leading eddy — at its qualitative temporal limit’s pace — appears to rotate solely counterclockwise.

This paper investigates the reduction of backscatter radar cross section (RCS) for a rectangular cavity embedded in the ground plane. The bottom of the cavity is coated by a thin, multilayered radar absorbing material (RAM) with possibly different permittivities. The objective is to minimize the backscatter RCS by the incidence of a plane wave over a single or a set of incident angles. By formulating the scattering problem as a Helmholtz equation with artificial boundary condition, the gradient with respect to the material permittivities is determined efficiently by the adjoint state method, which is integrated into a nonlinear optimization scheme. Numerical example shows the RCS may be significantly reduced.

An accurate cartesian method is devised to simulate incompressible viscous flows past an arbitrary moving body. The Navier-Stokes equations are spatially discretized onto a fixed Cartesian mesh. The body is taken into account via the ghost-cell method and the so-called penalty method, resulting in second-order accuracy in velocity. The accuracy and the efficiency of the solver are tested through two-dimensional reference simulations. To show the versatility of this scheme we simulate a three-dimensional self propelled jellyfish prototype.

We present a new numerical method to approximate the solutions of an Euler-Poisson model, which is inherent to astrophysical flows where gravity plays an important role. We propose a discretization of gravity which ensures adequate coupling of the Poisson and Euler equations, paying particular attention to the gravity source term involved in the latter equations. In order to approximate this source term, its discretization is introduced into the approximate Riemann solver used for the Euler equations. A relaxation scheme is involved and its robustness is established. The method has been implemented in the software HERACLES [29] and several numerical experiments involving gravitational flows for astrophysics highlight the scheme.

In this paper, a new numerical method for solving the decoupled forward-backward stochastic differential equations (FBSDEs) is proposed based on some specially derived reference equations. We rigorously analyze errors of the proposed method under general situations. Then we present error estimates for each of the specific cases when some classical numerical schemes for solving the forward SDE are taken in the method; in particular, we prove that the proposed method is second-order accurate if used together with the order-2.0 weak Taylor scheme for the SDE. Some examples are also given to numerically demonstrate the accuracy of the proposed method and verify the theoretical results.

The evolution of precipitates in stressed solids is modeled by coupling a quasi-steady diffusion equation and a linear elasticity equation with dynamic boundary conditions. The governing equations are solved numerically using a boundary integral method (BIM). A critical step in applying BIM is to develop fast algorithms to reduce the arithmetic operation count of matrix-vector multiplications. In this paper, we develop a fast adaptive treecode algorithm for the diffusion and elasticity problems in two dimensions (2D). We present a novel source dividing strategy to parallelize the treecode. Numerical results show that the speedup factor is nearly perfect up to a moderate number of processors. This approach of parallelization can be readily implemented in other treecodes using either uniform or non-uniform point distribution. We demonstrate the effectiveness of the treecode by computing the long-time evolution of a complicated microstructure in elastic media, which would be extremely difficult with a direct summation method due to CPU time constraint. The treecode speeds up computations dramatically while fulfilling the stringent precision requirement dictated by the spectrally accurate BIM.

This paper develops a high-order accurate gas-kinetic scheme in the framework of the finite volume method for the one- and two-dimensional flow simulations, which is an extension of the third-order accurate gas-kinetic scheme [Q.B. Li, K. Xu, and S. Fu, J. Comput. Phys., 229(2010), 6715-6731] and the second-order accurate gas-kinetic scheme [K. Xu, J. Comput. Phys., 171(2001), 289-335]. It is formed by two parts: quartic polynomial reconstruction of the macroscopic variables and fourth-order accurate flux evolution. The first part reconstructs a piecewise cell-center based quartic polynomial and a cell-vertex based quartic polynomial according to the “initial” cell average approximation of macroscopic variables to recover locally the non-equilibrium and equilibrium single particle velocity distribution functions around the cell interface. It is in view of the fact that all macroscopic variables become moments of a single particle velocity distribution function in the gas-kinetic theory. The generalized moment limiter is employed there to suppress the possible numerical oscillation. In the second part, the macroscopic flux at the cell interface is evolved in fourth-order accuracy by means of the simple particle transport mechanism in the microscopic level, i.e. free transport and the Bhatnagar-Gross-Krook (BGK) collisions. In other words, the fourth-order flux evolution is based on the solution (i.e. the particle velocity distribution function) of the BGK model for the Boltzmann equation. Several 1D and 2D test problems are numerically solved by using the proposed high-order accurate gas-kinetic scheme. By comparing with the exact solutions or the numerical solutions obtained the second-order or third-order accurate gas-kinetic scheme, the computations demonstrate that our scheme is effective and accurate for simulating invisid and viscous fluid flows, and the accuracy of the high-order GKS depends on the choice of the (numerical) collision time.

In this paper we develop a class of Implicit-Explicit Runge-Kutta schemes for solving the multi-scale semiconductor Boltzmann equation. The relevant scale which characterizes this kind of problems is the diffusive scaling. This means that, in the limit of zero mean free path, the system is governed by a drift-diffusion equation. Our aim is to develop a method which accurately works for the different regimes encountered in general semiconductor simulations: the kinetic, the intermediate and the diffusive one. Moreover, we want to overcome the restrictive time step conditions of standard time integration techniques when applied to the solution of this kind of phenomena without any deterioration in the accuracy. As a result, we obtain high order time and space discretization schemes which do not suffer from the usual parabolic stiffness in the diffusive limit. We show different numerical results which permit to appreciate the performances of the proposed schemes.

The carbuncle phenomenon has been regarded as a spurious solution produced by most of contact-preserving methods. The hybrid method of combining high resolution flux with more dissipative solver is an attractive attempt to cure this kind of non-physical phenomenon. In this paper, a matrix-based stability analysis for 2-D Euler equations is performed to explore the cause of instability of numerical schemes. By combining the Roe with HLL flux in different directions and different flux components, we give an interesting explanation to the linear numerical instability. Based on such analysis, some hybrid schemes are compared to illustrate different mechanisms in controlling shock instability. Numerical experiments are presented to verify our analysis results. The conclusion is that the scheme of restricting directly instability source is more stable than other hybrid schemes.

Numerical error caused by “ghost forces” in a quasicontinuum method is studied in the context of dynamic problems. The error in the discrete W1,∞ norm is analyzed for the time scale (ε) and the time scale (1) with ε being the lattice spacing.

In this study, we present a new numerical model for crystal growth in a vertical solidification system. This model takes into account the buoyancy induced convective flow and its effect on the crystal growth process. The evolution of the crystal growth interface is simulated using the phase-field method. A semi-implicit lattice kinetics solver based on the Boltzmann equation is employed to model the unsteady incompressible flow. This model is used to investigate the effect of furnace operational conditions on crystal growth interface profiles and growth velocities. For a simple case of macroscopic radial growth, the phase-field model is validated against an analytical solution. The numerical simulations reveal that for a certain set of temperature boundary conditions, the heat transport in the melt near the phase interface is diffusion dominant and advection is suppressed.

This paper presents a parallel Newton-Krylov-Schwarz method for the numerical simulation of unsteady flows at high Reynolds number around a high-speed train under crosswind. With a realistic train geometry, a realistic Reynolds number, and a realistic wind speed, this is a very challenging computational problem. Because of the limited parallel scalability, commercial CFD software is not suitable for supercomputers with a large number of processors. We develop a Newton-Krylov-Schwarz based fully implicit method, and the corresponding parallel software, for the 3D unsteady incompressible Navier-Stokes equations discretized with a stabilized finite element method on very fine unstructured meshes. We test the algorithm and software for flows passing a train modeled after China’s high-speed train CRH380B, and we also compare our results with results obtained from commercial CFD software. Our algorithm shows very good parallel scalability on a supercomputer with over one thousand processors.

Compressible flow past a circular cylinder at an inflow Reynolds number of 2 x 105 is numerically investigated by using a constrained large-eddy simulation (CLES) technique. Numerical simulation with adiabatic wall boundary condition and at a free-stream Mach number of 0.75 is conducted to validate and verify the performance of the present CLES method in predicting separated flows. Some typical and characteristic physical quantities, such as the drag coefficient, the root-mean-square lift fluctuations, the Strouhal number, the pressure and skin friction distributions around the cylinder, etc. are calculated and compared with previously reported experimental data, finer-grid large-eddy simulation (LES) data and those obtained in the present LES and detached-eddy simulation (DES) on coarse grids. It turns out that CLES is superior to DES in predicting such separated flow and that CLES can mimic the intricate shock wave dynamics quite well. Then, the effects of Mach number on the flow patterns and parameters such as the pressure, skin friction and drag coefficients, and the cylinder surface temperature are studied, with Mach number varying from 0.1 to 0.95. Non-monotonic behaviors of the pressure and skin friction distributions are observed with increasing Mach number and the minimum mean separation angle occurs at a subcritical Mach number of between 0.3 and 0.5. Additionally, the wall temperature effects on the thermodynamic and aerodynamic quantities are explored in a series of simulations using isothermal wall boundary conditions at three different wall temperatures. It is found that the flow separates earlier from the cylinder surface with a longer recirculation length in the wake and a higher pressure coefficient at the rear stagnation point for higher wall temperature. Moreover, the influences of different thermal wall boundary conditions on the flow field are gradually magnified from the front stagnation point to the rear stagnation point. Moreover, the influences of different thermal wall boundary conditions on the flow field are graduallymagnified from the front stagnation point to the rear stagnation point. It is inferred that the CLES approach in its current version is a useful and effective tool for simulating wall-bounded compressible turbulent flows with massive separations.