Efficient global resolvent analysis via the one-way Navier-Stokes equations. Part 1. Forced response

Resolvent analysis is a powerful tool for modeling and analyzing turbulent flows and in particular provides an approximation of coherent flow structures. Despite recent algorithmic advances, computing resolvent modes for flows with more than one inhomogeneous spatial coordinate remains computationally expensive. In this two-part paper, we show how efficient and accurate approximations of resolvent modes can be obtained using a well-posed spatial marching method for flows that contain a slowly varying direction. In this first part of the paper, we derive a well-posed and convergent one-way equation describing the downstream-traveling waves supported by the linearized Navier-Stokes equations. Integrating these equations, which requires significantly less CPU and memory resources than a direct solution of the linearized Navier-Stokes equations, approximates the action of the resolvent operator on a forcing vector. This capability is leveraged in part 2 of the paper to compute approximate resolvent modes. The method is validated and demonstrated using the examples of a simple acoustics problem and a supersonic turbulent jet.


Introduction
Resolvent analysis (also called input/output analysis or frequency response analysis) is a powerful and popular tool for studying linear energy amplification mechanism within the Navier-Stokes equations. The resolvent operator is derived from the linearized Navier-Stokes (LNS) equations and constitutes a transfer function between inputs and outputs of interest. It has been used to study the linear response of flows to external excitation (Trefethen et al., 1993;Farrell & Ioannou, 2001;Jovanović & Bamieh, 2005;Sipp et al., 2010) and to forcing from the nonlinear terms in the Navier-Stokes equations (McKeon & Sharma, 2010). In the latter context, the method can be derived by reorganizing the Navier-Stokes equations into terms that are linear and nonlinear with respect to perturbations to the turbulent mean flow. The singular value decomposition (SVD) of the resolvent operator associated with the linearized Navier-Stokes equations then identifies modes that are optimal in terms of their linear gain between the nonlinear terms and the perturbations to the mean. These optimal modes have been shown to provide a useful model of coherent structures within turbulent flows, and in particular provide an approximation of the space-time coherent structures educed from data using spectral proper orthogonal decomposition .
The computational resources required to compute resolvent modes depends critically on the number of spatial dimensions that must be numerically discretized. The linearized Navier-Stokes equations nominally contain three spatial dimensions, but the equations can be simplified by expanding the flow variables into Fourier modes in homogenous dimensions, i.e., directions in which the base flow (about which the equations are linearized) does not vary. This drastically reduces the size of the discretized operators that must be manipulated, decreasing the computational cost of the model. Methods in which all inhomogeneous dimensions are discretized are often called global methods (Theofilis, 2011).
Computing global resolvent modes for flows with more than one inhomogeneous direction remains a computationally intensive task. Whereas flows with one inhomogeneous direction can now be tackled on a lap top computer, high-memory workstations and large-scale clusters must be employed for flows with two and three inhomogeneous directions, respectively. These requirements limit the utility of global resolvent analysis for many flows of interest and mitigate some of their advantages compared to fully nonlinear numerical simulations.
Recent efforts to reduce the cost of resolvent analysis have focused on reducing the cost of computing the SVD of the resolvent operator. For example, Moarref et al. (2014) used a randomized singular value decomposition algorithm (RSVD) to compute resolvent modes for a turbulent channel flow (which has only one inhomogeneous direction) and reported that this reduced the cost of the calculations by a factor of two. Ribeiro et al. (2020) offered several improvements for the application of RSVD to resolvent analysis and achieved an order-of-magnitude speedup compared to standard SVD algorithms for a separated flow around an airfoil (which has two inhomogesous directions). Given the reduced cost of computing the SVD, the majority of the remaining cost in their algorithm is associated with computing the action of the resolvent operator on a vector. This is true also for time-stepping methods (Monokrousos et al., 2010;Martini et al., 2021;Farghadan et al., 2021), where the dominant cost is integrating direct and adjoint linear equation in the time domain used to apply the resolvent operator and its adjoint.
In this two-part paper, we develop a method that significantly reduces the cost of computing resolvent modes for flows that include a slowly varying spatial direction, i.e., a direction in which the mean flow is inhomogeneous but changes gradually. This common class of flows includes free-shear flows like mixing layers and jets as well as wall-bounded flows with spatially developing boundary layers on gradually changing objects like a flat plate, cone, airfoil, etc.
In part 1 of the paper, we show how the action of the resolvent operator on a forcing vector can be efficiently and accurately approximated for slowly varying flows using a well-posed spatial marching technique. In Part 2 (Rigas et al., 2021), we show how this capability can be used to compute the singular modes of the resolvent operator using iterative downstream and upstream marching.
Spatial marching methods are commonly applied to slowly varying flows in order to compute approximate eigenmodes or the downstream response to an initial disturbance introduced at some upstream location. The classical tool for these purposes is the parabolized stability equations (PSE), which constitute an ad hoc generalization of classical parallel flow stability theory (Bertolotti et al., 1992;Herbert, 1997). The basic idea of PSE is to separate the flow variables at each frequency into a slowly varying shape function and a rapidly varying wave-like component. Inserting this ansatz into the linearized Navier-Stokes equations leads to a modified set of equations for the shape function, which can be rapidly solved via spatial integration in the slowly varying direction, leading to an approximation of the downstream response to an initial disturbance. The initial disturbance is usually chosen to be a locally parallel eigenmode, in which case the PSE solution is interpreted as a weakly nonparallel eigenmode of the flow. Several authors have recently observed that for flows dominated by a single convective instability, the PSE solution appears to provide a reasonable approximation of the leading resolvent output mode, i.e., the left singular vector of the resolvent operator with the largest singular value (Jeun et al., 2016;Beneddine et al., 2016). However, there are several issues that limit the utility of PSE for approximating resolvent modes. First, PSE does not provide any information about the input resolvent modes or gains, both of which are critical for analysis and modeling. Second, by construction, PSE is incapable of computing sub-optimal modes, i.e., it can only approximate one mode per frequency. Third, PSE can accurately capture the influence of only a single instability mechanism. This limitation stems from the fact that, despite their name, the parabolized stability equations are, in fact, elliptic in the slowly varying direction due to the boundary value nature of the linearized Navier-Stokes equations (Li & Malik, 1997). As a result, the PSE spatial integration is mathematically ill-posed, and regularization methods are required to stabilize the spatial march. Recently, Towne et al. (2019) showed that these regularization methods contaminate the PSE solution, except in cases where the flow is dominated by a single instability mode at each frequency. As a result, PSE is not an appropriate tool for flows in which multiple modes are of interest, including multiple instability mechanisms, transient growth, or acoustics. Towne & Colonius (2015) introduced an alternative method for obtaining fast, approximate linear solutions for slowly varying flows that overcomes these issues by constructing well-posed spatial evolution equations that do not require detrimental PSE-like regularization. Using ideas originally developed for constructing high-order non-reflecting boundary conditions (Hagstrom & Warburton, 2004), the flow variables are decomposed into upstream and downstream propagating waves in the slowly varying direction. An approximate evolution equation is derived for the downstream waves, which can be solved via a well-posed spatial march. The method, which has come to be known as the one-way Navier-Stokes equations (OWNS), has been applied to both free-shear flows such as mixing layers (Towne & Colonius, 2013) and jets (Towne & Colonius, 2014;Rigas et al., 2017b) as well as wall-bounded flows (Rigas et al., 2017a;Kamal et al., 2020) and is typically more than an order-of-magnitude faster than global methods. A schematic comparison of global methods, PSE, and OWNS is provided in figure 1. Unfortunately, the original formulation of OWNS developed by Towne & Colonius (2015) cannot accommodate a forcing term on the linearized flow equations. Such a forcing term is fundamental to resolvent analysis, making this formulation unsuitable for computing resolvent modes.

Marching Methods
To enable the efficient computation of resolvent modes using spatial marching, we introduce a new variant of OWNS that naturally accommodates a forcing terms. The method is formulated in terms of a projection operator that splits the flow variables into upstream and downstream traveling components and which can be applied to the linearized Navier-Stokes equations to obtain evolution equations for each set of waves. Importantly, the projection operator can also be used to split arbitrary forcing terms into parts which influence downstream and upstream traveling waves, enabling its use for approximating resolvent modes.
To distinguish between the two variants of OWNS, we will refer to the original formulation as OWNS-O and the new formulation as OWNS-P, in recognition of their connections with outflow boundary conditions and a projection operator, respectively.
The remainder of this paper is organized as follows. The OWNS-P method is formulated and mathematically analyzed in § 2. It is then demonstrated in § 3 using two example problems -a simple acoustic wave propagation problem and a turbulent jet. Finally, the paper is summarized in § 4.
2 One-way Navier-Stokes equations -a projection approach

Problem setup
We begin with the compressible Navier-Stokes equations in an arbitrary orthogonal coordinate systems, written abstractly as ∂q ∂t = N (q).
(2.1) Equation (2.1) contains mass, momentum, and energy equations, and the state vector q(x, y, z, t) contains an appropriate set of variables, e.g., velocity components and two thermodynamic variables such as density and pressure. Applying the Reynolds decomposition q(x, y, z, t) =q(x, y, z) + q ′ (x, y, z, t) (2.2) and moving terms that are linear and nonlinear in the fluctuation q ′ to the left-and right-hand-sides of the equations, respectively, leads to an equation of the form The left-hand-side of (2.3) is the linearized Navier-Stokes equations, while the right-hand-side vector f (x, y, z, t) contains all remaining nonlinear terms, which within the context of resolvent analysis are interpreted as an external forcing on the linearized equations.
The relationship between the nonlinear term and the linear fluctuation can be expressed in the frequency domain asq = R Gf , (2.4) whereq andf are the Fourier transforms of q ′ and f , respectively, and is the global resolvent operator. Resolvent modes are given by the SVD of the resolvent operator The singular values, which appear within the diagonal positive-semi-definite matrix Σ, give the square root of the optimal gains between the input and output modes defined by the right and left singular vectors contained in the columns of the orthonormal matrices V and U , respectively. In practice, computing resolvent modes requires numerical discretization of (2.3) in all inhomogeneous directions. While the explicit construction of R G can usually be avoided (Jeun et al., 2016;Schmidt et al., 2018;Ribeiro et al., 2020), it is necessary to compute the action of the resolvent operator on an arbitrary forcing vector, i.e., to evaluate the discretized form of (2.4). In what follows, we will show how the slow variation of the mean flow present in many flows can be leveraged to efficiently and accurately approximate the action of the resolvent operator on a forcing vector using a spatial marching method. Then, in part 2 of the paper (Rigas et al., 2021), we will show how this capability can be used to computer optimal forces and responses, i.e., the singular modes of the approximate resolvent operator.

Identifying downstream-and upstream-traveling waves
The first critical step in developing a well-posed one-way equation is to identify parts of the solution that transfer energy in the positive and negative streamwise directions, which we call downstreamtraveling and upstream-traveling, respectively. To this end, we rewrite (2.3) as Here, x ∈ R is the slowly varying direction along which we will apply spatial marching, and y and z are additional, transverse spatial dimensions. In (2.7), we have isolated x-derivative terms arising from the convective terms in the Navier-Stokes equations (A ∂ ∂x ) and from streamwise viscous terms (C ∂ ∂x , D ∂ 2 ∂x 2 ); the linear operator B contains all other terms in the linearized Navier-Stokes equations.
To obtain a one-way equation, we wish to work in terms of a system including only first xderivatives (Towne & Colonius, 2015). This can be accomplished in one of several ways. The viscous terms can be parabolized by redefining the state vector to include both q ′ and ∂q ′ ∂x and writing (2.7) as an expanded system (of twice the orignal size) using this new state variable (Towne, 2016;Harris & Hack, 2020), analogous to the standard approach for solving quadratic eigenvalue problems (Tisseur & Meerbergen, 2001). Alternatively, the streamwise viscous terms can be moved to the right-hand-side of (2.7) and treated as a forcing term, which is later evaluated using the solution at the previous step in the spatial march (Kamal et al., 2020). Finally, following the standard boundary-layer approximation, the streamwise viscous terms can simply be neglected. We have found this simplification to be sufficient for all flows, including both free-shear and wallbounded flows, to which OWNS has been applied to date. Accordingly, we neglect C ∂ ∂x and D ∂ 2 in what follows.
Next, we discretize (2.7) in the transverse directions using a collocation method (such as finite differences) with N c collocation points. The semi-discrete approximation of (2.7) can then be written where q ′ (x, t), f (x, t) ∈ C N and A, B ∈ C N ×N are semi-discrete analogues of q, f , A and B, respectively. The entries of A consist of the values of A at the collocation points, while the matrix B contains discrete approximations of transverse derivatives contained within B as well as modifications required to enforce the desired transverse boundary conditions. The total size of the semi-discrete system is N = N q N c , where N q is the number of state variables in the Navier-Stokes equations (e.g., N q = 5 for three dimensional problems). Equation (2.8) is a one-dimensional strongly hyperbolic system since A is diagonalizable and has real eigenvalues. That is, there exists a transformation T (x) such that whereÃ is a diagonal matrix and the diagonal entries of the sub-matricesÃ contain the positive, negative, and zero eigenvalues of A, respectively. Here, N + , N − , and N 0 denote the number of positive, negative, and zero eigenvalues, respectively, and N = N + + N − + N 0 . The transformation T is known analytically since it is the discretization of the matrix that diagonalizes A.
To derive a one-way equation, it is convenient to work in terms of the characteristic variables of (2.8), which are defined in terms of the transformation T (x), (2.10) These characteristic variables can be split into three components associated with the positive, negative, a zero blocks ofÃ, For later use, we also define which contains only the characteristic variables associated with the nonzero block ofÃ, To be clear, the ± subscript here and throughout the paper does not indicate that we are choosing either the plus or minus characteristic, as in the typical usage of that symbol, but rather both component, as exemplified in (2.12) and (2.13).
In terms of the characteristic variables, (2.8) becomes (2.14) whereB = TBT −1 +ÃT dT −1 dx and f φ = T f . We ultimately wish to obtain the response to a forcing in the frequency domain. However, we proceed by applying a Laplace transform in time, rather than a Fourier transform, to (2.14), giving We will ultimately take η = 0 and set ω to a particular value to obtain the response to a forcing at that frequency, but keeping the possibility of non-zero η will help us distinguish between upstreamand downstream-travelling solutions of (2.14). Up to this point, we have made no approximation (aside from potentially discarding a small subset of viscous terms), and the action of the resolvent operator on a forcing vector could be computed by discretizing (2.15) in x, applying boundary conditions at the beginning and end of the x domain, solving forφ, and inverting the transformation (2.10) to obtainq. However, this involves solving a large system of equations, as discussed in § 2, which constitutes a large fraction of the cost of resolvent analysis. Instead, we will obtain an approximate solution of (2.15) via spatial integration. Directly integrating (2.15) is ill-posed (Li & Malik, 1997;Towne & Colonius, 2015;Towne et al., 2019), leading to exponential divergence of the solution. To circumvent this, we will derive a well-posed one-way equation that can be stably integrated.
To proceed, we isolate the x-derivatives within (2.15), giving (2.17) Here and throughout the paper, the subscripts of a matrix indicate its size, e.g., L 0± ∈ C N 0 ×(N + +N − ) . Equation (2.16) is a differential-algebraic equation (DAE) due to the zero left-hand-side of (2.16b), which occurs because of the zero eigenvalues of A contained in the zero matrixÃ 00 . These zero eigenvalues correspond to points in the base flow where the streamwise velocity is zero or exactly sonic (Towne & Colonius, 2015). The algebraic conditoins in (2.16b) function as a constraint on the allowable form ofφ. Assuming that L 00 is invertible, then (2.16) is an DAE of index 1 and the zero characteristic variableφ 0 is slaved to positive and negative characteristics variables aŝ To obtain a one-way equation, it is convenient to reduce (2.16) to an ordinary differential equation (ODE) forφ ± , which is accomplished by using (2.18) to eliminateφ 0 from (2.16a), leading to While these expressions formally contain L −1 00 , this inverse, and its potential detrimental effect on the sparsity of M, is in practice avoided by reversing the contraction of the system (from (2.16) to (2.19)) in the final set of one-way equations (see § 2.6 and Appendix C).
If N 0 = 0, which is the case, e.g., in subsoninc free-shear flows, then all matrices containing a zero index vanish,φ =φ ± , and (2.19) and (2.21) reduce to the simpler forms Since (2.8) is hyperbolic, the solutionφ ± of (2.19) consists of a summation of downstream-and upstream-traveling modes, i.e., waves that transfer energy in the positive and negative x-direction, respectively (waves that do not propagate were eliminated by the contraction of the system that lead to M). These downstream-and upstream-traveling components of the solution of (2.19) can be identified based on the eigenvalues and eigenvectors of M. Consider the eigen-expansion of the solutionφ ( 2.24) where each v k is an eigenvector of M with an associated eigenvalue iα k , and ψ k is an expansion coefficient defining the contribution of mode k to the solution. The well-posedness theory of Kreiss (1970), which can be thought of as an extension to x-dependent systems of Briggs' criteria (Briggs, 1964), provides a means to distinguish downstream-and upstream-traveling components of the solution: the mode associated with the eigenvalue The eigenvalue decomposition of M can be written where the columns of V + ∈ C N ×N + and V − ∈ C N ×N − and the rows of U + ∈ C N + ×N and U − ∈ C N − ×N contain the left and right eigenvectors associated with the downstream-and upstreamtraveling eigenvalues of M, respectively, which are contained in the diagonal matrices D ++ ∈ C N + ×N + and D −− ∈ C N − ×N − . The eigenvectors are normalized such that (2.28) Using this block matrix notation, (2.24) can be written aŝ and ψ + and ψ − are vectors of expansion coefficients for the downstream-and upstream-traveling modes, respectively. Therefore, the downstream-traveling part of the solution iŝ and the upstream-traveling part isφ

Exact projection operator
We define a projection operator that exactly splits the solutionφ ± into downstream-and upstream-traveling components at each (2.34b) Equation (2.34a) follows from (2.29) and (2.31): and (2.29) and (2.32) can be similarly manipulated to verify (2.34b). Using (2.33) and (2.28), it is straightforward to show that P is a projection operator, i.e., that PP = P.

One-way equation
We now obtain a one-way equation for the downstream-traveling component of the solutionφ ′ ± by applying the projection P to (2.19). This gives To obtain an evolution equation forφ ′ ± , we must move P inside the derivative. Using the chain rule, we have (2.37) Using that P and M commute (since they have the same eigenvectors by (2.27) and (2.33)), the first term on the right-hand-side of (2.36) can also be written in terms ofφ ′ ± , (2.38) Therefore, using (2.37) and (2.38), (2.36) can be written (2.39) Following the same steps, but with I − P replacing P, we similarly obtain (2.40) By neglecting dP dxφ ± in (2.39) and (2.40), we arrive at one-way equations for the downstreamand upstream-traveling components of the solution, When M is x-independent, dP dx = 0 and (2.41) and (2.42) exactly describe the evolution of downstream-and upstream-traveling waves, respectively. When M is x-dependent, dP dx = 0 and (2.41) and (2.42) are approximate. Insight into the nature of the approximation can be gained by solving both (2.39) and (2.40) for dP dxφ ± and equating the expressions, giving ( 2.43) Comparing (2.41) and (2.43) reveals that neglecting dP dxφ ± is equivalent to setting φ ′′ = 0 when calculating φ ′ . In other words, the one-way equation (2.41) neglects the influence of the upstreamtraveling waves on the evolution of the downstream-traveling waves. In the same way, comparing (2.42) and (2.43) reveals that neglecting dP dxφ ± is equivalent to setting φ ′ = 0 when calculating φ ′′ ; the one-way equation (2.42) neglects the influence of the downstream-traveling waves on the evolution of the upstream-traveling waves. As discussed in detail by Towne & Colonius (2015), it is reasonable to neglect the influence of upstream-traveling waves on the downstream-traveling waves (and vice versa) when M is slowly varying in x. Since M inherits its x-dependence from the mean flowq, the one-way equation will yield an accurate approximation of the downstream-or upstream-traveling response to the forcing for slowly varying flows.
Next, we will show that (2.41) and (2.42) are well-posed as one-way equations. This amounts to showing that their eigenvalues correspond to downstream-and upstream-traveling modes, respectively. Focusing first on (2.41), the relevant operator is Compared to the original elliptic operator M, the eigenvectors and downstream-traveling eigenvalues of the one-way operator PM are unchanged, but the upstream-traveling eigenvalues have been eliminated. Therefore, (2.41) is well-posed as a one-way equation and can be solved by integrating the equations in the positive x direction. The relevant operator for the well-posedness of (2.42) is The downstream-traveling eigenvalues have been eliminated, so (2.42) is well-posed as a one-way equation and can be solved by integrating the equations in the negative x direction. This projection-based paradigm for obtaining one-way equations and the projection operator defined in (2.33) were originally derived by Towne (2016) and was recently rederived by Harris & Hack (2020). While it is well posed, as we have shown, its computational efficiency is problematic; the eigen-decomposition of M is required at every x to construct the exact projection P, resulting in an intolerably high computational cost for large N due to the nominal O(N 3 ) scaling of the number of operations required to solve each eigenvalue problem. To obtain a practically useful one-way equation, we construct in the next section an approximation of P that can be efficiently computed.

Approximate projection operator
The following set of recursion equations approximate the action of P on an arbitrary vectorφ ± : Here, we have introduced a set of auxiliary variables {φ j ± : j = −N β , . . . , N β } and a set of complex scalar recursion parameters {β j + , β j − : j = 0, . . . , N β − 1}. N β is the order of the approximate projection.
In Appendix A, we show that the zero-indexed variable is the approximate projection ofφ ± , i.e., thatφ 0 ± = P N βφ ± ≈φ ′ ± . (2.47) The operator is the approximation of the exact projection P that is implicitly defined by the recursions (2.46). It is easy to verify that P N β P N β = P N β , so P N β is itself a projection. Furthermore, comparing (2.33) and (2.48), we see that P N β → P as R +− , R +− → 0. Therefore, the approximation converges if every entry of R +− and R +− converge toward zero as the order of the approximation increases. In appendix A, we show that the (n, m) entry of R +− and the (m, n) entry of R −+ are, respectively,

50)
α +,n is the n th downstream-traveling eigenvalue, α −,m is the m th upstream-traveling eigenvalue, and (w +− ) nm and (w −+ ) mn are scalar weights that do not depend on the recursion parameters.
Since the weights are fixed, the recursion parameters must be chosen such that F(α +,n )/F(α −,m ) goes to zero for all m, n.
A geometric interpretation of F is helpful for choosing parameters that accomplish this objective. Notice that the magnitude of each term in the product defining F is less than one for regions of the complex α plane that are closer to β j + than β j − (|α − β j + | < |α − β j − |) and greater than one for regions that are closer to β j − than β j . Therefore, F(α +,n ) is driven to zero by placing the β j + parameters near the downstream-traveling eigenvalues in the complex plane, and F(α −,m ) is driven to infinity by placing the β j − parameters near the upstream-traveling eigenvalues. This is exactly the same requirement for convergence as derived for the OWNS-O method by Towne & Colonius (2015). This has several important implications. First, as established by Towne & Colonius (2015), if α +,n = α −,m for all m,n, there always exist recursion parameters that make the approximation error arbitrarily small. Second, if the recursion parameters are well-placed, the convergence of the approximation is exponential. Third, any recursion parameters derived for OWNS-O can be used without modification for OWNS-P. A general strategy for choosing recursion parameters was outlined by Towne & Colonius (2015), and effective sets of recursion parameters have been developed for mixing layers (Towne & Colonius, 2014), jets (Towne & Colonius, 2013), subsonic boundary layers (Rigas et al., 2017a), and supersonic boundary layers (Kamal et al., 2020).
Finally, to be rigorous, we must show that (2.41) and (2.42) remain well-posed as one-way equations when P N β is used in place of P. This is confirmed in Appendix B.

Implementation
The recursion equations defining the approximate projection operator define a system of equations of the formφ ′ ± = P 3φ aux (2.51a) P 2φ aux = P 1φ± (2.51b) whereφ aux ∈ C Naux is a vector containing all of the auxiliary variables, the matrices P 1 ∈ C Naux×N and P 2 ∈ C Naux×Naux are defined by the recursion equations (2.46), P 3 ∈ C N ×Naux is a matrix that extracts the projected state from the auxiliary variables via (2.47), and N aux = 2N N β + N 0 . The structure of these matrices is exemplified in Appendix C. From (2.51) we see that applying the projection operator to a vectorφ ± to obtain the projected stateφ ′ ± requires the solution of an a linear system of size N aux ; the cost of this operation compared to those associated with a global solution strategy is discussed in the next section. We stress that in practice we never form the approximate projection operator P N β .
The approximate form of the one-way equation (2.41) can be expressed as a DAE input/output system The expanded state vector is (2.53) and the operators in (2.52) are and (2.55) The input to the system is the forcingf φ = Tf , while the output isφ ′ , the downstream-traveling component of the characteristic variable, from which the OWNS-P approximation of the physical state variable can be obtained using (2.10) asq ′ = T −1φ′ . The DAE (2.52) can be efficiently integrated in the positive x-direction given a specified value forφ ′ ± at the inlet of the domain; this value is physically a boundary condition for the global problem but functions as an initial condition for the spatial integration of the DAE.
An additional issue arises in the practical implementation of the method. Specifically, errors incurred during the numerical integration of (2.52) will not lie entirely in the downstream-traveling subspace. In other words, the numerical approximation ofφ ′ ± collects an error that projects onto the zero eigenvalues of P N β M, which is then propagated along during the march. This causes an accumulation of error that contaminates the solution. Fortunately, there is an easy fix: apply the projection operator to the solution after each step in the march.

Computational cost
The standard approach for obtaining the action of the resolvent operator on a forcing vector requires the solution of a linear system of equations of size N q N x N c , where N x and N c are the number of discretizations points in x and in all transverse directions, respectively, and N q is the number of state variables, e.g., N q = 5. In the following scaling estimates, we drop the dependence on N q since it is a constant for a given problem. Assuming the use of sparse discretization schemes, and direct solution via multi-frontal LU decomposition, the CPU cost (FLOPS) of solving the linear system is found, empirically, to scale as O(N a x N a c ), and the memory usage scales as O(N b x N b c ), where the factors 1 < a ≤ 3 and 1 < b ≤ 2 depend on the sparsity and structure of the matrix and the sophistication of the algorithm employed (Duff et al., 2017). In our global computations for 2D base flows corresponding to turbulent jets presented in part 2 (Rigas et al., 2021), we observed a ≈ 1.6 and b ≈ 1.2, whereas in some preliminary computations for 3D base flows, we observed a ≈ 2 and b ≈ 2. We have also implemented iterative solvers based on GMRES and Bi-CG-Stab, which decrease both a and b to near unity, but, without preconditioning, these were typically slower as the number of iterations required grew large.
The main cost for OWNS-P is solving the system of equations (2.51) of size N aux = 2N β N +N 0 = 2N β N q N c +N 0 if the DAE (2.52) is explicitly integrated or a system of size N ‡ = (2N β +1)N +N 0 = (2N β + 1)N q N c + N 0 if it is implicitly integrated. In both cases, the dominant factor in these size expressions is N c N β , and other terms are dropped. As the sparsity and structure of the OWNS matrices are analogous to those in the global solution, the FLOPS and memory of solving (2.51) . The additional factor of N x in the FLOPS scaling follows from the fact that an equation of this form must be solved at each step in the spatial march. Assuming N x > N β , OWNS-P represents a FLOPS and memory speedup by factors of N a−1 x /N a β and N b x /N b β , respectively.
Clearly OWNS-P will achieve significant savings in memory, and, for sufficiently large N x , significant saving in FLOPS, compared to global methods. These reductions allow problems that would otherwise require high-performance computing resources to be solved on a laptop. The advantage grows for 3D base flows where the factors a and b are larger.
Finally, comparing OWNS-P to OWNS-O, the set of recursion equations that must be solved within the OWNS-P method is roughly twice as large as those required for the OWNS-O method; the index of the auxiliary variables for an approximation of order N β span the range [−N β , N β ] and [0, N β ] for OWNS-P and OWNS-O, respectively. This is the price that must be paid to accommodate the inhomogeneous forcing term.

Examples
In this section, we present two example problems to validate and demonstrate the OWNS-P method. The starting point for both problems is the compressible Navier-Stokes equations, written here in terms of specific volume v, the velocity vector u, and pressure p, All variables have been appropriately non-dimensionalized by an ambient sound speed and density and a problem dependent length-scale. The fluid is approximated as a perfect gas with specific heat ratio γ and constant Reynolds number Re and Prandtl number P r. We have neglected viscous energy dissipation and assumed that the gradient of the dilatation is small. As an example, the linear operators in (2.7) obtained from this form of the Navier-Stokes equations using twodimensional Cartesian coordinates are reported in Appendix D.
For both example problems, the linearized equations are discretized in inhomegeneous transverse directions using fourth-order central finite differences with summation-by-parts boundary closure (Strand, 1994;Mattsson & Nordström, 2004). Far-field radiation boundary conditions are enforced at free transverse boundaries using a super-grid damping layer (Appelo & Colonius, 2009) truncated by Thompson characteristic conditions (Thompson, 1987). The numerical treatment of the xdirection is slightly different in each problem and is reported in the following subsections. Recursion parameters are selected using the strategies described in Towne & Colonius (2015).

Dipole forcing of a quiescent fluid
In this problem, a two-dimensional dipole force is used to excite waves in an inviscid quiescent fluid. The right-hand-side force that is applied to the energy equation is shown in Figure 2(a), and the proper forcing terms are applied to the other equations in order to produce a dipole response in the pressure field, which is shown in Figure 2(b). This is an exact solution of the Euler equations.
We use the OWNS-P method to obtain an approximate solution. The computational domain extends ten wavelengths in both x and y and is discretized using 200 equally spaced points in each The pressure-field obtained by integrating the one-way equations from left to right is shown in Figure 3(a). Clearly, the downstream-traveling waves are accurately captured. Similarly, the pressure-field obtained by integrating from right to left is shown in Figure 3(b). This time, the upstream-traveling waves are captured. Since the governing equations are x-independent in this problem, the full solution can be recovered by summing the downstream-and upstream-traveling solutions, which is shown in Figure 3(c).

Supersonic turbulent jet
Next, we demonstrate our method using the example of a turbulent jet. Resolvent analysis has been applied to jets by numerous authors (Garnaud et al., 2013;Jeun et al., 2016;Schmidt et al., 2018;Lesshafft et al., 2019), and resolvent modes have been shown to capture a rich set of physical phenomena. These include the Kelvin-Helmholtz instability, which leads to large-scale coherent wavepacket structures (Jordan & Colonius, 2013), the Orr mechanism (Tissot et al., 2017;Schmidt et al., 2018), the lift-up mechanism (Nogueira et al., 2019;Pickering et al., 2020), and trapped acoustic waves within the jet core (Tam & Hu, 1989;Towne et al., 2017b). The slow spread of the jet makes OWNS applicable, and the diverse set of physics embedded within the resolvent operator make it a challenging test case.
We will make comparisons between a standard global solution of the linearized Navier-Stokes equation (which we will call the LNS solution from here on out) and those obtained using OWNS-P (and OWNS-O and PSE when applicable). The LNS solution comprises the result of applying the the resolvent operator to a given forcing vector, giving a point of comparison for the approximation of this action provided by the OWNS-P method.
We consider the specific case of a jet with Mach number M = U j /c ∞ = 1.5, Reynolds number Re = ρ j U j D/µ = 1760000, and temperature ratio T j /T ∞ = 1, where the subscripts j and ∞ denote conditions at the jet nozzle exit and in the far-field, respectively. The mean flow about which the Navier-Stokes equations are linearized is obtained from a large-eddy simulation described by Brès et al. (2017). Within the linearized equation, we use a turbulent Reynolds number of Re T = 1760, three orders of magnitude less than the true Reynolds number. This choice is motivated by recent work showing that using an eddy-viscosity model or reduced effective Reynolds number improves both the near-field (Pickering et al., 2021a) and far-field (Pickering et al., 2021b) predictions in free-shear flows.
The linearized equations in cylindrical coordinates are discretized in the radial direction as described earlier, and boundary conditions at the polar axis are enforced using the approach of Mohseni et al. (2002). The physical portion of the domain extends to r/D = 6 and is discretized using 150 grid points with higher concentration around the shear layer at x/D = 0.5, while the damping layer contains an additional 50 points. Since the mean flow is axisymmetric, the azimuthal direction is homogeneous and can be decomposed into a Fourier series. In what follows, we focus on the axisymmetric mode, which is often of foremost interest in the study of jet aeroacoustics (Cavalieri et al., 2012).
The streamwise grid for each method extends from x/D = 0.5 to 30 with equidistant spacing of ∆x = 0.05 (with exception to the PSE method, which determines its own step size) and an additional sponge region included in the LNS computation at three diameters upstream and downstream. The OWNS-O and OWNS-P equations are integrated in the streamwise direction using the 2nd-order backward difference formula and a 5th Order Runge-Kutta Radau II scheme, respectively (Hairer & Wanner, 1971). N β = 13 recursion parameters are used for estimating the OWNS operators. In the following subsections, we show results for two frequencies, St = 0.26 and 0.52, which are close to the frequency of maximum acoustic radiation (Jordan & Colonius, 2013) and maximum near-field resolvent gain , respectively.

Near-nozzle Kelvin-Helmholtz forcing
First, we compute the downstream response of the linearized equations to a disturbance near the nozzle exit, and compare the OWNS-P solution to those obtained using PSE, OWNS-O, and LNS. Specifically, the local Kelvin-Helmholtz eigenfunction is specified at x/D = 0.5 as an initial condition in the PSE and OWNS marches and as a boundary condition for LNS. The solution can be interpreted as the response to a localized forcing at the upstream boundary. In other words, we are computing the result of applying the resolvent operator, approximated by each method, to a forcing vector that is localized at the upstream boundary and which specifically excites the Kelvin-Helmholtz instability.
Results of this test case are shown in figure 4, which shows the pressure field computed by each method for St = 0.26 and 0.52. The OWNS-P solution closely matches the LNS solution for both frequencies, demonstrating that the OWNS-P method provides a good approximation of the action of the resolvent operator on the Kelvin-Helmholtz forcing vector. The OWNS-P solution also matches the OWNS-O solution, which has been previously validated for a similar turbulent jet (Towne & Colonius, 2015). As shown by Sinha et al. (2014) and others, PSE also provides a reasonable solution for supersonic jets forced by a Kelvin-Helmholtz mode. PSE performs well in this case because its underlying assumptions are valid; we have artificially excited a single instability mode by using the Kelvin-Helmholtz mode as an initial condition. While the acoustic radiation is slightly damped compared to the other three solutions, espeically for St = 0.26, it is also reasonably captured because its wavelength is nearly the same as the Kelvin-Helmholtz mode for low-supersonic jets (Towne et al., 2019). However, PSE cannot handle more complex forcings such as those considered next.

Volumetric forcing using LES data
Second, we use both OWNS-P and the standard LNS approach to compute the response of the linearized equations to global forcing vectors extracted from LES data. The right-hand-side forcing term f in (2.3) is obtained by computing the two left-hand-side terms using the LES approximation of the nonlinear Navier-Stokes operator N . The first term can be computed as sinceq does not depend on time. The second term can be approximated as for some ǫ ≪ 1 (de Pando et al., 2012). We take ǫ = 10 −7 and note that the results are nearly independent of ǫ over a range spanning at least four orders of magnitude. Details of the procedure can be found in Towne (2016). After performing an azimuthal Fourier transform, an ensemble of realizations ofq andf are obtained from the LES data by segmenting q ′ and f info overlaping blocks, each containing N f snapshots of data. We use blocks of length N f = 256 with 80% overlap and employ the w C ∞ 4 window proposed by Martini et al. (2019) to minimize spectral leakage. This results in 189 realizations of the flow at discrete frequencies, separated by ∆St = 0.02604.
Each of theseq,f pairs approximately satisfy (2.4). That is, applying the resolvent operator tô f yieldsq. This relationship is approximate for the LES data for two reasons. First, by necessity, we used a finite length DFT in place of an infinite and continuous Fourier transform to obtainq and f , leading to aliasing and spectral leakage. In particular,q exhibits a relatively flat spectrum up to moderate frequencies, making it especially susceptible to aliasing (Towne et al., 2017a). Second, the forcing term f was defined implicitly using the LES operator as described above, whereas our resolvent operator is obtained by explicitly linearizing the Navier-Stokes equations and discretizing the linear equations using numerics different from those within the LES. Thus, we do not expect our LNS and OWNS-P approximation of the action of the resolvent operator onf to exactly reproduce the correspondingq obtained from the LES data. Nevertheless, thef vectors obtained from the LES data provide a physically motivated forcing that can be used to compare the action of the resolvent operators obtained via the OWNS-P march and the standard LNS solution. Figure 5 shows one realization of the forcing and response for St = 0.26 and 0.52. Specifically, we show the component of the forcing that is applied to the energy equation and the response of the pressure field; other components of the forcing and response are qualitatively similar (Towne, 2016). The forcing is rather incoherent and lacks readily visible large-scale structures. In contrast, the pressure field contains distinct wavepacket structures and an acoustic beam emanating from the near-field to the far-field. Despite the lack of structure within the forcing, its introduction as a volumetric forcing term to the LNS and OWNS-P operators produces a response that is quite similar to the LES pressure field (differences can be attributed to factors described earlier).
Of greater relevance to the current study, the OWNS-P response closely matches the LNS response. For both frequencies, the OWNS-P response accurately captures the relevant dynamics in the jet, including near-field structures such as Kelvin-Helmholtz wavepackets (similar to those observed in Figure 3) and Orr-type wavepackets (farther downstream, e.g., 15 < x/D, 20), as well as the angle and intensity of acoustic radiation to the far field. This substantial agreement between the LNS and OWNS-P responses indicates that the OWNS-P approximation of the action of the resolvent operator on this forcing vector accurately mimics that of the true resolvent operator. The ensemble of LES flow and forcing realizations can be used to compute the average error of the OWNS-P approximation over many potential forcing vectors. Figure 6 shows the power spectral density (PSD) of the LES, LNS, and OWNS-P pressure fields for St = 0.26 and 0.52, i.e., the squared-amplitude of the response averaged over the ensemble of realizations. In each case, the PSD is scaled by the maximum value of the PSD of the LES data. The PSD of the LNS and OWNS-P solutions accurately mimic that of the LES data. More importantly, the PSD computed from the ensemble of OWNS-P solutions closely matches the PSD computed from the LNS solutions. Both the envelope of the high-energy near-field region and intensity and directivity of the acoustic beam show good agreement for both frequencies.
The PSD of the error between the LNS and OWNS-P solutions, normalized by the maximum value of the PSD of the LNS solution, is shown in Figure 7. Comparing with figure 6, the average error is observed to be at least an order of magnitude smaller than the average solution in both the near-field hydrodynamic region and the acoustic field. Overall, these results indicate that the OWNS-P method provides an accurate approximation of the action of the resolvent operator for the ensemble of forcing realizations extracted from the LES data; we will see in Part 2 of the paper that it also provides an accurate approximation for the optimal forcing vectors that excite the largest response.

Conclusions
Computing resolvent modes for flows with multiple inhomogeneous directions remains a computationally intensive task. In this paper, we have laid the groundwork for a new method that reduces the cost of computing resolvent modes for flows that include a slowly varying direction. Specifically, we showed how the action of the resolvent operator on a forcing vector can be accurately approximated using a spatial marching method. This approach constitutes an extension of the one-way Navier-Stokes methodology introduced by Towne & Colonius (2015) to accommodate a forcing term on the linear equations, which is central to the concept of resolvent analysis. The method is based on an approximation of the projection operator that rigorously splits the solution vector into downstream-and upstream-traveling components. Applying this projection operator to the linearized Navier-Stokes equations yields a well-posed equation governing the downstream evolution of the flow, which can be stably integrated in the slowly varying direction. This projection-based method, which we call OWNS-P, uses the same recursion parameters and inherits the same convergence properties as the previous outflow-based OWNS methodology, OWNS-O, of Towne & Colonius (2015). Unlike the ubiquitous parabolized stability equations, the OWNS-P method is capable of capturing the complete downstream response of the flow rather than a single instability mechanism (Towne et al., 2019).
The OWNS-P method was demonstrated using an acoustic propagation problem and a turbulent jet. For the jet, we showed that the OWNS-P solution faithfully approximates the action of the resolvent operator on both a localized forcing near the nozzle exit designed to excite the Kelvin-Helmholtz instability and a distributed volumetric forcing extracted from LES data. In the latter case, the OWNS-P solution captured all three prominent downstream-traveling mechanisms present in the jet: Kelvin-Helmholtz, Orr, and acoustic waves.
In part 2 of the paper, we will show how the ability to efficiently approximate the action of the resolvent operator on a forcing vector via spatial marching can be leveraged to compute global resolvent modes at substantially reduced cost. This capability could enable computation of global resolvent modes for previously intractable three-dimensional turbulent flows that contain a slowly varying direction, including swept wings and other three-dimensional aerodynamic bodies and jets with non-circular nozzles, e.g., rectangular jets and other complex nozzle geometries typical of modern tactical aircraft.
The OWNS methodology is also applicable to other problems beyond computation of resolvent modes. For example, OWNS could be used to find the downstream receptivity to an incident disturbance, compute kernels for closed-loop flow control, or predict flow statistics based on stochastic forcing. Indeed, any problem to which PSE can be applied can also be addressed using OWNS, and the latter approach will provide greater accuracy for problems containing multiple relevant physical mechanisms, e.g., multiple instability modes, transient growth, or acoustics. The proper choice between the original OWNS-O method and the OWNS-P method developed in this paper is dictated by the need, or lack-there-of, for a forcing term: only the OWNS-P method can accommodate a forcing term, but at the cost of solving a system of recursion equations that is twice as large compared to the OWNS-O approach.
Since D is diagonal, each scalar component of ψ j can be treated separately, for k = 1, . . . , N . The intermediate auxiliary variables (j = −N β + 1, . . . , N β − 1) can then be easily eliminated, leaving The function F (α) is the same as in (2.50). Equation (A.4) can be written for all k as where, F ++ and F −− are diagonal matrices whose entries are the values of F(α) associated with each downstream-and upstream-traveling eigenvalue, respectively. To apply the termination conditions given by (2.46a) and (2.46e), it is necessary to write the lefthand-sides of (A.4) in terms ofφ where U +− ∈ C N + ×N − and so on. Then, (A.4) can be written as Applying equations (2.46a) and (2.46e) and eliminatingφ where P N β is given by (2.48).

Appendix B Well-posedness of the approximate one-way equation
To show that (2.41) and (2.42) remain well-posed as one-way equations when P N β is used in place of P, it suffices to show that P N β → P as η → ∞, since we know that the eigenvalues of PM and (I − P) M exhibit the correct behavior for well-posedness, as defined by (2.25) and (2.26), in this limit. Therefore, the task at hand is to show that R +− and R −+ go to zero as η → ∞. In this limit, M tends to the diagonal matrix −ηÃ −1 . Since it is diagonal, its eigenvector matrices are also diagonal and unitary; U ++ and U −− are appropriate sized identity matrices and U +− and U −+ are zero. Also, the eigenvalues of the asymptotic form of M approach infinity, which causes F (α k ) → 1 for every k since the recursion parameters are bounded as η → ∞ (Towne & Colonius, 2015). Consequently, F ++ and F −− become identity matrices. Putting this all together, we conclude from (A.8) that R +− and R −+ go to zero as η → ∞. Therefore, (2.41) and (2.42) are well-posed when P N β is used in place of P.
in terms of L andÃ. This can be achieved by introducing zero characteristic components to the auxiliary variables satisfyingφ j 0 = −L −1 00 L 0±φ j ± . (C.1) Then, (2.46), (2.47), and (C.1) can be written together in the form shown in (2.51). Here, we exhibit the structure of this system using the example N β = 2: The superscripts in (C.4) indicates which recursion parameter is subtracted, not a power of the matrix. The single plus, minus, and zero matrix subscripts indicate certain columns of a matrix, as in (2.27). For example, Q (−1) + contains the first N + columns of Q (−1) . These truncated blocks appear in (C.3) due to the terminal conditions (2.46a) and (2.46e) within the recursion equations, which are enforced by removing rows and columns corresponding toφ N β − andφ −N β + . Indeed, it is these terminal conditions, and their impact on P 2 , that force the recursions to be solved all at once as a coupled set.

Appendix D Linearized Navier-Stokes operators
Starting from (3.1), the linearized Navier-Stokes equations in two-dimensional Cartesian coordinates can be written in the form of (2.7) with , and q = {v, u x , u y , p} T . Second x-derivatives have been neglected as discussed in § 2.2. Similar expressions can be obtained for other coordinate systems, e.g., three-dimensional Cartesian or cylindrical coordinates, but are omitted for brevity.