April 16, 2014

April 20-22, 2009
Workshop on Computational Methods for Hyperbolic Problems
University of Waterloo, Department of Applied Mathematics



Timothy Barth, NASA Ames Research Center
Error Representation for Time Dependent Compressible Flow Calculations

For better or worse, our physical world is constantly evolving in time. Many important physical phenomena depend fundamentally on time either deterministically or through dynamical system behavior. These leads to a multitude of challenges associated with identifying, quantifying, and controlling numerical errors in complex time dependent numerical simulations. For example, in turbulent flow simulations it is well known that the control of pointwise solution errors quickly becomes intractable as the flow Reynolds number increases but the control of errors occurring in statistics and space-time averaged quantities may still be tractable [1].

In this presentation, we consider the representation and control of numerical solution errors in space-time FE methods using standard duality techniques as succinctly described in [2,3]. For computed output quantities of interest that are mathematically described as functionals, this technique exploits the precise relationship between functional errors and weighted combinations of computable element residuals. The associated non-perturbative theory provides the exact form of these residual weights and elucidates via Galerkin orthogonality why these weights depend not on the dual solution but rather on the difference between the dual solution and any projection of the dual solution into the primal approximation space. This latter property often leads to a non-intuitive space-time dependence of the local solution error on the accuracy of computed output functionals. Finally, we briefly mention recent work in adapting established error control techniques for ODEs to a space-time FEM setting. In this approach, the normed-error in each time slab is adaptively made smaller than some prescribed tolerance before the next time slab interval is solved. If the time slab error tolerances are suitably chosen, then error control after N time steps is achieved. The presentation of numerical results will hopefully accentuate the challenges and difficulties in space-time error representation as well as stimulate fruitful discussions addressing the feasibility of genuine error control for time dependent fluid flow problems.

[1] J. Hoffman and C. Johnson,"Adaptive Finite Element Methods in
Incompressible Fluid Flow", LNCSE, Vol. 25, Springer-Verlag Pub, 2002.

[2] K. Ericksson, D. Estep, P. Hansbo and C. Johnson,
Computational Differential Equations, Cambridge Press., 1996.

[3] R. Becker and R. Rannacher,"An Optimal Control Approach to
A-Posteriori Error Estimation in Finite Element Methods, Acta
Numerica, 2001.

Marsha Berger, New York University
Status Report on Embedded Boundary Grids

The Cartesian grid embedded boundary approach has attracted much interest in the last decade due to the ease of grid generation for
complicated geometries. This approach uses rectangular Cartesian meshes over most of the domain, with irregular (cut) cells only at the
intersection of the mesh with the boundary of a solid object. In this talk we first briefly survey our approach to embedded boundary computations, and describe what distinguishes it from level set or other immersed boundary approaches. We then discuss some of the
algorithmic issues that arise in this approach. In particular we concentrate on numerical discretizations that avoid loss of accuracy and
stability at irregular boundary cells. We end with some computations of 3D flows with collaborators at NASA Ames Research Center.

Martin Berzins, University of Utah
Nonlinear High-Order Space and Time Discretization Methods

The analysis of high order schemes based on ENO and WENO methods will be considered. It will be shown that the use of high-order provably data bounded polynomials, based on extensions of the work of [Berzins SIAM Review December 07], provides a way to develop positivity preserving (W)ENO methods of potentially very high orders. A series of numerical experiments are used to determine optimal orders for typical solution profiles. A similar idea is then used to consider variable-order positivity preserving methods for time integration, based on the variable-step variable-order backward differentiation code DASSL.

Hans De Sterck, University of Waterloo
Efficient solution of stationary Euler flows with critical points and shocks

It is well-known that stationary transonic solutions of the compressible Euler equations are hard to compute using the stationary form of the equations. Therefore, time marching methods with explicit or implicit time integration are normally employed. However, the computational complexity of the time marching approach is far from optimal, because convergence tends to be slow and tends to slow down even more as resolution increases. In this talk we explore the alternative of solving the stationary equations directly, which is a viable approach when the solution topology is known in advance. We first present a solution method for one-dimensional flows with critical points. The method is based on a dynamical systems formulation of the problem and uses adaptive integration combined with a two-by-two Newton shooting method. Example calculations show that the resulting method is fast and accurate. A sample application area for this method is the calculation of transonic hydrodynamic escape flows from extrasolar planets and the early Earth, and the method is also illustrated for quasi-one-dimensional nozzle flow and black hole accretion. The method can be extended easily to handle flows with shocks, using a Newton method applied to the Rankine-Hugoniot relations. Extension to flows with heat conduction is also discussed. The presentation will conclude with some thoughts on how the approach presented can be generalized to problems in higher dimensions.

Clinton P. T. Groth, University of Toronto
Numerical Solution of Continuum and Non-Equilibrium Flows Using Hyperbolic and Physically-Realizable Descriptions Based on Moment Closures

The numerical solution of continuum and non-equilibrium flows by using fully hyperbolic and realizable mathematical descriptions which follow from a hierarchy of moment closures is described. A somewhat novel hierarchy of both physically realizable and strictly hyperbolic moment closures is considered which follows from a slight modification to the more usual maximum-entropy closure hierarchies. An
arbitrary-Lagrangian-Eulerian (ALE) parallel finite-volume scheme with adaptive mesh refinement (AMR) and Riemann-solver based numerical flux function is also described for solving the system of partial differential equations arising from the closures on multi-block meshes with embedded and possibly moving boundaries. The purely hyperbolic nature of the moment equations is shown to be particularly insensitive and therefore well suited to discretizations involving grids with irregularities arising from embedded boundary and moving interface treatments. The capabilities of the proposed hyperbolic mathematical treatment and finite-volume scheme for predicting both continuum and non-equilibrium flows are explored by considering the a number of canonical problems, including one-dimensional shock structure flows, Couette flow, flat plate flows, flow past a circular cylinder, as well as transonic flow past a micro airfoil. The talk will conclude with a brief summary of findings and a discussion of the potential of the proposed hyperbolic descriptions for application to a wider class of flow problems.

Jae-Hun Jung, SUNY Buffalo
Spectral collocation methods for hyperbolic equations with singular sources

Hyperbolic equations with singular sources are found in many physical and engineering applications. The solution of such
equations contains the local jump discontinuity which is moving in time in general. Consequently its spectral approximation yields
the Gibbs phenomenon, that is, the solution is highly oscillatory and the error does not converge in the maximum error sense.
Regularizing singular sources reduces the Gibbs oscillations but the maximum error does not converge yet.
To minimize the Gibbs oscillations and improve convergence, we consider a consistent approach to the singular source terms in
the spectral approximation. Singular sources represented as the sum of a Dirac delta function and its derivative(s) are obtained
by the direct projection of the Heaviside functions on the collocation points. The direct projection is obtained by the
collocative derivative of the Heaviside functions and the singular source terms on the collocation points are highly oscillatory.
This approach, however, yields a convergent result for some hyperbolic equations with spectral accuracy due to some cancellations
of the Gibbs oscillations on the collocation points.
Several hyperbolic equations with singular sources are considered including the wave equations from the collision of two
black holes and a nonlinear Schr\"{o}dinger equation in homogeneous medium with defects. These equations contain the delta
function type sources which are either moving or static. For the moving sources, the given PDEs are also reduced into ODEs.
Numerical examples will be presented to verify the efficiency and accuracy of the proposed method.

Smadar Karni, University of Michigan
A Hybrid Scheme for the Baer-Nunziato Two-Phase Flow Model

The Baer-Nunziato two phase flow model describes flame propagation in gas-permeable reactive granular materials. We focus on the hydrodynamic part of the system and neglect the terms due to combustion processes, drag and heat transfer. The system is an averaged flow model, expressing conservation of mass, and momentum and energy balance of the gas and solid phases, plus an additional 'compaction' equation describing the evolution of porosity. The system is only conditionally hyperbolic, and may fail to have a complete set of eigenvectors. It is also in non-conservation form due to momentum and energy exchange between the phases. The presence of non-conservative terms has major consequences both theoretically and computationally. If the porosities are piecewise constant, the phases
'talk' to each other only across the porosity jump. Computing solutions of the Riemann problem rest on capturing the jump conditions across the porosity jump, but even for state-of-the-art numerical methods this may prove a difficult task.

We have formulated the Baer-Nunziato model in terms of the Riemann Invariants across the compaction wave, and propose a hybrid algorithm that uses the Riemann Invariants formulation across the compaction wave, and the conservative formulation away from the compaction wave. The talk will describe the hybrid scheme and present numerical results.

This is joint work with Gerardo Hernandez of the University of Michigan.

Barbara Keyfitz, Ohio State University
Linear and Nonlinear Stability of Shocks

Abstract: Shocks (nonlinear discontinuities) separating a region of hyperbolic states from a non-hyperbolic region can occur in one of two ways in conservation laws. Steady transonic flow is a well-known phenomenon, and transonic shocks appear to be both physically and mathematically stable. Unsteady non-hyperbolic systems, on the other hand, would seem to form ill-posed, even catastrophically ill-posed problems. In this talk, I will review some work that Karen Ames and I did that suggested a sense in which shocks in these ill-posed systems may enjoy a type of stability. Later work with Milton Lopes has confirmed our initial, linear results. These results can be compared to more recent work with a number of co-authors on construction of transonic shocks.

Lilia Krivodonova, University of Waterloo
A space and time adaptive discontinuous Galerkin method

Adaptive schemes provide computational advantages over uniform mesh computations by allocating resources in regions where they are most needed. The computational mesh can be refined ($h$-refinement) or the order of the scheme can be modified ($p$-refinement) to better resolve the solution. As a result, mesh sizes and local CFL numbers can vary greatly throughout the computational domain. This can be inefficient as a few small cells impose a restrictive time step on the whole mesh. We will present an adaptive high-order time integration scheme for solving partial differential equations with the method of lines and explicit Runge-Kutta integrators. A time step on each cell will be defined by the local CFL condition. That is, adaptivity in time will consist of taking a time step $dt$ on cells of size $h$, two $dt/2$ steps on cells of size $h/2$,and so on. The total number of function evaluations will be, therefore, reduced. The interface conditions will be imposed on interelemental boundaries so that the order of the underlying Runge-Kutta scheme is preserved. The novelty of this algorithm is its small stencil. This makes it suitable for unstructured grids, where multi-layer reconstructions are difficult. We will show how this
algorithm can be combined with adaptive mesh-refinement to achieve a significant speed-up for solution of transient flow problems.
Computational issues such as an efficient implementation of the scheme, integration of spatial and temporal refinement strategies, and limiting across the interfaces will be discussed.

Marc Laforest, Montreal Polytechnic
An adaptive version of the Illner-Rjasanow stochastic scheme for the Boltzmann equation

To model the dynamics of materials far from thermodynamic equilibrium, a kinetic model is often required, like the Boltzmann equation for dilute gases. Deterministic schemes for the solution of the Boltzmann equation are typically efficient and accurate but quite limited with respect to the collision operator that can be modeled. On the other hand, stochastic schemes, despite their slow convergence, are much simpler to implement and adapt for even the most complex collision kernels involving multiple complex species with long-range interactions. The key issue for stochastic schemes is to reduce the variance of the observed statistics while keeping the number of degrees of freedom under control.

The original Random Discrete Velocity Model (RDVM) of Illner and Rjasanow became the forerunner of what is now the most effective stochastic scheme for the Boltzmann equation, namely the Stochastic Weighted Particle Method of Rjasanow and Wagner. The RDVM scheme works by using a discrete set of velocities in each cell of a structured mesh. During the transport step, particles can be carried using any high-order transport scheme. During the collision step, the velocities are randomly subdivided into several small Discrete Velocity Models (DVM) which determine how these subsets interact. In this talk, we describe a new a posteriori error estimator for hierarchies of DVMs and describe how it can be used it to determine the size of the local discretization of velocity space in each cell. Using numerical experiments, we demonstrate the accuracy of the error estimator and it's efficiency at adapting the local discretization of velocity space.

Joint work with Kondo Assi

Emmanuel Lorin, UOIT
About the reservoir technique for hyperbolic conservation laws

The talk is devoted to the reservoir technique coupled with finite volume flux schemes for solving with low numerical diffusion nonlinear hyperbolic conservation laws. I will present some convergence results, simulations in 1d and 2d, and some results regarding the overall algorithmic complexity of the method"

Carl Ollivier-Gooch, University of British Columbia
High-Order Accurate Unstructured Mesh Finite-Volume Schemes for Computational Aerodynamics

High-order finite-volume methods for unstructured mesh CFD are reasonably mature as research tools, with demonstrable accuracy and efficiency benefits compared with second-order methods. This maturity comes as a result of recent advances in limiting, convergence acceleration, and adaptive techniques as applied to high-order methods, which will be the focus of this talk.

Steve Ruuth, Simon Frazer University
A Simple Technique for Solving Partial Differential Equations on Surfaces

Many applications require the solution of time-dependent partial differential equations (PDEs) on surfaces or more general manifolds. Methods for treating such problems include surface parameterization, methods on triangulated surfaces and embedding techniques. This talk describes an embedding approach based on the closest point representation of the surface and describes some of its advantages over other embedding methods. Noteworthy features of the method are its generality with respect to the underlying surface and its simplicity. In particular, the method requires only minimal changes to the corresponding three-dimensional codes to treat the evolution of PDEs on surfaces.

Chi-Wang Shu, Brown University
Superconvergence and Time Evolution of Discontinuous Galerkin Finite Element Solutions

In this talk, we study the convergence and time evolution of the error between the discontinuous Galerkin (DG) finite element solution
and the exact solution for conservation laws when upwind fluxes are used. We prove that if we apply piecewise polynomials to a one
dimensional scalar linear equation or a hyperbolic linear system, the DG solution will be superconvergent towards a particular projection of the exact solution. Thus, the error of the DG scheme will not grow for fine grids over a long time period (for $t$ up to $O( \frac{1}{\sqrt{h}} )$ where $h$ is the mesh size). We prove this result for general meshes and arbitrary polynomial degree $k$, and give numerical examples to demonstrate the superconvergence property, as well as the long time behavior of the error for more general cases including nonlinear equations and two-dimensional problems. Generalizations to local discontinuous Galerkin (LDG) method for convection diffusion equations will also be given. This is a joint work with Yingda Cheng.

Marek Stastna, University of Waterloo
Application of hyperbolic system methods to dispersive, nonlinear internal waves and porous media acoustics

I will begin by describing the manner in which numerical methods designed for hyperbolic systems improve the simulations of internal waves in an incompressible, density stratified fluid, especially in the high Reynolds number limit. I will pay particular attention to situations in which density overturns occur for which limiting has been termed "implicit Large Eddy Simulation (LES)". The relationship between methods can be turned around, and results for fully nonlinear internal waves can be used to improve the numerical simulation of classical, weakly nonlinear corrections to hyperbolic waves, such as the Korteweg de Vries equation. Finally I will consider porous media acoustics for which the corrections to hyperbolic systems are non-conservative and consist of absorbing terms. I will survey some of the techniques available for the numerical solution of such systems and demonstrate that classical, or Biot, frequency dependent absorption is non-causal.

Allen M. Tesdall, CUNY
A two-phase Stefan problem for the unsteady transonic small disturbance equations

Numerical solutions of weak shock reflection problems for the unsteady transonic small disturbance equations, the nonlinear wave system, and the compressible Euler equations at a set of parameter values for which regular reflection is impossible contain a complex structure. Instead of a mathematically inadmissible Mach reflection, as is apparently observed in experiments, the solutions contain a cascade of triple points and tiny supersonic patches behind a leading triple point. A centered expansion wave originates at each triple point. We call this structure of repeating supersonic patches and triple points ``Guderley Mach reflection,'' or GMR.

At the upstream side of each patch in GMR, a sonic line separates the patch from a region of subsonic flow. This sonic line can be considered a free boundary in the formulation of a free boundary problem, with the states on either side coupled through the boundary. As a step towards the goal of formulating this free boundary problem, we present a problem which retains its main features, but which is simpler. We choose the simplest model of weak shock reflection, the unsteady transonic small disturbance equations (UTSDE). We choose initial data which results in an expansion wave that reflects off a sonic line, similar to the reflection of a centered rarefaction off the sonic line in a single supersonic patch in GMR. In both GMR and our simpler problem for the UTSDE, the states on either side of the free boundary are coupled through the boundary, and the solutions in both the supersonic and subsonic regions are a priori unknown. These features make these problems analogous to two-phase Stefan problems for the heat equation. At the moment, we have linearized the simplified problem and solved it exactly. We have not yet formulated the free boundary problem, but we have solved the full nonlinear simplified problem numerically. We will present and describe our solutions and our solution method.

Tim Warburton, Rice University
Advances in Wave Propagation with the Discontinuous Galerkin Method

A range of important features relating to the practical application of discontinuous Galerkin (DG) method for wave propagation will be

Given the suitability of DG for solving Maxwell's equations and their ability to propagate waves over long distance, it is natural to seek
effective boundary treatments for artificial radiation boundary conditions. A new family of far field boundary conditions will be introduced which gracefully transmit propagating and evanescent components out of the domain. These conditions are specifically formulated with DG discretizations in mind, however they are also relevant for a range of numerical methods.

There is an Achilles heel to high order discontinuous Galerkin methods when applied to conservation laws. The methods are typically
constructed with polynomial field representations and unfortunatelythese suffer from excess maximum gradients near the edges of elements. I will describe a simple filtering process that allows us to reduce these anomalous gradients and provably yield a dramatic increase in the maximum allowable time step. Additional experiments with local time stepping methods will also be presented.

Finally, I will discuss the use of GPU hardware to accelerate computation for time-domain electromagnetics simulations.