Low-storage Runge-Kutta Research Articles

A split-explicit second order Runge–Kutta method for solving 3D hydrodynamic equations

Numerical models of marine hydrodynamics have to deal with processes exhibiting a wide range of timescales. These processes include fast external gravity waves and slower internal fully three-dimensional motions. In order to be both time-efficient and numerically stable, the time stepping scheme has to be chosen carefully to cope with the characteristic time scale of each phenomenon. An usual approach is to split the fast and slow dynamics into separate modes. The fast waves are modeled with a two-dimensional system through depth averaging while the other motions, where characteristic times are much longer, are dealt with in three dimensions. However, if the splitting is inexact, for instance in projecting the fields in a new 3D mesh, this procedure can lead to improper results regarding to the physical properties such as mass conservation and tracer consistency. In this work, a new split-explicit Runge–Kutta scheme is adapted and developed for the Discontinuous-Galerkin Finite Element method in order to obtain a new second-order time stepping, yielding to more accurate results. This method combines a three-stage low-storage Runge–Kutta for the slow processes and a two-stage low-storage for the fast ones. The 3D iterations are not affecting the surface elevation, hence an Arbitrary Lagrangian Eulerian implementation is straightforward. Water volume and tracers are conserved. A set of test cases for baroclinic flows as well as a laboratory application demonstrate the performance of the scheme. They suggest that the new scheme has little numerical diffusion.

Implementation of efficient low-storage techniques for 3-D seismic simulation using the curved grid finite-difference method

SUMMARYHigh-resolution 3-D seismic simulation imposes severe demands for computational memory, making low-storage seismic simulation particularly important. Due to its high-efficiency and low-storage, the half-precision floating-point 16-bit format (FP16) is widely used in heterogeneous computing platforms, such as Sunway series supercomputers and graphics processing unit (GPU) computing platforms. Furthermore, the low-storage Runge–Kutta (LSRK) technique requires lower memory resources compared with the classical Runge–Kutta. Therefore, FP16 and LSRK provide the possibility for low-storage seismic simulation. However, the orders of magnitude of the physical quantities (velocity, stress and Lamé constants) in the elastic wave equations are influenced by the P-wave and S-wave velocities and the densities of the elastic media. This results in a huge order of magnitude difference between the stored velocity and stress values, which exceed the range of the stored values of FP16. In this paper, we introduce three dimensionless constants, Cv, Cs and Cp, into elastic wave equations, and new elastic wave equations are derived. The three constants, Cv, Cs and Cp, keep the orders of magnitude of the velocity and stress at a similar level in the new elastic wave equations. Thus, the stored values of these variables in new equations remain within the range of the stored values of FP16. In addition, we introduce the use of the LSRK due to its low-storage characteristic. In this paper, based on the FP16 and LSRK low-storage techniques, we develop 3 optimized multi-GPU solvers for seismic simulation using the curved grid finite-difference method (CGFDM). Moreover, we perform a series of seismic simulations to verify the accuracy, stability, and validity of the optimized solver coupled with the two techniques. The verifications indicate that through maintaining the calculation accuracy, the computational efficiency of the solver is significantly optimized, and the memory usage is remarkably reduced. In particular, under the best conditions, the memory usage can be reduced to nearly 1/3 that of the original CGFDM solver.

Implementation of low-storage Runge-Kutta time integration schemes in scalable asynchronous partial differential equation solvers

The asynchronous computing method based on finite-difference schemes has shown promise in significantly improving the scalability of time-dependent partial differential equation (PDE) solvers by either relaxing data synchronization or avoiding communication between processing elements (PEs) on massively parallel machines. This method uses high-order accurate asynchrony-tolerant (AT) schemes coupled with appropriate time integration schemes to provide the desired accuracy required by high-fidelity solvers such as the direct numerical simulations for fluid flow. For time integration, Runge-Kutta (RK) schemes, particularly the low-storage implementation, are widely used due to their ability to provide good stability properties and be computationally efficient. However, the implementation of AT schemes with multi-stage RK schemes necessitates an over-decomposition of the spatial domain in a parallel setting, leading to increased message sizes for communication and redundant computations. In this paper, we propose a novel method to couple asynchrony-tolerant and low-storage explicit RK schemes in solving time-dependent PDEs that would result in a significant reduction in communications and relaxed synchronizations. We develop new asynchrony-tolerant schemes for ghost or buffer point updates that are necessary to maintain desired order of accuracy. The accuracy of this method is investigated, both theoretically and numerically, using simple one-dimensional linear model equations. Thereafter, we demonstrate the scalability of the proposed numerical method through three-dimensional simulations of decaying Burgers' turbulence, performed using two different asynchronous algorithms: communication-avoiding and synchronization-avoiding algorithms. The scalability studies up to 27,000 cores were found to yield speed-ups up to 6× compared to a baseline synchronous algorithm. Overall, the proposed approach shows the potential to improve the scalability of exascale PDE solvers significantly.

Development and verification of a high-speed compressible reactive flow solver in OpenFOAM

This paper proposes a newly developed density-based compressible solver called DCRFoam on a basis of OpenFOAM platform for simulating detonation phenomena with detailed chemical kinetics. It employs AUSM+-up scheme to determine the convective fluxes and a four steps low storage Runge-Kutta time integration method for time integration. A total variation diminishing (TVD) and symmetrical flux limiter sgva is implemented for interpolation process. Also, the dynamic mesh adaptation method is applied to refine or coarsen the computational mesh dynamically since the resolution requirements of supersonic combustion are localized. The numerical solver has been validated systematically with many numerical and experimental results reported in the literature, including those from non-reacting cases (Sod’s shock tube problem, two-dimensional shock/cylinder interaction) to reacting flow cases (the ignition sequence in a shock tube, oblique detonation wave (ODW) initiating process and detonation propagating in a bifurcated channel). The results obtained by current schemes are compared with original central upwind Kurganov-Noelle-Petrova (KNP) scheme in official OpenFOAM. The corresponding result demonstrates the robustness and accuracy of the current solver.

Explicit high-order exponential time integrator for discontinuous Galerkin solution of Maxwell's equations

In this paper, exponential propagator (EP) is introduced to the discontinuous Galerkin finite element method (DG-FEM) for solving transient electromagnetic problems. As a kind of explicit time-stepping method, the proposed EP has fourth order accuracy in time and requires only one storage unit for each of the degrees-of-freedom. It provides an attractive alternative to the high-order low-storage explicit Runge-Kutta (LSERK) approach, which is the most popular time-marching scheme currently used in high-fidelity simulation of Maxwell's equations with DG-FEM. When centered numerical flux is utilized, the discretized Maxwell's system in lossless region is symplectic structure-preserving (energy conserving), and thus very suitable for accurately computing the long-time propagation of waves in non-dissipative medium, or simulating the motion of charged particles in electromagnetic fields. Our method can also handle upwind numerical flux at the sacrifice of symplectic-conserving virtue for higher accuracy in spatial discretization and spurious inhibition. A major contribution of this work is that, the application range of this scheme is extended to lossy media such as conductors and particularly unsplit perfectly matched layers in the framework of DG-FEM. This greatly facilitates the treatment of open boundary conditions in scattering and/or radiation problems. Extensive numerical examples, including the analysis of photonic crystals, dielectric coated perfect conducting cylinder, propagation of waves in zero-index materials, and corrugated O-type Cherenkov oscillators etc., are presented to demonstrate its stability, accuracy, and performance. All the numerical experiments show that the presented method is superior to the LSERK scheme when centered numerical flux is adopted.

OpenSBLI: Automated code-generation for heterogeneous computing architectures applied to compressible fluid dynamics on structured grids

OpenSBLI is an open-source code-generation system for compressible fluid dynamics (CFD) on heterogeneous computing architectures. Written in Python, OpenSBLI is an explicit high-order finite-difference solver on structured curvilinear meshes. Shock-capturing is performed by a choice of high-order Weighted Essentially Non-Oscillatory (WENO) or Targeted Essentially Non-Oscillatory (TENO) schemes. OpenSBLI generates a complete CFD solver in the Oxford Parallel Structured (OPS) domain specific language. The OPS library is embedded in C code, enabling massively-parallel execution of the code on a variety of high-performance-computing architectures, including GPUs. The present paper presents a code base that has been completely rewritten from the earlier proof of concept Jacobs et al. (2017) [7], allowing shock capturing, coordinate transformations for complex geometries, and a wide range of boundary conditions, including solid walls with and without heat transfer. A suite of validation and verification cases are presented, plus demonstration of a large-scale Direct Numerical Simulation (DNS) of a transitional Shockwave Boundary Layer Interaction (SBLI). The code is shown to have good weak and strong scaling on multi-GPU clusters. We demonstrate that code-generation and domain specific languages are suitable for performing efficient large-scale simulations of complex fluid flows on emerging computing architectures. Program summaryProgram Title: OpenSBLI code-generation framework for compressible fluid dynamics on heterogeneous architecturesCPC Library link to program files:https://doi.org/10.17632/3sdb6hck2c.1Developer's repository link:https://github.com/opensbli/opensbliLicensing provisions: GPLv3Programming languages: Python, C/C++, OPS DSLNature of problem: The compressible 3D Navier-Stokes equations are solved via Implicit Large Eddy Simulation (ILES) or Direct Numerical Simulation (DNS).Solution method: OpenSBLI [1,2] is a Python-based code-generation system that uses symbolic algebra to generate a complete CFD solver in C/C++. The basic algorithm is a stencil-based finite-difference solver on structured curvilinear meshes. Shock-capturing is performed by a selection of high-order Weighted/Targeted Essentially Non-Oscillatory (WENO/TENO) schemes. Explicit low-storage Runge-Kutta schemes are used for time-advancement.Additional comments including restrictions and unusual features: The generated code is compliant with the Oxford Parallel Structured (OPS) [3] software library. OpenSBLI/OPS executables can be generated for the OpenMP, MPI, CUDA, OpenCL, and OpenACC parallel programming paradigms. Multi-GPU support is available via combinations of MPI with CUDA, OpenCL or OpenACC.

Efficient two-step Runge-Kutta methods for fluid dynamics simulations

Explicit two-step Runge-Kutta (TSRK) methods offer an efficient alternative to traditional explicit Low-Storage Runge-Kutta (LSRK) schemes for solving the Navier-Stokes equations. A special class of TSRK methods that reduce requirement compared to previous TSRK schemes are derived. Schemes of fourth, fifth and sixth order are implemented and tested.The new schemes are evaluated with two common test cases, a 2D cylinder and a 3D Taylor-Green vortex. The results are compared with classical time discretization strategies. Timings obtained in three different hardware configurations show that the new TSRK methods of order four are 25% faster than LSRK schemes of the same order.Fifth and sixth order TSRK methods are tested with the same 3D test case and the results are compared to LSRK algorithms. Results show TSRK schemes of fifth and sixth order are competitive compared to LSRK methods of the same orders, as LSRK methods are of second order for non linear differential systems.

Explicit two‐step Runge‐Kutta methods for computational fluid dynamics solvers

SummaryComputational fluid dynamics (CFD) has emerged as a successful tool for industry applications and basic science during the last decades. However, accurate solutions involving vortex propagation, and separated and turbulent flows, are still associated with high computing costs. In particular, large eddy simulations (LES) of complex geometries, such as a complete automobile, require several days on thousands of cores in order to obtain solutions with statistically relevant information. With an increase in the number of available cores, the number of degrees of freedom (DOF) per core can be reduced accordingly. When the number of DOF per core is below a certain threshold the total simulation time is not bounded by floating point operations (FLOPS), but by the time spend on communication between cores. To overcome this impediment we have identified and tested a class of two‐step Runge‐Kutta (TSRK) methods of high order with low number of stages, for time discretization of differential systems resulting from space discretization of weakly compressible Navier‐Stokes equations. These methods have not been used before in CFD simulations. The advantage of using these methods is reduction in communication times between cores. The numerical experiments indicate that the gains in computational performance of this new class of TSRK methods, as compared with classical Runge‐Kutta (RK) methods or low storage Runge‐Kutta (LSRK) schemes, are of the order of 25%, with no loss in accuracy.

An entropy–stable discontinuous Galerkin approximation for the incompressible Navier–Stokes equations with variable density and artificial compressibility

We present a provably stable discontinuous Galerkin spectral element method for the incompressible Navier–Stokes equations with artificial compressibility and variable density. Stability proofs, which include boundary conditions, that follow a continuous entropy analysis are provided. We define a mathematical entropy function that combines the traditional kinetic energy and an additional energy term for the artificial compressibility, and derive its associated entropy conservation law. The latter allows us to construct a provably stable split–form nodal Discontinuous Galerkin (DG) approximation that satisfies the summation–by–parts simultaneous–approximation–term (SBP–SAT) property. The scheme and the stability proof are presented for general curvilinear three–dimensional hexahedral meshes. We use the exact Riemann solver and the Bassi–Rebay 1 (BR1) scheme at the inter–element boundaries for inviscid and viscous fluxes respectively, and an explicit low storage Runge–Kutta RK3 scheme to integrate in time. We assess the accuracy and robustness of the method by solving a manufactured solution, the Kovasznay flow, a lid driven cavity, the inviscid Taylor–Green vortex, and the Rayleigh–Taylor instability.

Comparative Study of Standard Smagorinsky Model and Dynamic Smagorinsky Model in Large Eddy Simulation of Turbulent Channel Flow

This paper presents results of comparative study of large eddy simulation (LES) that is applied to a plane turbulent channel flow. The LES is performed by using a finite difference method of second order accuracy in space and a low-storage explicit Runge-Kutta method with third order accuracy in time. In the LES for subgrid-scale (SGS) modelling, Standard Smagorinsky Model (SSM) and Dynamic Smagorinsky Model (DSM) are used. Essential turbulence statistics from the two LES approaches are calculated and compared with those from direct numerical simulation (DNS) data. Comparing the results throughout the calculation domain, it has been found out that SSM performs better than DSM in the turbulent channel flow simulation. Flow structures in the computed flow field by the SSM and DSM are also discussed and compared through the contour plots and iso-surfaces.

Efficient SSP low-storage Runge–Kutta methods

In this paper we study the efficiency of Strong Stability Preserving (SSP) Runge–Kutta methods that can be implemented with a low number of registers using their Shu–Osher representation. SSP methods have been studied in the literature and stepsize restrictions that ensure numerical monotonicity have been found. However, for some problems, the observed stepsize restrictions are larger than the theoretical ones. Aiming at obtaining additional properties of the schemes that may explain their efficiency, in this paper we study the influence of the local error term in the observed stepsize restrictions. For this purpose, we consider the family of 5-stage third order SSP explicit Runge–Kutta methods, namely SSP(5,3), and the Buckley–Leverett equation. We deal with optimal SSP(5,3) schemes whose implementation requires at least 3 memory registers, and non-optimal 2-register SSP(5,3) schemes. The numerical experiments done show that small error constants improve the efficiency of the method in the sense that larger observed SSP coefficients are obtained.

A reduced-order discontinuous Galerkin method based on a Krylov subspace technique in nanophotonics

This paper is concerned with the design of a reduced-order model (ROM) based on a Krylov subspace technique for solving the time-domain Maxwell’s equations coupled to a Drude dispersion model, which are discretized in space by a discontinuous Galerkin (DG) method. An auxiliary differential equation (ADE) method is used to represent the constitutive relation for the dispersive medium. A semi-discrete DG scheme is formulated on an unstructured simplicial mesh, and is combined with a centered scheme for the definition of the numerical fluxes of the electric and magnetic fields on element interfaces. The ROM is established by projecting (Galerkin projection) the global semi-discrete DG scheme onto a low-dimensional Krylov subspace generated by an Arnoldi process. A low-storage Runge-Kutta (LSRK) time scheme is employed in the semi-discrete DG system and ROM. The overall goal is to reduce the computational time while maintaining an acceptable level of accuracy. We present numerical results on 2-D problems to show the effectiveness of the proposed method.

A low—storage Runge—Kutta OpenFOAM solver for compressible low—Mach number flows: aeroacoustic and thermo—fluid dynamic applications

A solver for compressible Navier–Stokes equations is presented in this paper. Low-storage RungeKutta schemes were adopted for time integration; on the other hand the finite volume approach available within OpenFOAM library has been adopted for space discretization. Kurganov-Noelle-Petrova approach was used for convective terms, while central schemes for diffusive ones. The aforementioned techniques were selected and tested in order to allow the possibility of solving a broad range of physical phenomena with particular emphasis to aeroacoustic and thermo-fluid dynamic problems. Indeed, that standard OpenFOAM solution techniques produce an unacceptable dissipation for acoustic phenomena computations. Non–reflective boundary treatment was also considered to avoid spurious numerical reflections. The reliability and the robustness of the solver is proved by computing several benchmarks. Lastly, the impact of the thermal boundary conditions on the sound propagation was analyzed.

Large eddy simulation of turbulent channel flow using differential equation wall model

A large eddy simulation (LES) of a plane turbulent channel flow is performed at a Reynolds number Re? = 590 based on the channel half width, ? and wall shear velocity, u? by approximating the near wall region using differential equation wall model (DEWM). The simulation is performed in a computational domain of 2?? x 2? x ??. The computational domain is discretized by staggered grid system with 32 x 30 x 32 grid points. In this domain the governing equations of LES are discretized spatially by second order finite difference formulation, and for temporal discretization the third order low-storage Runge-Kutta method is used. Essential turbulence statistics of the computed flow field based on this LES approach are calculated and compared with the available Direct Numerical Simulation (DNS) and LES data where no wall model was used. Comparing the results throughout the calculation domain we have found that the LES results based on DEWM show closer agreement with the DNS data, especially at the near wall region. That is, the LES approach based on DEWM can capture the effects of near wall structures more accurately. Flow structures in the computed flow field in the 3D turbulent channel have also been discussed and compared with LES data using no wall model.

A High-Order Discontinuous Galerkin Method for Coupled Wave Propagation in 1D Elastoplastic Heterogeneous Media

We propose a nodal high-order discontinuous Galerkin (DG) method for coupled wave propagation in heterogeneous elastoplastic soil columns. We solve the elastodynamic system written in velocity-strain formulation considering simultaneously the three components of motion. A 3D Iwan elastoplastic model for dry soils under dynamic loading is introduced in the DG method (DG-3C) based on centered fluxes and a low storage Runge–Kutta time scheme. Unlike the linear case, the nonlinear material behavior results in coupling effects between different components of motion. We focus on the choice of the numerical flux for nonlinear heterogeneous media and the classical centered flux is well adapted in terms of stability and direct implementation in 3D. Several numerical applications are studied using a soil column with realistic mechanical properties, and considering successively synthetic and borehole seismic recordings recorded at the Volvi test site (northern Greece).

Simulation of three-dimensional nanoscale light interaction with spatially dispersive metals using a high order curvilinear DGTD method

In this work, we present and study a flexible and accurate numerical solver in the context of three-dimensional computational nanophotonics. More precisely, we focus on the propagation of electromagnetic waves through metallic media described by a non-local dispersive model. For this model, we propose a discretization based on a high-order Discontinuous Galerkin time-domain method, along with a low-storage Runge–Kutta time scheme of order four. The semi-discrete stability of the scheme is analyzed for classical numerical fluxes, i.e. centered and upwind. Furthermore, the numerical treatment is enriched with an enhanced approximation of the geometry based on isoparametric curvilinear meshes. We finally assess our approach on several test cases, from academic to more physical ones.

Application of the Entropically Damped Artificial Compressibility model to direct numerical simulation of turbulent channel flow

The feasibility of the Entropically Damped Artificial Compressibility (EDAC) model for simulation of wall-bounded turbulence is examined. A favourable feature of EDAC is the purely parabolic type of governing equations resulting from the introduction of a damping term to the pressure evolution equation. This significantly reduces noise in the velocity divergence field when central schemes are applied for the spatial discretisation. The method of lines is used for the solution of the EDAC equations. A conservative finite difference method of purely 2nd- and mixed 4th/2nd-order is applied. The explicit 4th-order, six-stage low-storage Runge–Kutta method is used for time advancement. As expected, the mixed scheme is superior to the 2nd-order one in terms of both accuracy and computational efficiency.For the first time, the EDAC model is assessed for the direct numerical simulation of wall-bounded turbulent flow requiring a non-uniform mesh. As a particular test case we have chosen the channel flow at the friction Reynolds numbers Reτ=180 and 395. A very good agreement with the reference data is obtained, as documented by the chosen one-point velocity statistics: the mean and r.m.s. profiles and the budgets of the Reynolds stresses. The use of the explicit time discretisation and local (rather than compact) spatial schemes results in easily and efficiently parallelisable solution algorithm. We achieved very high computational performance on a desktop computer when solving the EDAC equations using a code dedicated for the graphics processing unit (GPU). The aspects of the GPU implementation are discussed.

Study of flow over object problems by a nodal discontinuous Galerkin-lattice Boltzmann method

The flow over object problems are studied by a nodal discontinuous Galerkin-lattice Boltzmann method (NDG-LBM) in this work. Different from the standard lattice Boltzmann method, the current method applies the nodal discontinuous Galerkin method into the streaming process in LBM to solve the resultant pure convection equation, in which the spatial discretization is completed on unstructured grids and the low-storage explicit Runge-Kutta scheme is used for time marching. The present method then overcomes the disadvantage of standard LBM for depending on the uniform meshes. Moreover, the collision process in the LBM is completed by using the multiple-relaxation-time scheme. After the validation of the NDG-LBM by simulating the lid-driven cavity flow, the simulations of flows over a fixed circular cylinder, a stationary airfoil and rotating-stationary cylinders are performed. Good agreement of present results with previous results is achieved, which indicates that the current NDG-LBM is accurate and effective for flow over object problems.

High-performance multi-GPU solver for describing nonlinear acoustic waves in homogeneous thermoviscous media

A double-precision numerical solver to describe the propagation of high-intensity ultrasound fluctuations using a novel finite-amplitude compressible acoustic model working in multiple processing units (GPUs) is presented. The present solver is based on a conservative hyperbolic formulation derived from a variational analysis of the compressible Navier–Stokes equations and is implemented using an explicit high-order finite difference strategy. In this work, a WENO–Z reconstruction scheme along with a high-order finite-difference stencil are used to approximate the contributions of convective and diffusive spatial operators, respectively. The spatial operators are then associated to a low–storage Runge–Kutta scheme to integrate the system explicitly in time. The present multi-GPU implementation aims to make the best use of every single GPU and gain optimal performance of the algorithm on the per-node basis. To assess the performance of the present solver, a typical mini-server computer with 4 Tesla K80 dual GPU accelerators is used. The results show that the present formulation scales linearly for large domain problems. Moreover, when compared to an OpenMP implementation running with an i7 processor of 4.2 GHz, this is outperformed by our MPI-GPU implementation by a factor of 99. In this work, the present multi-GPU solver is illustrated with a three-dimensional simulation of a highly-intense focused ultrasound propagation.

Petascale simulations of compressible flows with interfaces

We demonstrate a high throughput software for the efficient simulation of compressible multicomponent flow on high performance computing platforms. The discrete problem is represented on structured three-dimensional grids with non-uniform resolution. Discontinuous flow features are captured using a diffuse interface method. A distinguishing characteristic of the method is the proper treatment of the interface zone as a mixing region of liquid and gas. The governing equations are discretized by a Godunov-type finite volume method with explicit time stepping using a low-storage Runge-Kutta scheme. The presented flow solver Cubism-MPCF is based on our Cubism library which enables a highly optimized framework for the efficient treatment of stencil based problems on multicore architectures. The framework is general and not limited to applications in fluid dynamics. We validate our solver by classical benchmark examples. Furthermore, we examine a highly-resolved shock-induced bubble collapse and a cloud of O(103) collapsing bubbles, which demonstrate the high potential of the proposed framework and solver.