Abstract: Events in turbulent flows computed by direct numerical simulation (DNS) are often calibrated with
properties based on homogeneous isotropic turbulence, advanced by Kolmogorov, and given in Turbulence, U.
Frisch, Cambridge Univ. Press, UK (1995). However, these computational procedures are not calibrated using
numerical analyses in order to assess their strengths and weaknesses for DNS. This is with the exception in "A
critical assessment of simulations for transitional and turbulence flows- Sengupta, T.K., In Proc. of IUTAM
Symp. on Advances in Computation, Modeling and Control of Transitional and Turbulent Flows, 491-532,
(eds. Sengupta, Lele, Sreenivasan, Davidson), WSPC, Singapore (2016)", where such a calibration has been
advocated for numerical schemes using global spectral analysis (GSA) for the convection equation. In recent
times, due to growing computing power, simulations have been reported using pseudo-spectral methods, with
spatial discretization performed in Fourier spectral space and time-integration by multi-stage Runge-Kutta (RK)
methods. Here, we perform GSA of Fourier spectral methods for the first time with RK2 and other multi-
stage Runge-Kutta methods using the model linear convection and convection-diffusion equations. With the
help of GSA, various sources of numerical errors are quantified. The major limitations of the RK2-Fourier
spectral method are demonstrated for DNS and alternate choices are presented. We specify optimal parameters to
achieve the best possible accuracy for simulations. There is a one-to-one correspondence of numerical solutions
obtained by linear convection-diffusion equation and nonlinear Navier-Stokes equation with respect to numerical
parameters. This enables us to investigate the capabilities of the numerical methods for DNS.
Received: March 11, 2021. Revised: October 15, 2021. Accepted: December 16, 2021. Published: January 19, 2022.
Numerical Methods
Key-Words: Computer Simulation, Numerical Simulation, Computer Applications, Pseudo-Spectral Methods,
1 Introduction
Direct numerical simulations (DNS) of turbulent
ows have been attempted for decades, with one of
the early works being reported [19] for a Reynolds
number (based on Taylor’s microscale (
λ
)) of
R
λ
=35, using a 323computational periodic box.
This DNS of three-dimensional (3D) homogeneous
isotropic turbulence was compared with the pre-
dictions of the direct-interaction turbulence the-
ory. With increase in computing power, the
recorded R
λ
increased. Kaneda & Ishihara [15] re-
ported results showing a R
λ
=1200 using a 40963
periodic box for 3D homogeneous isotropic turbu-
lence simulations. These activities progressed to
a state where results are reported for R
λ
=1300
using a 122883periodic computational box for the
simulation of 3D homogeneous isotropic turbu-
lence [5]. A historic account and perspectives of
homogeneous isotropic turbulence are detailed in
[13]. The aforementioned simulations for homo-
geneous isotropic turbulence have been performed
using the incompressible Navier-Stokes equation
(INSE), given by,
∂
u
∂
t+u·
∂
u
∂
x=−∇P/
ρ
+
ν
∇2u+f,(1)
where uis the solenoidal velocity eld, Pis the
pressure,
ρ
is the density of the uid,
ν
is the
kinematic viscosity, and fis a forcing term im-
posed at a large length scale. Such a forcing term
being added to the INSE has been weakly justied
as necessary to ensure statistical stationarity of
the computed turbulent signal [12] for simulation
of homogeneous isotropic turbulence by the RK2-
Fourier spectral method, as in [5]. This method
was originally used in a DNS code [22] with the
added forcing term in INSE. The same version of
the code continues to be used by a large num-
ber of researchers [14, 4, 16, 6, 7, 37, 11, 10, 20]
Analysis of Pseudo-spectral Methods Used for Numerical Simulations
of Turbulence
TAPAN K. SENGUPTA1, VAJJALA K. SUMAN2,3, PRASANNABALAJI SUNDARAM2, ADITI
SENGUPTA1
1Department of Mechanical Engineering, Indian Institute of Technology (ISM) Dhanbad,
Dhanbad-826004, INDIA
2High Performance Computing Laboratory, Indian Institute of Technology Kanpur,
Kanpur-208 016, INDIA
3Computational & Theoretical Fluid Dynamics Division, CSIR-NAL,
Bangalore-560017, INDIA
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
9
Volume 10, 2022
among many others. This approach of using ad-
ditional forcing has been critiqued [2] for not only
this canonical homogeneous problem, but also for
so-called DNS of geophysical and engineering ows
[32, 3, 31]. For transitional ows, there are many
solutions of INSE which do not use such artices
and can be viewed truly as DNS [21, 26, 2].
However, it is noted that such simulations of-
ten display issues such as solution “blow-up” after
a nite time (especially in the limit of vanishing
ν
). Also, it has been observed that, “experimental
measurements in homogeneous turbulence at high
Reynolds numbers show the unequivocal presence
of a “bottleneck” eect [17]” that has encouraged
the authors in the reference to introduce hypervis-
cosity in the INSE as,
∂
u
∂
t+u·
∂
u
∂
x=−∇P/
ρ
+ (−1)h+1
ν
h∇2hu+f,(2)
where
ν
his the specied constant hyperviscosity
coecient, and fis the forcing function. The orig-
inal formulation of the INSE can be retrieved by
setting h=1and f=0. It was noted [17] that
the “(use of fis unnatural (as are most other ways
of forcing turbulence). However, we are primarily
concerned with bottleneck eects on energy spec-
tra, i.e., we investigate one particular characteris-
tic of small-scale turbulence. In this case, use of
a large-scale forcing (in order to analyze statisti-
cally stationary rather than decaying turbulence)
is justiable on the grounds that the details of the
forcing have little eect on the small-scale statis-
tics.”
Bracco & McWilliams [3] simulated the ho-
mogeneous stationary 2D turbulence problem by
solving the vorticity transport equation given by,
D
ω
Dt =D+F,(3)
where Fis, once again, a large scale forcing term,
and the other term is given by, D=
µ
∇−2
ω
+
ν
∇2
ω
. The rst term in this has been termed by
the authors as “physically articial hypoviscosity”.
These work on homogeneous, stationary turbu-
lence in 2D and 3D requiring “hypoviscosity” and/
or “hypervisocisity” is an indication that some
form of articiality is being introduced and one is
no longer truly solving the INSE. It is to be noted
[2] that the time-averaged, compensated spectral
density for (i) streamwise, (ii) wall-normal, and
(iii) spanwise velocity components plotted as a
function of streamwise wavenumber has demon-
strated the “bottleneck” eect noted in the ex-
periments [23] for the inhomogeneous ow over a
at plate. In these simulations [2], the INSE has
been solved in velocity-vorticity formulation by us-
ing high accuracy compact schemes that require
fourth diusion term for numerical stabilization.
There is another aspect of viscous action which
is obscured when one inspects periodic ow. It has
been demonstrated that the transport equation for
enstrophy (Ω1=
ω
i
ω
i) of inhomogeneous ow is
given by [30],
DΩ1
Dt =2
Re 1
2∇2Ω1−(∇
ω
)2+2
ω
i
ω
j
∂
ui
∂
xj
The rst set of terms on the right hand side
arise due to viscous diusion. Furthermore, it is
to be noted that the rst term (∇2Ω1) will drop
out for periodic problems. This has led to a misun-
derstanding about turbulent ows where diusion
term is negative denite for periodic problems, as
in homogeneous isotropic turbulence. The diu-
sion term is erroneously associated with dissipa-
tion, even when the ow is not periodic. The last
term on the right hand side of the above equation
arises due to vortex stretching, which can act as
a production term for enstrophy. However, this
term was identied [5] to additionally contribute
to self-attenuation in the presence of high non-
linearity during extreme events by creating length
scales smaller than the Kolmogorov scale. For Eu-
ler equation (in the limit of Re →∞or
ν
→0)
enstrophy was noted to grow unbounded in a -
nite time [5] citing an observation in Doering [9].
Interestingly, it was furthermore conjectured [5]
that such nite time blow-up “would correspond
to turbulent solutions of the INSE”. At the same
time, the authors state that “the question remains
open whether the non-linear amplication could
overcome viscous damping when the ow is su-
ciently turbulent”. During transition to turbulence
of inhomogeneous 2D ow past dierent bodies, it
has been shown [27] that after ows become fully
turbulent following transition, the velocity compo-
nents always achieve a limiting value, even when
vortex stretching is absent.
Having dwelt upon various physical aspects of
the ow, one notes that the role of numerical anal-
ysis of methods used to compute turbulent ows
has received less attention. One such analysis was
undertaken by solving the one-dimensional (1D)
convection equation, given as [25],
∂
u
∂
t+c
∂
u
∂
x=0 (4)
In this reference, the Fourier spectral-RK2 method
has been shown to be unconditionally unstable.
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
10
Volume 10, 2022
Here, we will elaborate on this aspect of oft-
implemented numerical methodology for comput-
ing turbulence by solving the convection-diusion
equation and associated properties. This equation
is more closely tied up with the INSE used to solve
turbulent ows, due to the added presence of the
diusion term.
The unknown in GSA is written in hybrid spec-
tral form as u(x,t) = Rˆ
U(k,t)eikxdk, where kis the
wavenumber in x-direction. If one starts solving
Eq. (4)subjected to the initial condition given by
u(xj,t=0) = u0
j=RU0(k)eikx jdk, then the numer-
ical solution after nth time-step will be obtained
as,
un
N,j=ZU0(k) [|Gj|]nei(kx j−n
φ
j)dk (5)
where Gj=ˆ
U(k,tn+∆t)
ˆ
U(k,tn)is the amplication fac-
tor and is in general, a complex quantity i.e.
Gj=Gr j +iGi j. Its modulus is given by |Gj|=
(G2
r j +G2
i j)1/2and the phase shift per time-step is
calculated from, tan
φ
j=−Gi j/Gr j. GSA has been
used to develop the error (e) propagation equation
for the convection equation given by [29],
∂
e
∂
t+c
∂
e
∂
x=−[1−cN
c]c
∂
¯uN
∂
x−ZVgN −cN
kZik0U0[|G|]neik0(x−cNt)dk0dk
−ZLn |G|
∆tU0[|G|]neik(x−cNt)dk
(6)
Here, the numerical error is given by, e(x,t) =
u(x,t)−¯uN(x,t), from which one can derive the
governing equation for e(x,t)[29]. The numerical
phase speed (cN) is obtained from
φ
j(phase shift
per time step), so that n
φ
j=kcNn∆t. The physi-
cal phase speed is cfor all wavenumbers, but cNis
noted to depend on k. Thus, the numerical solu-
tion is dispersive (in contrast to the non-dispersive
nature of 1D convection equation), and is rewrit-
ten as,
¯uN=ZU0(k) [|G|]t/∆teik(x−cNt)dk (7)
Having obtained the correct numerical phase
speed, the numerical group velocity at the jth-
node is expressed as, VgN
cj
=1
hNc
d
φ
j
dk . Equation
(6)clearly demonstrates that the accuracy of the
numerical solution is governed by the terms on
the right hand side. From this, we identify |G|,
cN/cand VgN /cas the main numerical parameters
for the convection equation, to determine the error
forcing. Such forcing arises due to nite discretiza-
tions, which dominates the round-o error as the
other contributor. Previously, looking at the lin-
ear nature of the governing Eq. (4), researchers
assumed that the error dynamics is also dictated
by this equation with the right hand side of Eq.
(6)going to zero.
It is important to understand the roles played
by these numerical parameters in ensuring accu-
racy of scientic computing. The concept of error
propagation was introduced in proper perspective
for convection equation [29] and the same has been
subsequently shown for the diusion equation [28]
and convection-diusion equation [33]. GSA has
also been used to demonstrate a linear dispersive
mechanism for numerical error growth [35]. Re-
cently, GSA has also been used to understand the
eect of dispersion and dissipation on both convec-
tion and diusion terms in the 2D linearized com-
pressible Navier-Stokes equation (LCNSE) for hy-
brid nite dierence/Fourier spectral scheme [35].
A major development was achieved [34], where
it was shown for the rst time that the convection-
diusion equation truly reects the dependence on
numerical parameters as it does for INSE. This
provides the framework to calibrate any numeri-
cal method for the solution of INSE (both homoge-
neous and inhomogeneous ows, with or without
periodicity) by checking the numerical parameters
to be used in the solution process. The process
involves calibration with the convection-diusion
equation which is given by,
∂
u
∂
t+c
∂
u
∂
x=
ν∂
2u
∂
x2(8)
As noted above, in numerical simulations
the physical phase speed (c) is altered to the
wavenumber-dependent numerical phase speed
(cN(k)) for convection equation. In Eq. (8), the
coecient of diusion (which is identical to the
kinematic viscosity in INSE) also changes from
constant
ν
to a wavenumber-dependent numeri-
cal coecient:
ν
N(k). It is customary [24] that,
the numerical parameters given by the Courant-
Friedrich-Lewy (CFL) number (Nc=c∆t/∆x) and
the Peclet number (Pe =
ν
∆t/(∆x)2) are impor-
tant in xing grid spacing (∆x) and time step (∆t)
for solving problems dominated by convection and
diusion.
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
11
Volume 10, 2022
2 Illustrative Simulations of
Convection- and
Convection-Di.usion Equation
First, we discretize the governing equation using
multi-stage Runge-Kutta time integration [24] and
Fourier spectral method for discretizing the con-
vection and diusion terms as,
∂
u
∂
x=Rik ˆ
Ueikxdk
and
∂
2u
∂
x2=R−k2ˆ
Ueikxdk. In the computations,
these derivatives are evaluated using fast Fourier
transform (FFT). From these discrete equations,
and using the denition of numerical factors, one
obtains the complex parameter Gas function of
k∆xand Ncfor the convection equation [29]; k∆x,
Ncand Pe for convection-diusion equation [33].
Having obtained |G|and
φ
, one can readily ob-
tain cN/cand Vg,N/cfor the convection equation.
For the convection-diusion equation, one would
also nd out
ν
N/
ν
, the non-dimensional dispersive
coecient of diusion [33, 34].
First, we demonstrate solution of convection
equation obtained by solving Eq. (6)using Fourier
spectral method for spatial discretization and two-
, three-, and four-stage Runge-Kutta method for
time integration. If one chooses a grid spacing
of ∆x, then the resolved maximum wavenumber
(kmax) is xed by the Nyquist criterion, kmax∆x=
π
.
The convection equation is solved in a domain,
0≤x≤10 discretized uniformly with 4096 points
for the propagation of a wave-packet. The initial
condition for the wave packet is given by,
u(x,0) = e−10(x−5)2sin(k0∆x)(9)
The wave-packet is characterized by the cen-
tral wavenumber given by, k0∆x=0.2442, a value
which is signicantly lower than the Nyquist limit.
This low value ensures that the error sources
due to phase error and dispersion error are sub-
dominant as compared to the error created by the
|G|term.
In Fig. 1, the propagation of the wave-packet
is shown after 30000 ∆t, where ∆t=4.884 ×10−4
is xed for Nc=0.1and c=0.5, with RK2 method
result shown in the top frame; the three-stage RK
(RK3) method result in the middle frame and four-
stage RK (RK4) method result depicted in the
bottom frame. As the exact solution is simply
translated to the right at the speed c, one can com-
pare it with the numerical solution. It is evident
that RK2-Fourier spectral method results are er-
roneous. As one is solving a periodic problem, the
round-o error magnies, as determined by |G|,
and that keeps convecting within the domain accu-
mulating in magnitude. However, no distinguish-
able errors are noted for the RK3-Fourier spectral
and RK4-Fourier spectral methods in Fig. 1for
this low value of CFL number, Nc=0.1.
In Fig. 2, the propagation of the same wave-
packet is shown after a time interval of 30000 ∆t,
for the same value of Nc=0.1with the three time
discretization methods used in Fig. 1, with Fourier
spectral method for spatial discretization, by solv-
ing the convection-diusion equation. Here, the
same time step is used, and the coecient of dif-
fusion is so chosen that the Peclet number is 0.01.
Unlike the non-dissipative, non-dispersive convec-
tion equation case in Fig. 1for which the physical
amplication rate is given by |Gphys|=1, here for
the convection-diusion equation case the physi-
cal amplication would indicate attenuation due
to diusion, such that |GPhys|<1. In this case,
the numerical amplication rate (|GNum|) will be
such that |GNum|/|GPhys|is equal to one for ac-
curacy, and |GNum|>1would indicate numerical
instability. In Fig. 2, the solutions are shown for
Nc=0.1and Pe =0.01 for the three Runge-Kutta
methods, with the same wave-packet starting from
the same initial locations. All the numerical solu-
tions coincide with the exact solution, implying
that the unstable RK2-Fourier spectral method
for convection equation has become stable for the
convection-diusion equation for the same value
of Nc=0.1. It has already been shown unambigu-
ously that the numerical behavior of convection-
diusion equation reects identically the behav-
ior of the numerical solution of INSE when the
numerical parameters are kept identical [34]. For
the case of RK3 and RK4 methods, the solutions
are error free for these parameter combinations
for both the convection- and convection-diusion
equations.
The conditional stability of the RK2-Fourier
spectral method for the solution of convection-
diusion equation, explains its use for the sim-
ulation of INSE, despite its failure to do so for
the convection equation. These aspects are fur-
ther explored by solving convection-diusion equa-
tion for a higher CFL number of Nc=0.5, while
keeping the Peclet number at the same level of
Pe =0.01. It is well known [18] that for DNS,
the diusion discretization places a stricter control
on the time step via the choice of Peclet number.
The solutions obtained by RK2, RK3 and RK4
methods are compared in Fig. 3for Nc=0.5and
Pe =0.01 and the simulations are run only for 101
time-steps. One notes erroneous solutions by the
RK2-Fourier spectral method, whereas RK3- and
RK4-Fourier spectral methods again reproduce er-
ror free propagation of the same wave-packet as
before. All of these results shown here indicate
the necessity to characterize the space-time dis-
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
12
Volume 10, 2022
x
u
0 2 4 6 8
-0.5
0
0.5
1
uex
uRK3
x
u
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5
-0.5
0
0.5
1
uex
uRK4
x
u
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5
-1
-0.5
0
0.5
1
uex
uRK2
Wave packet propagation (k0∆x=0.2442)
Nc=0.1, Pe=0, t=30,000∆t
Figure 1: Comparison of numerical solution with ex-
act solution at t=30,000∆tfor the convection of a
wavepacket using a 1D linear convection equation.
Here, Fourier-spectral method is adopted for spatial
derivative and RK2, RK3 and RK4 methods are em-
ployed for time integration. All simulations are per-
formed on a uniform grid with 4096 points and for a
CFL number Nc=0.1.
x
u
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5
-6E-08
-4E-08
-2E-08
0
2E-08
4E-08
6E-08
uex
uRK2
x
u
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5
-6E-08
-4E-08
-2E-08
0
2E-08
4E-08
6E-08
uex
uRK3
x
u
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5
-6E-08
-4E-08
-2E-08
0
2E-08
4E-08
6E-08
uex
uRK4
Wave packet propagation (k0∆x=0.2442)
Nc=0.1, Pe=0.01, t=30000∆t
Figure 2: Comparison of numerical solution with ex-
act solution at t=30,000∆tfor the convection of a
wavepacket using a 1D linear convection-diffusion
equation. Here, Fourier spectral method is adopted
for spatial derivative and RK2, RK3 and RK4 meth-
ods are employed for time integration. All simu-
lations are performed on a uniform grid with 4096
points, with a CFL number Nc=0.1 and Peclet num-
ber Pe =0.01.
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
13
Volume 10, 2022
x
u
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5
-1
-0.5
0
0.5
1
uex
uRK2
x
u
0 2 4 6 8
-0.5
0
0.5 uex
uRK3
x
u
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5
-0.5
0
0.5 uex
uRK4
Wave packet propagation (k0∆x=0.2442)
Nc=0.5, Pe=0.01, t=101∆t
Figure 3: Comparison of numerical solution with
exact solution at t=101∆tfor the convection of a
wavepacket using a 1D linear convection-diffusion
equation. Here, Fourier spectral method is adopted
for spatial derivative and RK2, RK3 and RK4 meth-
ods are employed for time integration. All simu-
lations are performed on a uniform grid with 4096
points, with a CFL number Nc=0.5 and Peclet num-
ber Pe =0.01.
cretization methods. Similar analysis based on the
convection equation has been reported in Fig. 1
in [25] for the RK2-Fourier spectral method. It is
shown that this method is unconditionally unsta-
ble for all CFL numbers for ows dominated by
convection. Only for very small values of k∆xone
observes neutral stability, which can explain why
researchers [1] have reported nite-time blow-up
for the solution of Euler equation. To understand
the implications of the reported results in Figs.
1to 3one must perform a thorough study of the
numerical properties of convection and convection-
diusion equations for the pseudo-spectral method
proposed in [22] using RK2 time integration. The
numerical evidences provided in Figs. 1to 3indi-
cate the shortcomings of the RK2-Fourier spectral
method. Next, we investigate the possibility of
overcoming these by using RK3 and/or RK4 time
integration methods.
3 Error Dynamics for Convection and
Convection-Di.usion Equations
It has been already noted in Eq. (6)that there are
three main sources driving the dynamics of error,
via the forcing by the phase error term, the disper-
sion error term and more importantly by numer-
ical instability and unphysical numerical stability
term, noted in Eq. (6)for the convection equation.
It is a common misconception that numerical sta-
bility is desirable [36], which makes many prac-
titioners complacent about this potential source
of error, whereas a numerically unstable method
is easily recognized due to solution blow-up. For
the present analysis of convection and convection-
diusion equation, numerical parameters are cho-
sen such that the error is solely due to numerical
instability and unphysical numerical stability. For
the convection equation, an error-free calculation
must have neutrally stable numerical amplication
factor, |GNum|=1. In contrast, for convection-
diusion equation the physical solution attenuates
such that |GPhys|<1and the numerical solution
should reect this, so that |GNum|/|GPhys|=1. The
procedures to calculate these quantities are given
for convection [29] and convection-diusion equa-
tions [33], and are provided here in Appendix A.
We use the GSA to map the property charts of
the method involving Fourier spectral method for
spatial discretization and RK2 method for time in-
tegration. We also report the properties at the in-
terior of the domain, as one is interested to analyze
the methods to solve 3D homogeneous isotropic
turbulence in a periodic box. In Fig. 4, RK2-
Fourier spectral method is considered for both the
canonical equations. For the one-dimensional con-
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
14
Volume 10, 2022
Nc
k∆x
0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3
Unstable
region
(|GNum|>1)
Pe=0.01
1D Convection-
Diffusion equation RK2-Spectral
|GNum|/|GPhys|-contours
Nc
k∆x
0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3
RK2-Spectral
|GNum|-contours
Unstable
region
(|GNum|>1)
1D Convection
equation
Nc
k∆x
0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3Unstable
region
(|GNum|>1)
1D Convection-
Diffusion equation RK2-Spectral
|GNum|/|GPhys|-contours Pe=0.1
Figure 4: Numerical amplification factor and ratio of
numerical amplification factor to the physical ampli-
fication factor for the RK2-Fourier spectral scheme
plotted in the (Nc,k∆x)-plane for linear convection
and convection-diffusion equations, respectively. For
the latter, contours are plotted for two representative
Peclet numbers- Pe =0.01 and 0.1. Dashed line rep-
resents the optimal Ncvalue.
vection equation, |G|-contours are shown plotted
in the (Nc,k∆x)-plane on the left side of Fig. 4.
It is noted that barring a small region very close
to the origin, the method is unconditionally un-
stable. This result was also reported [25], and one
notes the eects of numerical instability in Fig. 1,
in the top frame for which the choice of Nc=0.1
and k0∆x=0.2442 does not guarantee stability of
all the resolved scales by the Nyquist criterion.
The seeds of unstable signals are always present in
the round-o error and those wavenumber compo-
nents interact with the given wave-packet to create
multitudes of wave-packets noted for this periodic
problem. A non-periodic problem with dierent
inow and outow boundary conditions would al-
low some of the error-packets to leave the com-
putational domain. To emphasize the growth of
error that is inherent with the numerical method,
a periodic problem is considered so that the signal
is always inside the domain, along with the re-
circulating continuously evanescent error. For the
convection-diusion equation, the error dynamics
is additionally aected via a term involving
ν
N/
ν
and the error dynamics is given by,
et+cex−
ν
exx =Zkmax
−kmax
(
ν
N−
ν
)k2e−
ν
Nk2n∆tU0(k)eik(x−cNtn)dk
+ikcNe−
ν
Nk2n∆tU0(k)eik(x−cNtn)
kmax
−kmax
−Zkmax
−kmax VgN −cN
kZk
−kmax
ik0e−
ν
Nk02n∆tU0(k0)eik0(x−cNtn)dk0dk
−Zkmax
−kmax
ikc e−
ν
Nk2n∆tU0(k)eik(x−cNtn)dk
(10)
where
ν
N,cNand VgN denote the numerical diu-
sion coecient, numerical phase speed and numer-
ical group velocity, respectively. It is interesting
to note that the error evolution does not follow
the same dynamics as the governing equation. We
note that this behaviour was rst unequivocally
demonstrated for the linear convection equation
[29] thus dispelling the earlier held notion that the
error evolution always follows the same dynamics
as the governing linear equation. Further, it was
demonstrated in the same reference that the error
due to computing is solely due to three sources-
diusion/anti-diusion, phase and signal propa-
gation, with the errors from each source quanti-
ed by the error dynamics equation (given here by
Eq. (6)). This is a universal aspect of computing.
For the error dynamics equation of convection-
diusion equation, these three sources of error are
also noted. In this paper, we focus solely on the
error due to numerical amplication/diusion for
evaluating/analyzing the numerical errors for the
Fourier spectral method. Hence, mostly the con-
tours due to |GNum|/|GPhys|are shown for dierent
simulation cases in the following.
On the right hand side frames of Fig. 4,
we have shown the contours of |GNum|/|GPhys|in
(Nc,k∆x)-plane for Pe =0.01 and 0.1. One notices
that there is a nite range of Ncfor which there
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
15
Volume 10, 2022
Nc
k∆x
0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3Unstable
region
(|GNum|>1)
Pe=0.01
1D Convection-
Diffusion equation RK3-Spectral
|GNum|/|GPhys|-contours
0.9999
0.999
0.99
0.95
0.95
0.99
0.999999
Nc
k∆x
0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3
RK3-Spectral
|GNum|-contours
Unstable
region
(|GNum|>1)
1D Convection
equation
Nc
k∆x
0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3Unstable
region
(|GNum|>1)
1D Convection-
Diffusion equation RK3-Spectral
|GNum|/|GPhys|-contours Pe=0.1
Figure 5: Numerical amplification factor and ratio of
numerical amplification factor to the physical ampli-
fication factor for the RK3-Fourier spectral scheme
plotted in the (Nc,k∆x)-plane for linear convection
and convection-diffusion equations, respectively. For
the latter, contours are plotted for two representative
Peclet numbers- Pe =0.01 and 0.1. Dashed line rep-
resents the optimal Ncvalue for LES.
is no instability over the full resolved scale. The
instability is indicated whenever |GNum|>1. This
is true for both the Peclet numbers, but the con-
tour values in the stable region varies with k∆x.
Presented results here are for Nc=0.1and 0.5, for
which the quotient is almost equal to 1for only a
small range of k∆xwithin the Nyquist limit. To
nd such a value of Ncfor a chosen Pe with max-
imum range of wavenumbers within the Nyquist
limit, one requires to perform GSA.
In Figs. 5and 6,|GNum|/|GPhys|-contours are
plotted in (Nc,k∆x)-plane for Pe =0.0, 0.01 and 0.1
for RK3- and RK4-Fourier spectral methods. Un-
like the RK2-Fourier spectral method for the con-
vection equation, RK3 and RK4 time integration
methods do not exhibit unconditional numerical
instability, as shown in the left panel. However
to get neutral stability for this canonical equa-
tion, one will be forced to take small value of
Nc, with RK4 method signicantly better than
the RK3 time integration method allowing much
higher values of Nc. Such values of CFL number
will not display nite time blow-up of solution for
Euler equation. The right frames in Figs. 5and
6,|GNum|/|GPhys|-contours are plotted in (Nc,k∆x)-
plane for Pe =0.01 and 0.1 for simulating INSE.
Nc
k∆x
0.2 0.4 0.6 0.8 1 1.2 1.4
0
0.5
1
1.5
2
2.5
3Unstable
region
(|GNum|>1)
Pe=0.01
1D Convection-
Diffusion equation RK4-Spectral
|GNum|/|GPhys|-contours
0.999999
0.55
0.95
0.9999
0.999
0.99
0.99
0.75
0.91
0.75
Nc
k∆x
0.2 0.4 0.6 0.8 1 1.2 1.4
0
0.5
1
1.5
2
2.5
3
RK4-Spectral
|GNum|-contours
Unstable
region
(|GNum|>1)
1D Convection
equation
Nc
k∆x
0.2 0.4 0.6 0.8 1 1.2 1.4
0
0.5
1
1.5
2
2.5
3Unstable
region
(|GNum|>1)
1D Convection-
Diffusion equation RK4-Spectral
|GNum|/|GPhys|-contours Pe=0.1
Figure 6: Numerical amplification factor and ratio of
numerical amplification factor to the physical ampli-
fication factor for the RK4-Fourier spectral scheme
plotted in the (Nc,k∆x)-plane for linear convection
and convection-diffusion equations, respectively. For
the latter, contours are plotted for two representative
Peclet numbers- Pe =0.01 and 0.1. Dashed line rep-
resents the optimal Ncvalue for LES.
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
16
Volume 10, 2022
0.0 0.5 1.0 1.5 2.0 2.5 3.0
k
x
0
10 12
10 11
10 10
10 9
10 8
10 7
10 6
10 5
10 4
10 3
10 2
10 1
|1-|
GNum
|/|
GPhys
||
Solution of 1D CDE using
RK
2-Spectral method
Stable parameters
Nc=0.1,Pe=0.01
Nc=0.2,Pe=0.01
Nc=0.3,Pe=0.01
Figure 7: Determination of optimal Ncfor accurate
solution of 1D convection-diffusion equation using
RK2-Fourier spectral method for Pe =0.01. Accu-
racy is evaluated with respect to error in numerical
amplification factor (|1−(|Gnum|/|Gphys|)|) for all re-
solvable wavenumbers (k∆x).
For RK3 time integration method, performing ac-
curate simulation would require taking vanishingly
small values of Nc, for vanishingly small values of
Pe. It must however, be remembered that for no
parameter combinations, one will reproduce the
exact solution for the fully resolved Nyquist limit
of wavenumbers (k∆x=
π
). For the RK4 time in-
tegration method, the permissible values of Nc,
Pe for performing near-DNS simulation will be
relatively higher than that will be obtained by
the RK3 time integration method. Such imper-
fectly resolved computations can be viewed more
as an implicit large eddy simulation. Overall, com-
paring the |GNum|/|GPhys|-contours in Figs. 4to
6, it appears that RK4-Fourier spectral method
will provide the best results for Pe =0.01 and
Nc=0.243793, as described next.
In Figs. 4to 6, one can note the contri-
bution to error dynamics by the term involving
|GNum|/|GPhys|for specic values of Pe. On a cur-
sory observation, one notes that for a specic value
of Nc,|GNum|/|GPhys|value keeps changing with
k∆x. To quantify such variations and to obtain an
optimal value of Ncfor which the error will be least
due to this |GNum|/|GPhys|-term in Fig. 7, we have
plotted ||GNum|/|GPhys|−1|as a function of k∆xfor
the RK2-spectral method showing that the opti-
mal error forcing is obtained for Nc=0.2, which
is shown to be lower than that for Nc=0.1and
0.3. Having noted this optimal value of Nc=0.2,
one also notices that the forcing by this term in-
creases monotonically in nonlinear manner with
increase in wavenumber up to the Nyquist limit.
Such monotonic growth in error with wavenumber
is clearly apparent from Fig. 4, for these lower
0.0 0.5 1.0 1.5 2.0 2.5 3.0
k
x
0
10 12
10 11
10 10
10 9
10 8
10 7
10 6
10 5
10 4
10 3
|1-|
GNum
|/|
GPhys
||
1D CDE
Stable parameters
RK
2-Spectral,Nc=0.2,Pe=0.01
RK
3-Spectral,Nc=0.05,Pe=0.01
RK
4-Spectral,Nc=0.243793,Pe=0.01
Figure 8: Comparison of the accuracy of RK2, RK3
and RK4 time integration schemes in solving the
1D convection-diffusion equation for Pe =0.01 us-
ing Fourier spectral method. Accuracy is evaluated
with respect to error in numerical amplification factor
(|1−(|Gnum|/|Gphys|)|) for all resolvable wavenum-
bers (k∆x). Presented values of Nccorrespond to
the respective optimal values of the time integration
schemes.
values of CFL, which is also restricted to the right
due to numerical instability for the choice of Nc.
Even though the numerical instability region re-
duces with increasing Pe, the higher error with in-
crease in wavenumber is drastic.
Having noted the existence of optimal Ncvalue
for RK2-Fourier spectral method as equal to 0.2,
next we locate the same for RK3-Fourier spectral
and RK4-Fourier spectral methods. The variation
of contributions by the optimal Ncchoices to error
dynamics is compared in Fig. 8. In Fig. 5, the
variation of |GNum|/|GPhys|-contours in (Nc,k∆x)-
plane indicates that the error dynamics term de-
pendent on |GNum|/|GPhys|is more complex and
contribute higher values. For the RK3-Fourier
spectral method such an optimal value is noted
for Nc=0.05 and for RK4-Fourier spectral method
the optimal error is obtained for Nc=0.243793
for the choice of Pe =0.01 with the results de-
picted in Fig. 8. The quantitative variation of
error forcing for RK3-spectral method in this g-
ure was shown qualitatively in Fig. 5. But, the
RK4-Fourier spectral method shows monotonic
variation of |GNum|/|GPhys|with wavenumber, that
causes the least error among these three time in-
tegration methods.
4 Evaluating Schemes for High
Reynolds Number Computations
In this section, we analyze the performance of
RK2, RK3 and RK4 time integration schemes with
the Fourier spectral spatial discretization for high
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
17
Volume 10, 2022
Reynolds number ow computations. Such com-
putations imply very small values of Pe, as
ν
is
equal to 1
Re . Hence, GSA for low values of Pe
is warranted as it will reveal the performance in
terms of accuracy of the schemes, and in turn can
be used to pinpoint the causality between the error
and the chosen numerical parameters. Further-
more, this analysis will provide information to all
to make informed choices for the parameters to
attain a level of desired accuracy.
In Figs. 9and 10, contours of |GNum|/|GPhys|
are plotted for the RK2- and RK4-Fourier spec-
tral methods, respectively. In these gures, the
left frame corresponds to the properties for the
1D convection equation, whereas the right frames
are for the convection-diusion equation with low
values of Pe =0.001 and 0.0001. As noted ear-
lier, the RK2-Fourier spectral method is uncon-
ditionally unstable for the 1D convection equa-
tion, whereas it is conditionally stable for the
convection-diusion equation. From Fig. 9, one
notes the stable region of the RK2-Fourier spec-
tral scheme to decrease with decreasing Pe. For
Pe =0.001, numerical instability is noted for Nc≥
0.175, which is the critical limit i.e. Nccr =0.175.
This reduces to Nccr =0.096 for a reduced Pe =
0.0001. For a higher Pe =0.01, one notes Nccr =
0.33. Thus, Nccr decreases with decrease in Pe. For
the RK4-Fourier spectral method, Nccr =0.905,
0.905,0.92 for Pe =0.0001,0.001 and 0.01, re-
spectively. This, once again, shows that the region
of stability decreases with decrease in Pe. This is
a real characteristic of any numerical scheme for
which adding diusion (either numerical or physi-
cal) stabilizes it, as long as Pe value remains below
a critical value. It is noted that while the Nccr re-
duces to very low values for the RK2 scheme, a
marginal reduction is noted for the RK4 scheme.
This can be explained by the lowest bound on
Nccr =0being set for the convection equation
(Pe =0), for which the RK2 scheme is uncondi-
tionally unstable. But the RK4 scheme does not
suer from unconditional instability and the re-
duction in Nccr is minimal. For Pe 6=0, a stabilizing
eect overcomes the unconditional instability of
the RK2 scheme with Nccr increasing. Thus, RK2-
Fourier spectral method is unsuitable for comput-
ing high Reynolds number ows.
As noted earlier from the analysis of error evo-
lution and its propagation given by Eqn. (10)
for the convection-diusion equation, all numer-
ical schemes must satisfy zero diusion error
i.e. |GNum|/|GPhys|is equal to 1for the resolved
wavenumbers. Thus, this is a necessary condi-
tion to be satised by every scheme seeking to
obtain accurate results and therefore, qualify for
Nc
k∆x
0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3
RK2-Spectral
|GNum|-contours
Unstable
region
(|GNum|>1)
1D Convection
equation
Nc
k∆x
0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3
Unstable
region
(|GNum|>1)
Pe=0.0001
1D Convection-
Diffusion equation RK2-Spectral
|GNum|/|GPhys|-contours
Nc=0.12
k∆x=1.99
Nc
k∆x
0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3
Unstable
region
(|GNum|>1)
1D Convection-
Diffusion equation RK2-Spectral
|GNum|/|GPhys|-contours
Pe=0.001
Figure 9: Numerical amplification factor and ratio of
numerical amplification factor to the physical ampli-
fication factor for the RK2-Fourier spectral scheme
plotted in the (Nc,k∆x)-plane for linear convection
and convection-diffusion equations, respectively. For
the latter, contours are plotted for two representative
Peclet numbers- Pe =0.0001 and 0.001. Dashed line
represents the conditions at which numerical simula-
tions are performed. Blue coloured region indicates a
zone with zero diffusion error.
performing strict DNS. In Figs. 9and 10, the re-
gion with zero diusion error is shown marked by
blue colour. It is interesting to note that while
RK4 scheme is superior to RK2 due to the larger
extent of zero diusion error, strict DNS is im-
possible as one does not obtain zero error at all
resolved wavenumbers.
Next, the property charts obtained earlier are
corroborated by computing the propagation of
the wave-packet using RK2-, RK3- and RK4-
schemes with the Fourier pseudo-spectral method.
The computations are performed with uniformly
spaced grids having 4096 points in a domain x∼
[0,10], considered in section 2. The packet is cen-
tered around k∆x=0.22 and time step is chosen
as 5.86 ×10−4such that Nc=0.12 with c=0.5.
The diusion coecient is determined from Pe =
0.0001. Figure 11 shows the solutions plotted af-
ter 30,000∆tfor RK2-, RK3- and RK4-Fourier
pseudo-spectral methods, respectively. While
this parameter combination is noted to be sta-
ble for the RK3-and the RK4-methods, the RK2-
method suers a numerical instability due to anti-
diusion, as indicated by a dashed vertical line
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
18
Volume 10, 2022
0
0.9999
0.99
0.9
0.75
0.75
0.55
0.55
0.999999
Nc
k∆x
0.2 0.4 0.6 0.8 1 1.2 1.4
0
0.5
1
1.5
2
2.5
3
RK4-Spectral
|GNum|-contours
Unstable
region
(|GNum|>1)
1D Convection
equation
Nc
k∆x
0.2 0.4 0.6 0.8 1 1.2 1.4
0
0.5
1
1.5
2
2.5
3Unstable
region
(|GNum|>1)
1D Convection-
Diffusion equation RK4-Spectral
|GNum|/|GPhys|-contours
Pe=0.001
Nc
k∆x
0.2 0.4 0.6 0.8 1 1.2 1.4
0
0.5
1
1.5
2
2.5
3Unstable
region
(|GNum|>1)
Pe=0.0001
1D Convection-
Diffusion equation RK4-Spectral
|GNum|/|GPhys|-contours
Figure 10: Numerical amplification factor and ratio
of numerical amplification factor to the physical am-
plification factor for the RK4-Fourier spectral scheme
plotted in the (Nc,k∆x)-plane for linear convection
and convection-diffusion equations, respectively. For
the latter, contours are plotted for two representa-
tive Peclet numbers- Pe =0.0001 and 0.001. Blue
coloured region indicates a zone with zero diffusion
error.
x
u
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5
-200
-150
-100
-50
0
50
100
150
200
uex
uRK2
x
u
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5
-1
-0.5
0
0.5
1
uex
uRK3
x
u
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5
-1
-0.5
0
0.5
1
uex
uRK4
Wave packet propagation (k0∆x=0.22)
Nc=0.12, Pe=0.0001, t=30000∆t
Figure 11: Comparison of numerical solution with
exact solution for a time t=30,000∆tfor the con-
vection of a wavepacket using a 1D linear convection-
diffusion equation. Here, Fourier spectral method is
adopted for spatial derivative and RK2, RK3 and RK4
methods are employed for time integration. All sim-
ulations are performed on a uniform grid with 4096
points, with a CFL number Nc=0.12 and Peclet
number Pe =0.0001.
(blue) in Fig. 9. From this we note that instability
occurs for k∆x≥1.99. The computation validates
the observations from the analysis. Results for the
RK3- and RK4-methods show that the numeri-
cal solution follows the exact solution accurately,
whereas the solution obtained by RK2-method is
noted with unphysical wave-packets and high am-
plitudes attributed to anti-diusion.
One notes that although the wave-packet is
centered at k∆x=0.22, the computation with RK2
time integration shows instabilities. The explana-
tion for this behavior is noted from the analyti-
cal structure of the wavepacket which results in a
spectrum in the Fourier space containing contri-
butions from all wavenumbers upto the Nyquist
limit. These high wavenumber amplitudes, de-
spite having nearly negligible amplitudes initially,
are amplied by |Gnum|every time-step to cause
exponential growth leading to numerical instabil-
ity and eventual blow-up (not shown). We also
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
19
Volume 10, 2022
0.0 0.5 1.0 1.5 2.0 2.5 3.0
Wavenumber (k
x
)
10 19
10 16
10 13
10 10
10 7
10 4
10 1
102
|FFT(u)|
Nc
= 0.12,
Pe
= 0.0001
RK
2 Solution (
t
= 30, 000
t
)
RK
3 Solution (
t
= 30, 000
t
)
RK
4 Solution (
t
= 30, 000
t
)
Exact Solution (
t
= 30, 000
t
)
Figure 12: Comparison of FFT of the exact so-
lution with solution obtained from RK2, RK3 and
RK4 time integration schemes in solving the 1D lin-
ear convection-diffusion equation for Peclet number
0.0001 and Nc=0.12 using Fourier spectral method.
clear a misconception about the necessity of the
presence of relevant wavenumbers for the insta-
bility to occur. While it is intuitive to assume
that instabilities will only occur if the relevant
wavenumbers are present, practical computations
suer from random/round-o errors which excite
all scales and thus, the seeds for instability are cre-
ated. Therefore, it is not necessary for the relevant
wavenumbers to be present in the initial solution
for the instability to occur. Round-o errors in-
herent in any computing will trigger the unstable
scales, leading to eventual instabilities.
In Fig. 12, the FFT of the numerical solution at
t=30,000∆tis plotted for the numerical solutions
obtained with RK2-, RK3- and RK4-Fourier spec-
tral methods. In the same gure, the spectrum
of the exact solution is also shown for reference.
The scale selection for instability for RK2-Fourier
spectral method is evident by noting the increas-
ing amplitudes of the high wavenumber compo-
nents which corroborates very well with the verti-
cal dashed blue line in Fig. 9. It is interesting to
note for the RK2-Fourier spectral method that the
amplitudes around the central wavenumber of the
wave packet are accurately represented with re-
spect to the reference analytical solution whereas
amplitudes for k∆x≥2.2show a visible departure.
The solutions due to RK3- and RK4-Fourier spec-
tral methods show a very good comparison with
the reference solution. This behaviour is consis-
tent with GSA.
These results clearly demonstrate the utility of
GSA which enables an accurate characterization of
the behavior of numerical methods in solving prob-
lems involving convection and convection-diusion
processes. Furthermore, these analyses have a di-
rect bearing on the accuracy of Navier-Stokes sim-
ulations due to their one-to-one correspondence
with linear convection-diusion equation [34].
5 Summary and Conclusions
A detailed analysis and calibration is performed
using global spectral analysis (GSA) for the nu-
merical methods employing Fourier spectral and
multistage Runge-Kutta (RK) methods. In par-
ticular the RK2 time integration method is stud-
ied to assess the schemes for high accuracy simula-
tions, as it has been proposed as DNS over the last
ve decades [19, 12, 22, 5]. These pseudo-spectral
schemes are analyzed here for the linear convec-
tion and linear convection-diusion equations as
these serve as appropriate model equations for as-
sessing and calibrating numerical errors for Eu-
ler and Navier-Stokes equation governed ows, re-
spectively. In particular, the analysis based on
the latter is important due its one-to-one corre-
spondence with the incompressible Navier-Stokes
equation [34]. With the help of GSA, the schemes’
diusion, dispersion and signal propagation prop-
erties are accurately predicted. To our knowledge
such results have not been reported before.
Computations and numerical analysis pre-
sented here conrm the unconditional numerical
instability of RK2-Fourier spectral method in solv-
ing the 1D convection equation, as reported also
by [25]. Hence, this method is unsuitable for sim-
ulating inviscid problems such as the Euler equa-
tion which is often noted in the literature [1]. The
Fourier spectral- RK2 method displays stability
for the solution of convection-diusion equation
for Peclet value 0.01, shown in Fig. 2and anal-
ysed with results in Fig. 4. In comparison, RK3
or RK4 time integration schemes used with the
Fourier spectral method shows numerical stability
for both the convection and convection-diusion
equations, demonstrated by Figs. 1-3,5and 6.
Thus, choosing a higher order integration scheme
instead of RK2, leads to the elimination of numer-
ical instability.
Error dynamics analysis performed in section 3
help answer the question of which among the RK2,
RK3 and RK4 methods oer best accuracy for the
Fourier spectral method in solving the convection-
diusion equation. Although Eqn. (10)identies
all sources of numerical error, here only the er-
ror contributed by the schemes’ inability to repre-
sent diusion accurately is used as the necessary
condition for high accuracy computations. This
contribution is assessed by noting the dierence
between numerical and physical diusion, which
is also the dierence between numerical ampli-
cation and physical amplication factors, repre-
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
20
Volume 10, 2022
sented by |1−|GNum|/|GPhys||. From such analysis,
optimal Ncvalues are determined for minimum er-
ror in diusion. Comparing RK2, RK3 and RK4
time integration methods for a xed Peclet num-
ber for optimal Nc, the RK4 method shows supe-
rior accuracy in Fig. 8. For RK2-Fourier spec-
tral method the optimal Nc=0.2is for Pe =0.01,
whereas for the RK3- and RK4-Fourier spectral
methods optimal Ncvalues are 0.05 and 0.243793,
respectively. This conclusively demonstrates that
RK4 time integration oers the best accuracy for
Fourier spectral method among these time inte-
gration schemes.
More importantly, the analysis reveals that
strict DNS is impossible as all the schemes do not
display vanishing |1− |GNum|/|GPhys|| for the re-
solved wavenumbers for any choice of Nc. Thus,
the optimal parameter noted above are for under-
resolved DNS only. The best Fourier spectral-
RK4 method in Fig. 10, shows exact vanishing of
|1− |GNum|/|GPhys|| for a range of Ncas given by
the blue coloured region with the maximum limit
on accurately resolved wavenumbers (k∆x=0.8for
Pe=0.0001 and k∆x=0.6for Pe=0.001) noted at
lower Ncvalues.
GSA performed for the RK2-Fourier spectral
method for Peclet number Pe =0.01 in Fig. 13
shows the diusion, dispersion and signal propa-
gation errors for the method. The latter two quan-
tities are the sucient conditions used for evalut-
ing the accuracy of the method. For the optimum
Ncvalue determined using the necessary condition
for the RK2 method, a maximum of 5% and 15%
errors are noted for the phase speed and signal
propagation, respectively. Thus, one notes that
the sucient condition too is not satised at all
resolvable wavenumbers.
The reduced stability of RK2-Fourier spectral
method used for solving Navier-Stokes equations
at high Reynolds numbers by reducing the time
step (i.e. by decreasing Peclet number) is a
counter-intuitive result reported here for the rst
time. As Peclet number is decreased, the stability
region decreases rapidly as demonstrated in Fig.
9, rendering the method unsuitable unless other
non-physical xes are employed to suppress nu-
merical instability. The error due to dispersion
shown in Fig. 13 negates the use of RK2-Fourier
spectral method completely for DNS. The analy-
ses shows that higher order RK schemes can alle-
viate the instability problem of the RK2-Fourier
spectral method, but even then the higher order
time integration schemes are strictly incapable of
performing DNS.
Further evidence of the agreement between
GSA and numerical simulations is shown in Figs.
0.999999
1.00005
1
0.999
0.9975
1.005
1.05
Nc
k∆x
0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3
|GNum|/|GPhys|
Unstable region
(|GNum|>1)
1.00025
1.005
1.05
1.125
1.15
1.125
1.05
1.005
0.9
1
Nc
k∆x
0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3
cNum/c
0.99
1.0005
1.005
1.025
1
0.5
0.25
0
Nc
k∆x
0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3
νNum/ν
Unstable region
(|GNum|>1)
1
1.005
1.05
1.15
1.25
1.15
1.05
1
0.75
0.5
Nc
k∆x
0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3
Vg,Num/Vg,Phys
RK2-Fourier spectral method, Pe=0.01
Figure 13: GSA analysis of RK2-Fourier spectral
method for the 1D convection-diffusion equation
plotted in the (Nc,k∆x)-plane for Pe =0.01. Plotted
are the quantities of ratios of numerical amplifica-
tion to physical amplification factor (|GNum|/|GPhys|),
numerical diffusion to physical diffusion coefficient
(
ν
Num/
ν
), numerical phase speed to physical phase
speed (cNum/c) and numerical group velocity to phys-
ical group velocity (Vg,Num/Vg,Phys).
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
21
Volume 10, 2022
11 and 12, to detect instability and its scale se-
lection. This is demonstrated with an initial so-
lution containing unstable wavenumbers for early
time explosion of numerical solution. It is not nec-
essary to use such initial condition for instability,
as the round-o error inherent in computing, will
have the seed of the unstable wavenumber, capa-
ble of triggering the eventual instability in long
time solution.
The present work shows unequivocally that
RK3, RK4 time integration schemes are manda-
tory in using Fourier spectral method for spa-
tial discretization in order to ensure the necessary
condition of numerical stability for high Reynolds
number ow computations. For required accuracy
in DNS, one must thereafter investigate for the
sucient condition of ensuring negligible disper-
sion error, as given in Eqn. (10). In future, de-
tailed analysis considering errors due to phase and
signal propagation will be performed in order to
ne tune the simulation parameters for both the
necessary and sucient conditions.
Declaration of Interests
The authors report no conict of interest.
References:
[1] Beale J. T., Kato T., & Majda A., Remarks
on the breakdown of smooth solutions for the
3-D Euler equations, Commun. Math. Phys.,
94, 61-66 (1984)
[2] Bhaumik S., & Sengupta T. K., Precursor of
transition to turbulence: Spatiotemporal wave
front, Phys. Rev. E, 89(4), 043018 (2014)
[3] Bracco A., & McWilliams J. C., Reynolds-
number dependency in homogeneous, sta-
tionary two-dimensional turbulence, J. Fluid
Mech., 646, 517-526 (2010)
[4] Buaria D., Bodenschatz E., & Pumir A., Vor-
tex stretching and enstrophy production in
high Reynolds number turbulence, Phys. Rev.
Fluids, 5(10), 104602 (2020)
[5] Buaria D., Pumir A. & Bodenschatz E., Self-
attenuation of extreme events in Navier–Stokes
turbulence, Nat. Commun., 11, 5852 (2020)
[6] Buaria D., & Sreenivasan K. R., Dissipation
range of the energy spectrum in high Reynolds
number turbulence, Phys. Rev. Fluids, 5(9),
092601 (2020)
[7] Choi Y., Kim B., & Lee C., Alignment of ve-
locity and vorticity and the intermittent distri-
bution of helicity in isotropic turbulence, Phys.
Rev. E. 80(1), 017301 (2009)
[8] David Cl., Sagaut P., & Sengupta T., A lin-
ear dispersive mechanism for numerical error
growth: spurious caustics, Euro. J. Mech. -
B/Fluids, 28(1), 146-151 (2009)
[9] Doering C. R., The 3D Navier-Stokes problem,
Ann. Rev. Fluid Mech., 41, 109-128 (2009)
[10] Donzis D. A., & Yeung P. K., Resolution ef-
fects and scaling in numerical simulations of
passive scalar mixing in turbulence, Physica
D: Nonlinear Phenomena, 239(14), 1278-1287
(2010)
[11] Donzis D. A., Yeung P. K., & Sreenivasan
K. R., Dissipation and enstrophy in isotropic
turbulence: Resolution eects and scaling in
direct numerical simulations, Phys. Fluids.
20(4), 045108 (2008)
[12] Eswaran V., & Pope S. B., An examination of
forcing in direct numerical simulations of tur-
bulence, Computers & Fluids, 16(3), 257-278
(1988)
[13] Frisch U., Turbulence: The Legacy of A.
N. Kolmogorov, Cambridge Univ. Press, UK
(1995)
[14] Ishihara T., Gotoh T., & Kaneda Y., Study of
high–Reynolds number isotropic turbulence by
direct numerical simulation, Ann. Rev. Fluid
Mech., 41, 165-180 (2009)
[15] Ishihara T., Morishita K., Yokokawa M.,
Uno A. & Kaneda Y., Energy spectrum in
high-resolution direct numerical simulations of
turbulence, Phys. Rev. Fluids, 1(8), 082403
(2016)
[16] Kaneda Y., Ishihara T., Yokokawa M.,
Itakura K., & Uno A., Energy dissipation rate
and energy spectrum in high resolution direct
numerical simulations of turbulence in a peri-
odic box, Phys. Fluids, 15(2), L21-L24 (2003)
[17] Lamorgese A. G., Caughey D. A., & Pope
S. B., Direct numerical simulation of homo-
geneous turbulence with hyperviscosity, Phys.
Fluids, 17(1), 015106 (2005)
[18] Moin P., & Mahesh K., Direct numerical sim-
ulation: A tool in turbulence research, Ann.
Rev. Fluid Mech., 30, 539-578 (1998)
[19] Orszag S., & Patterson G., Numerical
Simulation of three-dimensional homogeneous
isotropic turbulence, Phys. Rev. Lett., 28(2),
76-79 (1972)
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
22
Volume 10, 2022
[20] Ranjan A., & Davidson P. A., DNS of a buoy-
ant turbulent cloud under rapid rotation, in
IUTAM Symp. Proc. Advances in Computa-
tion, Modeling and Control of Transitional and
Turbulent Flows, 452-460 (2016)
[21] Rist U., & Fasel H., Direct numerical sim-
ulation of controlled transition in a at-plate
boundary layer, J. Fluid Mech., 298, 211-248
(1995)
[22] Rogallo R. S., Numerical experiments in ho-
mogeneous turbulence, NASA Tech. Memo,
81315 (1981)
[23] Saddoughi S. G., & Veeravalli S. V., Local
isotropy in turbulent boundary layers at high
Reynolds number, J. Fluid Mech., 268, 333-372
(1994)
[24] Sengupta T. K., High Accuracy Computing
Methods: Fluid Flows and Wave Phenomena,
Cambridge Univ. Press, USA (2013)
[25] Sengupta T. K., A critical assessment of sim-
ulations for transitional and turbulent ows,
in IUTAM Symp. Proc. Advances in Compu-
tation, Modeling and Control of Transitional
and Turbulent Flows, 491-532 (2016)
[26] Sengupta T. K., & Bhaumik S., Onset of
turbulence from the receptivity stage of uid
ows, Phys. Rev. Lett., 107(15), 154501 (2011)
[27] Sengupta T. K., Bhaumik S., & Bhumkar
Y. G., Direct numerical simulation of two-
dimensional wall-bounded turbulent ows
from receptivity stage, Phys. Rev. E, 85(2),
026308 (2012)
[28] Sengupta T. K., & Bhole A., Error dynam-
ics of diusion equation: Eects of numeri-
cal diusion and dispersive diusion, J. Comp.
Phys., 266, 240-251 (2014)
[29] Sengupta T. K., Dipankar A., & Sagaut P.,
Error dynamics: Beyond von Neumann analy-
sis, J. Comp. Phys., 226(2), 1211-1218 (2007)
[30] Sengupta T. K., Singh H., Bhaumik S., &
Chowdhury R. R., Diusion in inhomogeneous
ows: Unique equilibrium state in an internal
ow, Computers & Fluids, 88, 440-451 (2013)
[31] Skote M., & Henningson D. S., Direct numer-
ical simulation of a separated turbulent bound-
ary layer, J. Fluid Mech., 471, 107-136 (2002)
[32] Smith L. M., & Yakhot V., Finite-size eects
in forced two-dimensional turbulence, J. Fluid
Mech., 274, 115-138 (1994)
[33] Suman V. K., Sengupta T. K., Jyothi Durga
Prasad C., Surya Mohan K., & Sanwalia D.,
Spectral analysis of nite dierence schemes
for convection diusion equation, Computers
& Fluids, 150, 95-114 (2017)
[34] Suman V. K., Sengupta T. K., & Mathur J.
S., Eects of numerical anti-diusion in closed
unsteady ows governed by two-dimensional
Navier-Stokes equation, Computers & Fluids,
201, 104479 (2020)
[35] Tan R., Ooi A., & Sandberg R. D., Two-
dimensional analysis of hybrid spectral/ nite
dierence schemes for linearized compressible
Navier–Stokes equations, J. Sci. Comp., 87(2),
42 (2021)
[36] von Neumann J., & Richtmyer R. D., On the
numerical solution of partial dierential equa-
tions of parabolic type, Los Alamos Rept., Se-
ries A, LA-657, 1-17 (1947)
[37] Yeung P. K., Donzis D. A., & Sreenivasan K.
R., Dissipation, enstrophy and pressure statis-
tics in turbulence simulations at high Reynolds
numbers, J. Fluid Mech., 700, 5-15 (2012)
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
23
Volume 10, 2022
A Appendix
The rst step in performing global spectral anal-
ysis for a numerical scheme is to obtain its com-
plex numerical amplication factor GNum. Here,
we provide GNum for RK2-, RK3- and RK4-
Fourier spectral methods for the linear convection-
diusion equation. The derivation is based on the
global spectral analysis for general nite dier-
ence schemes for the model linear convection and
convection-diusion equations provided in [24, 33].
The interested reader should consult these refer-
ences for complete details.
Following the approach in [33], the complex nu-
merical amplication factor for RK2-, RK3- and
RK4-Fourier spectral methods is given as
GNum,RK2=1−B+B2
2(11)
GNum,RK3=1−B+B2
2−B3
6(12)
GNum,RK4=1−B+B2
2−B3
6+B4
24 (13)
where B=Pe(k∆x)2+iNc(k∆x), with Ncand Pe de-
noting the CFL and Peclet numbers, respectively.
From the complex amplication factors given
above for the RK2-, RK3- and RK4-Fourier spec-
tral methods, the numerical diusion coecient
(
ν
Num), phase speed (cNum) and group velocity
(Vg,Num) are obtained as
ν
Num
ν
=−ln(|G|)
(k∆x)2Pe (14)
cNum
c=−1
(k∆x)Nc
tan−1ℑ(G)
ℜ(G)(15)
Vg,Num
Vg,Phys
=−1
Nc
d
d(k∆x)tan−1ℑ(G)
ℜ(G) (16)
where Gtakes GNum,RK2or GNum,RK3or GNum,RK4
depending on the scheme chosen for analysis.
Contribution of individual authors to
the creation of a scientic article
(ghostwriting policy)
Conceptualization and supervision: Tapan K Sen-
gupta
Analysis, investigation and visualization: Vajjala
K Suman and Prasannabalaji Sundaram.
Writing, reviewing and editing: Tapan K
Sengupta, Aditi Sengupta, Vajjala K Suman,
Prasannabalaji Sundaram
Sources of funding for research
presented in a scientic article or
scientic article itself
This study did not receive funding of any form.
Creative Commons Attribution License 4.0
(Attribution 4.0 International, CC BY 4.0)
This article is published under the terms of the Creative
Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.en_US
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.2
Tapan K. Sengupta, Vajjala K. Suman,
Prasannabalaji Sundaram, Aditi Sengupta
E-ISSN: 2415-1521
24
Volume 10, 2022