FC-based shock-dynamics solver with neural-network
localized artificial-viscosity assignment
Oscar P. Bruno
∗
Jan S. Hesthaven
†
Daniel V. Leibovici
∗
Abstract
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 previous shock capturing schemes
and artificial viscosity techniques, the combined FC-SDNN strategy effectively controls spuri-
ous 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 previ-
ous entropy 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.
Keywords:
Machine learning, Neural networks, Shocks, Artificial viscosity, Conservation laws,
Fourier continuation, Non-periodic domain, Spectral method
∗
Computing and Mathematical Sciences, Caltech, Pasadena, CA 91125, USA
†
Computational Mathematics and Simulation Science (MCSS), Ecole Polytechnique Federale de Lausanne, CH-
1015 Lausanne, Switzerland
1
arXiv:2111.01315v1 [math.NA] 2 Nov 2021
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 relies on use of
the FC-Gram Fourier Continuation method [1, 2, 5, 28] 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 [37]. The neural network ap-
proach [37] itself utilizes Fourier series to discretize the gas dynamics and related equations, 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 boundaries can-
not 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 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 underly-
ing FC-Gram spectral representations enable applicability to general non-periodic problems, and,
in view of the weak local viscosity assignments 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 methods, including low-order finite volume [25,26] and finite difference methods
equipped with slope limiters [25, 26], as well as higher order shock-capturing methods such as the
ENO [14] and WENO schemes [16, 27]. An efficient FC/WENO hybrid solver was proposed in [38].
The use of artificial viscosity as a computational device for conservation laws, on the other hand,
was first proposed in [36,44] and the subsequent contributions [9,21,22]. 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 [22] (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 were the fluid velocity varies rapidly. Reference [31] 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 [13] (EV) incorporates a nonlinear viscous “entropy-residual” term
which essentially vanishes 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 [20] (which additionally introduced a Hermite-based method to discretize the hyperbolic systems)
by modifying the EV viscosity term. Like the viscous term introduced by [44], the EV viscosity
assignments [13, 20] are themselves discontinuous in the vicinity of shocks, and, thus, their use may
2
introduce spurious oscillations. The C-method [32, 33, 35], 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 performance of classical shock capturing schemes [8, 34, 37, 42]. The ap-
proach [8, 34, 37] 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 [8, 34] and Fourier spectral schemes [37]. The
ML-based approach utilized in [42], in turn, modifies the finite volume coefficients utilized in the
WENO5-JS scheme by learning small perturbations of these coefficients leading to improved accu-
rate representations of functions at cell boundaries.
Like the strategy underlying the contribution [37], the FC-SDNN method proposed in this paper
relies on the occurrence of Gibbs oscillations for ML-based shock detection. Unlike the previous
approach, 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 enjoys spectral accuracy away from shocks (thus, delivering, in particular,
essentially vanishing dispersion in such regions; see Section 2.3 and Figure 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 [13]. (The modified version [20] 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 (con-
cerning 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. 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.
3
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), in-
cluding 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
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 [25]
a
=
√
γp
ρ
(9)
4
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 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
).
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 compu-
tational 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 related equations 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 derivatives 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 result-
ing 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 vis-
cosity 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.
5
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 [37] 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 [13] 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
) related to the local wave speed, and a normalization operator
N
=
N
[
e
](
x,t
) obtained from the function
η
.
In practice, reference [13] proposes
η
(
e
) =
u
2
2
,
ν
(
e
) =
a
u
2
2
and
C
(
e
) =
a
for the Linear Advection
equation (2),
η
(
e
) =
u
2
2
,
ν
(
e
) =
u
3
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 normalization operator,
reference [13] proposes
N
= 1 for 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 [13], 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)
6
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 the Einstein notation
{
(
∇
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
Fourier Continuation spatial approximation
The straightforward Fourier-based discretization of nonlinear conservation laws generally suffers
from crippling Gibbs oscillations resulting from two different sources: the physical flow discontinu-
ities, 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.
d = 5
1 -
C = 27
d = 5
Figure 1: Fourier Continuation of the non-periodic function
F
(
x
) =
x
on the interval [0
,
1]. With
reference to the text, the red triangles (resp. squares), represent the
d
= 5 left (resp. right)
matching points, while the blue circles represent the
C
= 27 continuation points.
The basic FC algorithm [1], called FC-Gram in view of its reliance on Gram polynomials for
near-boundary approximation, 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 continuation 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)
7
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
b
N/
2
c
.
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 Figure 1. (The
FC method can also be applied on the basis certain combinations of function values and derivatives
by constructing the continuation vector
F
c
, as described below in this section, on the basis of
e.g. the 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 paper we assume, for simplicity, that
N
+
C
is an odd integer, and, thus, the resulting series has bandwidth
M
=
N
+
C
−
1
2
; consideration of even
values of
N
+
C
would require a slight modification of the index range in (19).
To obtain the necessary discrete periodic extension
F
c
, the FC-Gram algorithm first produces two
polynomial interpolants, one per matching subinterval, using, as indicated above, a small number
d
of function values or a combination of function values and a derivative near each one of the endpoints
of the interval [0
,
1]. This approach gives rise to high-order interpolation of the function
F
over the
matching intervals [0
,
∆] and [1
−
∆
,
1]. The method for evaluation of the discrete periodic extension
is based on a representation of these two polynomials in a particular orthogonal polynomials basis
(the Gram polynomials), for each element of which the algorithm utilizes a precomputed smooth
function which blends the basis polynomial to the zero function over the distance
β
−
1 [1, 2].
Certain simple operations involving these “blending to zero” functions are then used, as indicated
in these references and as illustrated in Figure 1, to obtain smooth transitions-to-zero from the left-
most and right-most function values to the extension interval [1
,β
]. The values of this transition
function at the points
N/
(
N
−
1)
,
(
N
+ 1)
/
(
N
−
1)
,...,
(
N
+
C
−
1)
/
(
N
−
1) provide the necessary
C
additional point values from which, as mentioned above, the discrete extension
F
c
is obtained.
The continuation function
F
c
then easily results via an application of the FFT algorithm to the
function values
F
c
in the interval [0
,β
].
The discrete continuation procedure can be expressed in the matricial form
F
c
=
(
F
A
`
Q
T
F
`
+
A
r
Q
T
F
r
)
(20)
where the
d
-dimensional vectors
F
`
and
F
r
contain the point values of
F
at the first and last
d
discretization points in the interval [0
,
1], respectively; where
Q
is a
d
×
d
matrix, whose columns
contain the point values of the elements of the Gram polynomial bases on the left matching intervals;
and where
A
`
and
A
r
are
C
×
d
, matrices containing the
C
values of the blended-to-zero Gram
8
polynomials in the left and right Gram bases, respectively. These small matrices can be computed
once and stored on disc, and then read for use to produce FC expansions for functions
G
: [
a,b
]
→
R
defined on a given 1D interval [
a,b
], via re-scaling to the interval [0
,
1].
A minor modification of the procedure presented above suffices to produce a Fourier continuation
function on the basis of data points at the domain interior and a derivative at interval endpoints.
For example, given the vector
F
= (
F
0
,...,F
N
−
2
,F
′
N
−
1
)
T
, using an adequately modified version
̃
Q
of
the matrix
Q
, an FC series
F
c
(
x
) can be produced which matches the function values
F
0
,...,F
N
−
2
at
x
=
x
0
,...,x
N
−
2
, and whose derivative equals
F
′
N
−
1
at
x
N
−
1
. The matrix
̃
Q
is obtained by using
the matrix
Q
to obtain a value
F
N
−
1
such that the derivative
F
′
(
x
N
−
1
) equals the given value
F
′
N
−
1
.
Full details in this regard can be found in [2, Sec. 3.2].
Clearly, the approximation order of the Fourier Continuation method (whether derivative values
or function values are prescribed at endpoints) is restricted by the corresponding order
d
of the
Gram polynomial expansion, which, as detailed in various cases in Section 4, is selected as a small
integer:
d
= 2 or
d
= 5. The relatively low order of accuracy afforded by the
d
= 2 selection, which
must be used in some cases to ensure stability, is not a matter of consequence in the context of the
problems considered in the present paper, where high orders of accuracy are not expected from any
numerical method on account of shocks and other flow discontinuities. Importantly, even in this
context the FC method preserves one of the most significant numerical properties of Fourier series,
namely its extremely small numerical dispersion. In fact, with exception of the cyclic advection
example presented in Section 4.1.2, for which errors can accumulate on account of the spatio-
temporal periodicity, for all cases in Section 4 for which both the
d
= 2 and
d
= 5 simulations were
performed (which include those presented in Sections 4.1.1 (1D linear advection), 4.2.2 (2D Burgers
equation) and 4.3.1 (1D Euler equations), the lower and higher order results obtained were visually
indistinguishable.
The low dispersion character resulting from use of the FC method is demonstrated in Figure 6,
which displays solutions produced by means of two different methods, namely, the FC-based order-5
FC-SDNN algorithm (Section 3) and the order-6 centered finite-difference scheme (both of which
use the SSPRK-4 time discretization scheme), for a linear advection problem. The FC-SDNN
solution presented in the figure does not deteriorate even for long propagation times, thus illustrating
the essentially dispersion-free character of the FC-based approach. The finite-difference solution
included in the Figure, in turn, does exhibit clear degradation with time, owing to the dispersion
and diffusion effects associated with the underlying finite difference discretization.
3 FC-based time marching under neural network-controlled
artificial viscosity
3.1 Spatial grid functions and spatio-temporal FC-based differentiation
We consider in this work 1D problems on intervals
I
= [
ξ
`
,ξ
r
] as well as 2D problems on open
domains Ω contained in rectangular regions
I
×
J
, where
I
= [
ξ
`
,ξ
r
] and
J
= [
ξ
d
,ξ
u
] (
ξ
`
< ξ
r
,
ξ
d
< ξ
u
). Using a spatial meshsize
h
, the spatially discrete vectors of unknowns and certain related
flow quantities will be represented by means of scalar and vector grid functions defined on 1D or
2D discretization grids of the form
G
=
{
x
i
:
x
i
=
x
0
+
ih, i
= 0
,...,N
−
1
}
(
x
0
=
ξ
`
, x
N
−
1
=
ξ
r
)
,
9
and
G
=
Ω
∩
{
(
x
i
,y
j
) :
x
i
=
x
0
+
ih,y
j
=
y
0
+
jh,
0
≤
i
≤
N
1
−
1
,
0
≤
j
≤
N
2
−
1
}
,
respectively. Here
Ω denotes the closure of Ω,
x
0
=
ξ
`
,
x
N
1
−
1
=
ξ
r
,
y
0
=
ξ
d
and
y
N
2
−
1
=
ξ
u
. In
either case a function
b
:
G
→
R
q
will be called a “
q
-dimensional vector grid function”. Letting
I
=
{
(
i,j
)
∈{
0
,...,N
1
−
1
}×{
0
,...,N
2
−
1
}
: (
x
i
,y
j
)
∈
G
}
,
we will also write
b
(
x
i
) =
b
i
(0
≤
i
≤
N
−
1) and
b
(
x
i
,y
j
) =
b
ij
((
i,j
)
∈I
). The set of
q
-dimensional
vector grid functions defined on
G
will be denoted by
G
q
.
It is important to mention that, although the two-dimensional setting described above does
not impose any restrictions on the character of the domain Ω, for simplicity, the FC-SDNN solver
presented in this paper assumes that the boundary of
Ω is given by a union of horizontal and
vertical straight segments, each one of which runs along a Cartesian discretization line; see e.g.
the Mach 3 forward-facing step case considered in Figure 16a. Extensions to general domains
Ω, which could rely on either an embedded-boundary [5, 6, 28] approach, or an overlapping patch
boundary-conforming curvilinear discretization strategy [1, 2, 4], is left for future work.
A spatially-discrete but time-continuous version of the solution vector
e
(
x
,t
) considered in
Section 2 for 1D problems (resp. 2D problems) can be viewed as a time-dependent
q
-dimensional
vector grid function
e
i
=
e
i
(
t
) (resp.
e
ij
=
e
ij
(
t
)). Using
e
h
=
e
h
(
t
) to refer generically to the
1
D
and 2
D
time-dependent grid functions
e
i
and
e
ij
, the semidiscrete scheme for equation (11)
becomes
d
e
h
(
t
)
dt
=
L
[
e
h
(
t
)]
,
(21)
where
L
denotes a consistent discrete approximation of the spatial operator in (11).
The discrete time evolution of the problem, on the other hand, is produced, throughout this
paper, by means of the 4-th order strong stability preserving Runge-Kutta scheme (SSPRK-4) [12]—
which, while not providing high convergence orders for the non-smooth solutions considered in this
paper, does lead to low temporal dispersion and diffusion over smooth space-time regions of the
computational domain. The corresponding time step is selected adaptively at each time-step
t
=
t
n
according to the expression
∆
t
=
CFL
π
(
max
x
∈
Ω
|
S
[
e
](
x
,t
)))
|
h
+
max
x
∈
Ω
μ
[
e
](
x
,t
)
h
2
)
.
(22)
Here CFL is a constant parameter that must be selected for each problem considered (as illustrated
by the various selections utilized in Section 4), and
μ
[
e
](
x
,t
) and
S
=
S
[
e
](
x
,t
) denote the artificial
viscosity (equations (37) and (42)) and a
maximum wave speed bound
(MWSB) operator (which
must be appropriately selected for each equation; see Section 3.3) at the spatio-temporal point (
x
,t
).
(To avoid confusion we emphasize that equation (22) utilizes the
maximum
value for all
x
∈
Ω of
the selected bound
S
[
e
](
x
,t
) on the
maximum
wave speed.)
To obtain FC-based approximate derivatives of a function
F
:
K
→
R
defined on a one-
dimensional interval
K
= [
x
0
,x
N
−
1
], whose values (
F
0
,F
1
,...,F
N
−
1
)
T
are given on the uniform
mesh
{
x
0
,x
1
,...,x
N
−
1
}
, the interval
K
is re-scaled to [0
,
1] and the corresponding continuation
10
function
F
c
is obtained by means of the FC-Gram procedure described in Section 2.3. The ap-
proximate derivatives at all mesh points are then obtained by applying the IFFT algorithm to the
Fourier coefficients
(
ˆ
F
c
)
′
k
=
2
πik
β
ˆ
F
c
k
.
(23)
of the derivative of the series (19) and re-scaling back to the interval
K
.
All of the numerical derivatives needed to evaluate the spatial operator
L
[
e
h
(
t
)] are obtained
via repeated application of the 1D FC differentiation procedure described above. For a function
F
=
F
(
x,y
) defined on a two-dimensional domain Ω and whose values
F
ij
((
i,j
)
∈I
) are given on
a grid
G
of the type described above in this section, for example, partial derivatives with respect to
x
along the line
y
=
y
j
0
for a relevant value of
j
0
are obtained by differentiation of the FC expansion
obtained for the function values (
F
(
x
i
,y
j
0
))
i
for integers
i
such that (
i,j
0
)
∈I
. The
y
differentiation
process proceeds similarly. Mixed derivatives, finally, are produced by successive application of the
x
and
y
differentiation processes. Details concerning the filtered derivatives used in the proposed
scheme are provided in Section 3.4.
The boundary conditions of Dirichlet and Neumann considered in this paper are imposed as
part of the differentiation process described above. Dirichlet boundary conditions at time
t
n,ν
(
t
n
< t
n,ν
≤
t
n
+1
) corresponding to the
ν
-th SSPRK-4 stage (
ν
= 1
,...,
4) for the time-step starting
at
t
=
t
n
, are simply imposed by overwriting the boundary values of the unknown solution vector
e
h
obtained at time
t
=
t
n,ν
with the given boundary values at that time, prior to the evaluation of the
spatial derivatives needed for the subsequent SSPRK-4 stage. Neumann boundary conditions are
similarly enforced by constructing appropriate continuation vectors (Section 2.3) after each stage of
the SSPRK-4 scheme on the basis of the modified pre-computed matrix
̃
Q
mentioned in Section 2.3.
It is known that enforcement of the given physical boundary conditions at intermediate Runge-
Kutta stages, which is referred to as the “conventional method” in [7], may lead to a reduced
temporal order of accuracy at spatial points in a neighborhood of the boundary of the domain
boundary. This is not a significant concern in the context of this paper, where the global order of
accuracy is limited in view of the discontinuous character of the solutions considered. Alternative
approaches that preserve the order of accuracy for smooth solutions, such as those introduced
in [7,30], could also be used in conjunction with the proposed approach. Another alternative, under
which no boundary conditions are enforced at intermediate Runge-Kutta stages [19], can also be
utilized in our context, but we have found the conventional method leads to smoother solutions
near boundaries.
3.2 Neural network-induced smoothness-classification
3.2.1 Smoothness-classification operator and data pre-processing
The method described in the forthcoming Section 3.3 for determination of the artificial viscosity
values
μ
[
e
](
x
,t
) (cf. also Section 2.2) relies on the “degree of smoothness” of a certain function
Φ(
e
)(
x
,t
) (called the “proxy variable”) of the unknown solution vector
e
. In detail, following [37],
in this paper a proxy variable Φ(
e
) is used, which equals the velocity
u
, Φ(
e
) =
u
, (resp. the Mach
number, Φ(
e
) =
‖
u
‖
√
ρ
γp
) for equations (2) through (4) (resp. equations (5) and (6)). The degree of
smoothness of the function Φ(
e
) at a certain time
t
is characterized by a smoothness-classification
operator
τ
=
τ
[Φ(
e
)] that analyzes the oscillations in an FC expansion of Φ(
e
)—which is itself
obtained from the discrete numerical values
φ
= Φ(
e
h
), so that, in particular,
τ
[Φ(
e
)] = ̃
τ
[
φ
] for
11