FC-based shock-dynamics solver with neural-network localized artificial-viscosity assignment

Journal of Computational Physics: X 15 (2022) 100110

Contents lists available at

ScienceDirect

Journal

of

Computational

Physics:

X

www.elsevier.com/locate/jcpx

FC-based

shock-dynamics

solver

with

neural-network

localized

artificial-viscosity

assignment

Oscar

P. Bruno

a

,

∗

,

Jan

S. Hesthaven

b

,

Daniel

V. Leibovici

a

Computing

and

Mathematical

Sciences,

Caltech,

Pasadena,

CA

91125,

USA

b

Computational

Mathematics

and

Simulation

Science

(MCSS),

Ecole

Polytechnique

Federale

de

Lausanne,

CH-1015

Lausanne,

Switzerland

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Article

history:

Received

29

October

2021

Received

in

revised

form

1

June

2022

Accepted

6

June

2022

Available

online

9

June

2022

Keywords:

Machine

learning

Neural

networks

Shock

dynamics

Artificial

viscosity

Fourier

continuation

Non-periodic

domain

This

paper

presents

a

spectral

scheme

for

the

numerical

solution

of

nonlinear

conservation

laws

in

non-periodic

domains

under

arbitrary

boundary

conditions

.

The

approach

relies

on

the

use

of

the

Fourier

Continuation

(FC)

method

for

spectral

representation

of

non-periodic

functions

in

conjunction

with

smooth

localized

artificial

viscosity

assignments

produced

by

means

of

a

Shock-Detecting

Neural

Network

(SDNN).

Like

capturing

schemes

and

artificial

viscosity

techniques,

the

combined

FC-SDNN

strategy

effectively

controls

spurious

oscillations

in

the

proximity

of

discontinuities.

Thanks

to

its

use

of

a

localized

but

smooth

artificial

viscosity

term

,

whose

support

is

restricted

to

a

vicinity

of

flow-discontinuity

points,

the

algorithm

enjoys

spectral

accuracy

and

low

dissipation

away

from

flow

discontinuities,

and,

in

such

regions,

it

produces

smooth

numerical

solutions—

as

evidenced

by

an

essential

absence

of

spurious

oscillations

in

level

set

lines.

The

FC-

SDNN

viscosity

assignment,

which

does

not

require

use

of

problem-dependent

algorithmic

parameters,

induces

a

significantly

lower

overall

dissipation

than

other

methods,

including

the

Fourier-spectral

versions

of

the

viscosity

method.

The

character

of

the

proposed

algorithm

is

illustrated

with

a

variety

of

numerical

results

for

the

linear

advection,

Burgers

and

Euler

equations

in

one

and

two-dimensional

non-periodic

spatial

domains.

©

2022

The

Authors.

Published

by

Elsevier

Inc.

This

is

an

open

access

article

under

the

CC

BY-NC-ND

license

(

http://creativecommons.org/licenses/by-nc-nd/4.0/

).

1.

Introduction

This

paper

presents

a

new

“FC-SDNN”

spectral

scheme

for

the

numerical

solution

of

nonlinear

conservation

laws

under

arbitrary

boundary

conditions.

The

proposed

approaches

rely

on use

of

the

FC-Gram

Fourier

Continuation

method [

1

,

2

,

6

,

29

]

for

spectral

representation

of

non-periodic

functions

in

conjunction

with

localized

smooth

artificial

viscosity

assignments

prescribed

by

means

of

the

neural

network-based

shock-detection

method

proposed

in [

38

].

The

neural

network

ap-

proach [

38

]itself

utilizes

Fourier

series

to

discretize

the

gas

dynamics

and

it

eliminates

Gibbs

ringing

at

shock

positions

(which

are

determined

by

means

of

an

artificial

neural

network)

by

assigning

artificial

viscosity

over

a

small

number

of

discrete

points

in

a

close

vicinity

of

shocks.

The

use

of

the

classical

Fourier

spectral

method

in

that

contribution

restricts

the

method’s

applicability

to

periodic

problems

(so

that,

in

particular,

the

outer

computational

bound-

aries

cannot

be

physical

boundaries),

and

its

highly

localized

viscosity

assignments

give

rise

to

a

degree

of

non-smoothness,

resulting

in

certain

types

of

unphysical

oscillations

manifested

as

serrated

level-set

lines

in

the

flow

fields.

The

FC-SDNN

*

Corresponding

author.

E-mail

address:

obruno@caltech.edu

(O.P. Bruno).

https://doi.org/10.1016/j.jcpx.2022.100110

2590-0552/

©

2022

The

Authors.

Published

by

Elsevier

Inc.

This

is

an

open

access

article

under

the

CC

BY-NC-ND

license

(

http://creativecommons.org/licenses/by-nc-nd/4.0/

).

O.P. Bruno, J.S. Hesthaven and D.V. Leibovici

Journal of Computational Physics: X 15 (2022) 100110

method

presented

in

this

paper

addresses

these

challenges.

In

particular,

in

view

of

its

use

of

certain

newly-designed

smooth

viscosity

windows

introduced

in

Section

3.3

,

the

method

avoids

the

introduction

of

roughness

in

the

viscosity

assignments

and

thus

it

yields

smooth

flows

away

from

shocks

and

other

flow

discontinuities.

In

addition,

the

underlying

FC-Gram

spectral

representations

enable

applicability

to

general

non-periodic

problems,

and,

in

view

of

the

weak

local

viscosity

as-

signments

used,

it

gives

rise

to

sharp

resolution

of

shocks.

As

a

result,

and

as

demonstrated

in

this

paper

via

application

to

a

range

of

well

known

1D

and

2D

shock-wave

test

configurations,

the

overall

FC-SDNN

approach

yields

accurate

and

essentially

oscillation-free

solutions

for

general

non-periodic

problems.

The

computational

solution

of

systems

of

conservation

laws

has

been

tackled

by

means

of

a

variety

of

numerical

meth-

ods,

including

low-order

finite

volume [

26

,

27

] and

finite

difference

methods

equipped

with

slope

limiters [

26

,

27

],

as

well

as

higher

order

shock-capturing

methods

such

as

the

ENO [

15

] and

WENO

schemes [

17

,

28

].

An

efficient

FC/WENO

hybrid

solver

was

proposed

in [

39

].

The

use

of

artificial

viscosity

as

a

computational

device

for

conservation

laws,

on

the

other

hand,

was

first

proposed

in [

37

,

45

] and

the

subsequent

contributions [

10

,

22

,

23

].

The

viscous

terms

proposed

in

these

papers,

which

incorporate

derivatives

of

the

square

of

the

velocity

gradient,

may

induce

oscillations

in

the

vicinity

of

shocks [

23

](since

the

velocity

itself

is

not

smooth

in

such

regions),

and,

as

they

do

not

completely

vanish

away

from

the

discontinuities,

they

may

lead

to

significant

approximation

errors

in

regions

where

the

fluid

velocity

varies

rapidly.

Reference [

32

]proposes

the

use

of

a

shock-detecting

sensor

in

order

to

localize

the

support

of

the

artificial

viscosity,

which

is

used

in

the

context

of

a

Discontinuous

Galerkin

scheme.

The

entropy

viscosity

method [

14

](EV)

incorporates

a

nonlinear

viscous

“entropy-residual”

term

which

essentially

van-

ishes

away

from

discontinuities—in

view

of

the

fact

that

the

flow

is

isentropic

over

smooth

flow

regions—and

which

is

thus

used

to

limit

non-zero

viscosity

assignments

to

regions

near

flow

discontinuities,

including

both

shock

waves

and

contact

discontinuities.

This

method,

however,

relies

on

several

problem-dependent

algorithmic

parameters

that

require

tuning

for

every

application.

Additionally,

this

approach

gives

rise

to

a

significant

amount

of

dissipation

even

away

from

shocks,

in

particular

in

the

vicinity

of

contact

discontinuities

and

regions

containing

fast

spatial

variation

in

the

flow-field

variables.

Considerable

improvements

concerning

this

issue

were

obtained

in

[

21

] (which

additionally

introduced

a

Hermite-based

method

to

discretize

the

hyperbolic

systems)

by

modifying

the

EV

viscosity

term.

Like

the

viscous

term

introduced

by [

45

],

the

EV

viscosity

assignments [

14

,

21

]are

themselves

discontinuous

in

the

vicinity

of

shocks,

and,

thus,

their

use

may

in-

troduce

spurious

oscillations.

The

C-method [

33

,

34

,

36

],

which

augments

the

hyperbolic

system

with

an

additional

equation

used

to

determine

a

spatio-temporally

smooth

viscous

term,

relies,

like

the

EV

method,

on

use

of

several

problem

dependent

parameters

and

algorithmic

variations.

Recently,

significant

progress

was

made

by

incorporating

machine

learning-based

techniques

(ML)

to

enhance

the

per-

formance

of

classical

shock

capturing

schemes [

9

,

35

,

38

,

43

].

The

approach [

9

,

35

,

38

]utilizes

ML-based

methods

to

detect

discontinuities

which

are

then

smeared

by

means

of

shock-localized

artificial

viscosity

assignments

in

the

context

of

various

discretization

methods,

including

Discontinuous

Galerkin

schemes [

9

,

35

] and

Fourier

spectral

schemes [

38

].

The

ML-based

approach

utilized

in [

43

],

in

turn,

modifies

the

finite

volume

coefficients

utilized

in

the

WENO5-JS

scheme

by

learning

small

perturbations

of

these

coefficients

leading

to

improved

accurate

representations

of

functions

at

cell

boundaries.

Like

the

strategy

underlying

the

contribution [

38

],

the

FC-SDNN

method

proposed

in

this

paper

relies

on

the

occurrence

of

Gibbs

oscillations

for

ML-based

shock

detection.

Unlike

the

however,

the

present

method

assigns

a

smooth

(albeit

also

shock-localized)

viscous

term.

In

view

of

its

smooth

viscosity

assignments

this

procedure

effectively

eliminates

Gibbs

oscillations

while

avoiding

introduction

of

the

flow-field

roughness

that

is

often

evidenced

by

the

serrated

level

sets

produced

by

other

methods.

In

view

of

its

use

of

FC-based

Fourier

expansions,

further,

the

proposed

algorithm

en-

joys

spectral

accuracy

away

from

shocks

(thus,

delivering,

in

particular,

essentially

vanishing

dispersion

in

such

regions;

see

Section

2.3

and

Fig.

6

)

while

enabling

solution

under

general

(non-periodic)

boundary

conditions.

Unlike

other

techniques,

finally,

the

approach

does

not

require

use

of

problem-dependent

algorithmic

parameters.

The

capabilities

of

the

proposed

algorithm

are

illustrated

by

means

of

a

variety

of

numerical

results,

in

one

and

two-

dimensional

contexts,

for

the

Linear

Advection,

Burgers

and

Euler

equations.

In

order

to

provide

a

useful

reference

point,

this

paper

also

presents

an

FC-based

version,

termed

FC-EV,

of

the

EV

algorithm [

14

].

(The

modified

version [

21

]of

the

entropy

viscosity

approach,

which

was

also

tested

as

a

possible

reference

solver,

was

not

found

to

be

completely

reliable

in

our

FC-based

context,

since

it

occasionally

led

to

spurious

oscillations

in

shock

regions

as

grids

were

refined,

and

the

corresponding

results

were

therefore

not

included

in

this

paper.)

We

find

that

the

FC-SDNN

algorithm

generally

provides

significantly

more

accurate

numerical

approximations

than

the

FC-EV,

as

the

localized

artificial

viscosity

in

the

former

approach

induces

a

much

lower

dissipation

level

than

the

latter.

This

paper

is

organized

as

follows.

After

necessary

preliminaries

are

presented

in

Section

2

(concerning

the

hyperbolic

problems

under

consideration,

as

well

as

the

Fourier

Continuation

method,

and

including

basic

background

on

the

artificial-

viscosity

strategies

we

consider),

Section

3

describes

the

proposed

FC-SDNN

approach.

A

link

to

a

git

repository

containing

basic

Matlab

implementations

and

test

codes

for

the

1D

FC

algorithm

is

provided

at

the

end

of

Section

2.3

.

Section

4

then

demonstrates

the

algorithm’s

performance

for

a

variety

of

non-periodic

linear

and

nonlinear

hyperbolic

problems.

In

particular,

cases

are

considered

for

the

linear

advection

equation,

Burgers

equation

and

Euler’s

equations

in

one-dimensional

and

two-dimensional

rectangular

and

non-rectangular

spatial

domains,

including

cases

in

which

shock

waves

meet

smooth

and

non-smooth

physical

boundaries.

2

O.P. Bruno, J.S. Hesthaven and D.V. Leibovici

Journal of Computational Physics: X 15 (2022) 100110

2.

Preliminaries

2.1.

Conservation

laws

This

paper

proposes

novel

spectral

methodologies,

applicable

in

general

non-periodic

contexts

and

with

general

boundary

conditions,

for

the

numerical

solution

of

conservation-law

equations

of

the

form

∂

t

e

(

x

,

t

)

+∇·

f

(

e

(

x

,

t

))

=

0

(1)

on

a

bounded

domain

⊂

R

q

,

where

e

:

×[

0

,

T

]

→

R

r

and

f

:

R

r

→

R

r

×

R

q

denote

the

unknown

solution

vector

and

a

(smooth)

convective

flux,

respectively.

The

proposed

spectral

approaches

are

demonstrated

for

several

equations

of

the

form (

1

),

including

the

Linear

Advection

equation

∂

u

∂

t

+

a

∂

u

∂

x

=

0

(2)

with

a

constant

propagation

velocity

a

,

where

we

have

e

=

u

,

and

f

(

u

)

=

au

;

the

one- and

two-dimensional

scalar

Burgers

equations

∂

u

∂

t

+

1

2

∂

u

2

∂

x

=

0

(3)

and

∂

u

∂

t

+

1

2

∂

u

2

∂

x

+

1

2

∂

u

2

∂

y

=

0

,

(4)

for

each

of

which

we

have

e

=

u

and

f

(

u

)

=

u

2

;

as

well

as

the

one- and

two-dimensional

Euler

equations

∂

t

⎛

⎝

ρ

u

E

⎞

⎠

+

∂

x

⎛

⎜

⎝

ρ

u

ρ

u

2

+

p

u

(

E

+

p

)

⎞

⎟

⎠

=

0

(5)

and

∂

t

⎛

⎜

⎝

ρ

u

ρ

v

E

⎞

⎟

⎠

+

∂

x

⎛

⎜

⎝

ρ

u

ρ

u

2

+

p

ρ

uv

u

(

E

+

p

)

⎞

⎟

⎠

+

∂

y

⎛

⎜

⎝

ρ

v

ρ

uv

ρ

v

2

+

p

v

(

E

+

p

)

⎞

⎟

⎠

=

0(6)

with

E

=

p

γ

−

1

+

1

2

ρ

|

u

|

2

,

(7)

for

each

of

which

we

have

e

=

(

ρ

,

ρ

u

,

E

)

T

,

f

(

e

)

=

(

ρ

u

,

ρ

u

⊗

u

+

p

I

,

u

(

E

+

p

))

T

.

(8)

Here

I

denotes

the

identity

tensor,

a

⊗

b

=

(

a

i

b

j

)

denotes

the

tensor

product

of

the

vectors

a

=

(

a

i

)

and

b

=

(

b

j

)

,

and

ρ

,

u

,

E

and

p

denote

the

density,

velocity

vector,

total

energy

and

pressure,

respectively.

The

speed

of

sound [

26

]

a

=

γ

p

ρ

(9)

for

the

Euler

equations

plays

important

roles

in

the

various

artificial

viscosity

assignments

considered

in

this

paper

for

Euler

problems

in

both

1D

and

2D.

Remark

1.

As

an

example

concerning

notational

conventions,

note

that

in

the

case

of

the

2D

Euler

equations,

for

which

f

is

given

by (

8

),

∇·

f

(

e

)

can

be

viewed

as

a

three

coordinate

vector

whose

first,

second

and

third

coordinates

are

a

scalar,

a

vector

and

a

scalar,

respectively.

Using

a

super-index

notation

for

the

velocity

vector

u

(i.e.,

u

1

=

u

and

(

u

1

,

u

2

)

=

(

u

,

v

)

for

the

1D

and

2D

equations,

respectively),

together

with

the

Einstein

summation

convention,

these

three

components

are

respectively

given

by

∇·

(

ρ

u

)

=

∂

i

(

ρ

u

i

)

,

∇·

(

ρ

u

⊗

u

+

p

I

)

j

=

∂

i

(

ρ

u

j

u

i

+

p

)

and

∇·

((

E

+

p

)

u

)

=

∂

i

((

E

+

p

)

u

i

)

.

3

O.P. Bruno, J.S. Hesthaven and D.V. Leibovici

Journal of Computational Physics: X 15 (2022) 100110

2.2.

Artificial

viscosity

As

is

well

known,

the

shocks

and

other

flow

discontinuities

that

arise

in

the

context

of

nonlinear

conservation

laws

of

the

form (

1

)give

rise

to

a

number

of

challenges

from

the

point

of

view

of

computational

simulation.

In

particular,

in

the

framework

of

classical

finite

difference

methods

as

well

as

Fourier

spectral

methods,

such

discontinuities

are

associated

with

the

appearance

of

spurious

“Gibbs

oscillations”.

Artificial

viscosity

methods

aim

at

tackling

this

difficulty

by

considering,

instead

of

the

inviscid

equations (

1

),

certain

closely

which

include

viscous

terms

containing

second

order

spatial

derivatives.

Provided

the

viscous

terms

are

adequately

chosen

and

sufficiently

small,

the

resulting

solutions,

which

are

smooth

functions

on

account

of

viscosity,

approximate

well

the

desired

(discontinuous)

inviscid

solutions.

In

general

terms,

the

viscous

equations

are

obtained

by

adding

a

viscous

term

of

the

form

∇·

f

visc

[

e

]

to

the

right

hand

side

of

(

1

),

where

the

“viscous

flux”

operator

f

visc

[

e

]=

μ

[

e

]

D

[

e

]

,

(10)

(which,

for

a

given

vector-valued

function

e

(

x

,

t

)

,

produces

a

vector-valued

function

f

visc

[

e

]

(

x

,

t

)

defined

in

the

complete

computational

domain),

is

given

in

terms

of

a

certain

“viscosity”

operator

μ

[

e

]

(

x

,

t

)

(which

may

or

may

not

include

deriva-

tives

of

the

flow

variables

e

),

and

a

certain

matrix-valued

first

order

differential

operator

D

.

Once

such

a

viscous

term

is

included,

the

viscous

equation

∂

e

(

x

,

t

)

∂

t

+∇·

f

(

e

(

x

,

t

))

=∇·

f

visc

[

e

]

(

x

,

t

)

(11)

results.

Per

the

discussion

in

Section

1

,

this

paper

exploits

and

extends,

in

the

context

of

the

Fourier-Continuation

discretizations,

two

different

approaches

to

viscosity-regularization—each

one

resulting

from

a

corresponding

selection

of

the

operators

μ

and

D

.

One

of

these

approaches,

the

EV

method,

produces

a

viscosity

assignment

μ

[

e

]

(

x

,

t

)

on

the

basis

of

certain

differential

and

algebraic

operations

together

with

a

number

of

tunable

problem-dependent

parameters

that

are

specifically

designed

for

each

particular

equation

considered,

as

described

in

Section

2.2.2

.

The

resulting

viscosity

values

μ

[

e

]

(

x

,

t

)

are

highest

in

a

vicinity

of

discontinuity

regions

and

decrease

rapidly

away

from

such

regions.

The

neural-network

approach

introduced

in

Section

2.2.1

,

in

turn,

uses

machine

learning

methods

to

pinpoint

the

location

of

discontinuities,

and

then

produces

a

viscosity

function

whose

support

is

restricted

to

a

vicinity

of

such

discontinuity

locations.

As

a

significant

advantage,

the

neural-network

method,

which

does

not

require

use

of

tunable

parameters,

is

essentially

problem

independent,

and

it

can

use

a

single

pre-trained

neural

network

for

all

the

equations

considered.

Details

concerning

these

two

viscosity-assignment

methods

considered

are

provided

in

what

follows.

2.2.1.

Artificial

viscosity

via

shock-detecting

neural

network

(SDNN)

The

SDNN

approach

proposed

in

this

paper

is

based

on

the

neural-network

strategy

introduced

in [

38

]for

detection

of

discontinuities

on

the

basis

of

Gibbs

oscillations

in

Fourier

series,

together

with

a

novel

selection

of

the

operator

μ

in (

10

)

that

yields

spatially

localized

but

smooth

viscosity

assignments:

per

the

description

in

Section

3.3

,

the

FC-SDNN

viscosity

μ

[

e

]

(

x

,

t

)

is

a

smooth

function

that

vanishes

except

in

narrow

regions

around

flow

discontinuities.

The

differential

operator

D

,

in

turn,

is

simply

given

by

D

[

e

]

(

x

,

t

)

=∇

(

e

(

x

,

t

)),

(12)

where

the

gradient

is

computed

component-wise.

As

indicated

in

Section

1

,

the

smoothness

of

the

proposed

viscosity

assignments

is

inherited

by

the

resulting

flows

away

from

flow

discontinuities,

thus

helping

eliminate

the

serrated

level-set

lines

that

are

ubiquitous

in

the

flow

patterns

produced

by

other

methods.

2.2.2.

Entropy

viscosity

methodology

(EV)

The

operators

μ

and

D

employed

by

the

EV

approach [

14

]are

defined

in

terms

of

a

number

of

problem

dependent

functions,

vectors

and

operators.

Indeed,

starting

with

an

equation

dependent

entropy

pair

(

η

,

ν

)

where

η

is

a

scalar

function

and

ν

is

a

vector

of

the

same

dimensionality

as

the

velocity

vector,

the

EV

approach

utilizes

an

associated

scalar

entropy

residual

operator

R

EV

[

e

]

(

x

,

t

)

=

∂

η

(

e

(

x

,

t

))

∂

t

+∇·

ν

(

e

(

x

,

t

))

(13)

together

with

a

function

C

=

C

(

e

)

the

local

wave

speed,

and

a

normalization

operator

N

=

N

[

e

]

(

x

,

t

)

obtained

from

the

function

η

.

In

practice,

reference [

14

]proposes

η

(

e

)

=

u

2

,

ν

(

e

)

=

a

u

2

and

C

(

e

)

=

a

for

the

Linear

Advection

equation (

2

),

η

(

e

)

=

u

2

,

ν

(

e

)

=

u

3

and

C

(

e

)

=

u

for

the

1D

and

2D

Burgers

equations (

3

) and (

4

),

and

η

(

e

)

=

ρ

γ

−

1

log

(

p

/

ρ

γ

)

,

ν

(

e

)

=

u

ρ

γ

−

1

log

(

p

/

ρ

γ

)

and

C

(

e

)

=|

u

|

+

a

(where

a

denotes

the

speed

of

sound (

9

))

for

the

1D

and

2D

Euler

equations (

5

) and (

6

).

As

for

the

4

O.P. Bruno, J.S. Hesthaven and D.V. Leibovici

Journal of Computational Physics: X 15 (2022) 100110

normalization

operator,

reference [

14

]proposes

N

=

1for

the

Euler

equations

and

N

[

e

]

(

x

,

t

)

=|

η

(

e

)(

x

,

t

)

−

η

(

e

)(

t

)

|

for

the

Linear

advection

and

Burgers

equations,

where

η

(

e

)(

t

)

denotes

the

spatial

average

of

η

(

e

)

at

time

t

.

For

a

numerical

discretization

with

maximum

spatial

mesh

size

h

,

the

EV

viscosity

function

is

defined

by

μ

[

e

]

(

x

,

t

)

=

min

(

μ

max

[

e

]

(

t

),

μ

E

[

e

]

(

x

,

t

))

(14)

where

the

maximum

viscosity

μ

max

is

given

by

μ

max

[

e

]

(

t

)

=

c

max

h

max

x

∈

|

C

(

e

(

x

,

t

))

|

(15)

and

where

μ

E

[

e

]

(

x

,

t

)

=

c

E

h

2

|

R

EV

[

e

]

(

x

,

t

)

|

N

[

e

]

(

x

,

t

)

.

(16)

In

particular,

the

EV

viscosity

function

depends

on

two

parameters,

c

max

and

c

E

,

both

of

size

O

(

1

)

,

that,

following [

14

],

are

to

be

tuned

to

each

particular

problem.

Finally,

the

EV

differential

operator

D

for

the

Linear

Advection

and

Burgers

equations

is

defined

by

D

[

e

]

(

x

,

t

)

=∇

(

e

(

x

,

t

)),

(17)

while

for

the

Euler

equations

it

is

given

by

D

[

e

]

(

x

,

t

)

=

⎛

⎜

⎝

0

1

2

(

∇

u

+

(

∇

u

)

T

)

1

2

(

∇

u

+

(

∇

u

)

T

)

u

+

κ

∇

(

p

/

ρ

)

⎞

⎟

⎠

(18)

where,

using

once

again

the

Einstein

convention,

(

∇

u

+

(

∇

u

)

T

)

u

i

=

(∂

i

u

j

+

∂

j

u

i

)

u

j

,

and

where

κ

=

P

γ

−

1

μ

,

with

the

Prandtl

number

P

taken

to

equal

1.

2.3.

Spatial

approximation

via

Fourier

Continuation

The

straightforward

Fourier-based

discretization

of

nonlinear

conservation

laws

generally

suffers

from

crippling

Gibbs

oscillations

resulting

from

two

different

sources:

the

physical

flow

discontinuities,

on

one

hand,

and

the

overall

generic

non-

periodicity

of

the

flow

variables,

on

the

other.

Unlike

the

Gibbs

ringing

in

flow-discontinuity

regions,

the

ringing

induced

by

lack

of

periodicity

is

not

susceptible

to

treatment

via

artificial

viscosity

assignments

of

the

type

discussed

in

2.2

.

In

order

to

tackle

this

difficulty

we

resort

to

use

of

the

Fourier

Continuation

(FC)

method

for

equispaced-grid

spectral

approximation

of

non-periodic

functions.

The

basic

FC

algorithm [

1

],

called

FC-Gram

in

view

of

its

reliance

on

Gram

polynomials

for

near-boundary

approxima-

tion,

constructs

an

accurate

Fourier

approximation

of

a

given

generally

non-periodic

function

F

defined

on

a

given

one-

dimensional

interval—which,

for

definiteness,

is

assumed

in

this

section

to

equal

the

unit

interval

[

0

,

1

]

.

The

Fourier

contin-

uation

algorithm

relies

on

use

of

the

function

values

F

j

=

F

(

x

j

)

of

the

function

F

:[

0

,

1

]

→

R

at

N

points

x

j

=

jh

∈[

0

,

1

]

(

h

=

1

/(

N

−

1

)

to

produce

a

function

F

c

(

x

)

=

M

k

=−

M

ˆ

F

c

k

exp

(

2

π

ikx

/β)

(19)

which

is

defined

(and

periodic)

in

an

interval

[

0

,

β

]

that

strictly

contains

[

0

,

1

]

,

where

ˆ

F

c

k

denote

the

FC

coefficients

of

F

and

where,

as

detailed

below,

M

is

an

integer

that,

for

N

large,

is

close

to

(but

different

from)

the

integer

N

/

2

.

In

order

to

produce

the

FC

function

F

c

,

the

FC-Gram

algorithm

first

uses

two

subsets

of

the

function

values

in

the

vector

F

=

(

F

0

,

...,

F

N

−

1

)

T

(namely

the

function

values

at

the“matching

points”

{

x

0

,

..,

x

d

−

1

}

and

{

x

N

−

d

,

...,

x

N

−

1

}

located

in

small

matching

subintervals

[

0

,

]

and

[

1

−

,

1

]

of

length

=

(

d

−

1

)

h

near

the

left

and

right

ends

of

the

interval

[

0

,

1

]

,

where

d

is

a

small

integer

independent

of

N

),

to

produce,

at

first,

a

discrete

(but

“smooth”)

periodic

extension

vector

F

c

of

the

vector

F

.

Indeed,

using

the

matching

point

data,

the

FC-Gram

algorithm

produces

and

appends

a

number

C

of

continuation

function

values

in

the

interval

[

1

,

β

]

to

the

data

vector

F

,

so

that

the

extension

F

c

transitions

smoothly

from

F

N

−

1

back

to

F

0

,

as

depicted

in

Fig.

1

.

(The

FC

method

can

also

be

applied

on

the

basis

of

certain

combinations

of

function

values

and

derivatives

by

constructing

the

continuation

vector

F

c

,

as

described

below

in

this

section,

for

a

given

vector

F

=

(

F

0

,

...,

F

N

−

2

,

F

N

−

1

)

T

,

where

F

j

≈

F

(

x

j

)

for

1

≤

j

≤

N

−

2 and

where

F

N

−

1

≈

F

(

x

N

−

1

)

.

Such

a

procedure

enables

imposition

of

Neumann

boundary

conditions

in

the

context

of

the

FC

method.)

The

resulting

vector

F

c

can

be

viewed

as

a

discrete

set

of

values

of

a

smooth

and

periodic

function

which

can

be

used

to

produce

the

Fourier

continuation

function

F

c

via

an

application

of

the

FFT

algorithm.

The

function

F

c

provides

a

spectral

approximation

of

F

throughout

the

interval

[

0

,

1

]

which

does

not

suffer

from

either

Gibbs-ringing

or

the

associated

interval-wide

accuracy

degradation.

Throughout

this

5