Explicit Runge-Kutta Scheme Research Articles

Numerical Simulation of Transitional Flow Past a Wing Using Flux Reconstruction Method

In this paper, the capabilities of a high-order flux reconstruction scheme for the implicit large-eddy simulations of a transitional flow at a low Reynolds number of 60,000 around the SD7003 wing are examined. This is the case posed in the First and Second International Workshops on the high-order computational fluid dynamics method held in 2012 and 2013. The correction function used in the flux reconstruction scheme is the Radau polynomial, which is equivalent to the discontinuous Galerkin scheme. A time-accurate implicit lower/upper symmetric Gauss–Seidel solution algorithm for the application of the flux reconstruction scheme to complex unsteady flows is developed, and it is found to be able to produce comparable results to the explicit Runge–Kutta scheme while achieving better computational efficiency. Simulations are carried out at and 8 deg with polynomials of degree , 2, and 3 resulting in second-, third-, and fourth-order-accurate flux reconstruction schemes, respectively. Two structured hexahedral O-grid domains that differ in the grid resolution in a circumferential direction on the upper surface of the wing are considered, with a maximum of 2,000,000 degrees of freedom for the fourth-order () simulations on the finer domain. The results are validated by comparison with many reference data obtained from various high-order schemes using time-accurate explicit/implicit methods. The results agree reasonably well. The developed method with the flux reconstruction scheme can be a reliable tool of implicit large-eddy simulations for low-Reynolds-number flows.

An efficient conservative cut-cell method for rigid bodies interacting with viscous compressible flows

A Cartesian cut-cell method for viscous flows interacting with freely moving boundaries is presented. The method enables a sharp resolution of the embedded boundaries and strictly conserves mass, momentum, and energy. A new explicit Runge-Kutta scheme (PC-RK) is introduced by which the overall computational time is reduced by a factor of up to 2.5. The new scheme is a predictor-corrector type reformulation of a popular class of Runge-Kutta methods which substantially reduces the computational effort for tracking the moving boundaries and subsequently reinitializing the solver impairing neither stability nor accuracy. The structural motion is computed by an implicit scheme with good stability properties due to a strong-coupling strategy and the conservative discretization of the flow solver at the material interfaces. A new formulation for the treatment of small cut cells is proposed with high accuracy and robustness for arbitrary geometries based on a weighted Taylor-series approach solved via singular-value decomposition. The efficiency and the accuracy of the new method are demonstrated for several three-dimensional cases of laminar and turbulent particulate flow. It is shown that the new method remains fully conservative even for large displacements of the boundaries leading to a fast convergence of the fluid-solid coupling while spurious force oscillations inherent to this class of methods are effectively suppressed. The results substantiate the good stability and accuracy properties of the scheme even on relatively coarse meshes. An efficient Cartesian cut-cell method for freely moving solid boundaries is presented.The new predictor-corrector Runge-Kutta scheme (PC-RK) enables simulation speedups of up to 2.5.Mass, momentum, and energy are strictly conserved at the sharply resolved solid-fluid interfaces.The method exhibits very good stability and accuracy properties even for large displacements of the boundaries.Spurious force oscillations inherent to this class of methods are suppressed.

Nonlinear seismic response of a partially-filled rectangular liquid tank with a submerged block

The seismic response of partially-filled two-dimensional rigid rectangular liquid tanks with a bottom-mounted submerged block is numerically simulated. The Galerkin-weighted-residual based finite element method (FEM) is used for solving the governing Laplace equation with fully nonlinear free surface boundary conditions and also for velocity recovery. Based on the mixed Eulerian–Lagrangian (MEL) method, a fourth order explicit Runge–Kutta scheme is used for the time-stepping integration of free surface boundary conditions. A cubic-spline fitted regridding technique is used at every time step to eliminate possible numerical instabilities on account of Lagrangian node induced mesh distortion. An artificial surface damping term is used to mimic the viscosity induced damping. Three different earthquake motions characterized on the basis of low, intermediate and high frequency contents are used to study the effect of frequency content on the nonlinear dynamic response of this tank-liquid-submerged block system. The effect of the submerged block on the impulsive and convective response components of the hydrodynamic forces manifested in terms of base shear, overturning base moment and pressure distribution along the tank wall as well as the block wall has been quantified vis-a-vis frequency content of ground motions. It is observed that the convective response of this tank-liquid system is highly sensitive to the frequency content of the ground motion.

GPU IMPLEMENTATION OF A VISCOUS FLOW SOLVER ON UNSTRUCTURED GRIDS

Graphics processing units have gained popularities in scientific computing over past several years due to their outstanding parallel computing capability. Computational fluid dynamics applications involve large amounts of calculations, therefore a latest GPU card is preferable of which the peak computing performance and memory bandwidth are much better than a contemporary high-end CPU. We herein focus on the detailed implementation of our GPU targeting Reynolds-averaged Navier-Stokes equations solver based on finite-volume method. The solver employs a vertex-centered scheme on unstructured grids for the sake of being capable of handling complex topologies. Multiple optimizations are carried out to improve the memory accessing performance and kernel utilization. Both steady and unsteady flow simulation cases are carried out using explicit Runge-Kutta scheme. The solver with GPU acceleration in this paper is demonstrated to have competitive advantages over the CPU targeting one.

Multirate Technique for Explicit Discontinuous Galerkin Computations of Time-Domain Maxwell Equations on Complex Geometries

This paper presents a multirate technique to improve and optimize the time step of a second-order two-stage explicit Runge-Kutta (ERK) scheme. In this technique, the mesh elements are stored in different groups according to their stable time steps. These groups are sorted into two classes: 1) the bulk groups, where the two-stage ERK method is applied once or repeatedly; 2) the buffer groups that are used to accommodate the transition between two bulk groups. This technique is proposed for accelerating the explicit discontinuous Galerkin computations of the time-domain Maxwell equations. An application example on a human skull is proposed to show the efficiency of this technique and to simulate wave propagation on complex geometries.

Time consistent fluid structure interaction on collocated grids for incompressible flow

Abstract Consistent time integration on collocated grids for incompressible flow has been studied for static grids using the PISO method, in which the dependencies on time-step size and under-relaxation has been studied in detail. However, for moving grids a detailed description is still missing. Therefore, a step by step analysis of a time consistent fluid–structure interaction (FSI) method for incompressible flow on collocated grids is presented. The method consist of: face normal and area correction for moving grids, treatment of velocity boundary conditions for no-slip walls, time integration of structure equations and fluid force interpolation to structure. The basis of the method is the PISO method of the incompressible Navier–Stokes equations. Time consistency on static grids is shown first, after which time consistency on moving grids is described and analyzed. For moving grids consistent time integration is described in two ways: (1) constructing the face velocities from a normal and tangential component, and (2) correcting the face flux with a face normal and face area correction. For both descriptions the general formulation for the backward differencing schemes (BDF) are given and the correct behavior of the first, second and third order schemes is demonstrated by means of an academic test case (circular cavity flow). Additionally, the (force) coupling from the fluid to the structure is discussed in detail for combining a fourth order explicit Runge Kutta scheme for the structure with a BDF scheme for the fluid. Three interpolations for the fluid forces are shown, which result in either a first order or second order FSI scheme. Third order FSI is demonstrated when the third order BDF scheme is applied on both the structure and fluid equations. Also, under-relaxation for the fluid equations is considered and it is demonstrated that the order of the three BDF schemes are independent of the under-relaxation factor. Finally, the proposed method of time consistent FSI on collocated grids for incompressible flows is demonstrated by applying it to a three-dimensional flow over an elastic structure in a channel and the cylinder flap FSI benchmark case of Turek and Hron.

Dynamic earthquake rupture simulations on nonplanar faults embedded in 3D geometrically complex, heterogeneous elastic solids

Dynamic propagation of shear ruptures on a frictional interface in an elastic solid is a useful idealization of natural earthquakes. The conditions relating discontinuities in particle velocities across fault zones and tractions acting on the fault are often expressed as nonlinear friction laws. The corresponding initial boundary value problems are both numerically and computationally challenging. In addition, seismic waves generated by earthquake ruptures must be propagated for many wavelengths away from the fault. Therefore, reliable and efficient numerical simulations require both provably stable and high order accurate numerical methods.We present a high order accurate finite difference method for: a) enforcing nonlinear friction laws, in a consistent and provably stable manner, suitable for efficient explicit time integration; b) dynamic propagation of earthquake ruptures along nonplanar faults; and c) accurate propagation of seismic waves in heterogeneous media with free surface topography.We solve the first order form of the 3D elastic wave equation on a boundary-conforming curvilinear mesh, in terms of particle velocities and stresses that are collocated in space and time, using summation-by-parts (SBP) finite difference operators in space. Boundary and interface conditions are imposed weakly using penalties. By deriving semi-discrete energy estimates analogous to the continuous energy estimates we prove numerical stability. The finite difference stencils used in this paper are sixth order accurate in the interior and third order accurate close to the boundaries. However, the method is applicable to any spatial operator with a diagonal norm satisfying the SBP property. Time stepping is performed with a 4th order accurate explicit low storage Runge-Kutta scheme, thus yielding a globally fourth order accurate method in both space and time. We show numerical simulations on band limited self-similar fractal faults revealing the complexity of rupture dynamics on rough faults.

Mpdcm: A toolbox for massively parallel dynamic causal modeling

BackgroundDynamic causal modeling (DCM) for fMRI is an established method for Bayesian system identification and inference on effective brain connectivity. DCM relies on a biophysical model that links hidden neuronal activity to measurable BOLD signals. Currently, biophysical simulations from DCM constitute a serious computational hindrance. Here, we present Massively Parallel Dynamic Causal Modeling (mpdcm), a toolbox designed to address this bottleneck. New methodmpdcm delegates the generation of simulations from DCM's biophysical model to graphical processing units (GPUs). Simulations are generated in parallel by implementing a low storage explicit Runge–Kutta's scheme on a GPU architecture. mpdcm is publicly available under the GPLv3 license. ResultsWe found that mpdcm efficiently generates large number of simulations without compromising their accuracy. As applications of mpdcm, we suggest two computationally expensive sampling algorithms: thermodynamic integration and parallel tempering. Comparison with existing method(s)mpdcm is up to two orders of magnitude more efficient than the standard implementation in the software package SPM. Parallel tempering increases the mixing properties of the traditional Metropolis–Hastings algorithm at low computational cost given efficient, parallel simulations of a model. ConclusionsFuture applications of DCM will likely require increasingly large computational resources, for example, when the likelihood landscape of a model is multimodal, or when implementing sampling methods for multi-subject analysis. Due to the wide availability of GPUs, algorithmic advances can be readily available in the absence of access to large computer grids, or when there is a lack of expertise to implement algorithms in such grids.

Full-waveform inversion in three-dimensional PML-truncated elastic media

We are concerned with high-fidelity subsurface imaging of the soil, which commonly arises in geotechnical site characterization and geophysical explorations. Specifically, we attempt to image the spatial distribution of the Lamé parameters in semi-infinite, three-dimensional, arbitrarily heterogeneous formations, using surficial measurements of the soil’s response to probing elastic waves. We use the complete waveform response of the medium to drive the inverse problem , by using a partial-differential-equation (PDE)-constrained optimization approach, directly in the time-domain, to minimize the misfit between the observed response of the medium at select measurement locations, and a computed response corresponding to a trial distribution of the Lamé parameters. We discuss strategies that lend algorithmic robustness to our proposed inversion scheme. To limit the computational domain to the size of interest, we employ perfectly-matched-layers (PMLs). In order to resolve the forward problem, we use a recently developed hybrid finite element approach, where a displacement–stress formulation for the PML is coupled to a standard displacement-only formulation for the interior domain, thus leading to a computationally cost-efficient scheme. Time-integration is accomplished by using an explicit Runge–Kutta scheme, which is well-suited for large-scale problems on parallel computers. We verify the accuracy of the material gradients obtained via our proposed scheme, and report numerical results demonstrating successful reconstruction of the two Lamé parameters for both smooth and sharp profiles.

Phase Resolving Wave—Current Interactions with Improved Boussinesq-Type Equations

A set of fully nonlinear Boussinessq-type equations with improved linear and nonlinear properties is considered for wave—current interaction analysis. These phase-resolving equations are so that the highest order of the derivatives is three. We implement a new source function for the wave—current generation within the domain, which allows to generate a wide range of wave—current conditions. The set of equations is solved using a fourth order explicit numerical scheme which semi-discretizes the equations in space and then integrates in time using an explicit Runge—Kutta scheme. A novel treatment of the boundaries, which uses radiative boundary conditions for the current and damps the waves, is used to avoid boundary reflections. Several validation tests are presented to demonstrate the capabilities of the new model equations for wave—current interaction.

An Adaptive Fast Interface Tracking Method

An adaptive numerical scheme is developed for the propagation of an interface in a velocity field based on the fast interface tracking method proposed in [2]. A multiresolution stategy to represent the interface instead of point values, allows local grid refinement while controlling the approximation error on the interface. For time integration, we use an explicit Runge-Kutta scheme of second-order with a multiscale time step, which takes longer time steps for finer spatial scales. The implementation of the algorithm uses a dynamic tree data structure to represent data in the computer memory. We briefly review first the main algorithm, describe the essential data structures, highlight the adaptive scheme, and illustrate the computational efficiency by some numerical examples.

Adaptive Mesh Refinement for Hyperbolic Systems Based on Third-Order Compact WENO Reconstruction

In this paper we generalize to non-uniform grids of quad-tree type the Compact WENO reconstruction of Levy, Puppo and Russo (SIAM J. Sci. Comput., 2001), thus obtaining a truly two-dimensional non-oscillatory third order reconstruction with a very compact stencil and that does not involve mesh-dependent coefficients. This latter characteristic is quite valuable for its use in h-adaptive numerical schemes, since in such schemes the coefficients that depend on the disposition and sizes of the neighboring cells (and that are present in many existing WENO-like reconstructions) would need to be recomputed after every mesh adaption. In the second part of the paper we propose a third order h-adaptive scheme with the above-mentioned reconstruction, an explicit third order TVD Runge-Kutta scheme and the entropy production error indicator proposed by Puppo and Semplice (Commun. Comput. Phys., 2011). After devising some heuristics on the choice of the parameters controlling the mesh adaption, we demonstrate with many numerical tests that the scheme can compute numerical solution whose error decays as $\langle N\rangle^{-3}$, where $\langle N\rangle$ is the average number of cells used during the computation, even in the presence of shock waves, by making a very effective use of h-adaptivity and the proposed third order reconstruction.

An efficient time advancing strategy for energy-preserving simulations

Energy-conserving numerical methods are widely employed within the broad area of convection-dominated systems. Semi-discrete conservation of energy is usually obtained by adopting the so-called skew-symmetric splitting of the non-linear convective term, defined as a suitable average of the divergence and advective forms. Although generally allowing global conservation of kinetic energy, it has the drawback of being roughly twice as expensive as standard divergence or advective forms alone. In this paper, a general theoretical framework has been developed to derive an efficient time-advancement strategy in the context of explicit Runge–Kutta schemes. The novel technique retains the conservation properties of skew-symmetric-based discretizations at a reduced computational cost. It is found that optimal energy conservation can be achieved by properly constructed Runge–Kutta methods in which only divergence and advective forms for the convective term are used. As a consequence, a considerable improvement in computational efficiency over existing practices is achieved. The overall procedure has proved to be able to produce new schemes with a specified order of accuracy on both solution and energy. The effectiveness of the method as well as the asymptotic behavior of the schemes is demonstrated by numerical simulation of Burgers' equation.

Low-storage implicit/explicit Runge–Kutta schemes for the simulation of stiff high-dimensional ODE systems

Implicit/explicit (IMEX) Runge–Kutta (RK) schemes are effective for time-marching ODE systems with both stiff and nonstiff terms on the RHS; such schemes implement an (often A-stable or better) implicit RK scheme for the stiff part of the ODE, which is often linear, and, simultaneously, a (more convenient) explicit RK scheme for the nonstiff part of the ODE, which is often nonlinear. Low-storage RK schemes are especially effective for time-marching high-dimensional ODE discretizations of PDE systems on modern (cache-based) computational hardware, in which memory management is often the most significant computational bottleneck. In this paper, we develop and characterize eight new low-storage implicit/explicit RK schemes which have higher accuracy and better stability properties than the only low-storage implicit/explicit RK scheme available previously, the venerable second-order Crank–Nicolson/Runge–Kutta–Wray (CN/RKW3) algorithm that has dominated the DNS/LES literature for the last 25 years, while requiring similar storage (two, three, or four registers of length N) and comparable floating-point operations per timestep.

Linearly Implicit IMEX Runge--Kutta Methods for a Class of Degenerate Convection-Diffusion Problems

Multispecies kinematic flow models with strongly degenerate diffusive corrections give rise to systems of nonlinear convection-diffusion equations of arbitrary size. Applications of these systems include models of polydisperse sedimentation and multiclass traffic flow. Implicit-explicit (IMEX) Runge--Kutta (RK) methods are suitable for the solution of these convection-diffusion problems since the stability restrictions, coming from the explicitly treated convective part, are much less severe than those that would be deduced from an explicit treatment of the diffusive term. These schemes usually combine an explicit RK scheme for the time integration of the convective part with a diagonally implicit one for the diffusive part. In [R. Burger, P. Mulet, and L. M. Villada, SIAM J. Sci. Comput., 35 (2013), pp. B751--B777] a scheme of this type is proposed, where the nonlinear and nonsmooth systems of algebraic equations arising in the implicit treatment of the degenerate diffusive part are solved by smoothing o...

An Implicit Algorithm for High-Order DG/FV Schemes for Compressible Flows on 2D Arbitrary Grids

AbstractA Newton/LU-SGS (lower-upper symmetric Gauss-Seidel) iteration implicit method was developed to solve two-dimensional Euler and Navier-Stokes equations by the DG/FV hybrid schemes on arbitrary grids. The Newton iteration was employed to solve the nonlinear system, while the linear system was solved with LU-SGS iteration. The effect of several parameters in the implicit scheme, such as the CFL number, the Newton sub-iteration steps, and the update frequency of Jacobian matrix, was investigated to evaluate the performance of convergence history. Several typical test cases were simulated, and compared with the traditional explicit Runge-Kutta (RK) scheme. Firstly the Couette flow was tested to validate the order of accuracy of the present DG/FV hybrid schemes. Then a subsonic inviscid flow over a bump in a channel was simulated and the effect of parameters was investigated also. Finally, the implicit algorithm was applied to simulate a subsonic inviscid flow over a circular cylinder and the viscous flow in a square cavity. The numerical results demonstrated that the present implicit scheme can accelerate the convergence history efficiently. Choosing proper parameters would improve the efficiency of the implicit scheme. Moreover, in the same framework, the DG/FV hybrid schemes are more efficient than the same order DG schemes.

Achievable efficiency of numerical methods for simulations of solar surface convection

We investigate the achievable efficiency of both the time and the space discretisation methods used in Antares for mixed parabolic–hyperbolic problems. We show that the fifth order variant of WENO combined with a second order Runge–Kutta scheme is not only more accurate than standard first and second order schemes, but also more efficient taking the computation time into account. Then, we calculate the error decay rates of WENO with several explicit Runge–Kutta schemes for advective and diffusive problems with smooth and non-smooth initial conditions. With this data, we estimate the computational costs of three-dimensional simulations of stellar surface convection and show that SSPRK(3,2) is the most efficient scheme considered in this comparison.

Energy stable and large time-stepping methods for the Cahn–Hilliard equation

We present the numerical methods for the Cahn–Hilliard equation, which describes phase separation phenomenon. The goal of this paper is to construct high-order, energy stable and large time-stepping methods by using Eyre's convex splitting technique. The equation is discretized by using a fourth-order compact difference scheme in space and first-order, second-order or third-order implicit–explicit Runge–Kutta schemes in time. The energy stability for the first-order scheme is proved. Numerical experiments are given to demonstrate the performance of the proposed methods.

Time‐domain hybrid formulations for wave simulations in three‐dimensional PML‐truncated heterogeneous media

SUMMARYWe are concerned with the numerical simulation of wave motion in arbitrarily heterogeneous, elastic, perfectly‐matched‐layer‐(PML)‐truncated media. We extend in three dimensions a recently developed two‐dimensional formulation, by treating the PML via an unsplit‐field, but mixed‐field, displacement‐stress formulation, which is then coupled to a standard displacement‐only formulation for the interior domain, thus leading to a computationally cost‐efficient hybrid scheme. The hybrid treatment leads to, at most, third‐order in time semi‐discrete forms. The formulation is flexible enough to accommodate the standard PML, as well as the multi‐axial PML.We discuss several time‐marching schemes, which can be used à la carte, depending on the application: (a) an extended Newmark scheme for third‐order in time, either unsymmetric or fully symmetric semi‐discrete forms; (b) a standard implicit Newmark for the second‐order, unsymmetric semi‐discrete forms; and (c) an explicit Runge–Kutta scheme for a first‐order in time unsymmetric system. The latter is well‐suited for large‐scale problems on parallel architectures, while the second‐order treatment is particularly attractive for ready incorporation in existing codes written originally for finite domains.We compare the schemes and report numerical results demonstrating stability and efficacy of the proposed formulations. Copyright © 2014 John Wiley & Sons, Ltd.

Nonuniform-time-step explicit Runge–Kutta scheme for high-order finite difference method

Explicit Runge–Kutta method has been widely used for time-accurate simulations in aeroacoustics and aerodynamics, partly because of its strong stability properties. However, when dealing with problems involving irregular geometries or multi-scale phenomena where the grid is often refined in localized regions, a single global time step size dictated by the smallest grid size based on the explicit stability condition would inevitably result in excessive computational consumption. A local time stepping strategy that adopts time step size to the local stability requirement in different grid blocks is an effective way to reduce the cost of time integration and to increase the overall computational efficiency. In the present paper, a non-uniform time step (NUTS) explicit Runge–Kutta scheme, previously introduced in the framework of discontinuous Galerkin method, which makes the local time stepping strategy applicable to the explicit Runge–Kutta family without any interpolation and extrapolation in time, is further extended to high-order finite difference methods. The stability of non-uniform time step (NUTS) scheme in combination with high order difference schemes is studied. Numerical experiments are carried out using 1D and 2D acoustic problems for validation of the proposed approach. Computational cost reduction by the proposed algorithm is also discussed.