Exponential Time Differencing Algorithm for Debye Medium in FDTD

The exponential time differencing (ETD) formultion is proposed for Debye medium using finite-difference time-domain (FDTD) method. The schemes of the auxiliary differential equation (ADE) and ETD algorithms are presented. The algorithms are validated by comparing the simulation results with analytical values. The ETD method costs the same memory as the ADE method, and the numerical results indicate that the ETD technique is more accurate than the ADE method. This work is supported by the foundation under Grant No. Zl 0502.     

Operator Splitting Implicit Integration Factor Methods for Stiff Reaction-Diffusion-Advection Systems

For reaction-diffusion-advection equations, the stiffness from the reaction and diffusion terms often requires very restricted time step size, while the nonlinear advection term may lead to a sharp gradient in localized spatial regions. It is challenging to design numerical methods that can efficiently handle both difficulties. For reaction-diffusion systems with both stiff reaction and diffusion terms, implicit integration factor (IIF) method and its higher dimensional analog compact IIF (cIIF) serve as an efficient class of time-stepping methods, and their second order version is linearly unconditionally stable. For nonlinear hyperbolic equations, weighted essentially non-oscillatory (WENO) methods are a class of schemes with a uniformly high-order of accuracy in smooth regions of the solution, which can also resolve the sharp gradient in an accurate and essentially non-oscillatory fashion. In this paper, we couple IIF/cIIF with WENO methods using the operator splitting approach to solve reaction-diffusion-advection equations. In particular, we apply the IIF/cIIF method to the stiff reaction and diffusion terms and the WENO method to the advection term in two different splitting sequences. Calculation of local truncation error and direct numerical simulations for both splitting approaches show the second order accuracy of the splitting method, and linear stability analysis and direct comparison with other approaches reveals excellent efficiency and stability properties. Applications of the splitting approach to two biological systems demonstrate that the overall method is accurate and efficient, and the splitting sequence consisting of two reaction-diffusion steps is more desirable than the one consisting of two advection steps, because CWC exhibits better accuracy and stability.     

Finite Element $θ$-Schemes for the Acoustic Wave Equation

In this paper, we investigate the stability and convergence of a family of implicit finite difference schemes in time and Galerkin finite element methods in space for the numerical solution of the acoustic wave equation. The schemes cover the classical explicit second-order leapfrog scheme and the fourth-order accurate scheme in time obtained by the modified equation method. We derive general stability conditions for the family of implicit schemes covering some well-known CFL conditions. Optimal error estimates are obtained. For sufficiently smooth solutions, we demonstrate that the maximal error in the $L^2$-norm error over a finite time interval converges optimally as $\mathcal{O}(h^{p+1}+∆t^s)$, where $p$ denotes the polynomial degree, $s$=2 or 4, $h$ the mesh size, and $∆t$ the time step.     

Accuracy preserving limiter for the high-order accurate solution of the Euler equations

Higher-order finite-volume methods have been shown to be more efficient than second-order methods. However, no consensus has been reached on how to eliminate the oscillations caused by solution discontinuities. Essentially non-oscillatory (ENO) schemes provide a solution but are computationally expensive to implement and may not converge well for steady-state problems. This work studies the extension of limiters used for second-order methods to the higher-order case. Requirements for accuracy and efficient convergence are discussed. A new limiting procedure is proposed. Ringleb’s flow problem is used to demonstrate that nearly nominal orders of accuracy for schemes up to fourth-order can be achieved in smooth regions using the new limiter. Results for the fourth-order accurate solution of transonic flow demonstrates good convergence properties and significant qualitative improvement of the solution relative the second-order method. The new limiter can also be successfully applied to reduce the dissipation of second-order schemes with minimal sacrifices in convergence properties relative to existing approaches.     

Periodic Filtering as a subgrid-scale model for LES of laminar separation bubble flows

Laminar separation bubbles develop over many blades and airfoils at moderate angles of attack and Reynolds numbers ranging from 10⁴ to 10⁵. More accurate simulation tools are necessary to enable higher fidelity design optimisation for these airfoils and blades as well as to test flow control schemes. Following previous investigators, an equivalent problem is formulated by imposing suitable boundary conditions for flow over a flat plate which allows to use a high accuracy spectral solver. Large eddy simulation (LES) of such a flow were performed at drastically reduced resolution to assess the accuracy of several LES modelling approaches: the explicit dynamic Smagorinsky model, implicit LES, and the truncated Navier–Stokes approach (TNS). To mimic dissipation that occurs in implicit LES, the solution on a coarse mesh is filtered at every time step and two different filter strengths are used. In the TNS approach, the solution is filtered periodically, every few hundred time steps. The performance of each approach is evaluated against benchmark direct numerical simulation (DNS) data focusing on pressure and skin friction distributions, which are critical to airfoil designers. TNS results confirm that periodic filtering can act as an apt substitute for explicit subgrid-scale models, whereas filtering at every time step demonstrates the dependence of implicit LES on details of numerics.     

Alternating direction implicit schemes for the two-dimensional fractional sub-diffusion equation

New numerical techniques are presented for the solution of a two-dimensional anomalous sub-diffusion equation with time fractional derivative. In these methods, standard central difference approximation is used for the spatial discretization, and, for the time stepping, two new alternating direction implicit (ADI) schemes based on the L₁ approximation and backward Euler method are considered. The two ADI schemes are constructed by adding two different small terms, which are different from standard ADI methods. The solvability, unconditional stability and H¹ norm convergence are proved. Numerical results are presented to support our theoretical analysis and indicate the efficiency of both methods.     

A discontinuous Galerkin method for inviscid low Mach number flows

In this work we extend the high-order discontinuous Galerkin (DG) finite element method to inviscid low Mach number flows. The method here presented is designed to improve the accuracy and efficiency of the solution at low Mach numbers using both explicit and implicit schemes for the temporal discretization of the compressible Euler equations. The algorithm is based on a classical preconditioning technique that in general entails modifying both the instationary term of the governing equations and the dissipative term of the numerical flux function (full preconditioning approach). In the paper we show that full preconditioning is beneficial for explicit time integration while the implicit scheme turns out to be efficient and accurate using just the modified numerical flux function. Thus the implicit scheme could also be used for time accurate computations. The performance of the method is demonstrated by solving an inviscid flow past a NACA0012 airfoil at different low Mach numbers using various degrees of polynomial approximations. Computations with and without preconditioning are performed on different grid topologies to analyze the influence of the spatial discretization on the accuracy of the DG solutions at low Mach numbers.     

Locally implicit discontinuous Galerkin method for time domain electromagnetics

In the recent years, there has been an increasing interest in discontinuous Galerkin time domain (DGTD) methods for the solution of the unsteady Maxwell equations modeling electromagnetic wave propagation. One of the main features of DGTD methods is their ability to deal with unstructured meshes which are particularly well suited to the discretization of the geometrical details and heterogeneous media that characterize realistic propagation problems. Such DGTD methods most often rely on explicit time integration schemes and lead to block diagonal mass matrices. However, explicit DGTD methods are also constrained by a stability condition that can be very restrictive on highly refined meshes and when the local approximation relies on high order polynomial interpolation. An implicit time integration scheme is a natural way to obtain a time domain method which is unconditionally stable but at the expense of the inversion of a global linear system at each time step. A more viable approach consists of applying an implicit time integration scheme locally in the refined regions of the mesh while preserving an explicit time scheme in the complementary part, resulting in an hybrid explicit–implicit (or locally implicit) time integration strategy. In this paper, we report on our recent efforts towards the development of such a hybrid explicit–implicit DGTD method for solving the time domain Maxwell equations on unstructured simplicial meshes. Numerical experiments for 3D propagation problems in homogeneous and heterogeneous media illustrate the possibilities of the method for simulations involving locally refined meshes.     

Symplectic Euler Method for Nonlinear High Order Schrödinger Equation with a Trapped Term

In this paper, we establish a family of symplectic integrators for a class of high order Schrödinger equations with trapped terms. First, we find its symplectic structure and reduce it to a finite dimensional Hamilton system via spatial discretization. Then we apply the symplectic Euler method to the Hamiltonian system. It is demonstrated that the scheme not only preserves symplectic geometry structure of the original system, but also does not require to resolve coupled nonlinear algebraic equations which is different from the general implicit symplectic schemes. The linear stability of the symplectic Euler scheme and the errors of the numerical solutions are investigated. It shows that the semi-explicit scheme is conditionally stable, first order accurate in time and $2l^{th}$ order accuracy in space. Numerical tests suggest that the symplectic integrators are more effective than non-symplectic ones, such as backward Euler integrators.     

Partially implicit peer methods for the compressible Euler equations

When cut cells are used for the representation of orography in numerical weather prediction models this leads to very small cells. On one hand this results in very harsh time step restrictions for explicit methods due to the CFL criterion. On the other hand cut cells only appear in a small region of the domain. Therefore we consider a partially implicit method: In cut cells the Jacobian incorporates advection, diffusion and acoustics while in the full cells of the free atmosphere only the acoustic part is used, i.e. the method is linearly implicit in the cut cell regions and semi-explicit in the free regions and computes with time step sizes restricted only by the CFL condition in the free atmosphere. Furthermore we use a simplified Jacobian in the cut cell regions in order to save storage and gain computational efficiency. While the method retains the order independently of the Jacobian we present a linear stability theory which takes the effects of the simplifications of the Jacobian on stability into account. The presented method is as stable and accurate as the underlying split-explicit method but furthermore it can compute with cut cells with nearly no additional effort.     

A spectral FC solver for the compressible Navier–Stokes equations in general domains I: Explicit time-stepping

We present a Fourier continuation (FC) algorithm for the solution of the fully nonlinear compressible Navier–Stokes equations in general spatial domains. The new scheme is based on the recently introduced accelerated FC method, which enables use of highly accurate Fourier expansions as the main building block of general-domain PDE solvers. Previous FC-based PDE solvers are restricted to linear scalar equations with constant coefficients. The FC methodology presented in this text thus constitutes a significant generalization of the previous FC schemes, as it yields general-domain FC solvers for nonlinear systems of PDEs. While not restricted to periodic boundary conditions and therefore applicable to general boundary value problems on arbitrary domains, the proposed algorithm inherits many of the highly desirable properties arising from rapidly convergent Fourier expansions, including high-order convergence, essentially spectrally accurate dispersion relations, and much milder CFL constraints than those imposed by polynomial-based spectral methods—since, for example, the spectral radius of the FC first derivative grows linearly with the number of spatial discretization points. We demonstrate the accuracy and optimal parallel efficiency of the algorithm in a variety of scientific and engineering contexts relevant to fluid-dynamics and nonlinear acoustics.     

XTOR-2F: A fully implicit Newton–Krylov solver applied to nonlinear 3D extended MHD in tokamaks

XTOR-2F solves a set of extended magnetohydrodynamic (MHD) equations in toroidal tokamak geometry. In the original XTOR code, the time stepping is handled by a semi-implicit method 1, 2 and 3. Moderate changes were necessary to transform it into a fully implicit one using the NITSOL library with Newton–Krylov methods of solution for nonlinear system of equations [4]. After addressing the sensitive issue of preconditioning and time step tuning, the performances of the semi-implicit and the implicit methods are compared for the nonlinear simulation of an internal kink mode test case within the framework of resistive MHD including anisotropic thermal transport. A convergence study comparing the semi-implicit and the implicit schemes is presented. Our main conclusion is that on one hand the Newton–Krylov implicit method, when applied to basic one fluid MHD is more computationally costly than the semi-implicit one by a factor 3 for a given numerical accuracy. But on the other hand, the implicit method allows to address challenging issues beyond MHD. By testing the Newton–Krylov method with diamagnetic modifications on the dynamics of the internal kink, some numerical issues, to be addressed further, are emphasized.     

Optimized prefactored compact schemes

The numerical simulation of aeroacoustic phenomena requires high-order accurate numerical schemes with low dispersion and dissipation errors. In this paper we describe a strategy for developing high-order accurate prefactored compact schemes, requiring very small stencil support. These schemes require fewer boundary stencils and offer simpler boundary condition implementation than existing compact schemes. The prefactorization strategy splits the central implicit schemes into forward and backward biased operators. Using Fourier analysis, we show it is possible to select the coefficients of the biased operators such that their dispersion characteristics match those of the original central compact scheme and their numerical wavenumbers have equal and opposite imaginary components. This ensures that when the forward and backward stencils are added, the original central compact scheme is recovered. To extend the resolution characteristic of the schemes, an optimization strategy is employed in which formal order of accuracy is sacrificed in preference to enhanced resolution characteristics across the range of wavenumbers realizable on a given mesh. The resulting optimized schemes yield improved dispersion characteristics compared to the standard sixth- and eighth-order compact schemes making them more suitable for high-resolution numerical simulations in gas dynamics and computational aeroacoustics. The efficiency, accuracy and convergence characteristics of the new optimized prefactored compact schemes are demonstrated by their application to several test problems.     

ADI spectral collocation methods for parabolic problems

We discuss the Crank–Nicolson and Laplace modified alternating direction implicit Legendre and Chebyshev spectral collocation methods for a linear, variable coefficient, parabolic initial-boundary value problem on a rectangular domain with the solution subject to non-zero Dirichlet boundary conditions. The discretization of the problems by the above methods yields matrices which possess banded structures. This along with the use of fast Fourier transforms makes the cost of one step of each of the Chebyshev spectral collocation methods proportional, except for a logarithmic term, to the number of the unknowns. We present the convergence analysis for the Legendre spectral collocation methods in the special case of the heat equation. Using numerical tests, we demonstrate the second order accuracy in time of the Chebyshev spectral collocation methods for general linear variable coefficient parabolic problems.     

On the use of Gegenbauer reconstructions for shock wave propagation modeling

In therapeutic ultrasound, the presence of shock waves can be significant due to the use of high intensity beams, as well as due to shock formation during inertial cavitation. Although modeling of such strongly nonlinear waves can be carried out using spectral methods, such calculations are typically considered impractical, since accurate calculations often require hundreds or even thousands of harmonics to be considered, leading to prohibitive computational times. Instead, time-domain algorithms which generally utilize Godunov-type finite-difference schemes are commonly used. Although these time domain methods can accurately model steep shock wave fronts, unlike spectral methods they are inherently unsuitable for modeling realistic tissue dispersion relations. Motivated by the need for a more general model, the use of Gegenbauer reconstructions as a postprocess tool to resolve the band-limitations of the spectral methods are investigated. The present work focuses on eliminating the Gibbs phenomenon when representing a steep wave front using a limited number of harmonics. Both plane wave and axisymmetric 2D transducer problems will be presented to characterize the proposed method.     

A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous fluid in 2D

We present a new method for the evolution of inextensible vesicles immersed in a Stokesian fluid. We use a boundary integral formulation for the fluid that results in a set of nonlinear integro-differential equations for the vesicle dynamics. The motion of the vesicles is determined by balancing the non-local hydrodynamic forces with the elastic forces due to bending and tension. Numerical simulations of such vesicle motions are quite challenging. On one hand, explicit time-stepping schemes suffer from a severe stability constraint due to the stiffness related to high-order spatial derivatives and a milder constraint due to a transport-like stability condition. On the other hand, an implicit scheme can be expensive because it requires the solution of a set of nonlinear equations at each time step. We present two semi-implicit schemes that circumvent the severe stability constraints on the time step and whose computational cost per time step is comparable to that of an explicit scheme. We discretize the equations by using a spectral method in space, and a multistep third-order accurate scheme in time. We use the fast multipole method (FMM) to efficiently compute vesicle–vesicle interaction forces in a suspension with a large number of vesicles. We report results from numerical experiments that demonstrate the convergence and algorithmic complexity properties of our scheme.     

Finite-Difference Time-Domain Algorithm for Dispersive Media Based on Runge-Kutta Exponential Time Differencing Method

The electromagnetic propagation in dispersive media is modeled using finite difference time domain (FDTD) method based on the Runge-Kutta exponential time differencing (RKETD) method. The second-order RKETD-FDTD formulation is derived. The high accuracy and efficiency of the presented method is confirmed by computing the transmission and reflection coefficients for a nonmagnetized collision plasma slab in one dimension. The comparison of the numerical results of the RKETD and the exponential time differencing (ETD) algorithm with analytic values indicates that the RKETD is more accurate than the ETD algorithm.     

Upwind-biased FORCE schemes with applications to free-surface shallow flows

In this paper we develop numerical fluxes of the centred type for one-step schemes in conservative form for solving general systems of conservation laws in multiple-space dimensions on structured meshes. The proposed method is an extension of the multidimensional FORCE flux developed by Toro et al. (2009) [14]. Here we introduce upwind bias by modifying the shape of the staggered mesh of the original FORCE method. The upwind bias is evaluated using an estimate of the largest eigenvalue, which in any case is needed for selecting a time step. The resulting basic flux is first-order accurate and monotone. For the linear advection equation, the proposed UFORCE method reproduces exactly the upwind Godunov method. Extension to non-linear systems has been done empirically via the two-dimensional inviscid shallow water equations. Second order of accuracy in space and time on structured meshes is obtained in the framework of finite volume methods. The proposed method improves the accuracy of the solution for small Courant numbers and intermediate waves associated with linearly degenerate fields (contact discontinuities, shear waves and material interfaces). It achieves comparable accuracy to that of upwind methods with approximate Riemann solvers, though retaining the simplicity and efficiency of centred methods. The performance of the schemes is assessed on a suite of test problems for the two-dimensional shallow water equations.     

On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes

We construct uniformly high order accurate discontinuous Galerkin (DG) schemes which preserve positivity of density and pressure for Euler equations of compressible gas dynamics. The same framework also applies to high order accurate finite volume (e.g. essentially non-oscillatory (ENO) or weighted ENO (WENO)) schemes. Motivated by Perthame and Shu (1996) [20] and Zhang and Shu (2010) [26], a general framework, for arbitrary order of accuracy, is established to construct a positivity preserving limiter for the finite volume and DG methods with first order Euler forward time discretization solving one-dimensional compressible Euler equations. The limiter can be proven to maintain high order accuracy and is easy to implement. Strong stability preserving (SSP) high order time discretizations will keep the positivity property. Following the idea in Zhang and Shu (2010) [26], we extend this framework to higher dimensions on rectangular meshes in a straightforward way. Numerical tests for the third order DG method are reported to demonstrate the effectiveness of the methods.