Long-time Simulation Research Articles

Fokker–Planck equation and Feynman–Kac formula for multidimensional stochastic dynamical systems with Lévy noises and time-dependent coefficients

The aim of this paper is to establish a version of the Feynman–Kac formula for the time-dependent multidimensional nonlocal Fokker–Planck equation corresponding to a class of time-dependent stochastic differential equations driven by multiplicative symmetric (or asymmetric) α-stable Lévy noise. First the forward nonlocal Fokker–Planck equation is derived by the adjoint operator method, overcoming the challenges posed by time-dependent multidimensional nonlinear symmetric α-stable Lévy noise. Subsequently, the Feynman–Kac formula for the forward multidimensional time-dependent nonlocal Fokker–Planck equation is established by applying techniques for the backward nonlocal Fokker–Planck equations, which is associated with the backward stochastic differential equation driven by the multiplicative symmetric α-stable Lévy noise. Notably, in the case of asymmetric α-stable Lévy noise case, the presence of the characteristic function in the nonlocal operator adds complexity to the analysis. Using the Feynman–Kac formula, it is demonstrated that the solution of the forward nonlocal Fokker–Planck equation can be readily simulated through Monte Carlo approximation, especially in scenarios involving long-time simulation settings with large steps. These concepts are illustrated with intriguing examples, and the dynamic evolution of the probability density function corresponding to the stochastic SIS model and the stochastic model for the MeKS network (reflecting the interactions among the MecA complex, ComK and ComS) are investigated over an extended period.

Multi-objective optimization and long-time simulation of a multi-borehole ground heat exchanger system

Multi-objective optimization and CFD simulation are conducted to optimize the design of a multi-borehole ground heat exchanger (GHE) system and assess its long-time performance. The multi-objective optimization is performed to minimize the entropy generation number (EGN) and total cost rate by using various evolutionary algorithms, including NSGA-II, GDE-3, MOEA/D, PESA-II, SPEA-II, and SMPSO. NSGA-II and GDE-3 algorithms perform best in obtaining Pareto optimal solutions. Three prominent points on the NSGA-II Pareto frontier, representing the results of single-objective thermodynamic, single-objective economic, and multi-objective optimizations, are simulated in three dimensions over three months. The trends of EGN variations extracted from the transient CFD simulation agree well with those from the steady analytical model. The EGN obtained from multi-objective optimization is 58.8% lower than the EGN obtained using single-objective economic optimization and 1.9 times higher than that calculated from single-objective thermodynamic optimization. Likewise, the total cost rate obtained from multi-objective optimization is 64.4% lower than the value obtained from single-objective thermodynamic optimization and four times higher than that calculated using single-objective economic optimization. The proposed optimization approach can be reliably applied to improve the design of multi-borehole GHE systems.

A Poisson-bracket scheme for nonlinear shallow-water sloshing in an oscillating tank with irregular bottom surface

The mass-, energy-, and potential-enstrophy-conserving Poisson bracket numerical scheme introduced by Arakawa & Lamb (1981), and extended by Salmon (2004) and Stewart & Dellar (2016), is adapted to the problem of nonlinear shallow-water sloshing over a variable bottom surface in a rigid rectangular basin undergoing a prescribed coupled surge-sway motion. Adaptation to a finite domain requires a new approach to the boundary conditions at solid boundaries in the context of the Arakawa C grid. In this paper, the boundary condition at the vertical walls is taken to be vanishing of the normal derivative of the tangential velocity. This gives zero potential vorticity at the boundary, and is consistent with the material conservation of the potential vorticity. This condition, coupled to symmetric boundary conditions for the wave height arising from a reduced version of the evolution equations at boundaries, leads to an extension of the class of staggered C-grid Poisson-bracket schemes to interior flows with solid boundaries. The scheme is implemented, shown to preserve Casimirs, and applied to a range of problems in shallow water hydrodynamics including interaction of travelling hydraulic jumps at resonance, and X-type soliton interactions over a variable bottom surface. This structure-preserving scheme provides a robust building block for long-time simulation of floating ocean wave energy extractors.

An adaptive energy-based sequential method for training PINNs to solve gradient flow equations

Recent developments in numerous scientific computing fields have been made possible by physics-informed neural networks (PINNs), a class of machine learning techniques that incorporate physical knowledge into neural networks. However, PINNs still encounter difficulties in simulating a kind of partial differential equations (PDEs), i.e. gradient flow equations of complex dynamic systems. Through analysis from the neural tangent kernel's (NTK) point of view, we find that directly training PINNs on the large temporal window may lead to failure. To tackle this problem, we utilize one of the gradient flow equations' significant physical properties, the energy dissipation law, to propose an adaptive energy-based sequential training method for PINNs. Specifically, the gradient of free energy is employed to determine the adaptive temporal steps of the training window, which also enhances the training efficiency. Numerical experiments are demonstrated among gradient flow equations including Cahn-Hilliard, heat conduction, and Allen-Cahn equations. Our method not only outperforms PINNs but also preserves the energy dissipation law, which shows the ability to solve complex long-time simulation problems of gradient flow systems.

Extracting kinetic information from short-time trajectories: relaxation and disorder of lossy cavity polaritons

Abstract The emerging field of molecular cavity polaritons has stimulated a surge of experimental and theoretical activities and presents a unique opportunity to develop the many-body simulation methodology. This paper presents a numerical scheme for the extraction of key kinetic information of lossy cavity polaritons based on the transfer tensor method (TTM). Steady state, relaxation timescales, and oscillatory phenomena can all be deduced directly from a set of transfer tensors without the need for long-time simulation. Moreover, we generalize TTM to disordered systems by sampling dynamical maps and achieve fast convergence to disordered-averaged dynamics using a small set of realizations. Together, these techniques provide a toolbox for characterizing the interplay of cavity loss, disorder, and cooperativity in polariton relaxation and allow us to predict unusual dependences on the initial excitation state, photon decay rate, strength of disorder, and the type of cavity models. Thus, using the example of cavity polaritons, we have demonstrated significant potential in the use of the TTM toward both the efficient computation of long-time polariton dynamics and the extraction of crucial kinetic information about polariton relaxation from a small set of short-time trajectories.

Self-Organized Evolution of the Internal Transport Barrier in Ion-Temperature-Gradient Driven Gyrokinetic Turbulence.

Understanding the self-organization of the most promising internal transport barrier in fusion plasmas needs a long-time nonlinear gyrokinetic global simulation. The neighboring equilibrium update method is proposed, which solves the secularity problem in a perturbative simulation and speeds up the numerical computation by more than 10 times. It is found that the internal transport barrier emerges at the magnetic axis due to inward propagated turbulence avalanche, and its outward expansion is the catastrophe of self-organized structure induced by outward propagated avalanche.

Helium bubble evolution under cascade in bcc iron relevant to fusion conditions investigated by a novel coupling MD-OKMC method

The direct interaction between cascade and helium bubbles during the thermal spike phase plays an important role in the helium bubbles evolution in metals, which was not fully considered by the present long-time simulation methods. To address this issue, in this work we developed a coupling MD-OKMC method to investigate the helium bubble evolution under irradiation cascade in bcc iron. The short-time collision cascade was simulated by the molecular dynamics (MD), and the long-time helium bubble evolution was modeled by the object kinetic Monte Carlo (OKMC). The MD and OKMC simulations were alternated to imitate the gradual increase of both irradiation dose and helium concentration under experiment conditions. The impact of temperature on the helium bubble evolution under irradiation cascade in bcc iron was investigated and the results were compared with the data from neutron irradiation experiment. It's found that the results of the coupling MD-OKMC were consistent with the experiment results, and were better in agreement with the experiment than the traditional OKMC.

An integrative mathematical model of microwave ion thruster for mHz-frequency thrust noise analysis

Many high-precision space missions need thrusters to produce thrust with low noise to compensate for disturbances and ensure satellite platform stability. Microwave ion thruster is characterized with a wide thrust range and potential for these missions. A cost-effective and accurate mathematical model is crucial for mHz-frequency thrust noise analysis and feedback controller design. The Particle-In-Cell (PIC) and global models are two common simulation tools. The PIC model is characterized with high accuracy but huge computation cost, which is difficult to analyze long-time performance characteristics. Now, the global model is only used for the discharge chamber with low accuracy and cannot reflect ion extraction properties. In this paper, an integrative mathematical model is built for a 1-cm microwave ion thruster and can reflect ion beamlet divergence and impingement on the Accelerator Grid (AG). Simulation results show good agreement with experiments at 0.06 sccm. However, the model demonstrates worse consistency with experiments when the flux increases to 0.1 sccm, which may be because the influence of neutral gas on the Electron Cyclotron Resonance (ECR) is not considered in the model. A long-time (1000 s) simulation is conducted with this model under 35 μN. It takes 3 hrs, and the thrust noise reaches 1 μN/Hz0.5 at 1 mHz.

Effects of Dipping and Folded Structure on Gas Production from Hydrate-Bearing Sediments

Summary Due to a dipping and folded structure, hydrate-bearing sediments (HBS) have obvious fluctuation characteristics, and the internal temperature and pressure of HBS are unevenly distributed. Subsequently, gas and water production of natural gas hydrate (NGH) is affected. When using a numerical simulation method to predict effectively the productivity of HBS, it is necessary to establish a conceptual model that considers the formation fluctuation. However, few reported studies accurately describe the fluctuation characteristics of HBS in numerical simulation models. Therefore, the spatial evolution of gas production, water production, and seepage parameters of each model was compared by establishing the initial temperature and pressure model of each representative model pair, using the TOUGH + HYDRATE (T + H) code for a long-time simulation; the production process of gas and water and spatial evolution of seepage parameters of each model were compared; and then the spatial evolution of gas production, water production, and seepage parameters of the different dipping/folded HBS was obtained. The spatial evolution of water production and seepage parameters for different dipping/folded HBS is obtained. (a) The dipping and folded structure had an obvious influence on the spatial distribution of the initial temperature and pressure of HBS. (b) The limits of heat supply and seepage capacity of the fluctuating HBS gave lower gas production than in horizontal HBS. There should be more emphasis on heat supply conditions and the formation of secondary hydrates. (c) The additional pore water in fluctuating HBS is not conducive to the discharge of methane. Consequently, the development of improved water blocking measures is significant for the future large-scale production of NGH.

Solitary Wave Propagation of the Generalized Rosenau–Kawahara–RLW Equation in Shallow Water Theory with Surface Tension

This paper addresses a numerical approach for computing the solitary wave solutions of the generalized Rosenau–Kawahara–RLW model established by coupling the generalized Rosenau–Kawahara and Rosenau–RLW equations. The solution of this model is accomplished by using the finite difference approach and the upwind local radial basis functions-finite difference. Firstly, the PDE is transformed into a nonlinear ODEs system by means of the radial kernels. Secondly, a high-order ODE solver is implemented for discretizing the system of nonlinear ODEs. The main advantage of this technique is its lack of need for linearization. The global collocation techniques impose a significant computational cost, which arises from calculating the dense system of algebraic equations. The proposed technique estimates differential operators on every stencil. As a result, it produces sparse differentiation matrices and reduces the computational burden. Numerical experiments indicate that the method is precise and efficient for long-time simulation.

Redesigning OpenKMC for Multi-Component Trillion-Atom Simulations on the New Sunway Supercomputer

The atomic kinetic Monte Carlo method plays an important role in material simulations by connecting the microscale mechanism with macroscale evolution. However, the long-time simulation of multi-component materials is highly challenging because it demands significant computing resources. With the advent of exascale computing, ultra-high computing power can enable kinetic Monte Carlo (KMC) simulations. In this paper, we deeply optimize OpenKMC for the new-generation Sunway supercomputer. This includes optimizing the memory access for the SW39000 architecture, eliminating various redundant computations at growing scales, and proposing a communication strategy for heterogeneous platforms. In addition, we expanded OpenKMC's simulation for multi-component alloys. Finally, the acceleration framework can produces a <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"><tex-math notation="LaTeX">$37\times$</tex-math></inline-formula> performance enhancement on the Sunway platform. Furthermore, when powered by 10 million cores, our program can perform trillion-atom simulations of complex multi-component alloys with 85% parallel efficiency.

Reducing the low-wavenumber dispersion error by building the Lagrange dual problem with a powerful local restriction

Abstract The broadband finite-difference (FD) coefficients computed by a cost function have been widely applied to suppress of numerical dispersion. Under the same conditions, the FD coefficients with small low-wavenumber dispersion error will produce a more accurate numerical solution in the long-time seismic wave simulation. Thus, how to reduce the low-wavenumber dispersion error becomes a crucial problem. According to the research into the zero-point position at the dispersion curve for three types of cost function, we found that the more zero points concentrate on the low-wavenumber region, the less the dispersion error. Therefore, the concentration of zero points is a good way to reduce dispersion error, which can be implemented by modified wavenumbers of zero points. Then, we design a Lagrange dual problem with a restriction based on the modified wavenumbers. The requirements for constructing the Lagrange dual problem are the optimization function and restricted condition, where the former is based on the dispersion relation, and the latter comprises the modified wavenumbers. The solution of this optimization problem, calculated by the dual ascent method, affords a less low-wavenumber dispersion error than the solution yielded by the alternating direction method of multipliers (ADMM). The theoretical analysis and numerical modeling suggest that the proposed method is superior to the existing optimal FD coefficients in reducing numerical error accumulation in low-frequency simulation.

An Implementation Method of the Complex Frequency-Shifted Uniaxial/Multi-Axial PML Technique for Viscoelastic Seismic Wave Propagation

ABSTRACT A numerical implementation method of the complex frequency-shifted uniaxial/multi-axial PML (CFS UPML/MPML) is developed for modeling half-space viscoelastic wave propagation. A mixed CFS PML combining CFS UPML and CFS MPML is designed to improve PML’s performances in avoiding evanescent waves at grazing incidence and stabilizing long-time simulation of low-frequency viscoelastic wave propagation. Two numerical examples’ results showed that the developed method for the mixed CFS PML can implement elimination of evanescent waves at grazing incidence near the PML interface, long-time accurate simulation of low-frequency viscoelastic wave propagation and efficient absorption of the near-fault outgoing seismic waves in homogeneous earth medium.

Effects of favourable downstream pressure gradients on separated shock-wave/boundary-layer interactions

Shock-wave/boundary-layer interactions (SWBLI) strongly influence the aerodynamic behaviour of many aerospace, propulsion and transportation systems. They have therefore been studied intensively — with different methods, configurations, and test facilities. Regions of favourable pressure gradient are often placed downstream of the SWBLI region., e.g. to prevent wind-tunnel blockage. To be able to better understand their effect and compare cases from varying sources, a joint experimental–numerical study of a supersonic compression/expansion corner flow is carried out. A high-resolution long-time numerical simulation complements particle-image velocimetry measurements, both using the same configuration and conditions. Derived setups without the downstream expansion corner are also explored; both for the same compression-corner angle and for a ramp angle inducing a separation bubble of the same size as for the original compression/expansion configuration. Results demonstrate that the topology and very low-frequency dynamics of the separation bubble developing over the compression corner are altered by the downstream expansion corner. Differences are mainly found for the downstream motion of the shock/shrinking of the separated region and extend downstream of the reattachment point up to the expansion corner.

Linear schemes with several degrees of freedom for the transport equation and the long-time simulation accuracy

Abstract We consider linear schemes with several degrees of freedom for the transport equation on uniform meshes. For these schemes, the solution error is $O(h^p + th^q)$, where $p$ is equal to or greater by one than the order of the truncation error and $q \geqslant p$. We prove the existence of a mapping of smooth functions on the mesh space providing the $q$th order of the truncation error and deviating from the standard mapping ($L_2$-projection for example) by $O(h^p)$. In a one-dimensional case, this mapping can be found in the class of local mappings. In more dimensions, the existence of a local mapping with such properties is guaranteed only under additional assumptions.

High order WENO finite difference scheme with adaptive dual order ideal weights for hyperbolic conservation laws

An adaptive even order weighted essentially non-oscillatory polynomial reconstruction procedure (WPR) with Z-type weights (WENO-Za) is proposed to reduce the numerical dissipation and improve the resolution of the fine-scale structures in a long-time simulation of the hyperbolic conservation laws. The new spatial solver combines the even and odd order (even order plus one) WENO finite difference schemes together. The novel WENO-Za scheme adapts the ideal weights, instead of the nonlinear weights, between the even and odd order based on a nonlinear transition function, which is a function of the odd order global smoothness indicator and a user-adjustable parameter. The dissipation and dispersion behaviors of the WENO-Za scheme is illustrated by the approximate dispersion relation technique. For the Euler equations, the WPR procedure is performed on each positive and negative split flux variable in the characteristic space. The performance of the WENO-Za scheme, in terms of resolution, essentially non-oscillatory shock-capturing, and efficiency, are evaluated by solving several benchmark shock-tube problems, such as the extended shock-entropy interaction problem, the extended shock-density interaction problem, the two blast-waves interaction problem, the Riemann problem, and the forward-facing step problem. The results show that the WENO-Za scheme can 1) retain the fine-scale structures in the smooth regions better due to the reduction of numerical dissipation of the even order scheme in exchange for a slightly larger dispersive error in a long-time simulation and 2) capture localized strong gradients and shocks in an essentially non-oscillatory manner due to the strong sub-stencil biasing in the odd order scheme. The space-time evolutions of the transition functions demonstrate the robust adaptive nature of the WENO-Z4a5 scheme in the regions of shock, contact discontinuity, and rarefaction waves.

Conservative Continuous-Stage Stochastic Runge–Kutta Methods for Stochastic Differential Equations

In this paper, we develop a new class of conservative continuous-stage stochastic Runge–Kutta methods for solving stochastic differential equations with a conserved quantity. The order conditions of the continuous-stage stochastic Runge–Kutta methods are given based on the theory of stochastic B-series and multicolored rooted tree. Sufficient conditions for the continuous-stage stochastic Runge–Kutta methods preserving the conserved quantity of stochastic differential equations are derived in terms of the coefficients. Conservative continuous-stage stochastic Runge–Kutta methods of mean square convergence order 1 for general stochastic differential equations, as well as conservative continuous-stage stochastic Runge–Kutta methods of high order for single integrand stochastic differential equations, are constructed. Numerical experiments are performed to verify the conservative property and the accuracy of the proposed methods in the longtime simulation.

Numerical Study on Viscous Fingering Using Electric Fields in a Hele-Shaw Cell

We investigate the nonlinear dynamics of a moving interface in a Hele-Shaw cell subject to an in-plane applied electric field. We develop a spectrally accurate boundary integral method where a coupled integral equation system is formulated. Although the stiffness due to the high order spatial derivatives can be removed, the long-time simulation is still expensive since the evolving velocity of the interface drops dramatically as the interface expands. We remove this physically imposed stiffness by employing a rescaling scheme, which accelerates the slow dynamics and reduces the computational cost. Our nonlinear results reveal that positive currents restrain finger ramification and promote overall stabilization of patterns. On the other hand, negative currents make the interface more unstable and lead to the formation of thin tail structures connecting the fingers and a small inner region. When no flux is injected, and a negative current is utilized, the interface tends to approach the origin and break up into several drops. We investigate the temporal evolution of the smallest distance between the interface and the origin and find that it obeys an algebraic law $\displaystyle (t_*-t)^b$, where $t_*$ is the estimated pinch-off time.

Influence of Convection Term on Temperature Field during Soil Freezing

In recent years, the numerical model of frozen soil has mainly focused on the water-heat equation of unsaturated soil, which is of great significance to predicting engineering stability in frozen soil areas. In these numerical models, the change of freezing temperature is usually ignored, and 0 °C is often used as the freezing condition. In addition, most equations only consider the effect of latent heat released during water freezing on the frozen soil and ignore the effect of high-temperature water on the heat transfer of frozen soil when unfrozen water migrates to the frozen zone. Therefore, there will be deviations under long-time simulation. At the same time, due to excessive attention to the moisture field and neglect of the selection of the temperature field, there is no clear conclusion on when to choose the heat transfer equation with the convection term. The equation in this paper considered the change of freezing temperature under different initial saturation conditions and the convection effect during moisture migration. Through COMSOL 5.5 software simulation, they were combined with the experiment to verify. Three different hydrothermal equations were selected to analyze the effects of latent heat and the effects of convection on the temperature field under different initial saturation conditions. The results show that the convection term plays an essential role in the heat transfer equation for unsaturated soil with high initial saturation. Additionally, the frost heave occurs mainly above the ice front interface. This study provides a reference for when to choose the heat transfer equation with convective terms and can provide help for the construction and prediction of frozen soil in the future.

A new particle shifting technique for SPH methods based on Voronoi diagram and volume compensation

Particle shifting technique (PST) in smoothed particle hydrodynamics (SPH) is a useful strategy to render particle distribution isotropic for higher numerical accuracy and stability. However, the non-isotropic particle distribution near the free surfaces and volume non-conservation still often occur, which will degenerate accuracy in the long-time simulation with violent flows. In this study, the two major concerns are solved by proposing a new PST. The Voronoi cell is constructed for the particles with incomplete kernel support, and the difference between the Voronoi cell centroid and the corresponding particle location is used to determine the shifting vector, which ensures a very isotropic particle distribution and that no clustering or nonphysical gap occurs in these areas. For the free-surface particles, the shifting vector is determined according to the volume change and thus the volume conservation is ensured, where the volume is calculated conveniently with the merit of Voronoi diagram. The present method is compared with the other two advanced PST methods, and a set of challenging cases with violent flows is simulated to validate its superior properties in terms of maintaining isotropic particle distribution and volume conservation. Moreover, the efficiency of the present method is almost the same as those of other advanced PSTs. Overall, the new PST method presents a strong potential in the long-time simulation of violent flows in coastal engineering.