Third-order Runge-Kutta Method Research Articles

LAPS: An MPI-parallelized 3D pseudo-spectral Hall-MHD simulation code incorporating the expanding box model

Numerical simulations have been an increasingly important tool in space physics. Here, we introduce an open-source three-dimensional compressible Hall-Magnetohydrodynamic (MHD) simulation code LAPS (UCLA-Pseudo-Spectral, https://github.com/chenshihelio/LAPS). The code adopts a pseudo-spectral method based on Fourier Transform to evaluate spatial derivatives, and third-order explicit Runge-Kutta method for time advancement. It is parallelized using Message-Passing-Interface (MPI) with a “pencil” parallelization strategy and has very high scalability. The Expanding-Box-Model is implemented to incorporate spherical expansion effects of the solar wind. We carry out test simulations based on four classic (Hall)-MHD processes, namely, 1) incompressible Hall-MHD waves, 2) incompressible tearing mode instability, 3) Orszag-Tang vortex, and 4) parametric decay instability. The test results agree perfectly with theory predictions and results of previous studies. Given all its features, LAPS is a powerful tool for large-scale simulations of solar wind turbulence as well as other MHD and Hall-MHD processes happening in space.

Numerical study of shock wave boundary layer interaction flows using a novel turbulence model

This study proposes a novel partial average fluctuation velocity (PAFV) turbulence based on the PAFV and the average strain rate tensor, which is used to calculate the shock wave boundary layer interaction flows, including cascade flow and transonic compressor flow. In order to illustrate the adequacy, the calculation results of several other turbulence models were also compared with PAFV, including k-ω, shear stress transport, SA-neg (negative Spalart-Allmaras), and stress-BSL (ω-based Reynolds stress model). The suitability of the PAFV turbulence model was initially tested by calculating the transonic flow in the L030-4 cascade. Subsequently, numerical computations were performed for the 3D (three-dimensional) National Aeronautics and Space Administration Rotor 67 transonic compressor rotor. In the numerical simulation, the flow control equations were transformed in a general curvilinear coordinate system for precision enhancement, and the spatial discretization of control equations was executed using the finite difference method. The convection term was processed using the Steger–Warming flux vector splitting method, the diffusion term was discretized using a second-order central difference scheme, and the time derivative term was discretized using the third-order Runge–Kutta method. A parallel computing strategy with multi-block partitioning was employed to speed up the calculations and facilitate faster convergence alongside a local time step method. All numerical computation procedures were written in Python, revealing the following: (1) The PAFV turbulence model aptly simulates the transonic flow within the cascade channel. The shock wave pattern in the cascade channel aligns well with the experimental schlieren images. (2) During the Rotor 67 transonic compressor rotor calculation, the flow–pressure ratio, flow efficiency, and inlet and outlet total temperatures and pressures align well with the experimental data. Furthermore, the resolution of the flow field structure related to the interaction between the tip leakage vortex and channel shock wave is high. (3) The turbulence model uses PAFV as the turbulence velocity scale, thereby suppressing the fluctuation velocity and turbulence viscosity near shock wave areas, preventing excessive turbulence viscosity. Moreover, PAFV inherently exhibits transport properties and anisotropic characteristics, making it well-suited for handling transonic compressor flows.

Numerical simulation and fast method for the 0D-1D multi-scale coupled model and its application in ischemic brain tissue blood flow problems.

As living standards rise, more and more people are paying attention to their own health, especially issues such as cerebral thrombosis, cerebral infarction, and other cerebral blood flow problems. An accurate simulation of blood flow within cerebral vessels has emerged as a crucial area of research. In this study, we focus on microcirculatory blood flow in ischemic brain tissue and employ a 0D-1D geometric multi-scale coupled model to characterize this process. Given the intricate nature of human cerebral vessels, we apply a numerical method combining the finite element method and the third-order Runge-Kutta method to resolve the coupled model. To enhance computational efficiency, we introduce a fast method based on the reduced-order extrapolation algorithm. Our numerical example underscores the stability of the method and convergence accuracy to , while significantly improving the accuracy and efficiency of blood flow simulation, making the mechanism analysis more accurate. Additionally, we present examples detailing variations and distribution of intracranial pressure and blood flow in ischemic brain tissue throughout a cardiac cycle. Both reduced vascular compliance and vascular stenosis can have adverse effects on intracranial cerebral pressure and blood flow, leading to insufficient local oxygen supply and negative effects on brain function. Meanwhile, there will also be corresponding changes in volume flow and pulsatile blood pressure.

Iterated Crank–Nicolson Method for Peridynamic Models

In this paper, we explore the iterated Crank–Nicolson (ICN) algorithm for the one-dimensional peridynamic model. The peridynamic equation of motion is an integro-differential equation that governs structural deformations such as fractures. The ICN method was originally developed for hyperbolic advection equations. In peridynamics, we apply the ICN algorithm for temporal discretization and the midpoint quadrature method for spatial integration. Several numerical tests are carried out to evaluate the performance of the ICN method. In general, the ICN method demonstrates second-order accuracy, consistent with the Störmer–Verlet (SV) method. When the weight is 1/3, the ICN method behaves as a third-order Runge–Kutta method and maintains strong stability-preserving (SSP) properties for linear problems. Regarding energy conservation, the ICN algorithm maintains at least second-order accuracy, making it superior to the SV method, which converges linearly. Furthermore, selecting a weight of 0.25 results in fourth-order superconvergent energy variation for the ICN method. In this case, the ICN method exhibits energy variation similar to that of the fourth-order Runge–Kutta method but operates approximately 20% faster. Higher-order convergence for energy can also be achieved by increasing the number of iterations in the ICN method.

High-order accurate well-balanced energy stable adaptive moving mesh finite difference schemes for the shallow water equations with non-flat bottom topography

This paper proposes high-order accurate well-balanced (WB) energy stable (ES) adaptive moving mesh finite difference schemes for the shallow water equations (SWEs) with non-flat bottom topography. To enable the construction of the ES schemes on moving meshes, a reformulation of the SWEs is introduced, with the bottom topography as an additional conservative variable that evolves in time. The corresponding energy inequality is derived based on a modified energy function, then the reformulated SWEs and energy inequality are transformed into curvilinear coordinates. A two-point energy conservative (EC) flux is constructed, and high-order EC schemes based on such a flux are proved to be WB that they preserve the lake at rest. Then high-order ES schemes are derived by adding suitable dissipation terms to the EC schemes, which are newly designed to maintain the WB and ES properties simultaneously. The adaptive moving mesh strategy is performed by iteratively solving the Euler-Lagrangian equations of a mesh adaptation functional. The fully-discrete schemes are obtained by using the explicit strong-stability preserving third-order Runge-Kutta method. Several numerical tests are conducted to validate the accuracy, WB and ES properties, shock-capturing ability, and high efficiency of the schemes.

Acoustic characteristics of supersonic planar impinging jets

The present work aims to study the tonal and broadband noise associated with a supersonic planar jet, of Mach number 1.28 and Reynolds number 50000, impinging on a flat plate. Five different cases are considered corresponding to different plate distances and angles of impingement. The near-field noise is analyzed using Large Eddy Simulation (LES). For this, we employ a low dispersive and dissipative Finite Volume solver using AUSM for inviscid fluxes and third-order Runge–Kutta method for temporal discretization. The far-field noise is calculated using a hybrid approach in which the acoustic sources and mean flow are provided by the fluid solver while the far-field acoustic pressure is calculated by an acoustic solver. It solves Acoustic Perturbation Equations using a spectral/hp element method with a Discontinuous Galerkin projection. The present work highlights many important acoustics studies of impinging jets. Firstly, the far-field sound with the effect of non-uniform base fields is obtained which is important as the far-field acoustic pressure in most of the literature is reported either without any base fields or by the application of uniform base fields. Secondly, a correlation between the different modes of oscillations of the impinging jet and the production of impinging tones is studied using Modal Decomposition techniques. This analysis also sheds some light on the number of concurrent cycles for the different nozzle-to-plate distances. Lastly, the relation between the phase lag term and the distance between the plate and the nozzle exit is highlighted using a vortex tracking method. These observations can serve to provide a better understanding of the sound generation process.

Leveraging Google's Tensor Processing Units for tsunami-risk mitigation planning in the Pacific Northwest and beyond

Abstract. Tsunami-risk mitigation planning has particular importance for communities like those of the Pacific Northwest, where coastlines are extremely dynamic and a seismically active subduction zone looms large. The challenge does not stop here for risk managers: mitigation options have multiplied since communities have realized the viability and benefits of nature-based solutions. To identify suitable mitigation options for their community, risk managers need the ability to rapidly evaluate several different options through fast and accessible tsunami models, but they may lack high-performance computing infrastructure. The goal of this work is to leverage Google's Tensor Processing Unit (TPU), a high-performance hardware device accessible via the Google Cloud framework, to enable the rapid evaluation of different tsunami-risk mitigation strategies available to all communities. We establish a starting point through a numerical solver of the nonlinear shallow-water equations that uses a fifth-order weighted essentially non-oscillatory method with the Lax–Friedrichs flux splitting and a total variation diminishing third-order Runge–Kutta method for time discretization. We verify numerical solutions through several analytical solutions and benchmarks, reproduce several findings about one particular tsunami-risk mitigation strategy, and model tsunami runup at Crescent City, California whose topography comes from a high-resolution digital elevation model. The direct measurements of the simulation's performance, energy usage, and ease of execution show that our code could be a first step towards a community-based, user-friendly virtual laboratory that can be run by a minimally trained user on the cloud thanks to the ease of use of the Google Cloud platform.

Research on Aerodynamic Characteristics of Crescent Iced Conductor Based on S-A Finite Element Turbulence Model

Galloping is a typical wind-induced phenomenon in iced conductors, which can have serious impacts on the safe and stable operation of power systems. The aerodynamic characteristics of an iced conductor are the key factor in the study of galloping, which can be determined mainly by the numerical simulation of flow past an iced conductor. Based on the Reynolds-averaged Navier-Stokes (RANS) equations closed by the Spalart-Allmaras (S-A) turbulence model, the third-order Runge-Kutta method along the uniform streamline and Galerkin method are used for temporal and spatial discretization, respectively. The convection and diffusion terms in the discretization scheme are treated semi-implicitly, and the finite element scheme based on the S-A turbulence model is presented and used to numerically simulate flow past a crescent iced conductor. We systematically investigate the effects of icing thickness, wind speed, and wind attack angle on aerodynamic coefficients and flow patterns. Based on the experimental results, the effectiveness of the present algorithm is verified. Using the streamline diagram and pressure distribution diagram of the crescent-shaped iced conductor, the mechanism for the sharp peak of the lift coefficient is explored. Combined with the galloping mechanism of Den Hartog and Nigol, the galloping instability zone of the crescent-shaped iced conductor is analyzed.

ANALYSIS OF THE SPRUCE BUDWORM MODEL USING THE HEUN METHOD AND THIRD-ORDER RUNGE-KUTTA

This study discusses the analysis of the Spruce Budworm model using numerical methods, namely the Heun method and the Third Order Runge-Kutta method. The purpose of this study is to determine the numerical results of the Heun method and the Third Order Runge-Kutta method on the cypress caterpillar model and to determine the comparison of errors from the two methods, namely the Heun method and the Third Order Runge-Kutta method in analyzing the Spruce Budworm model. The results of the study using the Heun method for the initial conditions at years, for , the result obtained is and . For the calculation result of the Spruce Budworm model using the third-order Runge-Kutta method, the result obtained is and . while the error of Caterpillar Density with the third-order Runge-Kutta method is bigger than the Heun method, while the error of Branch Surface Area and the error of Reserve Food error using the Heun method bigger than using the third-order Runge-Kutta method.

Toward the Optimal Design and FPGA Implementation of Spiking Neural Networks.

The performance of a biologically plausible spiking neural network (SNN) largely depends on the model parameters and neural dynamics. This article proposes a parameter optimization scheme for improving the performance of a biologically plausible SNN and a parallel on-field-programmable gate array (FPGA) online learning neuromorphic platform for the digital implementation based on two numerical methods, namely, the Euler and third-order Runge-Kutta (RK3) methods. The optimization scheme explores the impact of biological time constants on information transmission in the SNN and improves the convergence rate of the SNN on digit recognition with a suitable choice of the time constants. The parallel digital implementation leads to a significant speedup over software simulation on a general-purpose CPU. The parallel implementation with the Euler method enables around 180× ( 20× ) training (inference) speedup over a Pytorch-based SNN simulation on CPU. Moreover, compared with previous work, our parallel implementation shows more than 300× ( 240× ) improvement on speed and 180× ( 250× ) reduction in energy consumption for training (inference). In addition, due to the high-order accuracy, the RK3 method is demonstrated to gain 2× training speedup over the Euler method, which makes it suitable for online training in real-time applications.

Third-order Paired Explicit Runge-Kutta schemes for stiff systems of equations

The ability to advance locally-stiff systems of equations in time depends on accurate and efficient temporal schemes. Recently, a new family of Paired Explicit Runge-Kutta (P-ERK) methods has been proposed. This approach allows different Runge-Kutta schemes with different numbers of active stages to be assigned based on local stiffness criteria. Whereas the original P-ERK formulation was only second-order accurate, in this paper we propose a new third-order family. We present the general formulation for these schemes, and then optimize them for the high-order discontinuous Galerkin method recovered via the flux reconstruction approach. We then verify that these schemes obtain their designed third-order accuracy for non-linear systems of equations. Performance results show that they achieve speedup factors up to four relative to a classical third-order Runge-Kutta methods for laminar and turbulent flow over an SD7003 airfoil, and turbulent flow over a tandem sphere configuration. Based on these results, this new family of third-order P-ERK schemes provides an appealing approach for accurate and efficient solution of locally-stiff systems of equations.

Effects of Divergence Inlet on Kerosene/Air Rotating Detonation Engines

In this paper, an improved three-dimensional conservation element and solution element method with chemical source terms is developed and applied to investigate the effects of different divergence inlet patterns on vapor kerosene/air rotating detonation engines. A two-step kerosene/air chemistry model is adopted and solved by a third-order Runge–Kutta method. A series of cases with varying the slot inlet width, in other words, 6, 4, and 2 mm, and divergence section length, in other words, 4, 6, and 10 mm, are simulated. The simulation results indicate the supersonic injection generated by the divergence section is helpful to decrease the disturbance from the downstream flow. Affected by the divergence section, the detonation front is divided into two sections, in other words, the main detonation front and the banded detonation front. For the main detonation front, the detonation strength on the outer wall is significantly stronger than that on the inner wall, while it is not for the banded detonation front. The total pressure loss caused by the inlet normal shock significantly increases with the decrease of slot inlet width. A declining deflagration surface will appear when a small slot inlet width of 2 mm is used. The rotating detonation engine (RDE) performance analysis shows that the slot inlet width has a significant impact on the thrust, specific impulse, total pressure gain, and outlet Mach number. A bigger slot inlet area ratio contributes to a higher total pressure gain. The simulation results can get insight into rotating detonation behaviors and guide the development of the kerosene/air RDE.

Formulative Visualization of Numerical Methods for Solving Non-Linear Ordinary Differential Equations

Many physical problems in the real world are frequently modeled by ordinary differential equations (ODEs). Real-life problems are usually non-linear, numerical methods are therefore needed to approximate their solution. We consider different numerical methods viz., Explicit (Forward) and Implicit (Backward) Euler method, Classical second-order Runge-Kutta (RK2) method (Heun’s method or Improved Euler method), Third-order Runge-Kutta (RK3) method, Fourth-order Runge-Kutta (RK4) method, and Butcher fifth-order Runge-Kutta (BRK5) method which are popular classical iteration methods of approximating solutions of ODEs. Moreover, an intuitive explanation of those methods is also be presented, comparing among them and also with exact solutions with necessary visualizations. Finally, we analyze the error and accuracy of these methods with the help of suitable mathematical programming software.

The modified high-order Haar wavelet scheme with Runge–Kutta method in the generalized Burgers–Fisher equation and the generalized Burgers–Huxley equation

In this work, we focus on the modified high-order Haar wavelet numerical method, which introduces the third-order Runge–Kutta method in the time layer to improve the original numerical format. We apply the above scheme to two types of strong nonlinear solitary wave differential equations named as the generalized Burgers–Fisher equation and the generalized Burgers–Huxley equation. Numerical experiments verify the correctness of the scheme, which improves the speed of convergence while ensuring stability. We also compare the CPU time, and conclude that our scheme has high efficiency. Compared with the traditional wavelets method, the numerical results reflect the superiority of our format.

A third-order weighted essentially non-oscillatory scheme in optimal control problems governed by nonlinear hyperbolic conservation laws

The weighted essentially non-oscillatory (WENO) methods are popular and effective spatial discretization methods for nonlinear hyperbolic partial differential equations. Although these methods are formally first-order accurate when a shock is present, they still have uniform high-order accuracy right up to the shock location. In this paper, we propose a novel third-order numerical method for solving optimal control problems subject to scalar nonlinear hyperbolic conservation laws. It is based on the first-disretize-then-optimize approach and combines a discrete adjoint WENO scheme of third order with the classical strong stability preserving three-stage third-order Runge–Kutta method SSPRK3. We analyze its approximation properties and apply it to optimal control problems of tracking-type with non-smooth target states. Comparisons to common first-order methods such as the Lax–Friedrichs and Engquist–Osher method show its great potential to achieve a higher accuracy along with good resolution around discontinuities.

Application of Numerical Methods for the Analysis of Damped Parallel RLC Circuit

A sudden application of sources results in time-varying currents and voltages in the circuit known as transients. This phenomenon occurs frequently during switching. A simple circuit constituting a resistor, an inductor, and a capacitor is termed an RLC circuit. It may be in parallel or series configuration or both. Different values of damping factors determine the different nature of the transient response. We applied different numerical solution methods such as explicit (forward) Euler method, third-order Runge-Kutta (RK3) method, and Butcher's fifth-order Runge-Kutta (BRK5) method to approximate the solution of second-order differential equation with initial value problem (IVP). We thoroughly compared the numerical solutions obtained by these methods with the necessary visualization and analysis of error. We also examined the superiority of these methods over one another and the appropriateness of numerical methods for different damping conditions is explored. With high accuracy of the approximation and thorough analysis of the observation, we found Butcher's fifth-order Runge-Kutta (BRK5) method to be the best numerical technique. Regarding the different values of damping factors, we considered the further possibility of discussion and analysis of this iterative method.

Matrix-free TriGlobal adjoint stability analysis of compressible Navier–Stokes equations

A numerical method for TriGlobal (i.e., fully three-dimensional) adjoint stability analysis for compressible flows was developed and is presented in this paper. The developed method solves the adjoint stability problem using a matrix-free method based on Krylov-Schur method and a time-stepping approach. Because of the low memory (RAM) requirement of the matrix-free approach, the developed method can analyze fully three-dimensional flows that are difficult to analyze with conventional matrix-based methods. To perform time-stepping on the adjoint variables, the adjoint equations, including appropriate boundary conditions for the compressible Navier–Stokes equations, were derived. The equations were discretized using a finite compact difference method, and time integration was conducted using the three-step third-order Runge–Kutta method. A flow field around a two-dimensional square cylinder and a cubic cavity flow were analyzed using the developed method, and it was confirmed that the method reproduces the dominant instabilities reported in the literature. In addition, in the square cylinder flow analysis, the receptivity and sensitivity regions of the secondary wake mode, which corresponds to far wake instability, were clarified. Finally, the TriGlobal direct and adjoint stabilities of compressible flows over a finite width cavity were analyzed for the first time, and it was proven that large-scale adjoint stability analysis can be performed with the developed method. The results also show that instability phenomena similar to those obtained with BiGlobal stability analysis appear, but sidewall effects exist.

Research on GPU parallel algorithm for direct numerical solution of two-dimensional compressible flows

In this paper, a novel parallel algorithm is proposed to solve the problems of heavy computation and long simulation time in the field of compressible flows. In this algorithm, a third-order upwind scheme and a fourth-order central difference scheme are employed, with a third-order Runge-Kutta method for time stepping. Considering the powerful floating-point computing ability of the Graphics Processing Unit (GPU), this paper establishes the algorithm on the basis of GPU. Moreover, the direct numerical simulation method is adopted in this algorithm to improve the solution accuracy of the simulation results. To further enhance the efficiency of the algorithm, several optimization strategies are explored in the design of the algorithm as well. Both accuracy and feasibility of the algorithm are verified by a classical two-dimensional example. Compared with solving this example on the Central Processing Unit platform, the experimental results demonstrate that the maximum speedup ratio achieved by our approach is 18.03 times.

Sharp-Interface Immersed-Boundary Method for Compressible Flows with Shock–Particle Interaction

A sharp-interface immersed-boundary method that does not smear the boundaries, is easy to implement, and is straightforward in handling the Neumann condition as compared to previous methods is developed for compressible flows to simulate the interaction of solid particles with shocks. The inviscid and viscous fluxes of compressible flow equations in curvilinear coordinates are discretized with a third-order weighted essentially nonoscillatory (WENO) scheme and a central scheme, respectively. The equations are advanced in time using a third-order Runge–Kutta method. The sharp interface at the immersed boundaries is maintained by reconstructing the flow variables along the normal direction to the boundary. The WENO discretization is reverted to a biased essentially nonoscillatory scheme near the immersed boundaries to avoid using the nodes that are inside the immersed boundary. The method is validated against experimental measurements and shown to be between second- and third-order-accurate in the presence of immersed boundaries. The method is applied to simulate shock-tube experiments involving the interaction of a moving normal shock with a stationary cylinder as well as a cylindrical and a spherical particle accelerating by a shock. The numerical results capture all of the shock features observed in the experiments and show great agreement with the measurements and previous benchmark solutions. The results show that the acceleration of a sphere due to an incident shock highly depends on the density ratio of the sphere to the incoming fluid.

Evaluation of breakup models for liquid side jets in supersonic cross flows

In this study, numerical simulations of the interaction between supersonic crossflow and a side liquid jet were performed. The processes of the liquid jet breakup, flow separation and shock waves resulting from the interaction of the gas flows and liquid droplets were predicted. In the numerical simulation, a two-way coupling Eulerian-Lagrangian method was used. The two-dimensional unsteady compressible Navier-Stokes Equations were solved based on the MUSCL type (Monotone Upstream-Centred Scheme for Conservation Laws) spatial differencing coupling with Roe-type flux-difference splitting and a third-order accurate Runge-Kutta method for time-splitting. The atomization of liquid droplets was simulated using the combined KH-RT breakup model. The modified breakup model considered both the interface instability and the sudden deformation of a liquid spray induced by the compressible flows. Shock waves, recirculation zones and breakup processes of droplet particles were captured. The predicted droplet distributions at various operating conditions were discussed. Satisfied agreement between simulations and measurements is shown in the current work.