NASA/TM–2014–218175

First-Order Hyperbolic System

Method for Time-Dependent

Advection-Diffusion Problems

Alireza Mazaheri

Langley Research Center, Hampton, Virginia

Hiroaki Nishikawa

National Institute of Aerospace, Hampton, Virginia

March 2014
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. 

It collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI 

program provides access to the NASA Aeronautics 

and Space Database and its public interface, the 

NASA Technical Report Server, thus providing one 

of the largest collections of aeronautical and space 

science STI in the world. Results are published in 

both non-NASA channels and by NASA in the 

NASA STI Report Series, which includes the 

following report types: 

 TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

Programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but 

having less stringent limitations on manuscript 

length and extent of graphic presentations. 

 TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

 CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

 CONFERENCE PUBLICATION. 

Collected papers from scientific and 

technical conferences, symposia, seminars, 

or other meetings sponsored or cosponsored by NASA. 

 SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, 

often concerned with subjects having 

substantial public interest. 

 TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission. 

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and feeds, 

providing information desk and personal search 

support, and enabling data exchange services. 

For more information about the NASA STI 

program, see the following: 

 Access the NASA STI program home page 

at http://www.sti.nasa.gov

 E-mail your question to help@sti.nasa.gov

 Fax your question to the NASA STI 

Information Desk at 443-757-5803 

 Phone the NASA STI Information Desk at 

443-757-5802 

 Write to: 

 STI Information Desk 

 NASA Center for AeroSpace Information 

 7115 Standard Drive 

 Hanover, MD 21076-1320 
NASA/TM–2014–218175

First-Order Hyperbolic System

Method for Time-Dependent

Advection-Diffusion Problems

Alireza Mazaheri

Langley Research Center, Hampton, Virginia

Hiroaki Nishikawa

National Institute of Aerospace, Hampton, Virginia

National Aeronautics and

Space Administration

Langley Research Center

Hampton, Virbinia 23681-2199

March 2014
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 

443-757-5802 

The use of trademarks or names of manufacturers in this report is for accurate reporting and does not 

constitute an official endorsement, either expressed or implied, of such products or manufacturers by the 

National Aeronautics and Space Administration. 
Abstract

A time-dependent extension of the first-order hyperbolic system method [J.

Comput. Phys., 227 (2007) 315-352] for advection-diffusion problems is introduced. Diffusive/viscous terms are written and discretized as a hyperbolic

system, which recovers the original equation in the steady state. The resulting

scheme offers advantages over traditional schemes: a dramatic simplification in

the discretization, high-order accuracy in the solution gradients, and orders-ofmagnitude convergence acceleration. The hyperbolic advection-diffusion system is discretized by the second-order upwind residual-distribution scheme in

a unified manner, and the system of implicit-residual-equations is solved by

Newton’s method over every physical time step. The numerical results are

presented for linear and nonlinear advection-diffusion problems, demonstrating solutions and gradients produced to the same order of accuracy, with rapid

convergence over each physical time step, typically less than five Newton iterations.

1
Contents

1 Introduction 3

2 Hyperbolic Advection-Diffusion System 4

3 Discretization, Steady Solver, and Boundary Condition 6

3.1 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3.2 Steady Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3.3 Boundary Condition . . . . . . . . . . . . . . . . . . . . . . . 10

4 Extension to Time-Dependent Problems 11

5 Extension to Nonlinear Advection-Diffusion System 13

6 Results 17

6.1 Steady Linear Advection-Diffusion (High Re Applications) . . 17

6.2 Unsteady Linear Advection-Diffusion . . . . . . . . . . . . . . 19

6.3 Unsteady Nonlinear Advection-Diffusion . . . . . . . . . . . . 23

6.4 Accuracy Verification . . . . . . . . . . . . . . . . . . . . . . . 25

7 Conclusions 28

2
1 Introduction

In this paper, we introduce an unsteady hyperbolic advection-diffusion scheme.

The first-order hyperbolic system method was initially proposed in Ref. [1]

as a radical approach to steady diffusion problems: discretize an equivalent

first-order hyperbolic system to solve the diffusion equation. It has a number

of attractive features such as the discretization of high-order derivatives by

methods for hyperbolic systems, equal order of accuracy for the solution and

the gradients (viscous/heat fluxes), and orders of magnitude acceleration in

steady convergence. These advantages have been demonstrated for the diffusion equation [1] and the advection-diffusion equation [2] by the second-order

residual-distribution method, and for the compressible Navier-Stokes equations [3] by the second-order finite-volume method, and for the advectiondiffusion equation by first, second and third order finite-volume methods [4,5].

These lines of development, however, have been restricted exclusively to steady

state problems. Towards enabling accurate practical time-dependent computations by the first-order hyperbolic system method, this paper presents the first

study on the unsteady extension of the first-order hyperbolic system method.

We demonstrate that time-accurate schemes constructed via the hyperbolic

method produce high-order accurate solutions and derivatives at every time

step and allow the construction of a highly efficient inner solver for implicit

time stepping.

Time-dependent computations are possible by implicit time-integration

methods such as the second-order backward-difference scheme, which is widely

used for practical computations. To perform the implicit time integration, it

is required to solve a system of global time-dependent residual equations over

each time step. We employ a steady solver developed by the first-order hyperbolic system method to solve the system efficiently and accurately. The

spatial discretization is performed by the residual-distribution method as in

Refs. [1, 2]. Taking advantage of the compactness of the residual-distribution

schemes, we develop a fully-implicit steady solver with exact linearization, i.e.,

Newton’s method, for a second-order upwind scheme. The strategy taken here

may be thought of as a dual-time formulation [6]. However, we introduce

a pseudo-time just for convenience of discretization: the first-order system

is hyperbolic in the pseudo-time and thus will be discretized by the upwind

residual-distribution scheme. Once the spatial discretization is complete, the

pseudo-time step is taken to be infinity to define the system of time-dependent

residual equations, which is then solved by Newton’s method. In this paper,

we consider one-dimensional linear and nonlinear advection-diffusion equations. These examples are sufficient to illustrate the general methodology

for enabling time-dependent computations by the hyperbolic method; this is

the main goal of this paper. The introduction of this unsteady hyperbolic

advection-diffusion scheme along with its great advantages is served as a foun3
dation for the extension of the scheme to higher-dimensions and more complex

systems.

It is emphasized also that there are real practical applications that onedimensional analyses are sufficient and routinely being used in industries (e.g.,

NASA). For instance, in-depth material thermal response calculations of thermal protection systems of atmospheric entry vehicles, where advection-diffusion

equation is solved through the material, is a prime example of such applications [7–9]. Another example is the experimental aeroheating data reduction

that are obtained with global phosphor thermography techniques [10, 11]. Efficient and accurate one-dimensional schemes as presented in this paper will

significantly improve these analyses and could potentially make such calculations part of the routine aerothermodynamic environment predictions. The

present paper serves, therefore, also as an important foundation for the extension to more practical one-dimensional problems.

There exist other methods for enabling time-accurate computations by

the residual-distribution method: the Crank-Nicolson method [12], the spacetime method [13–15], the backward-difference method [16,17], and the explicit

Runge-Kutta-type method [18]. These methods are equally applicable to the

construction of time-accurate schemes in the hyperbolic method. In this paper, we employ the backward-difference method for simplicity but without losing practical applicability. The main difference between our time-integration

method and those in Refs. [16, 17] lies in the method for solving the system of residual equations arising from the implicit time integration scheme:

we employ Newton’s method whereas Ref. [16] employs the dual-time stepping method with a pseudo-time marching iteration and Ref. [17] employs the

multigrid method with a point-implicit iteration.

The paper is organized as follows. In the next section, the hyperbolic

advection-diffusion system is described. In Section 3, a compact second-order

residual-distribution scheme, a steady solver, and the boundary conditions are

discussed. In Section 4, the construction of time-accurate schemes is given.

In Section 5, the time-accurate scheme is extended to a nonlinear advectiondiffusion equation. In Section 6, numerical results are presented. Finally,

Section 7 concludes the study with remarks.

2 Hyperbolic Advection-Diffusion System

Consider the one-dimensional (1-D) advection-diffusion equation,

∂tu + a ∂xu = ν ∂xxu, (1)

where a and ν are both taken to be positive constant. We will follow the

procedure described in Ref. [2] and rewrite the above equation with a first4
order advection-diffusion system:

∂tu = −a ∂xu + ν ∂xp, (2)

∂tp = (∂xu − p)/Tr, (3)

where the relaxation time, Tr > 0, is arbitrary. Towards the steady state,

the variable p approaches the solution gradient and hence the above equation

becomes identical to the original advection-diffusion equation in the steady

state. In the vector form, the first-order advection-diffusion system can be

written as

∂tU + A∂xU = S, (4)

where

U =



u

p



, A =



a −ν

−1/Tr 0



, S =



0

−p/Tr



. (5)

The above system is hyperbolic in time because it has the following real eigenvalues for any positive Tr:

λ1 =

a

2



1 −

r

1 +

4ν

a

2Tr



, λ2 =

a

2



1 + r1 +

4ν

a

2Tr



, (6)

with linearly independent eigenvectors [1].

Note that the above hyperbolic formulation is equivalent to the original

Eq. (1) only in the steady state. The time derivatives in the hyperbolic system, at least ∂tp, should therefore be taken as pseudo-time derivatives. They

serve mainly as a guide for discretization by making the whole system hyperbolic in time for which various discretization techniques are available, e.g.,

upwinding. The benefits are, however, much more than just the convenience in

discretization. First, the hyperbolic discretization typically results in a strong

coupling between the two variables and achieves the equal order of accuracy

for u and p on arbitrary grids. Second, the resulting explicit numerical scheme

is stable with an O(h) time step through the diffusion limit or O(1/h) condition number in the residual Jacobian for implicit solvers. It has been shown

to yield O(1/h) acceleration in convergence over traditional methods for the

diffusion equation [1,4], for the advection-diffusion equation [2,5], and for the

compressible Navier-Stokes equations [3]. Unsteady computations are possible

by implicit time stepping, which is the main subject of this paper.

In the next section, we first fully describe the advection-diffusion equation

discretization process, the implicit steady state formulation, and the boundary

condition implementation in the steady state limit. We then extend the scheme

to time-accurate and nonlinear advection-diffusion equation in the sections

that follow.

5
i-1 i i+1

hL hR

Figure 1: Schematic of grid spacing for a 1-D grid.

3 Discretization, Steady Solver, and Boundary Condition

3.1 Discretization

To discretize the first-order hyperbolic advection-diffusion system, we employ

Residual-Distribution (RD) method. The method consists of two steps; 1)

residual evaluation over the cells (or elements), and 2) distribution of the

residuals to the nodes.

Consider a one-dimensional domain discretized with N nodes that are distributed arbitrarily over the domain of interest with the solution, u, and the

solution gradient, p, data stored at the nodes denoted by xi, i = 1, 2, 3, ..., N

(Fig.1). We define the cell-residual ΦC by integrating the spatial part of the

Eq. (4) over the cell, C, defined by the nodes i and i + 1:

Φ

C =

Z xi+1

xi

(−A∂xU + S)dx =







−a(ui+1 − ui) + ν(pi+1 − pi)

1

Tr



ui+1 − ui −

pi+1 + pi

2

(xi+1 − xi)

 .

(7)

Note that the source term integration has been evaluated by the trapezoidal

rule to ensure second-order accuracy. It should be noted also that the flux balance terms e.g., ux, has been integrated exactly, which is not possible in higher

dimensions and a high-order quadrature would be required to achieve the equal

order of accuracy for the solution and the gradients in the residual-distribution

method as described in Ref. [2]. We again remark that the relaxation time,

Tr, is arbitrary because in the steady state limit p = ux for any Tr. We complete the evaluation of the cell residual with the definition of Tr. Using the

dimensional analysis [2], this is defined as

Tr =

Lr

max(λ1, λ2)

=

Lr

a

2



1 + r

1 +

4ν

a

2Tr

, (8)

where Lr is a length scale that may be optimized using Fourier analysis to

enhance the convergence [19]. The effect of the optimized length scale remains

to be investigated. Here we use the value recommended in Ref. [1, 4, 5]:

Lr =

1

2π

. (9)

6
We now solve for Tr from Eq. (8):

Tr =

Lr

a + ν/Lr

. (10)

Substituting back the relaxation time into Eq. (6), we express the eigenvalues

as

λ1 = −ν/Lr, λ2 = a + ν/Lr. (11)

At this point the residual evaluation is completed. We now define the

distribution of the residuals. We achieve this by first diagonalizing the matrix

A with the following right and left eigenvectors:

R =



−λ1Tr −λ2Tr

1 1 

, L =

aLr + ν

Lr(aLr + 2ν)



1 −Lr

−1 (νLr)/(aLr + ν)



.(12)

We then expand the matrix A to two separate parts that are distinguished by

its corresponding wave speeds (eigenvalues):

A = RΛL = λ1B

− + λ2B+, (13)

where

Λ = 

λ1 0

0 λ2



, (14)

B

− =

1

aLr + 2ν



ν νLr

(aLr + ν)/Lr aLr + ν



, (15)

B

+ =

1

aLr + 2ν



aLr + ν −νLr

−(aLr + ν)/Lr ν



. (16)

The matrices B+ and B− project the solution, respectively, onto the leftand right-running waves. Hence, this is upwinding. We also note that the

distribution matrices have the following property

B

− + B+ = I, (17)

which is required for conservation. We can now distribute the cell residual

ΦC to the nodes using the projection matrices B+ and B− as described in

Fig. 2. The distribution is done this way because the left running wave, λ2,

for example, is contributing to the B+ and thus the cell residual on the left cell

is weighted with the B+. After the distribution step, we obtain the following

semi-discrete scheme in pseudo-time:

dUi

dτ =

1

hi

(B

+ΦL + B−ΦR

), (18)

where ΦL and ΦR denote the cell-residuals over the left and right cells of the

node i, respectively, and hiis the dual volume (see Fig. 1) defined by

hi =

hL + hR

2

, hL = xi − xi−1, hR = xi+1 − xi. (19)

7
i-1 i i+1

B+!i"1 B!"i B+!i

i+2

B!"i+1

Figure 2: Residual distribution to the nodes.

3.2 Steady Solver

Equation (18) can be solved with either an explicit or an implicit formulation. We choose to focus on the implicit formulation for a better transition

to the next sections, where we extend this scheme to general nonlinear timedependent advection-diffusion problems. Therefore, we drop the pseudo-time

derivative from (18) and solve the resulting system of steady residual equations,

0 = Res, (20)

by a fully implicit solver. We remark that in some complex problems the

pseudo time derivative should be kept for stability purposes. In our onedimensional problems, however, we could choose an infinite pseudo time step

without encountering any stability issue.

The implicit formulation is defined by

Uk+1 = Uk + ∆Uk, (21)

where U = (u1, p1, u2, p2, . . . , uN , pN ) and k is the iteration counter. The

correction ∆Uk = Uk+1 − Ukis determined as the solution to the linear

system:

∂Res

∂U

∆Uk = −Resk, (22)

where Reskis the steady residual vector evaluated by Uk. The Jacobian

matrix is sparse, having three 2×2 blocks in each row for all interior nodes

and two blocks for boundary nodes. The i-th interior pair of row of the linear

system is given by

Ji−1∆U

k

i−1 + Ji∆U

k

i + Ji+1∆U

k

i+1 = −

1

hi

(B

+ΦL + B−ΦR

)

k

, (23)

8
where ∆U

k

i−1 = (∆u

k

i−1

, ∆p

k

i−1

), ∆U

k

i = (∆u

k

i

, ∆p

k

i

), ∆U

k

i+1 = (∆u

k

i+1, ∆p

k

i+1),

Ji−1 =

1

hi

B

+ ∂Φ

L

∂Ui−1

, (24)

Ji =

1

hi



B

+ ∂Φ

L

∂Ui

+ B

− ∂Φ

R

∂Ui



, (25)

Ji+1 =

1

hi

B

− ∂Φ

R

∂Ui+1

. (26)

and the derivatives of the cell-residuals are given by

∂Φ

R

∂Ui

=







a −ν

−

1

Tr

−

hR

2Tr





 , (27)

∂Φ

L

∂Ui

=







−a ν

1

Tr

−

hL

2Tr





 . (28)

The linear system may be solved efficiently by Thomas’ algorithm, which

is an O(N) method. In this paper, however, we employ the Gauss-Seidel (GS)

relaxation, which is also an O(N) method for the discretization arising from

hyperbolic advection-diffusion system and extends straightforwardly to more

complex systems, irregular grids, and higher dimensions, whereas Thomas’

algorithm is not applicable in higher dimensions. The GS relaxation scheme

may be applied equation by equation (decoupled relaxation) or node by node

(collective relaxation) updating the solution set at a node simultaneouously

[20]. In this work, the collective GS relaxation is employed for robustness; the

decoupled GS relaxation was found to give better convergence in most cases

but also found to be unstable in some cases whereas the collective relaxation

encountered no such problems. In actual computations, we do not fully solve

but relax the linear system to reduce the linear residual by two or three orders

of magnitude.

Note that the Jacobian is exact and thus the implicit algorithm is Newton’s method. It is one of the advantages of the residual-distribution method

that second-order accuracy is achieved within a compact stencil, thus allowing

the exact linearization with a sparse Jacobian matrix. The same is true for

extension to higher dimensions. For linear problems, the method is essentially

a direct method. That is, if we solve the linear system exactly, then we get

the solution at the first Newton iteration.

9
1 2

B!"1 B+!1 B!"2

3 N

N-1

B+!N "2 B!"N !1 B+!N "1

N-2

Figure 3: Schematic of the Dirichlet and the Neumann boundary conditions

implementation in the hyperbolic advection-diffusion scheme.

3.3 Boundary Condition

In this section we describe the boundary condition implementation for our firstorder hyperbolic system scheme. Because in our hyperbolic advection-diffusion

system the second equation in the hyperbolic system solves for the solution

gradient, p, the Neumann boundary condition is implemented exactly in the

same way as the Dirichlet boundary condition. Therefore, we only show the

formulation for the u variable and the same can be repeated for the p variable.

Note that the discrete problem has a unique solution for linear problems with

two values given at the left and right boundaries (this is consistent with the

differential equation allowing two boundary conditions): 2(N−1) cell-residuals

for 2N − 2 unknowns, which is the case in both of the following boundary

conditions.

Dirichlet or Neumann Boundary Condition

For Dirichlet or Neumann boundary condition type schematically shown in

Fig. 3, the discretized equation for the boundary condition imposed on u takes

the following form at the first node:

J

∗

1∆U

k

1 + J

∗

2∆U

k

2 = −

1

h1

(B

−ΦR

)

k

, (29)

where for the imposed boundary condition on the u variable, the Jacobian

matrices are

J

∗

1 =

1

h1

B

−







1 0

∂Φ

R(2)

∂u1

∂Φ

R(2)

∂p1





 , (30)

J

∗

2 =

1

h2

B

−







0 0

∂Φ

R(2)

∂u2

∂Φ

R(2)

∂p2





 , (31)

where ΦR(2) denotes the second component of ΦR corresponding to the second

equation (the p equation) of the first-order hyperbolic system. The same

process is repeated for imposed boundary condition on the last node (for either

u or p.)

10
N 1 2

B+!N "1 B!

"1

N-1

Figure 4: Schematic of the periodic boundary condition implementation in the

hyperbolic advection-diffusion scheme.

Periodic Boundary Condition

To implement a periodic boundary condition we set UN = U1 in Eq. (23) and

obtain the following discretized equation and applied it on both two boundaries

(Fig. 4):

JN−1∆U

k

N−1 + J1∆U

k

1 + J2∆U

k

2 = −

(B+Φ

L + B−ΦR)k

h1

, (32)

where L denotes the cell defined by the nodes N − 1 and N, and R denotes

the cell defined by the nodes 1 and 2.

4 Extension to Time-Dependent Problems

We extend the steady state formulation by modifying the first-order hyperbolic

system given by Eq. (2) and rewriting it as a dual-time step system:

∂τu = −aux + ν∂xp −

α

∆t

u + s, (33)

∂τp = (∂xu − p)/Tr, (34)

where t is the physical time, ∆t is the physical time step, τ is the pseudo-time,

and s and α both depend on the physical time-step discretization scheme.

Adopting first- or second-order Backward-Differencing-Formulation (BDF1 or

BDF2, respectively) for the physical time integration, the α and s(x) terms

are:

α = 1, s =

u

n

∆t

: BDF1,

α =

2∆t

n−1 + ∆tn

∆t

n−1 + ∆tn

, s =

∆t

n−1 + ∆tn

∆t

n−1∆tn

u

n −

∆t

n−1

∆t

n(∆tn−1 + ∆tn)

u

n−1

: BDF2,

(35)

where a nonuniform time step formulation is used in the BDF2 formulation.

The BDF1 is only used at the very beginning of the simulation with an initial

time interval of ∆t

0

(<< ∆t

n

). After the first time interval, BDF2 with a

uniform time step, ∆t, is used. This procedure will ensure a second-order

11
accurate result through all times. The use of BDF1 is necessary because

explicit time stepping is not possible with the hyperbolic system method.

We solve the above system over each physical time interval to obtain the

steady state solution in the pseudo-time space. In the steady state in τ , we

have p = ∂xu and therefore recover the consistent implicit time discretization

scheme for the original time-dependent advection-diffusion equation.

In the dual-time step system, the physical time derivative term acts as a

source term to the first equation. This system thus shares the same eigenvalues

and eigenvectors as the steady hyperbolic advection-diffusion system. We will

see this more clearly by writing the time-dependent hyperbolic system in the

vector form:

∂U

∂τ + A

∂U

∂x = DU + S, (36)

where

A =



a −ν

−1/Tr 0



, D =



−α/∆t 0

0 −1/Tr



, S =



s(x)

0



. (37)

Because of the similarity in wave definitions in both time-dependent and steady

hyperbolic systems, the distribution matrices B− and B+ are the same for both

of these systems. Once the discretization is obtained, we drop the pseudo-time

derivatives and solve the resulting system of discrete equations by Newton’s

method as described in Section 3.2. To develop a time-dependent scheme, we

therefore only need to add the source term effect in the cell-residual.

The cell-residual ΦC of the time-dependent hyperbolic advection-diffusion

system is

Φ

C =

Z xi+1

xi

(−AUx + DU + S)dx

=









−a(ui+1 − ui) + ν(pi+1 − pi) − (xi+1 − xi)

α

∆t

(ui+1 + ui)/2

1

Tr



ui+1 − ui −

pi+1 + pi

2

(xi+1 − xi)











k+1

+

"

(xi+1 − xi)(si+1 + si)/2

0

#n−1,n

, (38)

where k and n are the Newton iteration counter and the physical time index,

respectively. Note that the second term of the Eq. (38), which is computed at

the two previous physical time steps, is constant during the Newton iteration,

and thus will not contribute to the Jacobian. The residual Jacobians needed

in the Newton solver are exactly in the same form as Eq. (26), but the derivatives of the cell-residuals now include the contribution from the physical time

12
derivative:

∂ΦR

∂Ui

=







a − hR

α

2∆t

−ν

−

1

Tr

−

hR

2Tr





 , (39)

∂ΦL

∂Ui

=







−a − hL

α

2∆t

ν

1

Tr

−

hL

2Tr





 . (40)

The hyperbolic advection-diffusion scheme is now time-accurate. At each

physical time level n, we solve the pseudo steady problem by Newton’s method

with the current solution at n as the initial solution.

5 Extension to Nonlinear Advection-Diffusion

System

Consider the generalized time-dependent nonlinear advection-diffusion equation:

∂tu + ∂xf = ∂x(νux), (41)

where f is a nonlinear function of u and ν = ν(u) is a function of u. We write

the above equation using the dual-time stepping formulation in the following

first-order hyperbolic system as

∂τu = −∂xf + ∂xp −

α

∆t

u + s, (42)

Tr

ν(u)

∂τp = ∂xu − p/ν(u). (43)

This is a preconditioned conservative form introduced in Ref. [3] to extend the

hyperbolic method to nonlinear equations. Note that the variable p will be

the diffusive flux in the pseudo steady state, not the gradient. In the vector

form, the preconditioned conservative system is written as

P

−1Uτ = −Fx + S, (44)

where

P

−1 =



1 0

0 Tr/ν(u)



, F =



f − p

−u



, S =

 −α

∆t

u + s(x)

−p/ν(u)



, (45)

and P is the preconditioning matrix. The main role of the preconditioning

matrix here is not to optimize the condition number of the differential system,

but rather to simplify the discretization. In particular, in this form, there is

13
no need to differentiate the diffusion coefficient with respect to the solution to

derive the flux Jacobian. Multiplying both sides by P from the left, the flux

term becomes

PFx = P

∂F

∂U

Ux = AUx (46)

where

A =



a −1

−ν/Tr 0



. (47)

Here, we have introduced a = ∂f/∂u, which is not a constant but a function of u. Note that the preconditioned Jacobian is formally equivalent to

the Jacobian of the linear time dependent hyperbolic advection-diffusion system. Hence, the eigenvalues and eigenvectors of this system are identical to

those of the linear equations. It greatly simplifies the construction of the distribution matrices as described below. This is one of the advantages of the

preconditioned conservative formulation.

To define the cell-residual, we first integrate the right hand side of Eq. (44)

to get the unpreconditioned residual ΨC:

Ψ

C =

Z xi+1

xi

(−Fx + S)dx

= −(Fi+1 − Fi) + Si+1 + Si

2

(xi+1 − xi). (48)

For the advective flux that is quadratic in the solution, it is well known that

a conservative linearization exists:

Ψ

C = −

∂F

∂U

(Ui+1 − Ui) + Si+1 + Si

2

(xi+1 − xi), (49)

where ∂F

∂U

is the analytical Jacobian evaluated at the arithmetic average of the

solution, (Ui+1 + Ui)/2 and it satisfies

Fi+1 − Fi =

∂F

∂U

(Ui+1 − Ui), (50)

exactly. The preconditioned cell-residual, which is to be distributed, is then

given by

Φ

C = PΨ

C

, (51)

where the matrix P is evaluated at the arithmetic average of the solution so

that we have

Φ

C = −A(Ui+1 − Ui) + Si+1 + Si

2

(xi+1 − xi). (52)

14
The distribution matrices can then be constructed in the same way as in the

linear case based on A:

B

−

C =

1

aLr + 2ν

"

ν νLr

(aLr + ν)/Lr aLr + ν

#

, (53)

B

+

C =

1

aLr + 2ν

"

aLr + ν −νLr

−(aLr + ν)/Lr ν

#

. (54)

where

a ̄ =

ai + ai+1

2

, ν ̄ =

νi + νi+1

2

. (55)

The subscript C indicates that these matrices are defined over the cell C.

Note that the distribution matrices are not globally constant in general for

nonlinear equations. The construction of the distribution matrices based on

the conservative linearization is important for proper upwinding. If different

averages are used, the distribution may not be strictly upwind. If the target

equation does not allow the conservative linearization, the construction proposed in Ref. [21] may be employed to design a proper upwind distribution

matrix. Note also that we evaluate ν in the source term by the average value

over the cell to simplify the second component of the cell-residual as follows:

Φ

C

i =









−(fi+1 − fi) + (pi+1 − pi) − (xi+1 − xi)

α

∆t

(ui+1 + ui)/2

1

Tr

[ν(ui+1 − ui) − p(xi+1 − xi)]









k+1

+

"

(xi+1 − xi)(si+1 + si)/2

0

#n−1,n

, (56)

where

Tr =

Lr

a ̄ + ν/Lr

. (57)

Observe that we have p = ν(ui+1 − ui)/(xi+1 − xi) in the pseudo steady state.

The cell-residual is then distributed to the nodes in all cells to yield the residual

equation at node i:

0 =

1

hi

(B

+

L Φ

L + B

−

R Φ

R

), (58)

where the pseudo-time derivatives have been dropped.

15
We note that the Jacobian of this nonlinear advection-diffusion system

is slightly different than that of the linear problem. The Jacobian blocks,

however, are function of derivative of the distribution matrices B− and B+:

Ji−1 =

1

hi

∂



B

+

L Φ

L



∂Ui−1

, (59)

Ji =

1

hi

 

∂



B

+

L Φ

L



∂Ui

+

∂



B

−

R Φ

R



∂Ui

!

, (60)

Ji+1 =

1

hi

∂



B

−

R Φ

R



∂Ui+1

. (61)

In the special case a/ν = constant, which is the only case considered in this paper, the distribution matrices reduce entirely to constant matrices. Therefore,

we can entirely ignore the derivatives of B

−

i

and B

+

i

in deriving the Jacobian. As a result, the Jacobians are given formally as in Eq. (26). This makes

the implementation of both linear and nonlinear systems extremely easy and

straightforward. Effects of nonlinearity arise only in the derivatives of the left

and right cell-residuals:

∂Φ

R

∂Ui

=









ai − hR

α

2∆t

−1

ui+1 − ui

Tr

∂ν

∂ui

+

ν(ui+1 − ui) − p hR

Lr



∂a

∂ui

+

1

Lr

∂ν

∂ui



−

ν ̄

Tr

−

hR

2Tr









k

,

(62)

∂Φ

L

∂Ui

=









−ai − hL

α

2∆t

1

ui − ui−1

Tr

∂ν

∂ui

+

ν(ui − ui−1) − p hL

Lr



∂a

∂ui

+

1

Lr

∂ν

∂ui



+

ν ̄

T r

−

hL

2Tr









k

,

(63)

where the average values indicated by the over bar are understood to be taken

over the corresponding cell. Note that these derivatives depend on the solution

at the current iteration, and therefore the superscript k has been introduced.

The same applies to the distribution matrices also. Obviously, the nonlinear

system approaches the linear system in the event that both the advection and

the diffusion coefficients are constant. Also, the boundary conditions can be

implemented as described in Section 3.3 for the linear case. Therefore, our

hyperbolic nonlinear time dependent advection-diffusion system scheme can

be implemented for linear, nonlinear, steady and unsteady problems without

loss of accuracy. For simplify the discussion we only showed Jacobian for nonlinear cases with a/ν = constant but the scheme works for general nonlinear

problems.

16
6 Results

In this section we present the results in three categories: 1) Steady advectiondiffusion equation for high Reynolds (or Peclet) number applications, 2) Unsteady linear advection-diffusion, and 3) Unsteady nonlinear advection-diffusion

problems. We then present the order of accuracy results in the last subsection. Note that time-accurate computations are started by BDF1 in the first

time step with extremely small time step (e.g. t = 10−8), and then by BDF2

thereafter with much larger time step (e.g. t = 0.001 − 0.5) for all unsteady

cases. This will ensure that the order of accuracy stays second order through

all times. We remark that explicit time stepping is not available for timeaccurate computations with the hyperbolic system method.

6.1 Steady Linear Advection-Diffusion (High Re Applications)

Consider the advection-diffusion equation in x ∈ (0, 1) with u(0) = 0 and

u(1) = 1 boundary conditions:

∂tu + a∂xu = ν∂xxu + s(x), (64)

where

s(x) = π

Re

[a cos(πx) + πν sin(πx)], Re = a/ν. (65)

This is a boundary layer problem with a non-trivial steady state solution in

the diffusion limit as a result of the source term addition [2] . This equation

develops a very narrow boundary layer near the right boundary (x = 1) where

the advection becomes dominant. The exact steady state solution to this

problem is given by

u

exact(x) = e

−Re − e(x−1)Re

e

−Re − 1

+

1

Re

sin(πx). (66)

We chose Re values ranging from 1 to 106 and solved the equation on nonuniform grid sizes up to 30000 nodes. We remark that the high-Re cases required extremely fine grids to meet the well-known requirement on the mesh

Reynolds-number [2]. If desired, the computations can be performed on substantially coarser grids with more aggressive and customized grid stretching

as well as with some nonlinear mechanism to prevent numerical oscillations.

However, we simply refined the grid to meet the mesh Reynolds-number requirement because our method is powerful enough to solve the problem very

efficiently (i.e., less than 10 Newton iterations) even on such dense grids. The

ability to efficiently solve the problem on highly refined grids is a great advantage of these schemes. Shown in Fig. 5 are the comparisons between the

17
0 0.2 0.4 0.6 0.8 1 −0.2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

x

u(x)

Re=1.0

Re=10.0

Re=103

ï105

Exact

Hyp−ADE

0.8 0.85 0.9 0.95 1

10−9

10−8

10−7

10−6

10−5

10−4

10−3

10−2

x

log(u(x))

Re=104

Re=105

Re=103

Exact

Hyp−ADE

Figure 5: Residual-Distribution scheme on hyperbolic advection-diffusion system for boundary layer equation.

computational and the exact solutions. In this comparison, the exact analytic

solution is represented with + (plus) symbols while the numerical results are

shown with solid lines (denoted Hyp-ADE). The exact solutions are plotted at

nodes, and therefore they also show the computational grids. The solutions

were obtained with the Newton-GS method as described in Section 3.2. For

this problems, within each Newton iteration, the GS relaxation were conducted

until three orders of magnitude reduction is achieved in the linear solver. The

computations were continued until the residuals of both the solution and the

solution gradient were reduced by eight orders of magnitude. Table 1 gives

the convergence data obtained using the implicit formulation of the hyperbolic advection-diffusion scheme. The linear dependency of the iterations on

Table 1: Boundary layer problem (Convergence criteria: Residuals < 10−8.)

log10 Re Number of nodes GS relaxations/Newton iteration Newton iteration

0 10 17 7

0 100 324 7

0 200 631 7

0 300 967 7

1 300 549 7

2 300 131 7

3 300 28 7

4 300 38 6

5 3000 60 6

6 30000 50 5

the grid size is demonstrated at Re = 1, which is a consequence of solving the

advection-diffusion equation as a hyperbolic system for all Reynolds numbers

18
through the diffusion limit. This is remarkable because the linear convergence

is retained for any irregular grid in any dimensions (N is approximately the

number of nodes in each coordinate direction in two and three dimensions)

as demonstrated in Refs. [1–5]. In conventional methods, e.g., such as the

alternating direction implicit method with Thomas’ algorithm, the linear convergence can be achieved only on Cartesian grids (one-dimensional grid must

exist in each coordinate direction). The O(N) convergence on general problems and arbitrary grid systems makes the current scheme extremely attractive

because it leads to orders of magnitude faster convergence in comparison with

conventional methods whose convergence is typically O(N2).

Table 1 shows for the grid of 300 nodes that the number of GS relaxations

decreases significantly as Re increases. In Ref. [2], a formula for Lr was derived

that minimizes the condition number of the hyperbolic advection-diffusion

system and a nearly uniform convergence was obtained for an explicit pseudotime marching scheme. We tested the formula in our method, but found that it

did not change the significant decrease although it did reduce the number of GS

relaxations approximately by 10 %. A different approach may be necessary to

derive a formula to achieve Re-independent convergence of the GS relaxation

in the implicit solver. Further study on Re-independent convergence is left

as future work. Nevertheless, it should be noted that the number of Newton

iterations is nearly Re-independent and very small: only 5 to 7 iterations to

reduce the residual by eight orders of magnitude. Considering the fact that the

cost of one GS relaxation is significantly cheaper than one Newton iteration,

we find that the developed solver is still a very powerful steady solver. It thus

serves as a very efficient solver for the unsteady residual equations as shown

in the following sections.

Finally, we remark that the high-Re cases required extremely fine grids to

meet the well-known requirement on the mesh Reynolds-number as described

in Ref. [2]. If desired, the computations can be performed on substantially

coarser grids with more aggressive grid stretching. However, we simply refined

the grid to meet the mesh Reynolds-number requirement because our method

is powerful enough to solve the problem very efficiently (i.e., 5 to 7 Newton

iterations) even on such dense grids. The ability to efficiently solve the problem

on highly refined grids is a great advantage.

6.2 Unsteady Linear Advection-Diffusion

We present unsteady computations for the two types of boundary conditions;

Dirichlet (or Neumann) and periodic. We again note that for hyperbolic system scheme the implementation of Dirichlet and Neumann boundary conditions are similar.

19
Consider the time-dependent advection-diffusion equation in x ∈ (0, 1)

∂tu + a∂xu = ν∂xxu, (67)

with the following time-dependent boundary condition (Dirichlet):

u(0, t) = 0, (68)

u(1, t) = U cos(ωt), (69)

where U is the amplitude of the oscillation and ω is the frequency of the oscillation on the right boundary. This problem has the following exact analytic

solution:

u

exact(x, t) = Real 

e

λ1x − eλ2x

e

λ1 − eλ2

Ueiωt, λ1,2 =

a ±

√

a

2 + 4iων

2ν

, (70)

where i =

√

−1. We solved this time-dependent advection-diffusion equation

with the following first-order hyperbolic system:

∂τu + a∂xu = ν∂xp −

3

2∆t

u +

4u

n − un−1

2∆t

, (71)

∂τp = (∂xu − p)/Tr. (72)

The solution evolution was started with the exact initial solution given as

u(x, 0) = Real 

e

λ1x − eλ2x

e

λ1 − e

λ2

U



, (73)

and was advanced in physical time after the time-dependent residual was reduced by two orders of magnitude at each physical time step. We tested the

scheme on a number of nonuniform grid systems ranging from N = 100...10000

nodes. Shown in Fig. 6 are the solution and the solution gradient obtained

on the coarsest mesh (N=100) for U = 1, ν = a = 1, and ω = 7π/2. We

show the results on the coarsest mesh just to demonstrate the accuracy of the

scheme. Similar results were obtained on finer grid systems. The results are

then compared with the analytic solution. Here the numerical results (denoted

Hyp-ADE) are over plotted with the exact solutions, meaning excellent agreement. We will show the residual values in the accuracy verification section.

We also obtained convergence data for different grid systems over 100 time

steps. These are tabulated in Table 2, which shows a linear increase in number

of iterations with increase of the grid system. This is, again, a consequence of

solving the advection-diffusion equation as a hyperbolic system, where there

are only first-derivatives and therefore the number of iterations is O(N), not

O(N2). It is also shown in Table 2 that only four Newton iterations are

required for all grids considered.

20
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ï2

ï1.5

ï1

ï0.5

0

0.5

1

1.5

2

x

u(x,t)

Exact

HypïADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ï2

ï1.5

ï1

ï0.5

0

0.5

1

1.5

2

x

u(x,t)

Exact

HypïADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ï2

ï1.5

ï1

ï0.5

0

0.5

1

1.5

2

x

p(x,t)

Exact

HypïADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ï2

ï1.5

ï1

ï0.5

0

0.5

1

1.5

2

x

p(x,t)

Exact

HypïADE

Figure 6: Time-dependent advection-diffusion with oscillatory boundary condition u(0, t) = 0, u(1, t) = cos( 7π

2

t) on N = 100 irregular nodes. Time step =

0.01s. Left figures at t = 0.1s, right figures at t = 0.4s.

Table 2: Hyperbolic time-dependent linear advection-diffusion problem with

BDF2. Average data over 100 time steps are given. (Convergence criteria:

Residuals < 10−2)

Number of nodes GS relaxations/Newton iteration Newton iterations

100 127 4

300 370 4

500 607 4

1000 1201 4

21
For the periodic boundary condition demonstration, consider the timedependent advection-diffusion equation in x ∈ (0, 1) given as in Eq. (67) with

the following initial condition:

u(x, t = 0) = sin(κx), (74)

where κ is an arbitrary constant. The exact solution to this problem with a

periodic boundary condition is:

u

exact(x, t) = e−κ

2νt sin(κ(x − at)). (75)

We solved this problem with the same first-order hyperbolic advection-diffusions

equation given as in Eqs. (71) and (72). For each physical time, we reduced

the residuals by two orders of magnitude before advancing in time. During

each time step, we also relaxed the linear system using GS relaxations until

two orders of magnitude reduction in the linear system residuals was achieved.

Shown in Fig. 7 are the solution and solution gradient on N = 25 uniform grid

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ï2

ï1.5

ï1

ï0.5

0

0.5

1

1.5

2

x

u(x,t)

Exact

HypïADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ï2

ï1.5

ï1

ï0.5

0

0.5

1

1.5

2

x

u(x,t)

Exact

HypïADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ï2

ï1.5

ï1

ï0.5

0

0.5

1

1.5

2

x

u(x,t)

Exact

HypïADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ï10

ï8

ï6

ï4

ï2

0

2

4

6

8

10

x

p(x,t)

Exact

HypïADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ï10

ï8

ï6

ï4

ï2

0

2

4

6

8

10

x

p(x,t)

Exact

HypïADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ï10

ï8

ï6

ï4

ï2

0

2

4

6

8

10

x

p(x,t)

Exact

HypïADE

Figure 7: Time-dependent linear advection-diffusion problem(a = 1, ν = 0.03)

with periodic boundary condition on N = 25 uniform nodes. Time step =

0.01s. From left to right, solutions at t = 0.1s, 0.5s, 1.0s.

for a = 1, ν = 0.03, and κ = 2π. The results (Hyp-ADE) are over plotted

with the exact solution, indicating excellent accuracy on such a coarse grid.

Similar results were also obtained for larger grid systems.

We examined the convergence rate of this problem on several grid systems.

Given in Table 3, are the average numbers of GS relaxations per Newton iteration obtained over 100 time steps. Clearly the convergence rate is of O(N), not

O(N2) as typical for numerical methods for the advection-diffusion equation.

Observe that for most grid systems only one or two Newton iterations were

22
Table 3: Hyperbolic time-dependent linear advection-diffusion problem with

periodic BC (a = 1, ν = 0.03) Average data over 100 time steps are

given.(Convergence criteria: Residuals < 10−2)

Number of nodes GS relaxations/Newton iteration Newton iterations

25 4 3

100 11 3

300 33 2

500 55 2

1000 116 2

sufficient to obtain accurate solutions. We also note that we used ∆t = 0.01–

sec for all the grid systems; the BDF schemes are unconditionally stable. The

time step is orders-of-magnitude larger than that required for conventional

explicit schemes, which is limited by O(h

2

). Of course, conventional implicit

schemes also allow unconditionally large time steps, but it requires O(N2)

convergence in an iterative linear solver and potentially a much larger number of outer iterations as well if the exact linearization is not possible and

Newton’s method cannot be constructed. Hence, the method developed here

has two major advantages over conventional methods: the exact linearization

(Newton’s method) and O(N) iterative convergence in the linear solver. The

latter advantage can be potentially huge with increase of the grid system as

the speed-up factor is O(N) and thus grows for finer grids. For example, for

N = 1000, our scheme is at least four orders of magnitude faster than the

conventional schemes, which is remarkable.

6.3 Unsteady Nonlinear Advection-Diffusion

Consider the time-dependent nonlinear viscous Burgers equation with a timedependent source term:

∂tu + ∂xf = ∂x (ν ux) + g(x, t), x ∈ (0, 1), (76)

where f = u

2/2, ν = u, and

g(x, t) = u

e

t +

1

2

((u

e

)

2

)x − (u

e

x

)

2 − ueue

xx. (77)

The source term has been generated by the following function:

u

e

(x, t) = Real 

sinh(x

p

iω/ν)

sinh(piω/ν)

Ueiωt!+ C, C > 1, (78)

23
so that it is the exact solution to Eq. (76). The boundary conditions are

u(0, t) = C, (79)

u(1, t) = C + U cos(ωt), (80)

where ω is the frequency of the oscillation on the right boundary, and U is the

amplitude of the oscillation. The constant C must be greater than 1 to keep

the diffusion coefficient positive.

We solved this time-dependent nonlinear advection-diffusion equation with

the following equivalent first-order hyperbolic system:

∂τu + ∂x



u

2



= ∂xp −

3

2∆t

u +

4u

n − un−1

2∆t

+ g(x, t), (81)

Tr

ν

∂τp = (∂xu − p/ν). (82)

Shown in Fig. 8 are the results obtained on N = 100 nonuniform nodes for

U = 1, ω = 7π/2, and C = 2. The scheme produces very accurate time

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0

0.5

1

1.5

2

2.5

3

3.5

4

x

u(x,t)

Exact

HypïADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0

0.5

1

1.5

2

2.5

3

3.5

4

x

u(x,t)

Exact

Hyp−ADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0

0.5

1

1.5

2

2.5

3

3.5

4

x

u(x,t)

Exact

Hyp−ADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ï5

ï4

ï3

ï2

ï1

0

1

2

3

4

5

x

p(x,t)

Exact

HypïADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −5

−4

−3

−2

−1

0

1

2

3

4

5

x

p(x,t)

Exact

Hyp−ADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −5

−4

−3

−2

−1

0

1

2

3

4

5

x

p(x,t)

Exact

Hyp−ADE

Figure 8: Time-dependent nonlinear advection-diffusion (a = ν = u) with

oscillatory boundary condition u(0, t) = 2, u(1, t) = 2 + cos( 7π

2

t) on N = 100

nonuniform nodes. Time step = 0.01s. From left to right, solutions at t =

0.1s, 2.5s, 10.0s.

evolutions of the solution (u) and the gradient (p/ν) on nonuniform grids (see

also Fig. 9.)

The O(N) convergence rate was once again achieved for the time-dependent

nonlinear hyperbolic advection-diffusion system. This is given in Table 4,

where the average number of iterations were obtained over 1000 time steps

24
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

t

u(1,t)

Exact

HypïADE

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 ï5

ï4

ï3

ï2

ï1

0

1

2

3

4

5

t

p(1,t)

Exact

HypïADE

Figure 9: Time-dependent nonlinear advection-diffusion (a = ν = u) with

oscillatory boundary condition u(0, t) = 2, u(1, t) = 2 + cos( 7π

2

t) on N = 100

nonuniform nodes. Time step = 0.01s

(over 17 periods). Note also that the method converged at four Newton iterations for all grids considered. Again, a constant time step (∆t = 0.01–sec)

was used for all the grid systems. And the method is substantially more efficient here also than conventional explicit/implicit schemes as discussed for

the linear case in Section 6.2.

Table 4: Hyperbolic time-dependent nonlinear advection-diffusion problem

with BDF2. Average data over 100 time steps are given.(Convergence criteria:

Residuals < 10−2)

Number of nodes GS relaxations/Newton iteration Newton iterations

100 88 4

300 263 4

500 440 4

1000 882 4

6.4 Accuracy Verification

To demonstrate the overall order of accuracy of our time-dependent hyperbolic

advection-diffusion scheme, we computed the L1 =

PN

i=1(U

exact

i −Ui)/N as well

as L∞ = M ax(U

exact

i − Ui) norms of errors on a series of nonuniform grids of

N nodes, and time step increments ∆t. For each case, we terminate the solver

when the residuals were dropped by two orders of magnitude. We remark

that in some practical applications two orders of magnitude reduction in the

residuals may not be sufficient to achieve the appropriate order of accuracy.

But for the cases considered here, terminating at even nine orders of magnitude

25
resulted in similar order of accuracy. The spatial accuracy was obtained by

comparing several nonuniform grid systems on a fixed ∆t, while the temporal

accuracy was obtained by comparing various ∆t on a fixed nonuniform grid

system. We also included boundary nodes in both L1 and L∞ norms. For

unsteady cases, the L1 and L∞ norms were evaluated at t = 1.0–sec.

−3.5 −3 −2.5 −2 −1.5 −1 −0.5

−6

−5

−4

−3

−2

−1

0

Slope 2

∆t = 0.001

log10 h

log10 L

∞

u

p(=u

x

)

−3.5 −3 −2.5 −2 −1.5 −1 −0.5 0

−5

−4

−3

−2

−1

0

1

2

Slope 2

hMax. = 0.045

log10 dt

log10 L

∞

u

p(=u

x

)

Figure 10: L∞ error obtained for time-dependent linear hyperbolic advectiondiffusion system scheme. Left: Spatial accuracy, Right: Temporal accuracy.

Figures 10 and 11 show the L∞ error convergence results for both linear

and nonlinear time-dependent hyperbolic advection diffusion system, where h

is the representative mesh spacing defined by h = /(N − 1). For discussion

purposes, we only present the accuracy plots for Re = 1 results; similar results

were obtained for other Reynolds numbers. These results show that our scheme

is second-order accurate for all the variables and the gradients (note that again

that our flux integral is exact) at all the grid nodes including the boundary

nodes (see also Tables 5, 6, 7, and 8). We note that it is natural for the error

convergence in the temporal space to deteriorate when the O(∆t

2

) << O(h

2

)

as well as when ∆t is too large for the time discretization to be in asymptotic

range; i.e. a lower order of accuracy for coarse discretization in time is typical.

These limits are also included in the temporal accuracy plot for completeness.

Table 5: Spatial accuracy for the linear advection-diffusion problem with oscillatory boundary condition on several nonuniform grid systems (∆t = 0.001s).

Number of nodes L1 error of u L∞ error of u Order L1 error of p L∞ error of p Order

25 9.65E-04 3.44E-03 6.58E-03 1.71E-02

50 2.28E-05 8.98E-04 1.94 1.08E-03 4.22E-03 2.02

100 6.10E-05 2.13E-04 2.08 3.08E-04 9.95E-04 2.08

200 1.59E-06 5.04E-05 2.08 8.76E-05 2.35E-04 2.08

300 7.58E-06 2.06E-05 2.21 4.87E-05 9.53E-05 2.23

26
−3 −2.5 −2 −1.5 −1 −0.5 0

−6

−5

−4

−3

−2

−1

0

Slope 2

∆t = 0.001

log10 h

log10 L

∞

u

p(=u

x

)

−4 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0

−6

−5

−4

−3

−2

−1

0

1

2

Slope 2

hMax. = 0.045

log10 dt

log10 L

∞

u

p(=u

x

)

Figure 11: L∞ error obtained for time-dependent nonlinear hyperbolic

advection-diffusion system scheme. Left: Spatial accuracy, Right: Temporal

accuracy.

Table 6: Temporal accuracy for the linear advection-diffusion problem with

oscillatory boundary condition on a nonuniform grid (N=100) with hMax. =

0.045.

∆t L1 error of u L∞ error of u Order L1 error of p L∞ error of p Order

0.500 3.01E-01 6.63E-01 2.46E+00 4.09E+00

0.100 2.32E-02 7.55E-02 1.35 1.25E-01 2.32E-01 1.78

0.050 4.37E-03 1.72E-02 2.13 3.57E-02 6.28E-02 1.88

0.025 9.66E-04 3.16E-03 2.44 1.16E-02 2.26E-02 1.47

0.010 1.48E-04 2.99E-04 2.57 2.21E-03 4.63E-03 1.73

0.001 6.02E-05 2.13E-04 0.15 3.05E-04 9.95E-04 0.67

Table 7: Spatial accuracy for the nonlinear advection-diffusion problem with

oscillatory boundary condition on several nonuniform grid systems (∆t =

0.001s).

Number of nodes L1 error of u L∞ error of u Order L1 error of p L∞ error of p Order

10 1.83E-03 8.41E-03 1.20E-02 4.46E-02

25 4.36E-04 1.50E-03 1.88 2.79E-03 7.29E-03 1.98

50 9.00E-05 3.74E-04 2.00 5.48E-04 1.78E-03 2.03

100 2.80E-05 9.39E-05 1.99 1.41E-04 4.36E-04 2.03

200 8.19E-06 2.33E-05 2.01 4.81E-05 1.07E-04 2.03

27
Table 8: Temporal accuracy for the nonlinear advection-diffusion problem

with oscillatory boundary condition on a nonuniform grid system (N=100)

with hMax. = 0.045.

∆t L1 error of u L∞ error of u Order L1 error of p L∞ error of p Order

0.200 6.46E-02 1.55E-01 4.97E-01 6.83E-01

0.100 1.52E-02 4.67E-02 1.73 9.00E-02 1.54E-01 2.15

0.050 2.03E-03 7.40E-03 2.66 1.99E-02 3.45E-02 2.16

0.020 3.52E-04 7.32E-04 2.52 5.01E-03 1.03E-02 1.32

0.010 1.34E-04 3.20E-04 1.19 1.70E-03 3.06E-03 1.75

0.005 6.68E-05 1.41E-04 1.18 6.20E-04 9.90E-04 1.63

0.001 2.80E-05 9.39E-05 0.25 1.41E-04 4.36E-04 0.51

7 Conclusions

In this paper, we have extended the first-order hyperbolic system method to

time-accurate computations. The time integration scheme is implicit with the

second-order backward-difference formula, and the resulting system of timedependent residual equations is solved by a steady solver developed based on

the hyperbolic method. The steady solver is Newton’s method constructed

from a compact upwind residual-distribution scheme for linear and nonlinear

advection-diffusion equations. The construction of the numerical scheme was

greatly simplified for the nonlinear equation by the preconditioned conservative

formulation. We have demonstrated that the developed time-accurate schemes

are capable, at every physical time step, of producing second-order accurate

solution and gradient on highly refined nonuniform grids by less than five

Newton iterations with the linear convergence of the block GS relaxations for

all test cases considered. Thus, the notable features of the hyperbolic method

have been shown to carry over to time-dependent problems.

The developed methodology based on the second-order backward Euler

time integration can be extended to practical one-dimensional problems [7–11],

bringing significant improvements in both accuracy and efficiency. Extension

to higher-order accuracy is also possible by improving the residual evaluation

and higher-order backward-difference formulas. It is relatively simple in one

dimension because the integration of the flux balance term is exact in one

dimension as mentioned in Section 3.1. In higher dimensions, although the

construction of time-accurate schemes is as simple as presented in the paper,

it is known that a high-order quadrature would be required for the flux balance

integration in order to achieve the equal order of spatial accuracy in the solution and the gradient [2]. Although not considered in the present paper, the

construction of non-oscillatory schemes is an important area of development

for problems with discontinuities such as shock waves.

28
Acknowledgments

The second author acknowledges the support by the U.S. Army Research Office under the contract/grant number W911NF-12-1-0154 with Dr. Frederick

Ferguson as the program manager.

References

1. H. Nishikawa. A first-order system approach for diffusion equation. I:

Second-order residual-distribution schemes. J. Comput. Phys., 227:315–

352, 2007.

2. H. Nishikawa. A first-order system approach for diffusion equation. II:

Unification of advection and diffusion. J. Comput. Phys., 229:3989–4016,

2010.

3. H. Nishikawa. New-generation hyperbolic Navier-Stokes schemes: O(1/h)

speed-up and accurate viscous/heat fluxes. In Proc. of 20th AIAA Computational Fluid Dynamics Conference, AIAA Paper 2011-3043, Honolulu,

Hawaii, 2011.

4. H. Nishikawa. First, second, and third order finite-volume schemes for

diffusion. In Proc. of 51st AIAA Aerospace Sciences Meeting, AIAA Paper

2011-1125, Grapevine, Texas, January 2013.

5. H. Nishikawa. First, second, and third order finite-volume schemes for

advection-diffusion. In Proc. of 21st AIAA Computational Fluid Dynamics

Conference, AIAA Paper 2011-2568, San Diego, California, June 2013.

6. A. Jameson. Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings. AIAA Paper 91-1596,

1991.

7. Y.-K. Chen and F. S. Milos. Ablation and thermal response program

for spacecraft heatshield analysis. J. Spacecraft and Rockets, 36:475–583,

1999.

8. M. Mahzari, R. D. Braun, T. White, and D. Bose. Preliminary analysis

of the mars science laboratory entry aerothermodynamic environment and

thermal protection system performance. In Proc. of 51st AIAA Aerospace

Sciences Meeting including the New Horizons Forum and Aerospace, AIAA

Paper 2013-0185, Grapevine, Texas, January 2013.

9. D. Bose, T. White, M. Mahzari, and K. Edquist. A reconstruction of

aerothermal environment and thermal protection system response of the

29
mars science laboratory entry vehicle. In Proc. of American Astronomical Society 221st Meeting, AAS 2013-311-0712, Long Beach, California,

January 2013.

10. N. R. Merski. Reduction and analysis of phosphor thermography data with

the IHEAT software package. In Proc. of 36th AIAA Aerospace Sciences

Meeting and Exhibits, AIAA Paper 98-0712, Reno, Nevada, January 1998.

11. T. J. Horvath and s. A. Berry and b. R. Hollis and d. S. Liechty and h. H.

Hamilton and n. R. Merski. J. Spacecraft and Rockets, 38:635–645, 2001.

12. A. Ferrante and H. Deconinck. Solution of the unsteady euler equations using residual distribution and flux corrected transport. In Technical Report

VKI-PR 97-08. 1997.

13. A. Csik, M. Ricchiuto, H. Deconninck, and S. Poedts. Space-time residual distribution schemes for hyperbolic conservation laws. In Proc. of

15th AIAA Computational Fluid Dynamics Conference, AIAA Paper 2001-

2617, Anaheim, 2001.

14. R. Abgrall and M. Mezine. Construction of second order accurate monotone and stable residual distribution schemes for unsteady flow problems.

J. Comput. Phys., 188:16–55, 2002.

15. A. Bonfiglioli and R. Paciorri. A mass-matrix formulation of unsteady

fluctuation splitting schemes consistent with Roe’s parameter vector. Int.

J. Comput. Fluid Dyn., 27, 2013.

16. G. Rossiello, P. De Palma, G. Pascazio, and M. Napolitano. Third-orderaccurate fluctuation-splitting schemes for unsteady hyperbolic problems.

J. Comput. Phys., 222:332–352, 2007.

17. D. Caraeni and L. Fuchs. Compact third-order multidimensional upwind

scheme for Navier-Stokes simulations. Appl. Math., 15:373–401, 2002.

18. M. Ricchiuto and R. Abgrall. Explicit Runge-Kutta residual distribution

schemes for time dependent problems: Second order case. J. Comput.

Phys., 229:5653–5691, 2010.

19. M. Ricchiuto and H. Nishikawa. Dual-time hyperbolic diffusion systems.

Unpublished, 2012.

20. U. Trottenberg, C. W. Oosterlee, and A. Sch ̈uller. Multigrid. Academic

Press, 2000.

21. A. Csik, M. Ricchiuto, and H. Deconinck. A conservative formulation

of the multidimensional upwind residual distribution schemes for general

nonlinear conservation laws. J. Comput. Phys., 179:286–312, 2002.

30
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Memorandum

 4. TITLE AND SUBTITLE

First-Order Hyperbolic System Method for Time-Dependent 

Advection-Diffusion Problems 

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Mazaheri, Alireza; Nishikawa, Hiroaki

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-20378

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 64

Availability: NASA CASI (443) 757-5802

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

A time-dependent extension of the first-order hyperbolic system method for advection-diffusion problems is introduced. 

Diffusive/viscous terms are written and discretized as a hyperbolic system, which recovers the original equation in the steady 

state. The resulting scheme offers advantages over traditional schemes: a dramatic simplification in the discretization, 

high-order accuracy in the solution gradients, and orders-of-magnitude convergence acceleration. The hyperbolic 

advection-diffusion system is discretized by the second-order upwind residual-distribution scheme in a unified manner, and the 

system of implicit-residual-equations is solved by Newton's method over every physical time step. The numerical results are 

presented for linear and nonlinear advection-diffusion problems, demonstrating solutions and gradients produced to the same 

order of accuracy, with rapid convergence over each physical time step, typically less than five Newton iterations.

15. SUBJECT TERMS

Finite Volume; Higher-Order Scheme; Nonlinear; Residual Distribution; Unsteady

18. NUMBER

 OF 

 PAGES

35

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

747797.02.09.30.04 

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/TM-2014-218175

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 03 - 2014