The superconvergent property of the Discontinuous Galerkin (DG) method for linear hyperbolic systems of partial differential equations in one dimension is explained by relating the DG method to a particular continuous method, whose accuracy depends in part on a local analysis, and in part on information transferred from upwind elements.

In this paper we generalize a new type of compact Hermite weighted essentially non-oscillatory (HWENO) limiter for the Runge-Kutta discontinuous Galerkin (RKDG) method, which was recently developed in [38] for structured meshes, to two dimensional unstructured meshes. The main idea of this HWENO limiter is to reconstruct the new polynomial by the usage of the entire polynomials of the DG solution from the target cell and its neighboring cells in a least squares fashion [11] while maintaining the conservative property, then use the classical WENO methodology to form a convex combination of these reconstructed polynomials based on the smoothness indicators and associated nonlinear weights. The main advantage of this new HWENO limiter is the robustness for very strong shocks and simplicity in implementation especially for the unstructured meshes considered in this paper, since only information from the target cell and its immediate neighbors is needed. Numerical results for both scalar and system equations are provided to test and verify the good performance of this new limiter.

We propose a general framework of computing interfacial structures between two modulated phases. Specifically we propose to use a computational box consisting of two half spaces, each occupied by a modulated phase with given position and orientation. The boundary conditions and basis functions are chosen to be commensurate with the bulk phases. We observe that the ordered nature of modulated structures stabilizes the interface, which enables us to obtain optimal interfacial structures by searching local minima of the free energy landscape. The framework is applied to the Landau-Brazovskii model to investigate interfaces between modulated phases with different relative positions and orientations. Several types of novel complex interfacial structures emerge from the calculations.

A Lagrangian surface hopping algorithm is implemented to study the two dimensional massless Dirac equation for Graphene with an electrostatic potential, in the semiclassical regime. In this problem, the crossing of the energy levels of the system at Dirac points requires a particular treatment in the algorithm in order to describe the quantum transition—characterized by the Landau-Zener probability— between different energy levels. We first derive the Landau-Zener probability for the underlying problem, then incorporate it into the surface hopping algorithm. We also show that different asymptotic models for this problem derived in [O. Morandi, F. Schurrer, J. Phys. A:Math. Theor. 44 (2011) 265301]may give different transition probabilities. We conduct numerical experiments to compare the solutions to the Dirac equation, the surface hopping algorithm, and the asymptotic models of [O. Morandi, F. Schurrer, J. Phys. A: Math. Theor. 44 (2011) 265301].

An explicit spectrally accurate order-adaptive Hermite-Taylor method for the Schrödinger equation is developed. Numerical experiments illustrating the properties of the method are presented. The method, which is able to use very coarse grids while still retaining high accuracy, compares favorably to an existing exponential integrator – high order summation-by-parts finite difference method.

This paper is concerned with numerical approximations of a nonlocal heat equation define on an infinite domain. Two classes of artificial boundary conditions (ABCs) are designed, namely, nonlocal analog Dirichlet-to-Neumann-type ABCs (global in time) and high-order Padé approximate ABCs (local in time). These ABCs reformulate the original problem into an initial-boundary-value (IBV) problem on a bounded domain. For the global ABCs, we adopt a fast evolution to enhance computational efficiency and reduce memory storage. High order fully discrete schemes, both second-order in time and space, are given to discretize two reduced problems. Extensive numerical experiments are carried out to show the accuracy and efficiency of the proposed methods.

We propose here a class of numerical schemes for the approximation of weak solutions to nonlinear hyperbolic systems in nonconservative form—the notion of solution being understood in the sense of Dal Maso, LeFloch, and Murat (DLM). The proposed numerical method falls within LeFloch-Mishra's framework of schemes with well-controlled dissipation (WCD), recently introduced for dealing with small-scale dependent shocks. We design WCD schemes which are consistent with a given nonconservative system at arbitrarily high-order and then analyze their linear stability. We then investigate several nonconservative hyperbolic models arising in complex fluid dynamics, and we numerically demonstrate the convergence of our schemes toward physically meaningful weak solutions.

A parallel, high-order direct Discontinuous Galerkin (DDG) method has been developed for solving the three dimensional compressible Navier-Stokes equations on 3D hybrid grids. The most distinguishing and attractive feature of DDG method lies in its simplicity in formulation and efficiency in computational cost. The formulation of the DDG discretization for 3D Navier-Stokes equations is detailed studied and the definition of characteristic length is also carefully examined and evaluated based on 3D hybrid grids. Accuracy studies are performed to numerically verify the order of accuracy using flow problems with analytical solutions. The capability in handling curved boundary geometry is also demonstrated. Furthermore, an SPMD (single program, multiple data) programming paradigm based on MPI is proposed to achieve parallelism. The numerical results obtained indicate that the DDG method can achieve the designed order of accuracy and is able to deliver comparable results as the widely used BR2 scheme, clearly demonstrating that the DDG method provides an attractive alternative for solving the 3D compressible Navier-Stokes equations.

Determining the drag of a flowover a rough surface is a guiding example for the need to take geometric micro-scale effects into account when computing a macroscale quantity. A well-known strategy to avoid a prohibitively expensive numerical resolution of micro-scale structures is to capture the micro-scale effects through some effective boundary conditions posed for a problem on a (virtually) smooth domain. The central objective of this paper is to develop a numerical scheme for accurately capturing the micro-scale effects at essentially the cost of twice solving a problem on a (piecewise) smooth domain at affordable resolution. Here and throughout the paper “smooth” means the absence of any micro-scale roughness. Our derivation is based on a “conceptual recipe” formulated first in a simplified setting of boundary value problems under the assumption of sufficient local regularity to permit asymptotic expansions in terms of the micro-scale parameter.

The proposed multiscale model relies then on an upscaling strategy similar in spirit to previous works by Achdou et al. [1], Jäger and Mikelic [29, 31], Friedmann et al. [24, 25], for incompressible fluids. Extensions to compressible fluids, although with several noteworthy distinctions regarding e.g. the “micro-scale size” relative to boundary layer thickness or the systematic treatment of different boundary conditions, are discussed in Deolmi et al. [16,17]. For proof of concept the general strategy is applied to the compressible Navier-Stokes equations to investigate steady, laminar, subsonic flow over a flat plate with partially embedded isotropic and anisotropic periodic roughness imposing adiabatic and isothermal wall conditions, respectively. The results are compared with high resolution direct simulations on a fully resolved rough domain.

The computational work and storage of numerically solving the time fractional PDEs are generally huge for the traditional direct methods since they require total memory and work, where NT and NS represent the total number of time steps and grid points in space, respectively. To overcome this difficulty, we present an efficient algorithm for the evaluation of the Caputo fractional derivative of order α∈(0,1). The algorithm is based on an efficient sum-of-exponentials (SOE) approximation for the kernel t–1–α on the interval [Δt, T] with a uniform absolute error ε. We give the theoretical analysis to show that the number of exponentials Nexp needed is of order for T≫1 or for TH1 for fixed accuracy ε. The resulting algorithm requires only storage and work when numerically solving the time fractional PDEs. Furthermore, we also give the stability and error analysis of the new scheme, and present several numerical examples to demonstrate the performance of our scheme.

Lipid vesicles appear ubiquitously in biological systems. Understanding how the mechanical and intermolecular interactions deform vesicle membranes is a fundamental question in biophysics. In this article we develop a fast algorithm to compute the surface configurations of lipid vesicles by introducing surface harmonic functions to approximate themembrane surface. This parameterization allows an analytical computation of the membrane curvature energy and its gradient for the efficient minimization of the curvature energy using a nonlinear conjugate gradient method. Our approach drastically reduces the degrees of freedom for approximating the membrane surfaces compared to the previously developed finite element and finite difference methods. Vesicle deformations with a reduced volume larger than 0.65 can be well approximated by using as small as 49 surface harmonic functions. The method thus has a great potential to reduce the computational expense of tracking multiple vesicles which deform for their interaction with external fields.

In this paper, we present an efficient rescaling scheme for computing the long-time dynamics of expanding interfaces. The idea is to design an adaptive time-space mapping such that in the new time scale, the interfaces evolves logarithmically fast at early growth stage and exponentially fast at later times. The new spatial scale guarantees the conservation of the area/volume enclosed by the interface. Compared with the original rescaling method in [J. Comput. Phys. 225(1) (2007) 554–567], this adaptive scheme dramatically improves the slow evolution at early times when the size of the interface is small. Our results show that the original three-week computation in [J. Comput. Phys. 225(1) (2007) 554–567] can be reproduced in about one day using the adaptive scheme. We then present the largest and most complicated Hele-Shaw simulation up to date.

Offline computation is an essential component in most multiscale model reduction techniques. However, there are multiscale problems in which offline procedure is insufficient to give accurate representations of solutions, due to the fact that offline computations are typically performed locally and global information is missing in these offline information. To tackle this difficulty, we develop an online local adaptivity technique for local multiscale model reduction problems. We design new online basis functions within Discontinuous Galerkin method based on local residuals and some optimally estimates. The resulting basis functions are able to capture the solution efficiently and accurately, and are added to the approximation iteratively. Moreover, we show that the iterative procedure is convergent with a rate independent of physical scales if the initial space is chosen carefully. Our analysis also gives a guideline on how to choose the initial space. We present some numerical examples to show the performance of the proposed method.

How to effectively simulate the interaction between fluid and solid elements of different sizes remains to be challenging. The discrete element method (DEM) has been used to deal with the interactions between solid elements of various shapes and sizes, while the material point method (MPM) has been developed to handle the multiphase (solid-liquid-gas) interactions involving failure evolution. A combined MPM-DEM procedure is proposed to take advantage of both methods so that the interaction between solid elements and fluid particles in a container could be better simulated. In the proposed procedure, large solid elements are discretized by the DEM, while the fluid motion is computed using the MPM. The contact forces between solid elements and rigid walls are calculated using the DEM. The interaction between solid elements and fluid particles are calculated via an interfacial scheme within the MPM framework. With a focus on the boundary condition effect, the proposed procedure is illustrated by representative examples, which demonstrates its potential for a certain type of engineering problems.

Plasma waves with frequencies close to the particular gyrofrequencies of the charged particles in the plasma lose energy due to cyclotron damping. We briefly discuss the gyro-resonance of low frequency plasma waves and ions particularly with regard to particle-in-cell (PiC) simulations. A setup is outlined which uses artificially excited waves in the damped regime of the wave mode's dispersion relation to track the damping of the wave's electromagnetic fields. Extracting the damping rate directly fromthe field data in real or Fourier space is an intricate and non-trivial task. We therefore present a simple method of obtaining the damping rate Γ from the simulation data. This method is described in detail, focusing on a step-by-step explanation of the course of actions. In a first application to a test simulation we find that the damping rates obtained from this simulation generally are in good agreement with theoretical predictions. We then compare the results of one-, two- and three-dimensional simulation setups and simulations with different physical parameter sets.

A novel mesh deformation technique is developed based on the Delaunay graph mapping method and the inverse distance weighting (IDW) interpolation. The algorithm maintains the advantages of the efficiency of Delaunay graph mapping mesh deformation while it also possesses the ability of better controlling the near surface mesh quality. The Delaunay graph is used to divide the mesh domain into a number of sub-domains. On each sub-domain, the inverse distance weighting interpolation is applied, resulting in a similar efficiency as compared to the fast Delaunay graph mapping method. The paper will show how the near-wall mesh quality is controlled and improved by the new method

We propose and investigate a novel solution strategy to efficiently and accurately compute approximate solutions to semilinear optimal control problems, focusing on the optimal control of phase field formulations of geometric evolution laws. The optimal control of geometric evolution laws arises in a number of applications in fields including material science, image processing, tumour growth and cell motility. Despite this, many open problems remain in the analysis and approximation of such problems. In the current work we focus on a phase field formulation of the optimal control problem, hence exploiting the well developed mathematical theory for the optimal control of semilinear parabolic partial differential equations. Approximation of the resulting optimal control problemis computationally challenging, requiring massive amounts of computational time and memory storage. The main focus of this work is to propose, derive, implement and test an efficient solution method for such problems. The solver for the discretised partial differential equations is based upon a geometric multigrid method incorporating advanced techniques to deal with the nonlinearities in the problem and utilising adaptive mesh refinement. An in-house two-grid solution strategy for the forward and adjoint problems, that significantly reduces memory requirements and CPU time, is proposed and investigated computationally. Furthermore, parallelisation as well as an adaptive-step gradient update for the control are employed to further improve efficiency. Along with a detailed description of our proposed solution method together with its implementation we present a number of computational results that demonstrate and evaluate our algorithms with respect to accuracy and efficiency. A highlight of the present work is simulation results on the optimal control of phase field formulations of geometric evolution laws in 3-D which would be computationally infeasible without the solution strategies proposed in the present work.

In this paper we consider a discontinuous Galerkin discretization of the ideal magnetohydrodynamics (MHD) equations on unstructured meshes, and the divergence free constraint (∇·B=0) of its magnetic field B. We first present two approaches for maintaining the divergence free constraint, namely the approach of a locally divergence free projection inspired by locally divergence free elements [19], and another approach of the divergence cleaning technique given by Dedner et al. [15]. By combining these two approaches we obtain a efficient method at the almost same numerical cost. Finally, numerical experiments are performed to show the capacity and efficiency of the scheme.

Moment models are often used for the solution of kinetic equations such as the Boltzmann equation. Unfortunately, standard models like Grad's equations are not hyperbolic and can lead to nonphysical solutions. Newly derived moment models like the Hyperbolic Moment Equations and the Quadrature-Based Moment Equations yield globally hyperbolic equations but are given in partially conservative form that cannot be written as a conservative system.

In this paper we investigate the applicability of different dedicated numerical schemes to solve the partially conservative model equations. Caused by the non-conservative type of equation we obtain differences in the numerical solutions, but due to the structure of the moment systems we show that these effects are very small for standard simulation cases. After successful identification of useful numerical settings we show a convergence study for a shock tube problem and compare the results to a discrete velocity solution. The results are in good agreement with the reference solution and we see convergence considering an increasing number of moments.

The lattice Boltzmann method is employed to simulate the steady flow in a two-dimensional lid-driven semi-elliptical cavity. Reynolds number (Re) and vertical-to-horizontal semi-axis ratio (D) are in the range of 500-5000 and 0.1-4, respectively. The effects of Re and D on the vortex structure and pressure field are investigated, and the evolutionary features of the vortex structure with Re and D are analyzed in detail. Simulation results show that the vortex structure and its evolutionary features significantly depend on Re and D. The steady flow is characterized by one vortex in the semi-elliptical cavity when both Re and D are small. As Re increases, the appearance of the vortex structure becomes more complex. When D is less than 1, increasing D makes the large vortexes more round, and the evolution of the vortexes with D becomes more complex with increasing Re. When D is greater than 1, the steady flow consists of a series of large vortexes which superimpose on each other. As Re and D increase, the number of the large vortexes increases. Additionally, a small vortex in the upper-left corner of the semi-elliptical cavity appears at a large Re and its size increases slowly as Re increases. The highest pressures appear in the upper-right corner and the pressure changes drastically in the upper-right region of the cavity. The total pressure differences in the semi-elliptical cavity with a fixed D decrease with increasing Re. In the region of themain vortex, the pressure contours nearly coincide with the streamlines, especially for the cavity flow with a large Re.