Half Time Step Research Articles

Asymptotic-preserving gyrokinetic implicit particle-orbit integrator for arbitrary electromagnetic fields

We extend the asymptotic preserving and energy conserving time integrator for charged-particle motion developed in Ricketson and Chacón (2020) [36] to include finite Larmor-radius (FLR) effects in the presence of electric-field length-scales comparable to the particle gyro-radius (the gyro-kinetic limit). We introduce two modifications to the earlier scheme. The first is the explicit gyro-averaging of the electric field at the half time-step, along with an analogous modification to the current deposition, which we show preserves total energy conservation in implicit PIC schemes. The number of gyrophase samples is chosen adaptively, ensuring proper averaging for large timesteps and the recovery of full-orbit dynamics in the small time-step limit. The second modification is an alternating large and small time-step strategy that ensures the particle trajectory samples gyrophases evenly. We show that this strategy relaxes the time-step restrictions on the scheme, allowing even larger speed-ups than previously achievable. We demonstrate the new method with several single-particle motion tests in a variety of electromagnetic field configurations featuring gyro-scale variation in the electric field. The results demonstrate the advertised ability to capture FLR effects accurately even when significantly stepping over the gyration time-scale.

Novel barostat implementation for molecular dynamics.

We propose a novel implementation of the extended-dynamics equations for isothermal-isobaric ensemble in molecular dynamics, as the Martyna-Tobias-Klein thermostat and barostat. This method is suitable for systems with constraints and the Verlet-family integrators. Instead of iterations or the Trotter-expansion-based methods, both velocities and box sizes (scaling of bond lengths) are predicted. The algorithm begins with force calculation, requiring neither quarter nor half time steps, and necessitating iterations only inside SHAKE. Several tests demonstrate that the quality is comparable to other implementations. It is found that the formula relating the extended barostat mass to the characteristic time of volume fluctuations is inaccurate for condensed systems, which has consequences for the parameter setup. Emphasis is also put on the verification of the precise isothermal-isobaric ensemble and finite-size effects.

An improved partitioned time stepping method based on modified characteristic FEM for the evolutionary dual-porosity-Navier–Stokes model

In this paper, we study an improved partitioned time stepping method based on the modified characteristic finite element method for the evolutionary dual-porosity-Navier–Stokes coupling model with the Beavers–Joseph interface condition. By referring to the idea of partitioned time stepping, we lead into the half-time steps to separate the original coupling problem into three subproblems, implement two-level decoupling to simplify the complexity of solving coupling systems. Before that, the matrix and microfracture equations on the half-time steps in the dual-porosity media are calculated first, with the aim of replacing the initial values of the decoupling equations. Besides, for the nonlinear convection term in Navier–Stokes equations, the modified characteristic finite element method is chosen to avoid the problem of low computational efficiency. The convergence of the decoupling algorithm is rigorously derived by using the idea of mathematical induction. Finally, the reliability of the theoretical analysis and the rationality of the scheme are verified by numerical experiments.

Novel Schemes of No-Slip Boundary Conditions for the Discrete Unified Gas Kinetic Scheme Based on the Moment Constraints.

The boundary conditions are crucial for numerical methods. This study aims to contribute to this growing area of research by exploring boundary conditions for the discrete unified gas kinetic scheme (DUGKS). The importance and originality of this study are that it assesses and validates the novel schemes of the bounce back (BB), non-equilibrium bounce back (NEBB), and Moment-based boundary conditions for the DUGKS, which translate boundary conditions into constraints on the transformed distribution functions at a half time step based on the moment constraints. A theoretical assessment shows that both present NEBB and Moment-based schemes for the DUGKS can implement a no-slip condition at the wall boundary without slip error. The present schemes are validated by numerical simulations of Couette flow, Poiseuille flow, Lid-driven cavity flow, dipole-wall collision, and Rayleigh-Taylor instability. The present schemes of second-order accuracy are more accurate than the original schemes. Both present NEBB and Moment-based schemes are more accurate than the present BB scheme in most cases and have higher computational efficiency than the present BB scheme in the simulation of Couette flow at high Re. The present Moment-based scheme is more accurate than the present BB, NEBB schemes, and reference schemes in the simulation of Poiseuille flow and dipole-wall collision, compared to the analytical solution and reference data. Good agreement with reference data in the numerical simulation of Rayleigh-Taylor instability shows that they are also of use to the multiphase flow. The present Moment-based scheme is more competitive in boundary conditions for the DUGKS.

New decoupled method for the evolutionary dual-porosity-Stokes model with Beavers-Joseph interface conditions

In this paper, we present and analyze a decoupled method for the nonstationary dual-porosity-Stokes coupling problem with the Beavers-Joseph interface conditions. Instead of solving the whole coupling problem on each time step, the proposed scheme allows the coupling problem to be divided into four subproblems in a non-iterative manner, which improves the computational efficiency. Specifically, a Stokes subproblem on a half-time step is calculated first, by using its solution as the interface information, we solve uncoupled matrix pressure equation, microfracture pressure equation and Stokes equation separately. The stability analysis and error estimation of the decoupled algorithm are presented. They show that if the rescaling factor is small enough, our scheme is stable on a bounded time interval and optimal convergent. The reliability of the theoretical analysis and the feasibility of the scheme are verified by the numerical experiments.

New Stable, Explicit, Shifted-Hopscotch Algorithms for the Heat Equation

Our goal was to find more effective numerical algorithms to solve the heat or diffusion equation. We created new five-stage algorithms by shifting the time of the odd cells in the well-known odd-even hopscotch algorithm by a half time step and applied different formulas in different stages. First, we tested 105 = 100,000 different algorithm combinations in case of small systems with random parameters, and then examined the competitiveness of the best algorithms by testing them in case of large systems against popular solvers. These tests helped us find the top five combinations, and showed that these new methods are, indeed, effective since quite accurate and reliable results were obtained in a very short time. After this, we verified these five methods by reproducing a recently found non-conventional analytical solution of the heat equation, then we demonstrated that the methods worked for nonlinear problems by solving Fisher’s equation. We analytically proved that the methods had second-order accuracy, and also showed that one of the five methods was positivity preserving and the others also had good stability properties.

On the numerical solution of Burgers-Fisher equation by the Strang Splitting Method

In this paper, we apply the Strang splitting method to approximate the solution of Burgers-Fisher equation. This method offers a second-order accurate approximation to evolution equation by combining two non-commuting operators. This accuracy is achieved through a symmetric decomposition in which one operator is applied twice for half timesteps, and the other operator is applied once for a full timestep. The stability criteria is derived using the invariant region. The numerical results obtained are compared with the exact solution. The numerical error shows that the exact and the numerical solution are agrees with each other with a good accuracy.

A Two-Stage Fourth-Order Multimoment Global Shallow-Water Model on the Cubed Sphere

AbstractA new multimoment global shallow-water model on the cubed sphere is proposed by adopting a two-stage fourth-order Runge–Kutta time integration. Through calculating the values of predicted variables at half time step t = tn + (1/2)Δt by a second-order formulation, a fourth-order scheme can be derived using only two stages within one time step. This time integration method is implemented in our multimoment global shallow-water model to build and validate a new and more efficient numerical integration framework for dynamical cores. As the key task, the numerical formulation for evaluating the derivatives in time has been developed through the Cauchy–Kowalewski procedure and the spatial discretization of the multimoment finite-volume method, which ensures fourth-order accuracy in both time and space. Several major benchmark tests are used to verify the proposed numerical framework in comparison with the existing four-stage fourth-order Runge–Kutta method, which is based on the method of lines framework. The two-stage fourth-order scheme saves about 30% of the computational cost in comparison with the four-stage Runge–Kutta scheme for global advection and shallow-water models. The proposed two-stage fourth-order framework offers a new option to develop high-performance time marching strategy of practical significance in dynamical cores for atmospheric and oceanic models.

A simplified discrete unified gas kinetic scheme for incompressible flow

The discrete unified gas kinetic scheme (DUGKS) is a new finite volume (FV) scheme for continuum and rarefied flows, which combines the benefits of both the lattice Boltzmann method and UGKS. By the reconstruction of the gas distribution function using particle velocity characteristic lines, the flux contains more detailed information of fluid flow and more concrete physical nature. In this work, a simplified DUGKS is proposed with the reconstruction stage on a whole time step instead of a half time step in the original DUGKS. Using the temporal/spatial integral Boltzmann Bhatnagar–Gross–Krook equation, the auxiliary distribution function with the inclusion of the collision effect is adopted. The macroscopic and mesoscopic fluxes of the cell on the next time step are predicted by the reconstruction of the auxiliary distribution function at interfaces along particle velocity characteristic lines. According to the conservation law, the macroscopic variables of the cell on the next time step can be updated through its flux, which is a moment of the predicted mesoscopic flux at cell interfaces. The equilibrium distribution function on the next time step can also be updated. The gas distribution function is updated by the FV scheme through its predicted mesoscopic flux in a time step. Compared with the original DUGKS, the computational process of the proposed method is more concise because of the omission of half time step flux calculation. The numerical time step is only limited by the Courant–Friedrichs–Lewy condition, and a relatively good stability has been preserved. Several test cases, including the Couette flow, lid-driven cavity flow, laminar flows over a flat plate, a circular cylinder, and an airfoil, and microcavity flow cases, are conducted to validate the present scheme. The observed numerical simulation results reasonably agree with the reported results.

Thermodynamical Extension of a Symplectic Numerical Scheme with Half Space and Time Shifts Demonstrated on Rheological Waves in Solids.

On the example of the Poynting–Thomson–Zener rheological model for solids, which exhibits both dissipation and wave propagation, with nonlinear dispersion relation, we introduce and investigate a finite difference numerical scheme. Our goal is to demonstrate its properties and to ease the computations in later applications for continuum thermodynamical problems. The key element is the positioning of the discretized quantities with shifts by half space and time steps with respect to each other. The arrangement is chosen according to the spacetime properties of the quantities and of the equations governing them. Numerical stability, dissipative error, and dispersive error are analyzed in detail. With the best settings found, the scheme is capable of making precise and fast predictions. Finally, the proposed scheme is compared to a commercial finite element software, COMSOL, which demonstrates essential differences even on the simplest—elastic—level of modeling.

The Numerical Solution of Heat Equation by the Fourth-Order Iterative Alternating Decomposition Explicit Method with MPI

The numerical solution of the heat equation in one space dimension is obtained using the Fourth-Order Iterative Alternating Decomposition Explicit Method (4-IADE) on a parallel platform with Message Passing Interface (MPI). Here, a higher fourth-order Crank-Nicolson type scheme is used in the approximation which gives rise to a Penta diagonal matrix in the solution of the system at each time level. The method employs a splitting strategy which is applied alternately at each half time step. The method is shown to be computationally stable and appropriate parameters chosen to accelerate convergence. The accuracy of the method is comparable to that of existing well known methods. Results obtained by this method for several different problems were compared with the exact solution and agreed closely with those obtained by other finite-difference methods with correlation between speedup and efficiency.

Optimal Temperature Evaluation in Molecular Dynamics Simulations with a Large Time Step.

In molecular dynamics (MD) simulations, an accurate evaluation of temperature is essential for controlling temperature as well as pressure in the isothermal-isobaric conditions. According to the Tolman's equipartition theorem, all motions of all particles should share a single temperature. However, conventional temperature estimation from kinetic energy does not include Hessian terms properly, and thereby, the equipartition theorem is not satisfied with a large time step. In this paper, we show how to evaluate temperature the most accurately without increasing computational cost. We define two kinds of kinetic energies, evaluated at full- and half-time steps that underestimate or overestimate temperature, respectively. A combination of these two kinetic energies provides an optimal instantaneous temperature up to the third order of the time step. The method is tested for a one-dimensional harmonic oscillator, pure water molecules, a Bovine pancreatic trypsin inhibitor (BPTI) protein in water molecules, and a hydrated 1,2-dispalmitoyl- sn-phosphatidylcholine (DPPC) lipid bilayer in water molecules. In all tests, the optimal temperature estimator fulfills the equipartition theorem better than existing methods and reproduces well the usual physical properties for time steps up to and including 5 fs.

Symmetrized operator split schemes for force and source modeling in cascaded lattice Boltzmann methods for flow and scalar transport.

Operator split forcing schemes exploiting a symmetrization principle, i.e., Strang splitting, for cascaded lattice Boltzmann (LB) methods in two- and three-dimensions for fluid flows with impressed local forces are presented. Analogous scheme for the passive scalar transport represented by a convection-diffusion equation with a source term in a novel cascaded LB formulation is also derived. They are based on symmetric applications of the split solutions of the changes on the scalar field or fluid momentum due to the sources or forces over half time steps before and after the collision step. The latter step is effectively represented in terms of the post-collision change of moments at zeroth and first orders, respectively, to represent the effect of the sources on the scalar transport and forces on the fluid flow. Such symmetrized operator split cascaded LB schemes are consistent with the second-order Strang splitting and naturally avoid any discrete effects due to forces or sources by appropriately projecting their effects for higher-order moments. All the force or source implementation steps are performed only in the moment space and they do not require formulations as extra terms and their additional transformations to the velocity space. These result in particularly simpler and efficient schemes to incorporate forces or sources in the cascaded LB methods unlike those considered previously. Numerical study for various benchmark problems in 2D and 3D for fluid flow problems with body forces and scalar transport with sources demonstrate the validity and accuracy, as well as the second-order convergence rate of the symmetrized operator split forcing or source schemes for the cascaded LB methods.

Kinetic energy definition in velocity Verlet integration for accurate pressure evaluation.

In molecular dynamics (MD) simulations, a proper definition of kinetic energy is essential for controlling pressure as well as temperature in the isothermal-isobaric condition. The virial theorem provides an equation that connects the average kinetic energy with the product of particle coordinate and force. In this paper, we show that the theorem is satisfied in MD simulations with a larger time step and holonomic constraints of bonds, only when a proper definition of kinetic energy is used. We provide a novel definition of kinetic energy, which is calculated from velocities at the half-time steps (t - Δt/2 and t + Δt/2) in the velocity Verlet integration method. MD simulations of a 1,2-dispalmitoyl-sn-phosphatidylcholine (DPPC) lipid bilayer and a water box using the kinetic energy definition could reproduce the physical properties in the isothermal-isobaric condition properly. We also develop a multiple time step (MTS) integration scheme with the kinetic energy definition. MD simulations with the MTS integration for the DPPC and water box systems provided the same quantities as the velocity Verlet integration method, even when the thermostat and barostat are updated less frequently.

A three-step Boris integrator for Lorentz force equation of charged particles

This paper provides a numerical procedure for integrating the equations of motion for electrically charged particles in magnetic fields with higher accuracy than the conventional Boris integrator. It is confirmed by both of numerical test and theoretical analysis that the proposed three-step integrator has the same accuracy as the standard (two-step) Boris integrator with a half time step.

Improved quaternion-based integration scheme for rigid body motion

Rotation quaternions are frequently used for describing the orientation of non-spherical rigid bodies. Their compact representation by four numbers and disappearance of numerical problems, such as gimbal lock, are reasons for using them. We describe an improvement of a predictor–corrector direct multiplication (PCDM) method for numerically integrating the rigid body equations of motion with rotation quaternions. The method only uses quaternions to describe the orientation, so no rotation matrices are needed in the implementation. A predictor–corrector approach is used to update the quaternions each time step, such that no renormalization is needed at the end of the time step. The PCDM method suggested by Zhao and Van Wachem is improved such that forces and torques are calculated at the correct time using position and orientation information at that same time. This is achieved by using a leapfrog approach in the improved version, in which the linear and angular velocities and rotation quaternions are defined at half time steps, while whole time step information of these quantities is calculated as part of the improved integration scheme. The improved PCDM scheme is compared with the original implementation for rotational kinetic energy conservation, accuracy of object orientation and angular velocity, and rate of convergence for different time steps. With the modifications that we propose, the improved method has a true second-order rate of convergence, without the need for explicit renormalization of the quaternions. Furthermore, the method is applicable to problems with position and velocity dependent torques, while still only a single force/torque evaluation is needed per time step. For objects experiencing torque, the improved PCDM method performs better than the original method, now showing a true second-order rate of convergence, and much smaller errors in the prediction of object orientation and angular velocity while still requiring only a single torque evaluation per time step.

Pseudo-Spectral simulation of 1D nonlinear propagation in heterogeneous elastic media

In this work, a 1D Pseudo-Spectral Time Domain (PSTD) algorithm has been developed for solving elastic wave equation in nonlinear heterogeneous solids using FFTs for calculation of the spatial differential operator on staggered grid. The solver uses a staggered fourth order Adams–Bashforth method, by which stress and particle velocity are updated at alternating half time steps, to integrate forward in time. To circumvent wraparound inherent to FFT-based pseudo-spectral simulation, Convolution Perfectly Matched Layer (CPML) boundary condition has been used to eliminate implementation problems linked in classical PML to the introduction in nonlinear elasticity of a time dependent bulk modulus. Different kinds of nonlinear elastic models (quadratic and cubic nonlinearity, Nazarov hysteretic nonlinearity, bi-modular nonlinearity, PM-Space nonlinearity) have been implemented. The present study will focus on the comparison of nonlinear signature (harmonics generation, shock, frequency shift and attenuation) of these different kinds of nonlinearity for rod resonance, shock wave generation. These results are expected to be useful in helping to determine the predominant nonlinear mechanism in a specific experiment.

Simulation of 2-D Dam Break Using Improved Incompressible Smoothed Particle Hydrodynamics Based on Projection Method

This paper deals with numerical modeling of water flow which is generated by the break of a dam. The problem is solved by applying a new Incompressible Smoothed Particle Hydrodynamics (ISPH) algorithm based on projection method. The proposed ISPH model has two steps. In the first step, the incompressibility of fluid is maintained regarding to the changes of intermediate and initial particles densities at the first half-time step (stability step). In the second step, by computing the divergence of the intermediate secondary velocity at the second half-time step (accuracy step), the incompressibility is satisfied completely. In fact, by using this method both stability and accuracy are increased. The simulation illustrates the formation and subsequent propagation of a wave over the horizontal plane. It is shown that the model predictions compare well with experimental data.

On particle movers in cylindrical geometry for Particle-In-Cell simulations

Three movers for the orbit integration of charged particles in a given electromagnetic field are analyzed and compared in cylindrical geometry. The classic Boris mover, which is of leap-frog type with position and velocities staggered by half time step, is connected to a second order Strang operator splitting integrator. In general the Boris mover is about 20% faster than the Strang splitting mover without sacrificing much in terms of accuracy. Furthermore, the Boris mover is second order accurate only for a very specific choice of the initial half step needed by the algorithm to get started. Unlike the case in Cartesian geometry, where any initial half step which is at least first order accurate does not compromise the second order accuracy of the method, in cylindrical geometry any attempt to use a more accurate initial half step does in fact decrease the accuracy of the scheme to first order. Through the connection with the Strang operator splitting integrator, this counter-intuitive behavior is explained by the fact that the error in the half step velocities of the Boris mover is proportional to the time step of the simulation.For the case of a uniform and static magnetic field, we discuss the leap-frog cyclotronic mover, cylindrical analogue of the cyclotronic mover of Ref. [L. Patacchini and I. Hutchinson, J. Comput. Phys. 228 (7), 2009], where the step involving acceleration due to inertial forces is combined with the acceleration due to the magnetic part of the Lorentz force. The advantage of a cyclotronic mover is that the gyration of a charged particle in a magnetic field is treated analytically and therefore only the dynamics associated with the electric field needs to be resolved.

COLLOCATED SIBC-FDTD METHOD FOR COATED CONDUCTORS AT OBLIQUE INCIDENCE

A collocated surface impedance boundary condition (SIBC) { flnite difierence time domain (FDTD) method is developed for conductors coated with lossy dielectric coatings at oblique incidence. The method is based on the collocated electric and magnetic fleld components on the planar interface between two media, and rational approximation for tangent function of surface impedance formulation is adopted. In contrast to the traditional SIBC-FDTD implementation which is approximated with the magnetic fleld component on the boundary located at half-cell distance from the interface and half time step earlier in time, the collocation approach is more accurate for both magnitude and phase of re∞ection coe-cient. By the comparison with exact results, the proposed model is numerically verifled in the frequency domain for both parallel polarization plane wave and vertical polarization plane wave at varying oblique angles of incidence.