12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

THE ROLE OF AN OPTIMAL MODELING BASISIN VARIATIONAL RESOLVENT ANALYSIS

Benedikt Barthel
Graduate Aerospace Laboratories
California Institute of Technology
1200 E California Blvd, Pasadena CA, 91125
bbarthel@caltech.edu

Salvador Gomez
Graduate Aerospace Laboratories
California Institute of Technology
1200 E California Blvd, Pasadena CA, 91125

Beverley J. McKeon
Graduate Aerospace Laboratories
California Institute of Technology
1200 E California Blvd, Pasadena CA, 91125

ABSTRACTIn this definition introduced by Barthel et al. (2022), theresolvent response modes are defined as stationary points of
the operator norm of the linearized dynamics. The variational
formulation discussed here constitutes an extension fromthe min-max principle principle which concerns only the
optimal resolvent mode, to what Barthel et al. (2022) coined“variational resolvent analysis” (VRA) which concurrently
defines all resolvent modes.

This work is concerned with the variational formulationof resolvent analysis (VRA) introduced by Barthel et al.
(2022). In this framework, resolvent modes are approximatedby projecting onto a lower dimensional modeling basis. There
are no inherent restrictions on this basis, but in general theefficacy of the model relies on a priori knowledge of thespatial support of the modes of interest. We apply the VRA
framework to a a (1D) Ginzburg-Landau model as well asa (2D) streamwise developing mean flow to investigate theeffect of the boundary conditions of the modeling basis on
the accuracy of the VRA approximation. For the latter we
employ a modeling basis of local (1D), and thus streamwise
periodic, resolvent modes. We find that for weak to moderatestreamwise development of the background flow the periodic
basis functions provide accurate approximations of the global
modes. However, as the streamwise dependence of the
background flow grows, the VRA approximation degrades.

1This variational definition is based on the solutions of
the Euler-Lagrange equations associated with the constrained
variation of the operator norm of the linearized dynamics.
Critically, this definition does not involve the inversion of
any operator, which is useful from both a theoretical andpractical sense. The inversion of large matrices is both costly
and obscures the intuitive interpretation of the underlying
linear differential operator. While in general the resulting
Euler-Lagrange equations remain difficult to solve exactly,
this variational formulation allows for the approximation ofresolvent modes as an expansion in any convenient basis, forexample the much cheaper one-dimensional resolvent basis
in a two- or three-dimensional problem, an analytical basis
such as that described by Dawson & McKeon (2019) or adata-driven one. Further, it requires only the eigenvalue
decomposition of a matrix of reduced size. In this sense
the freedom to choose a convenient modeling basis is both astrength and a weakness, as it allows for greater flexibility,
but places the burden of choosing an efficient basis on theuser. In order to provide any practical advantage, the modelingbasis should be small compared to the dimension of the fullsystem, however it must be large enough to (reasonably) span
the modes being modeled. While in many cases physicalmechanisms such as localization about a critical layer orsymmetries such as periodicity within a domain can guide thechoice of modeling basis, this may not always be the case. In
this paper we expand on the results of Barthel et al. (2022)
and investigate in more detail the limitations of the variational
approximation of global resolvent modes when the boundary
conditions of the modeling basis do not match those of thelinear operator under investigation.

Introduction

Resolvent analysis (RA) can be used to give insight
into the forced response of a linearized dynamical system.
This concept was introduced by Trefethen et al. (1993) and
Jovanović & Bamieh (2005) who considered the stability
and amplification of linearly stable flows to external forcing.
These ideas were later applied to turbulent flows by McKeon& Sharma (2010) who interpreted the nonlinear term in the
Navier-Stokes equations (NSE) as a forcing to the linearized
system. The conceptual framework of RA is inspired bycontrol theory (CT), such that the resolvent operator, theinverse of the linearized operator, is interpreted as a transfer
function from the forcing to the response. A singular value
decomposition (SVD) of the discretized resolvent operator
provides two distinct orthonormal bases (left and right singularmodes) for both the response and the forcing, ordered bya set of gains (singular values) which quantify the linear
amplification of the system.
In this work we take an alternative approach and propose
an equivalent definition based on the calculus of variations.

112th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
2

A variational definition of resolvent modesof amplitudes a j . The Euler-Lagrange equations then take theform of an eigenvalue problem

As described by Barthel et al. (2022), the resolvent
response modes may be defined as the stationary points,q∗ , of the operator norm of L under the condition that theargument q∗ satisfies some norm constraint. More explicitly,
the resolvent modes of the linear operator H are defined as thestationary points of the functional
J = ∥Lψ∥2a(1)∥ψ∥2b = 1

(2)

Ma − σ −2 Qa = 0,

where M, Q ∈ Cr×r . The eigenvectors a contain the amplitudes
a j which optimally approximate the resolvent response modes
and the σ are the approximate singular values. Throughout
the paper we use r to refer to the size of the reduced system(7) and n to refer to the size of the original system, and forlarge systems we achieve a model reduction, r/n, of up to twoorders of magnitude. In this work we use analytical or 1D
local resolvent modes as a modeling basis. However, othertypes of basis are possible, for example see Towne et al. (2015)
for a discussion of using data as a modeling basis for global
resolvent modes in the context of turbulent jets.
While the freedom to choose a modeling basis offers
flexibility in the application of the VRA method, it alsopresents inherent limitation. Namely one may not know apriori what basis is optimal, or how many basis elements are
required for satisfactory convergence. This is of particularconcern when analyzing flows with non-periodic and nonhomogeneous boundary conditions. For example, Barthel
et al. (2022) found that when the boundary conditions of
the modeling basis q j differed significantly from the natural
boundary conditions of the linear operator L, the accuracyof the VRA approximation deteriorated. Those authors usedstreamwise periodic local resolvent modes as a modelingbasis for global modes for a streamwise developing boundary
layer. They found that when the global modes underinvestigation were localized close to the wall, and thus not
significantly affected by the streamwise development of themean flow, the streamwise periodic modeling basis led toaccurate reconstructions. However, when the global modes
extended into the wake region, and thus subject to thestrong streamwise development of the mean flow, the VRA
approximation deteriorated.

subject to the constraint

where the norms ∥x||a and ∥x||b need not be the same. For
simplicity, each norm ∥x||i will be associated with the innerproduct ⟨a, b⟩i = a∗ Qi b, where Qi is a positive definite matrix.
The method of Lagrange multipliers allows us to formulate aconstrained variational problem and define a Lagrangian
L (ψ) = ψ H LH Qa Lψ − σ −2 ψ H Qb ψ.

(3)

The resolvent response modes of H = L−1 are then defined as
the solutions to the Euler-Lagrange equations given by
LH Qa Lψ j = σ −2
j Qb ψ j .

(4)

The resolvent forcing modes are recovered throughφ j = σ j Lψ j .

(5)

We would like to highlight a specific advantage of the
proposed variational formulation in Equations 1 and 2. We
note that the norms ∥x∥a and ∥x∥b norms need not only not
be the same, they need not be quadratic. In many physicalproblems the norms of interest are not quadratic, such as incompressible flows where the temperature contributes a linear
term in the total energy. This variational formulation presents astraightforward way to define a ‘resolvent basis’ in such cases.
The SVD based definition is inherently limited to quadratic
norms since in a practical sense it is based on the Euclidean
2-norm. However, in this paper we focus on quadratic norms
to allow for the comparison of our methods to classical SVD
based algorithms.

3

4

Ginzburg-Landau Model

To investigate the effects of using a modeling basiswith suboptimal boundary conditions we apply the VRA to
the simple case of the Ginzburg-Landau (GL) equation, inparticular, the case studied in Bagheri et al. (2009). Here, wedenote the linear operator
A = −iω + L = −iω − ν∂x + γ∂x x + (µ0 − c2µ ) +µ2 2
x , (8)2

and the resolvent operator, H = A −1 for x ∈ R, assuming
compact support. The coefficients ν, γ, µ0 , and cµ are chosen
to be the same as in Bagheri et al. (2009). Here the frequency,
ω is set to 2 and µ2 is either −0.01 or 1. Finally, the standard
inner product over R is chosen such that

Estimation of Global Modes

In general, the Euler-Lagrange equations (4) are bothanalytically intractable and computationally intensive forcomplex flows with multiple non-homogeneous spatial
dimensions. However, the variational definition provides aconvenient way to estimate resolvent modes as an expansion
in any convenient known basis: q j (x) with ( j = 1...r). The
global resolvent response modes may then be expanded as
ψ = a jq j

(7)⟨a, b⟩ =Z ∞−∞a∗ (x)b(x)dx.

(9)

As noted by Cossu et al. (2009), the GL equation has
many of the hallmarks of the problems studied in wallbounded shear flows, namely the convective nonnormality andnonparallel effects, paramaterized by the coefficients ν and µ2 ,respectively. By varying µ2 between −0.01 and −1 we can(6)which upon insertion into (3) transforms the continuous vector
field q ∈ C∞ into a discrete field a ∈ Cr , where a is the vector

212th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

3

model the effect of increased nonparallel terms in the operator
on the choice of modelling basis.
This operator is discretized to create a matrix
representation, with A and H being the discretized versionsof A and H , respectively. Differentiation is performed
with fourth-order central difference schemes using N = 1920
equispaced grid points for x ∈ [−L, L] with L = 50. Dirichletboundary conditions are enforced in A. Integration is
performed with a trapezoidal-rule integration scheme, whichis enforced through the diagonal, positive definite matrix W.
As such, the inner product is discretized such that⟨a, b⟩ = a∗ Wb.

<j

2
1
0

10

15

20

j
Figure 1: Comparison of singular values using 23 basis

elements for the GL system with µ2 = −.01. The
blue triangles and red dots are the VRA approximatedsingular values using a Gaussian basis and Fourier basis,
respectively. The black open circles are the singularvalues computed using the SVD.

(11)

jAi j

0.2

where r denotes the number of basis elements. Here, σr =2L
2L
.675r and xr,h = r j − L, so that the width and location of
the Gaussians depend on the number of basis elements. The
second set of basis elements are chosen to be Fourier modes
such that
12π
bFj (x; r) = √ exp(i j x).
L
L

5

(10)

Although A is sparse, the required matrix operations forresolvent analysis scale with O(N 3 ).
To investigate the role of the boundary conditions, we will
illustrate the VRA with two sets of modeling basis elements.
The first set of basis elements will be created with Gaussians,
centered at different locations. The basis is defined as(x − xr, j )2
1
bGj (x; r) = p
exp(
),
2
2σr2
2πσr

0

0

-0.2

jAi j

0.2
0

-0.2

(12)

-50

0

50

-50

x

0

50

x

(a) Response modes
Note that the Fourier basis does not have compact support inall the basis elements, unlike the Gaussian basis.
In Figure 1, the singular values from the VRA
approximation are compared with the SVD based modes using
23 basis elements with µ2 = −0.01. The blue triangles denotethe singular values from the Gaussian basis, the red circles
denote the singular values from the Fourier basis, and the black
circles are the singular values from the SVD. Note that theconvergence is best for the leading singular values and drops
off for the higher order modes, which is consistent with Barthel
et al. (2022).
In Figure 2, the response modes and the forcing modes areplotted to compare the approximation of the resolvent modes
using 23 basis elements for µ2 = −0.01. The colors are thesame as in Figure 1 and here only the leading four response andforcing modes are plotted. The forcing modes are computed
using Equation 5. The error in the forcing modes is that thematrix A is not a directional amplifier in the direction of theleading forcing modes. Any error in the approximation of
the response modes becomes amplified when computing theforcing modes.
Finally, the convergence of the singular values, responsemodes, and forcing modes for µ2 = −0.01 and µ2 = −1 isplotted in Figures 3 and 4, respectively. The error in the
singular values is defined as eσ ,i = |σi − σ̃i |, where σ̃i denotes
the singular value approximated from the VRA in either basis.
The error in the forcing modes and response modes is definedas eφ ,i = ||φ i − φ̃ i || and eψ,i = ||ψ i − ψ̃ i ||, respectively. Foreφ ,i and eψ,i , the norm denotes the induced norm from theinner product and the tilde denotes the approximate mode

j?i j

0.2
0

-0.2

j?i j

0.2
0

-0.2
-50

0

50

-50

x

0

50

x

(b) Forcing modes
Figure 2: Comparison of response modes and forcing

modes with 23 basis elements for the GL system with
µ2 = −0.01. The solid lines are the VRA approximatedmodes using the Gaussian basis (blue) and the Fourier
basis (red) while the dotted black lines are the SVD
modes. For each subplot j = 1 − 4 from top left.

from the VRA approximation. In this case, the number ofbasis elements is varied, and generally, convergence increases
with increasing number of basis elements. The closed circles
denote the error with the Fourier basis, while the open triangles
denote the error with the Gaussian basis. For larger number ofbasis elements, r, the Gaussian basis converges faster than the
Fourier basis. The reason for this quick convergence is because312th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
the Fourier basis is not compact whereas the resolvent modesof the Ginzburg Landau operator are. The broad support of theFourier modes means that a large number of them are required
to resolve the resolvent modes, which have compact support.
This is similar to resolving a localized bump with a Fourier
series. The Gaussian basis does not have this limitation and
it is better able to resolve the modes. Comparing Figures 3
and 4, it is clear that the convergence in both bases is slower
with the increase in µ2 . This decrease in convergence rate is
likely due to the fact that the term µ22 x2 in Equation 8 acts as apotential well and increasing the magnitude of the coefficient
makes the term stronger, effectively making the potential wellthinner thereby making the resolvent modes narrower. Thesenarrow modes pose problems for the Gaussian basis functions
for small r, since their widths depends on r, but ultimately thisbasis can resolve these modes as shown in Figure 4. While
the Fourier basis can resolve these modes, this basis is not
as efficient for µ2 = −1 by the fact that the Fourier series iseffectively resolving a more localized signal in space.
Finally, in Figure 3, the convergence of the forcing modes
is slower than the convergence of the response modes because
A amplifies the error in the approximation. Although notshown here, this is also the case for the forcing modes ofµ2 = −1

Error in <
100

Error

10!1
10!2
10!3
10!4

Error in A

Error

10!1
10!2
10!3
10!4
10!5

Error in ?

Error

10!1

5

10!2

Next we consider resolvent modes computed abouta streamwise developing mean flow, the configuration
investigated in Barthel et al. (2022). In this case the reference
2D resolvent modes are computed using the NS operator
linearized about the 2-dimensional, 2-component (2D/2C)
mean flow, Ū(x, y) modelled with the Monkewitz composite
profile (Monkewitz et al., 2007) with inlet friction Reynolds
number Reτ = 1200. The modes are parameterized by the
spanwise wavenumber kz and the temporal frequency, ω,which are defined with respect to inlet δ99 and free stream
velocity.
The global reference modes are computed using
Chebyshev collocation in the wall normal direction andFourier differentiation in the streamwise direction. The natural
boundary conditions relevant to the fluctuations of the stateq = [u, v, w, p]T about the spatially developing mean velocity
are: u(x, 0) = 0, vy (x, y) = 0, and uy (x, ymax ) = 0. An artificial
sponge layer is applied to damp any artificial reflections due tothe boundary conditions as in Ran et al. (2017, 2019), which
also has the effect of damping the modes in this region andrelaxes the periodic constraint on the modes. This is in contrast
to the modeling basis of 1D resolvent modes, computed using
ikx x , which
the inlet mean velocity, q j (x, y) = ψ 1D
j (y; kx , kz , c)eare periodic in the streamwise direction and have support
across the entire domain. The reference modes are computed
using an LU decomposition and Arnoldi Method which isapplied as in Sipp & Marquet (2013) and Schmidt et al. (2018)
to compute the SVD of the resolvent by solving linear systems,
as opposed to computing the matrix inverse. See Barthel et al.
(2022) for a more details of both the modeling basis as well as
the numerical strategies used to compute the reference modes.
We show the results for three wavenumber - frequency
combinations, [kz , ω] = [301, 5.6], [75, 2.8], and [2π, 1.9]. heratio of the number basis elements to the degrees of freedom
in the fully discretized system is 1/288, 5/648, and 5/432,
respectively, a reduction by two orders of magnitude. These
modes exhibit support over an increasing wall-normal extent
and are thus increasingly affected by the wake development

10!3
10!4
10!5
20

40

60

80

100

120

r
Figure 3: Error in the singular values, response modes,

and forcing modes using r Gaussian basis elements(open triangles) or r Fourier basis elements (closed
circles) for the GL system with µ2 = −.01. The red,
blue, green, and purple denote the error in the 1st, 2nd,
3rd, and 4th mode, respectively.

10!1

Error

10!2
10!3
10!4
10!5
10!6
10!7
20

40

60

80

100

Streamwise Developing Boundary Layer

120r
Figure 4: Error in the response modes using r basis

elements, with the same markers as figure 3, for the GL
system with µ2 = −1.

4

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
is sufficient mismatch between the boundary conditions of themodeling basis and the linear operator whose resolvent modes
are being modeled. We illustrated this phenomenon throughthe analysis of a Ginzburg-Landau model and found thatmodeling basis which satisfy the proper boundary conditions
lead to drastically improved convergence properties. We
also applied the VRA approach to the NSE governing thefluctuations about a streamwise developing boundary layer
flow. In this case the VRA approximation using a streamwise
periodic modeling basis is very effective in modeling near wallmodes, but deteriorates as the global modes being modeled
extend further into the wake region of the boundary layer.
We acknowledge funding from ONR grants numbers:
N00014-17-1-2307 and N0014–17–1-3022.

of the mean velocity. This allows us to to illustrate how theincreasing streamwise development of the mean flow affects
the accuracy of the VRA reconstruction using the streamwise
periodic modeling basis. We will refer to the results from these
combinations as short, medium, and tall, respectively.
The VRA approximation of the first two resolvent
response modes is compared to the SVD-based reference inFigure 5. The short modes, in Figures 5a and 5b, are localizednear the wall extending only up to approximately y+ = 30. In
this region, the flow is essentially parallel which makes theFourier (in x) modeling basis a nearly optimal choice to model
these short modes and results in a VRA approximation that is
almost indistinguishable from the reference mode.
The medium modes, shown in Figures 5c and 5d, extend
further, up to about y+ = 100, these modes begin to exhibit
signs of non-periodicity in the streamwise direction. In thiscase the Fourier modeling basis is no longer capturing theexact boundary conditions of the reference modes. However,
the error in the VRA approximation is largely confined tothe streamwise boundaries of the domain. In the bulk of thedomain the VRA approximation captures the qualitative, andquantitative features of the reference. It is interesting to notethat the support of the VRA approximation is generally shiftedslightly upstream relative to the reference. The VRA modes
tend to be roughly centered in the domain, while the reference
is increasingly concentrated towards the downstream side as
the streamwise development strengthens.
Finally, we plot the tall modes in Figures 5e and 5f. These
modes have support all the way to approximately y+ = 500
and are thus strongly affected by the streamwise developmentof the mean flow. Here the downstream bias of the referencemode is even more pronounced, with the bulk of the modes
amplitude being in the downstream half of the domain. The
VRA approximation largely fails to capture this asymmetry,displaying support throughout the domain. Despite this, wenote that the VRA approximation does capture the amplitude
and the fundamental streamwise wave number observed in thereference mode.
The VRA approximation recovers the near wall modes
more faithfully than the taller modes. Part of this can be
attributed to the increase in nonparallel effects of the mean
as distance from the wall is increased. This can serve as
an analogy to what was observed in the case of the GL
operator with increase in |µ2 | since this parameter increased
the nonparallel terms in the GL operator. On the other hand,the effect of the distance from the wall does not constrain thespatial support of the modes like in the GL case. Nonetheless,the choice of modelling basis for this operator is suboptimal
because of the significant streamwise variation in the mean andthe presence of nonreflecting boundary conditions that limit
the periodicity assumed in the basis elements. In particular,
this can be gleaned by noting that the true global modes in
Figures 5e and 5f abruptly decay to 0 at around x = 32 becauseof the sponge. This abrupt decay behaves like a Heavisidefunction, which requires a lot of Fourier modes to resolve.
On the other hand, for the near wall modes, these modes had
already begun to decay due to viscous dissipation downstream
so the effect of the sponge setting the modes to 0 is minimal.

6

REFERENCES
Bagheri, Shervin, Henningson, Dan S, Hoepffner, J & Schmid,Peter J 2009 Input-output analysis and control design
applied to a linear model of spatially developing flows.
Applied Mechanics Reviews 62 (2).
Barthel, Benedikt, Gomez, Salvador & McKeon, Beverley J.
2022 Variational formulation of resolvent analysis. PhysicalReview Fluids 7 (1), 013905.
Cossu, Carlo, Pujals, Gregory & Depardon, Sebastien 2009
Optimal transient growth and very large–scale structures inturbulent boundary layers. Journal of Fluid Mechanics 619,
79–94.
Dawson, Scott T. M. & McKeon, Beverley J. 2019 On
the shape of resolvent modes in wall-bounded turbulence.
Journal of Fluid Mechanics 877, 682–716.
Jovanović, Mihailo R. & Bamieh, Bassam 2005
Componentwise energy amplification in channel flows.
Journal of Fluid Mechanics 534, 145–183.
McKeon, Beverley. J. & Sharma, Ati. S. 2010 A criticallayer framework for turbulent pipe flow. Journal of Fluid
Mechanics 658, 336–382.
Monkewitz, Peter A, Chauhan, Kapil A & Nagib, Hassan M
2007 Self-consistent high-reynolds-number asymptotics forzero-pressure-gradient turbulent boundary layers. Physics
of Fluids 19 (11), 115101.
Ran, Wei, Zare, Armin, Hack, M. J. Philipp & Jovanović,Mihailo R. 2019 Stochastic receptivity analysis of boundary
layer flow. Physical Review Fluids 4 (9), 093901.
Ran, Wei, Zare, Armin, Nichols, Joseph W. & Jovanovic,
Mihailo R. 2017 The effect of sponge layers on global
stability analysis of Blasius boundary layer flow. In 47th
AIAA Fluid Dynamics Conference. American Institute of
Aeronautics and Astronautics.
Schmidt, Oliver T., Towne, Aaron, Rigas, Georgios, Colonius,
Tim & Brès, Guillaume A. 2018 Spectral analysis of jet
turbulence. Journal of Fluid Mechanics 855, 953–982.
Sipp, Denis & Marquet, Olivier 2013 Characterization ofnoise amplifiers with global singular modes: the case of
the leading-edge flat-plate boundary layer. Theoretical andComputational Fluid Dynamics 27 (5), 617–635.
Towne, A., Colonius, T., Jordan, P., Cavalieri, A. V. & Brès,
G. A. 2015 .
Trefethen, Lloyd N., Trefethen, Anne E., Reddy, Satish C. &Driscoll, Tobin A. 1993 Hydrodynamic Stability Without
Eigenvalues. Science 261 (5121), 578–584.

Summary

In this work we investigated the impact of the modelingbasis on the efficacy of the variational resolvent analysis
modeling of global resolvent modes introduced by Barthel
et al. (2022). This approach is susceptible to errors when there

512th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

45

y+

y+

45
15

45

y+

y+

45

15

15
0

1

2

3

4

5

15

6

0

1

2

3

x

y+

95

6

95

30

30

95

95

30

30
0

5

10

15

20

0

5

10

x

15

20

x

(c)

(d)

900

y+

900
500

500

100

100

900

900

y+

y+

5

(b)

y+

y+

y+

(a)

y+

4

x

500
100

500
100

0

5

10

15

20

25

30

0

5

10

15

x

20

25

30

x

(e)

(f)

Figure 5: First two resolvent response modes (ψ j ). Real part of the streamwise component u, j = 1 (left column), andj = 2 (right column). [kz , ω] = [301, 5.6] (a, b), [75, 2.8] (c, d), [2π, 1.9] (e, f). Top panels: true global modes, bottom
panels: VRA model. lower x-axis represents outer units x.

6