Weak solutions for partial differential equations with source terms: Application to the shallow water equations

Weak solutions of problems with m equations with source terms are proposed using an augmented Riemann solver defined by m + 1 states instead of increasing the number of involved equations. These weak solutions use propagating jump discontinuities connecting the m + 1 states to approximate the Riemann solution. The average of the propagated waves in the computational cell leads to a reinterpretation of the Roe’s approach and in the upwind treatment of the source term of Vázquez-Cendón. It is derived that the numerical scheme can not be formulated evaluating the physical flux function at the position of the initial discontinuities, as usually done in the homogeneous case. Positivity requirements over the values of the intermediate states are the only way to control the global stability of the method. Also it is shown that the definition of well-balanced equilibrium in trivial cases is not sufficient to provide correct results: it is necessary to provide discrete evaluations of the source term that ensure energy dissipating solutions when demanded. The one and two dimensional shallow water equations with source terms due to the bottom topography and friction are presented as case study. The stability region is shown to differ from the one defined for the case without source terms, and it can be derived that the appearance of negative values of the thickness of the water layer in the proximity of the wet/dry front is a particular case, of the wet/wet fronts. The consequence is a severe reduction in the magnitude of the allowable time step size if compared with the one obtained for the homogeneous case. Starting from this result, 1D and 2D numerical schemes are developed for both quadrilateral and triangular grids, enforcing conservation and positivity over the solution, allowing computationally efficient simulations by means of a reconstruction technique for the inner states of the weak solution that allows a recovery of the time step size.     

Adiabatic perturbations for compactons under dissipation and numerically-induced dissipation

Compacton propagation under dissipation shows amplitude damping and the generation of tails. The numerical simulation of compactons by means of dissipative schemes also show the same behaviors. The truncation error terms of a numerical method can be considered as a perturbation of the original partial differential equation and perturbation methods can be applied to its analysis. For dissipative schemes, or when artificial dissipation is added, the adiabatic perturbation method yields evolution equations for the amplitude loss in the numerical solution and the amplitude of the numerically-induced tails. In this paper, such methods are applied to the K(2,2)$K(2\text{,}2)$ Rosenau–Hyman equation, showing a very good agreement between perturbative and numerical results.     

A finite element variational multiscale method for incompressible flows based on two local gauss integrations

In this article, we present a finite element variational multiscale (VMS) method for incompressible flows based on two local Gauss integrations, and compare it with common VMS method which is defined by a low order finite element space L_h$L_h$ on the same grid as X_h$X_h$ for the velocity deformation tensor and a stabilization parameter α$\alpha$. The best algorithmic feature of our method is using two local Gauss integrations to replace projection operator. We theoretically discuss the relationship between our method and common VMS method for the Taylor–Hood elements, and show that the nonlinear system derived from our method by finite element discretization is much smaller than that of common VMS method computationally.     

An energy-conserving Galerkin scheme for a class of nonlinear dispersive equations

A Galerkin scheme is presented for a class of conservative nonlinear dispersive equations, such as the Camassa–Holm equation and the regularized long wave equation. The scheme has two advantageous features: first, it is conservative in that it keeps the discrete analogue of the continuous energy conservation property in the original equations; second, it can be formulated only with cheap H¹$H^1$-elements even if the original equations include third derivative u_xxx$u_{\mathit{xxx}}$. Numerical experiments confirm the stability and effectiveness of the proposed scheme.     

Local derivative post-processing for the discontinuous Galerkin method

Obtaining accurate approximations for derivatives is important for many scientific applications in such areas as fluid mechanics and chemistry as well as in visualization applications. In this paper we discuss techniques for computing accurate approximations of high-order derivatives for discontinuous Galerkin solutions to hyperbolic equations related to these areas. In previous work, improvement in the accuracy of the numerical solution using discontinuous Galerkin methods was obtained through post-processing by convolution with a suitably defined kernel. This post-processing technique was able to improve the order of accuracy of the approximation to the solution of time-dependent symmetric linear hyperbolic partial differential equations from order k+1$k+1$ to order 2k+1$2k+1$ over a uniform mesh; this was extended to include one-sided post-processing as well as post-processing over non-uniform meshes. In this paper, we address the issue of improving the accuracy of approximations to derivatives of the solution by using the method introduced by Thomée [19]. It consists in simply taking the α$\alpha$th-derivative of the convolution of the solution with a sufficiently smooth kernel. The order of convergence of the approximation is then independent   of the order of the derivative, |α|$|\alpha |$. We also discuss an efficient way of computing the approximation which does not involve differentiation but the application of simple finite differencing. Our results show that the above-mentioned approximations to the α$\alpha$th-derivative of the exact solution of linear, multidimensional symmetric hyperbolic systems obtained by the discontinuous Galerkin method with polynomials of degree k$k$ converge with order 2k+1$2k+1$ regardless of the order |α|$|\alpha |$ of the derivative.     

A bubble-stabilized least-squares finite element method for steady MHD duct flow problems at high Hartmann numbers

In this paper we devise a stabilized least-squares finite element method using the residual-free bubbles for solving the governing equations of steady magnetohydrodynamic duct flow. We convert the original system of second-order partial differential equations into a first-order system formulation by introducing two additional variables. Then the least-squares finite element method using C⁰$C^0$ linear elements enriched with the residual-free bubble functions for all unknowns is applied to obtain approximations to the first-order system. The most advantageous features of this approach are that the resulting linear system is symmetric and positive definite, and it is capable of resolving high gradients near the layer regions without refining the mesh. Thus, this approach is possible to obtain approximations consistent with the physical configuration of the problem even for high values of the Hartmann number. Before incoorperating the bubble functions into the global problem, we apply the Galerkin least-squares method to approximate the bubble functions that are exact solutions of the corresponding local problems on elements. Therefore, we indeed introduce a two-level finite element method consisting of a mesh for discretization and a submesh for approximating the computations of the residual-free bubble functions. Numerical results confirming theoretical findings are presented for several examples including the Shercliff problem.     

A high-order low-dispersion symmetry-preserving finite-volume method for compressible flow on curvilinear grids

A new high-order finite-volume method is presented that preserves the skew symmetry of convection for the compressible flow equations. The method is intended for Large-Eddy Simulations (LES) of compressible turbulent flows, in particular in the context of hybrid RANS–LES computations. The method is fourth-order accurate and has low numerical dissipation and dispersion. Due to the finite-volume approach, mass, momentum, and total energy are locally conserved. Furthermore, the skew-symmetry preservation implies that kinetic energy, sound-velocity, and internal energy are all locally conserved by convection as well. The method is unique in that all these properties hold on non-uniform, curvilinear, structured grids. Due to the conservation of kinetic energy, there is no spurious production or dissipation of kinetic energy stemming from the discretization of convection. This enhances the numerical stability and reduces the possible interference of numerical errors with the subgrid-scale model. By minimizing the numerical dispersion, the numerical errors are reduced by an order of magnitude compared to a standard fourth-order finite-volume method.     

Tensor decomposition in electronic structure calculations on 3D Cartesian grids

In this paper, we investigate a novel approach based on the combination of Tucker-type and canonical tensor decomposition techniques for the efficient numerical approximation of functions and operators in electronic structure calculations. In particular, we study applicability of tensor approximations for the numerical solution of Hartree–Fock and Kohn–Sham equations on 3D Cartesian grids. We show that the orthogonal Tucker-type tensor approximation of electron density and Hartree potential of simple molecules leads to low tensor rank representations. This enables an efficient tensor-product convolution scheme for the computation of the Hartree potential using a collocation-type approximation via piecewise constant basis functions on a uniform n×n×n$n\times n\times n$ grid. Combined with the Richardson extrapolation, our approach exhibits O(h³)$O(h^3)$ convergence in the grid-size h=O(n^-1)$h=O(n^{-1})$. Moreover, this requires O(3rn+r³)$O(3\mathit{rn}+r^3)$ storage, where r   denotes the Tucker rank of the electron density with r=O(logn)$r=O(\log n)$, almost uniformly in n  . For example, calculations of the Coulomb matrix and the Hartree–Fock energy for the CH₄${\text{CH}}_4$ molecule, with a pseudopotential on the C atom, achieved accuracies of the order of 10^-6${10}^{-6}$ hartree with a grid-size n of several hundreds. Since the tensor-product convolution in 3D is performed via 1D convolution transforms, our scheme markedly outperforms the 3D-FFT in both the computing time and storage requirements.     

A transformation-free HOC scheme for incompressible viscous flows past an impulsively started circular cylinder

In this paper, we present a higher order compact scheme for the unsteady two-dimensional (2D) Navier–Stokes equations on nonuniform polar grids specifically designed for the incompressible viscous flows past a circular cylinder. The scheme is second order accurate in time and at least third order accurate in space. The scheme very efficiently computes both unsteady and time-marching steady-state flow for a wide range of Reynolds numbers (Re)$(\mathit{Re})$ ranging from 10 to 9500 for the impulsively started cylinder. The robustness of the scheme is highlighted when it accurately captures the vortex shedding for moderate Re   represented by the von Kármán street and the so called α$\alpha$ and β$\beta$-phenomena for higher Re. Comparisons are made with established numerical and experimental results and excellent agreement is found in all the cases, both qualitatively and quantitatively.     

The basic bundle gerbe on unitary groups

We consider the construction of the basic bundle gerbe on SU(n)$SU\left(n\right)$ introduced by Meinrenken and show that it extends to a range of groups with unitary actions on a Hilbert space including U(n)$U\left(n\right)$ and _Up(H)$U_p\left(H\right)$, the Banach Lie group of unitaries differing from the identity by an element of a Schatten ideal. In all these cases we give an explicit connection and curving on the basic bundle gerbe and calculate the real Dixmier–Douady class. Extensive use is made of the holomorphic functional calculus for operators on a Hilbert space.     

Zeta regularized determinant of the Laplacian for classes of spherical space forms

We derive the spectral zeta function in terms of certain Dirichlet series for a variety of spherical space forms _MG$M_G$. Extending results in [C. Nash, D. O’Connor, Determinants of Laplacians on lens spaces, J. Math. Phys. 36 (3) (1995) 1462–1505] the zeta-regularized determinant of the Laplacian on _MG$M_G$ is obtained explicitly from these formulas. In particular, our method applies to manifolds of dimension higher than 3 and it includes the case where G$G$ arises from the dihedral group of order 2m$2m$. As a crucial ingredient in our analysis we determine the dimension of eigenspaces of the Laplacian in form of some combinatorial quantities for various infinite classes of manifolds from the explicit form of the generating function in [A. Ikeda, On the spectrum of a Riemannian manifold of positive constant curvature, Osaka J. Math. 17 (1980) 75–93].     

Very-high-order weno schemes

We study weno(2r − 1) reconstruction [D.S. Balsara, C.W. Shu, Monotonicity prserving weno schemes with increasingly high-order of accuracy, J. Comput. Phys. 160 (2000) 405–452], with the mapping (wenom) procedure of the nonlinear weights [A.K. Henrick, T.D. Aslam, J.M. Powers, Mapped weighted-essentially-non-oscillatory schemes: achieving optimal order near critical points, J. Comput. Phys. 207 (2005) 542–567], which we extend up to weno17 (r=9)$(r=9)$. We find by numerical experiment that these procedures are essentially nonoscillatory without any stringent cfl limitation (cfl∈[0.8,1])$(\text{<small-caps>cfl</small-caps>}\in [0.8\text{,}1])$, for scalar hyperbolic problems (both linear and scalar conservation laws), provided that the exponent p_β$p_{\beta }$ in the definition of the Jiang–Shu [G.S. Jiang, C.W. Shu, Efficient implementation of weighted eno schemes, J. Comput. Phys. 126 (1996) 202–228] nonlinear weights be taken as p_β=r$p_{\beta }=r$, as originally proposed by Liu et al. [X.D. Liu, S. Osher, T. Chan, Weighted essentially nonoscillatory schemes, J. Comput. Phys. 115 (1994) 200–212], instead of p_β=2$p_{\beta }=2$ (this is valid both for weno and wenom reconstructions), although the optimal value of the exponent is probably p_β(r)∈[2,r]$p_{\beta }(r)\in [2\text{,}r]$. Then, we apply the family of very-high-order wenom_pβ=r${\text{<small-caps>wenom</small-caps>}}_{p_{\beta }=r}$ reconstructions to the Euler equations of gasdynamics, by combining local characteristic decomposition [A. Harten, B. Engquist, S. Osher, S.R. Chakravarthy, Uniformly high-order accurate essentially nonoscillatory schemes iii, J. Comput. Phys. 71 (1987) 231–303], with recursive-order-reduction (ror) aiming at aleviating the problems induced by the nonlinear interactions of characteristic fields within the stencil. The proposed ror algorithm, which generalizes the algorithm of Titarev and Toro [V.A. Titarev, E.F. Toro, Finite-volume weno schemes for 3-D conservation laws, J. Comput. Phys. 201 (2004) 238–260], is free of adjustable parameters, and the corresponding rorwenom_pβ=r${\text{<small-caps>rorwenom</small-caps>}}_{p_{\beta }=r}$ schemes are essentially nonoscillatory, as Δx→0$\mathrm{\Delta }x\to 0$, up to r=9$r=9$, for all of the test-cases studied. Finally, the unsplit linewise 2-D extension of the schemes is evaluated for several test-cases.     

Maximum entropy algorithm with inexact upper entropy bound based on Fup basis functions with compact support

The maximum entropy (MaxEnt) principle is a versatile tool for statistical inference of the probability density function (pdf) from its moments as a least-biased estimation among all other possible pdf’s. It maximizes Shannon entropy, satisfying the moment constraints. Thus, the MaxEnt algorithm transforms the original constrained optimization problem to the unconstrained dual optimization problem using Lagrangian multipliers. The Classic Moment Problem (CMP) uses algebraic power moments, causing typical conventional numerical methods to fail for higher-order moments (m>5–10)$(m>5''–''10)$ due to different sensitivities of Lagrangian multipliers and unbalanced nonlinearities. Classic MaxEnt algorithms overcome these difficulties by using orthogonal polynomials, which enable roughly the same sensitivity for all Lagrangian multipliers. In this paper, we employ an idea based on different principles, using Fup_n${\mathit{Fup}}_n$ basis functions with compact support, which can exactly describe algebraic polynomials, but only if the Fup order-n   is greater than or equal to the polynomial’s order. Our algorithm solves the CMP with respect to the moments of only low order Fup₂${\mathit{Fup}}_2$ basis functions, finding a Fup₂${\mathit{Fup}}_2$ optimal pdf with better balanced Lagrangian multipliers. The algorithm is numerically very efficient due to localized properties of Fup₂${\mathit{Fup}}_2$ basis functions implying a weaker dependence between Lagrangian multipliers and faster convergence. Only consequences are an iterative scheme of the algorithm where power moments are a sum of Fup₂${\mathit{Fup}}_2$ and residual moments and an inexact entropy upper bound. However, due to small residual moments, the algorithm converges very quickly as demonstrated on two continuous pdf examples – the beta distribution and a bi-modal pdf, and two discontinuous pdf examples – the step and double Dirac pdf. Finally, these pdf examples present that Fup MaxEnt algorithm yields smaller entropy value than classic MaxEnt algorithm, but differences are very small for all practical engineering purposes.     

A positive MUSCL-Hancock scheme for ideal magnetohydrodynamics

We present a highly robust second order accurate scheme for the Euler equations and the ideal MHD equations. The scheme is of predictor–corrector type, with a MUSCL scheme following as a special case. The crucial ingredients are an entropy stable approximate Riemann solver and a new spatial reconstruction that ensures positivity of mass density and pressure. For multidimensional MHD, a new discrete form of the Powell source terms is vital to ensure the stability properties. The numerical examples show that the scheme has superior stability compared to standard schemes, while maintaining accuracy. In particular, the method can handle very low values of pressure (i.e. low plasma β$\beta$ or high Mach numbers) and low mass densities.     

A general strategy for the optimization of Runge–Kutta schemes for wave propagation phenomena

We analyze optimized explicit Runge–Kutta schemes (RK) for computational aeroacoustics, and wave propagation phenomena in general. Exploiting the analysis developed in [S. Pirozzoli, Performance analysis and optimization of finite-difference schemes for wave propagation problems, J. Comput. Phys. 222 (2007) 809–831], we rigorously evaluate the performance of several time integration schemes in terms of appropriate error and cost metrics, and provide a general strategy to design Runge–Kutta methods tailored for specific applications. We present families of optimized second- and third-order Runge–Kutta schemes with up to seven stages, and describe their implementation in the framework of Williamson’s 2N$2N$-storage formulation [J.H. Williamson, Low-storage Runge–Kutta schemes, J. Comput. Phys. 35 (1980) 48–56]. Numerical simulations of the 1D linear advection equation and of the 2D linearized Euler equations are performed to demonstrate the validity of the theory and to quantify the improvement provided by optimized schemes.     

A first-order system approach for diffusion equation. II: Unification of advection and diffusion

In this paper, we unify advection and diffusion into a single hyperbolic system by extending the first-order system approach introduced for the diffusion equation [J. Comput. Phys., 227 (2007) 315–352] to the advection–diffusion equation. Specifically, we construct a unified hyperbolic advection–diffusion system by expressing the diffusion term as a first-order hyperbolic system and simply adding the advection term to it. Naturally then, we develop upwind schemes for this entire   system; there is thus no need to develop two different schemes, i.e., advection and diffusion schemes. We show that numerical schemes constructed in this way can be automatically uniformly accurate, allow O(h)$O(h)$ time step, and compute the solution gradients (viscous stresses/heat fluxes for the Navier–Stokes equations) simultaneously to the same order of accuracy as the main variable, for all Reynolds numbers. We present numerical results for boundary-layer type problems on non-uniform grids in one dimension and irregular triangular grids in two dimensions to demonstrate various remarkable advantages of the proposed approach. In particular, we show that the schemes solving the first-order advection–diffusion system give a tremendous speed-up in CPU time over traditional scalar schemes despite the additional cost of carrying extra variables and solving equations for them. We conclude the paper with discussions on further developments to come.     

An implicit high-order hybridizable discontinuous Galerkin method for nonlinear convection–diffusion equations

In this paper, we present hybridizable discontinuous Galerkin methods for the numerical solution of steady and time-dependent nonlinear convection–diffusion equations. The methods are devised by expressing the approximate scalar variable and corresponding flux in terms of an approximate trace of the scalar variable and then explicitly enforcing the jump condition of the numerical fluxes across the element boundary. Applying the Newton–Raphson procedure and the hybridization technique, we obtain a global equation system solely in terms of the approximate trace of the scalar variable at every Newton iteration. The high number of globally coupled degrees of freedom in the discontinuous Galerkin approximation is therefore significantly reduced. We then extend the method to time-dependent problems by approximating the time derivative by means of backward difference formulae. When the time-marching method is (p+1)$(p+1)$th order accurate and when polynomials of degree p?0$p?0$ are used to represent the scalar variable, each component of the flux and the approximate trace, we observe that the approximations for the scalar variable and the flux converge with the optimal order of p+1$p+1$ in the L²$L^2$-norm. Finally, we apply element-by-element postprocessing schemes to obtain new approximations of the flux and the scalar variable. The new approximate flux, which has a continuous interelement normal component, is shown to converge with order p+1$p+1$ in the L²$L^2$-norm. The new approximate scalar variable is shown to converge with order p+2$p+2$ in the L²$L^2$-norm. The postprocessing is performed at the element level and is thus much less expensive than the solution procedure. For the time-dependent case, the postprocessing does not need to be applied at each time step but only at the times for which an enhanced solution is required. Extensive numerical results are provided to demonstrate the performance of the present method.     

Finite difference approximations for a fractional advection diffusion problem

The use of the conventional advection diffusion equation in many physical situations has been questioned by many investigators in recent years and alternative diffusion models have been proposed. Fractional space derivatives are used to model anomalous diffusion or dispersion, where a particle plume spreads at a rate inconsistent with the classical Brownian motion model. When a fractional derivative replaces the second derivative in a diffusion or dispersion model, it leads to enhanced diffusion, also called superdiffusion. We consider a one-dimensional advection–diffusion model, where the usual second-order derivative gives place to a fractional derivative of order α$\alpha$, with 1<α?2$1<\alpha ?2$. We derive explicit finite difference schemes which can be seen as generalizations of already existing schemes in the literature for the advection–diffusion equation. We present the order of accuracy of the schemes and in order to show its convergence we prove they are stable under certain conditions. In the end we present a test problem.