Adjoint sensitivity analysis of chaotic dynamical systems with non-intrusive least squares shadowing

Accepted Manuscript Adjoint sensitivity analysis of chaotic dynamical systems with non-intrusive least squares shadowing

Patrick J. Blonigan

PII: DOI: Reference:

S0021-9991(17)30573-9 http://dx.doi.org/10.1016/j.jcp.2017.08.002 YJCPH 7504

To appear in:

Journal of Computational Physics

Received date: Revised date: Accepted date:

20 April 2017 1 July 2017 2 August 2017

Please cite this article in press as: P.J. Blonigan, Adjoint sensitivity analysis of chaotic dynamical systems with non-intrusive least squares shadowing, J. Comput. Phys. (2017), http://dx.doi.org/10.1016/j.jcp.2017.08.002

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Highlights • • • •

A discrete adjoint non-intrusive least squares shadowing (NILSS) is presented. The NILSS approach is closely related to multiple shooting shadowing (MSS). Adjoint NILSS prevents exponential growth in time of the adjoint field. Adjoint NILSS is demonstrated on a simulation of wall-bounded turbulent flow.

Adjoint Sensitivity Analysis of Chaotic Dynamical Systems with Non-Intrusive Least Squares Shadowing Patrick J. Blonigan NASA Ames Research Center, Moffett Field, CA 94035, United States

Abstract This paper presents a discrete adjoint version of the recently developed non-intrusive least squares shadowing (NILSS) algorithm, which circumvents the instability that conventional adjoint methods encounter for chaotic systems. The NILSS approach involves solving a smaller minimization problem than other shadowing approaches and can be implemented with only minor modifications to preexisting tangent and adjoint solvers. Adjoint NILSS is demonstrated on a small chaotic ODE, a one-dimensional scalar PDE, and a direct numerical simulation (DNS) of the minimal flow unit, a turbulent channel flow on a small spatial domain. This is the first application of an adjoint shadowing-based algorithm to a three-dimensional turbulent flow. Keywords: Sensitivity Analysis, Adjoint, Chaos, Shadowing

1. Introduction Adjoint-based sensitivity analysis is a powerful tool for engineers and scientists who use simulations to study physical systems ranging from aerospace vehicles to the climate system of the earth. Using the adjoint enables one to efficiently compute the sensitivity of an objective function to many parameters, which makes adjoint-based sensitivity analysis invaluable for gradient-based design optimization [1], error estimation [2, 3], flow control [4], and uncertainty quantification [5]. Conventional adjoint- and tangent-based sensitivity analysis approaches are mature technologies for steady and some unsteady simulations. However, these approaches break down when applied to timeaveraged objective functions of chaotic dynamical systems [6]. A number of approaches have been proposed to mitigate the issues faced by conventional sensitivity analysis approaches, including the ensemble adjoint approach [6, 7], Fokker-Planck approaches [8, 9], fluctuation-dissipation theorem approaches [10, 11, 12], and shadowing-based approaches including least squares shadowing (LSS) [13, 14]. All of these approaches are quite computationally costly and can be difficult to implement for a given simulation. Multiple shooting shadowing (MSS) is a modification of LSS which reduces the size of the minimization problem solved by LSS [15]. Additionally, MSS can be built using existing primal, tangent, and adjoint solvers. The recently developed “Non-Intrusive LSS” (NILSS) is a modification of MSS [16, 17, 18]. Like MSS, it can be built using existing primal, tangent, and adjoint solvers. Additionally, it solves a relatively small minimization problem after running the tangent solver a number of times, unlike MSS, which requires repeated runs of the tangent and adjoint solvers while solving the minimization problem. This decoupled minimization problem, along with the use of existing tangent and adjoint solvers makes NILSS “Nonintrusive”. This paper presents a discrete adjoint version of NILSS 1 . Additionally, the first application of the shadowing-based adjoint sensitivity analysis approach to a three-dimensional scale-resolving turbulent flow Email address: [email protected] (Patrick J. Blonigan) discrete adjoint approach is pursued because it can be verified to machine precision using the existing tangent NILSS approach 1A

Preprint submitted to Journal of Computational Physics

August 4, 2017

simulation is demonstrated. This paper starts with an overview of shadowing-based sensitivity analysis, followed by derivations of tangent and then adjoint NILSS. Adjoint NILSS is then demonstrated on the Lorenz 63 equations, the Kuramoto-Sivashinsky (K-S) equation, and finally, a direct numerical simulation (DNS) of the minimal flow unit, a turbulent channel flow with the smallest domain that still captures the key near-wall flow structures needed to sustain a turbulent flow [19]. 2. Shadowing-based Sensitivity Analysis for Chaotic Dynamical Systems Suppose one wants to carry out sensitivity analysis on the following nonlinear ODE du (1) = f (u; s), u(T0 ) = u0 dt where u(t) is an n × 1 vector with the state at time t and s is some system parameter. Equation (1) can represent any simulation from the simple three degree of freedom (DoF) Lorenz 63 system to a large DNS with millions of DoFs. Many objective functions of interest to engineers and scientists are time averages of some function J(u(t); s) that depends on the state u(t) ¯ = 1 J(s) ΔT



T1

J(u(t); s) dt

ΔT = T1 − T0 ,

(2)

T0

¯ could be the lift of a stalled airfoil as a function of the geometric parameters that define For example, J(s) ¯ to a parameter s by differentiating equation (2) to the airfoil shape. One can compute sensitivities of J(s) obtain   T1 1 ∂J ∂J  dJ¯ = dt, (3) v(t) +  ds ΔT T0 ∂u t ∂s where v(t) ≡

∂u ∂s

is the solution to the tangent equation   ∂f  ∂f  dv v + , = dt ∂u t ∂s t

v(T0 ) = 0,

(4)

where it is assumed the initial condition of equation (1), u0 , is independent of s. ¯ then it is much more efficient to use an adjoint If there are many parameters s but only one objective J, approach than a tangent one. In this case   T1 1 ∂J ∂f  dJ¯ = dt, (5) w(t) + ds ΔT T0 ∂s t ∂s where w(t) is the solution to the following adjoint equation −

T  dw 1 ∂J  ∂f  w + , = dt ∂u t ΔT ∂u t

w(T1 ) = 0

(6)

Note that the adjoint is solved backwards in time from a terminal condition w(T1 ). Equation (5), shows that only one adjoint solution w(t) is required to compute sensitivities to different parameters s. Conventional sensitivity analysis approaches break down when equation (1) exhibits chaotic dynamics. This is because a chaotic system has at least one positive Lyapunov exponent [20, 21]. The Lyapunov exponents Λi are the growth rates of perturbations to the state u(t) of a dynamical system governed by equation (1). Since state perturbations are well approximated by the tangent equation (4), the Lyapunov exponents are the growth or decay rates of modes of the tangent equation called Lyapunov covariant vectors ψ 1 (u), ψ 2 (u), ..., ψ i (u). If the tangent equation (4) has a ∂f /∂s parallel to the ith covariant vector ψ i (u), or the initial condition v(0) is equal to ψ i (u), the tangent solution v(t) will grow or decay at rate Λi . If any of the Lyapunov exponents are positive, then the tangent solution will likely grow exponentially with time. 2

This means that as the time horizon T0 to T1 in equations (2) and (3) is expanded, the sensitivity also grow exponentially. The Lyapunov covariant vectors are governed by the following equation [22]  ∂f  d i ψ i (u(t)) − Λi ψ i (u(t)). ψ (u(t)) = dt ∂u u(t)

dJ¯ ds

will

(7)

Like equation (4), equation (7) is a tangent equation. Equation (7) shows that the covariant vectors ψ(u(t)) vary with the state u and therefore vary with time, so the perturbation that triggers growth or decay at the rate Λi depends on t. The adjoint solution w(t) also grows or decays backwards in time at the same rates, Λ i , as the tangent solution v(t) grows or decays forwards in time. Each Lyapunov exponent Λ i has a corresponding adjoint Lyapunov covariant vectors, which are governed by the following adjoint equation [22] ∗ d ˆi ∂f  ˆi − ψ (u(t)) = (8) ψ (u(t)) − Λi ψˆi (u(t)). dt ∂u t i ˆi where ψˆi (u(t)) is the ith adjoint covariant vector. If ∂J ∂u |t aligns with ψ (u(t)) corresponding to any Λ > 0 dJ¯ the adjoint will grow exponentially backwards in time [6]. This means that sensitivities ds computed using equation (5) will also grow exponentially as the time horizon T0 to T1 is expanded. Shadowing-based approaches compute special tangent or adjoint solutions that do not grow exponentially with time and avoid the issues encountered by conventional sensitivity analyses. To explain the basic idea of shadowing, we first recall that sensitivity analysis depends on a linearization derived from the difference between a solution u(t; s) and a perturbed solution u(t; s + δs),

u(t; s + δs) − u(t; s) . δs The conventional approach uses the same initial conditions for u(t; s) and u(t; s + δs), v(t) = lim

δs→0

u(0; s + δs) = u(0; s) = u0 .

(9)

(10)

A shadowing-based approach selects a different initial condition for the perturbed solution, u(0; s) = u0 ,

u(0; s + δs) = uS (0, s + δs).

(11)

S

This choice of initial condition makes u(t; s + δs) = u (t; s + δs), the shadow trajectory. The shadow trajectory remains close to u(t; s) for all time, so uS (t; s + δs) − u(t; s)2 is bounded for all time. This means that the tangent v(t) will not grow exponentially if it is defined with u(t; s + δs) = uS (t; s + δs) in equation (9) [13]. It has been found that the adjoint of this linearization will also remain bounded [13]. Equation (11) requires that the system governed by equation (1) is ergodic. An ergodic system is a system with long-time behavior that does not depend on initial conditions [20]. This means that a timeaveraged objective function J¯ computed for an ergodic system is independent of the initial condition u(0), a property that holds for a wide range of physical systems of interest to engineers and scientists. For an ergodic system, the existance of a shadow trajectory follows from the shadowing lemma [ 23], which states that Consider a solution u to du = f (u; s) dt If this system has a hyperbolic strange attractor and if some system parameter s is slightly perturbed by some δs: For any δ > 0 there exists ε > 0, such that for every u that satisfies du/dt − f (u; s + δs) < ε, there exists a true solution uS and a time transformation τ (t), such that uS (τ (t)) − u(t) < δ, |1 − dτ /dt| < δ and duS /dτ − f (uS ; s + δs) = 0. Note that  ·  is the L2 norm in phase space. 3

Figure 1: Diagram showing the checkpoints and time segments used in MSS and NILSS. The times t+ refers to a time infinitesimally greater than time t.

The shadowing lemma holds for hyperbolic strange attractors, a special class of strange attractors which have a number of properties including ergodicity and that the tangent solution vector can be decomposed into stable and unstable subspaces[24, 25]. In other words, the Lyapunov covariant vectors span phase space at all points on the surface of a hyperbolic attractor. Although hyperbolic attractors are rare in practice, Gallavotti and Cohen’s chaotic hypothesis conjectures that large chaotic systems behave more like hyperbolic systems than non-hyperbolic systems, which gives good reason to believe that shadowing-based approaches will work for larger chaotic systems like scale-resolving turbulent flow simulations [ 26, 27, 28]. Shadowing-based sensitivity analysis approaches draw heavily from continuation methods, which are used by the dynamical systems community to find homoclinic and heteroclinic orbits of dynamical systems [29, 30]. Shadowing has also been used for noise reduction [31]. The application of these ideas to sensitivity analysis was first proposed by Wang [32] and was successfully applied to the Lorenz 63 system. However, the initial implementation required computing all Lyapunov exponents and covariant vectors, the current algorithm for which is infeasible for larger systems [22]. To avoid computing Lyapunov covariant vectors, Wang et al. introduced LSS [13], which computes the tangent solution corresponding to the shadow trajectory by solving a least-squares minimization problem. For hyperbolic strange attractors, it can be shown that shadowing-based approaches will compute accurate sensitivities with respect to system parameters [14, 33]. LSS has been applied to a number of chaotic systems, most of which are quasi-hyperbolic but not hyperbolic. These include small model problems, for instance the Lorenz 63 system, a few chaotic maps, the non-chaotic Van der Pol oscillator, and a 4 DoF model of a chaotic aeroelastic oscillator [14, 13]. LSS has computed accurate sensitivities in many of these cases, but in some cases exhibited an ergodicity-breaking error [34]. Studying ever larger chaotic systems will show if Gallavotti and Cohen’s chaotic hypothesis holds and if the ergodicity-breaking error is limited to smaller systems and pathological cases. LSS has been applied to only a few larger systems, most notably a spectral method simulation of Homogeneous Isotropic Turbulence [35] and chaotic two-dimensional flow around an airfoil [36]. In all of these cases LSS computed accurate gradients. These results are promising, but the computations were large, time consuming, and complicated to implement. The large cost of LSS arrises from solving the minimization problem [13] 1 min v(t),η(t) 2



T1

T0

2

2

2

v(t) + a η(t) dt

s.t.

  ∂f  dv ∂f  = v+ + ηf (u(t); s), dt ∂u t ∂s t

T0 < t < T 1

(12)

where a is a weighting parameter for the norm of the time dilation term η in the minimization statement 2 . Since the constraint equation is applied for T0 < t < T1 , the discretization of the equation (12) can be reduced to m n-dimensional constraint equations, where m is the number of time steps used to discretize T0 < t < T1 and n is the number of DoFs in equation (1). It can be shown that the discrete Karush-KuhnTucker (KKT) equations for this minimization problem are a mn × mn linear system [13]. For a n = 105 DoF simulation simulated for 10 4 time steps, this system is 10 9 × 109 . The multiple shooting shadowing (MSS) approach was introduced to reduce the size of the minimization problem solved for LSS. It does this by solving a modified version of equation (12) that only minimizes v(t) at discrete checkpoints in time 2 Note that η = 1 − dτ /dt, where τ (t) is a time transformation like the one defined in the shadowing lemma. For more details consult [13]

4

1 T [vi ] vi vi 2 i=0   ∂f  dv ∂f  = s.t. v + + ηf (u(t); s), dt ∂u t ∂s t K

min

(13) ti−1 < t < ti ,

1≤i≤K

(14)

s.t. f (u(t); s), v(t) = 0

(15)

s.t. v(t− i )

(16)

=

v(t+ i )

≡ vi

where index i is the checkpoint index, as shown in figure 1, η(t) is the time dilation term and vi is an n × 1 vector with the tangent initial condition for time segment i + 1. The first constraint equation (14) is the tangent equation corresponding to the governing equation (1). The second constraint equation (15) is required because the contribution of η to the MSS minimization is zero, i.e. a = 0 in the notation of equation (12). Equation (15) ensures that v(t) does not grow linearly with time due to phase shifts. It does this by fixing η(t) so that v(t) is always perpendicular to f (u(t); s), the direction of perturbtions that result in phase shifts. This formulation results in a more straightforward implementation of MSS [ 15]. The final constraint equation (16) ensures that v(t) is continuous at each checkpoint ti . The time between checkpoints should be chosen so that the tangent solution does not grow too large over each segment [15]. For example, one conservative choice of time segment size is ti − ti−1 ≈ 1/Λ1 , where Λ1 is the largest Lyapunov exponent. This choice of ti − ti−1 limits the growth of the magnitude of the tangent v(t) to v(ti ) ∼ O(e1 )v(ti−1 ) ≈ O(2.718)v(ti−1 ) Limiting the growth of v(t) on a each segment helps prevent the MSS KKT system from becoming poorly conditioned [15]. MSS decreases the size of the KKT system for the minimization problem from mn × mn to Kn × Kn. Additionally, MSS can be implemented using a conventional primal, tangent, and adjoint solver [ 15], making it less intrusive than the original LSS algorithms. 3. Non-Intrusive Least Squares Shadowing Like MSS, NILSS divides the time horizon T0 to T1 into K time segments as shown in figure 1 but it uses a decomposition of the tangent solution v(t) to reduce the size of the minimization problem solved [16, 17, 18]. In other words, NILSS is a reduced-order MSS approach (see Appendix B). NILSS reduces the cost of MSS by computing the solutions to p tangent equations with orthogonal initial conditions and expressing the MSS optimization problem as a minimization over a linear combination of these p tangent solutions. It has been shown that the number of extra tangent equations to solve should be equal to the number of positive Lyapunov exponents and that p << n for large scale-resolving turbulent flow simulations [37, 38, 39]. This means that the solution to the NILSS minimization problem is found by solving a Kp × Kp linear system which is smaller than the Kn × Kn problem solved by MSS. The following section explains the tangent and adjoint NILSS algorithms in detail. 3.1. Tangent NILSS NILSS seeks to reduce the cost of solving the MSS minimization problem (13) by considering the following decomposition of the tangent solution for each time segment v(t) = vi∗ (t) + Wi (t)αi , where

vi∗ (t)

is the solution of the tangent equation   ∂f  ∗ ∂f  dvi∗ v + + ηvi∗ f, = dt ∂u t i ∂s t 5

ti−1 < t < ti

ti−1 < t < ti

(17)

and Wi (t) is a n × p matrix whose pth column is Wip (t), the solution to the following tangent equation at time t on segment i  ∂f  p dWip W + ηWip f ti−1 < t < ti = dt ∂u t i Finally, αi is a p × 1 vector with the weights for each column of W (t) for the ith time segment. The matrix of initial conditions Wi (ti−1 ) can be any rank p matrix, and v ∗ (t0 ) = 0. The choice of each v ∗ (ti ) is related to the LSS continuity constraint (16) and will be discussed below. 3.1.1. Karush-Kuhn-Tucker System

Figure 2: Representative schematic of tangent NILSS, showing tangent norms versus time for Wi (t) (blue), vi∗ (t) (red), and v(t) (green). In this case p = 3, so there are there homogeneous solutions. NILSS seeks the linear combination of the blue homogeneous tangent solutions Wi (t)αi that cancel out the exponential growth of the red inhomogeneous tangent v ∗ (t) and obtain the green curve, a bounded, continuous shadowing tangent v(t).

Figure 2 shows that the minimization problem of equations (13)-(16) can be recast using equation (17) as finding the optimal linear combination of homogeneous solutions and the inhomogeneous solution on all segments so that their norm is minimized and the solution satisfies the constraints (14)-(16). To solve equations (13)-(15) the Karush-Kuhn-Tucker (KKT) equations are formed. To find these for NILSS, start by applying the decomposition of equation (17) to the continuity constraint (16) of the LSS minimization statement − + + ∗ vi∗ (t− i ) + Wi (ti )αi = vi+1 (ti ) + Wi+1 (ti )αi+1

To further simplify this expression, take the QR-decomposition of Wi (t− i ) Qi Ri = Wi (t− i ) and set Wi+1 (t+ i ) = Qi to obtain + ∗ vi∗ (t− i ) + Qi Ri αi = vi+1 (ti ) + Qi αi+1

(18)

To ensure continuity of v(t) at ti , the following must hold ∗ T ∗ − vi+1 (t+ i ) = (I − Qi Qi )vi (ti ) ∗ (t+ This means that vi+1 i ) contains any ∗ (t+ basis Qi . Therefore QTi vi+1 i ) = 0 T equation with Qi

(19)

components of v(t) orthogonal to the span of the n × p orthonormal and equation (18) can be simplified by left-multiplying the entire QTi vi∗ (t− i ) + Ri αi = αi+1

(20)

T − βiT Ri = 0, αiT + βi−1

(21)

The dual equation is then

6

as shown in Appendix ⎡ −I ⎢ −I ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ R1 −I ⎢ ⎢ R2 ⎢ ⎢ ⎣

A. Equations (20) and (21) can be combined to form the reduced KKT system ⎤⎛ ⎞ ⎛ ⎞ 0 α1 R1T ⎥ ⎜ α2 ⎟ ⎜ ⎟ 0 −I R2T ⎥⎜ ⎟ ⎜ ⎟ .. ⎥⎜ ⎟ ⎜ ⎟ .. .. .. .. ⎥⎜ ⎟ ⎜ ⎟ . . . . . ⎥ ⎜ ⎟ ⎜ ⎟ T ⎥⎜ ⎟ ⎜ ⎟ 0 α −I RK −I K ⎥⎜ ⎟ ⎜ ⎟ ⎜ αK+1 ⎟ = ⎜ ⎟ 0 −I −I ⎥ (22) ⎥⎜ ⎟ ⎜ ⎟ ⎥ ⎜ β1 ⎟ ⎜ −QT1 v1∗ (t− ) ⎟ 1 ⎥⎜ ⎟ ⎜ ⎟ ⎥ ⎜ β2 ⎟ ⎜ −QT v ∗ (t− ) ⎟ −I 2 2 2 ⎥⎜ ⎟ ⎜ ⎟ ⎥⎜ ⎟ ⎜ ⎟ .. .. .. .. ⎦ ⎝ ⎠ ⎝ ⎠ . . . . RK

−I

∗ (t− −QTK vK K)

βK

Note that if the QR-decomposition in equation (18) is omitted and all Qi are set equal to the identity matrix, the KKT system equation (22) is equivalent to the one derived for MSS, as shown in Appendix B. Equation (22) also be written in the following short-hand form     −I B T  α  0 =− (23) β b B where ⎡ ⎢ ⎢ B=⎢ ⎣

R1

−I R2

⎤ −I .. .

⎥ ⎥ ⎥, ⎦

..

. RK

−I

⎛ ⎜ ⎜ ⎜ α=⎜ ⎜ ⎝

⎛

⎞

α1 α2 .. .

⎜ ⎜ ⎜ β=⎜ ⎜ ⎝ βK−1 βK

⎟ ⎟ ⎟ ⎟, ⎟ ⎠

αK αK+1

β1 β2 .. .

⎞ ⎟ ⎟ ⎟ ⎟, ⎟ ⎠

⎛ ⎜ ⎜ b=⎜ ⎝

QT1 v1∗ (t− 1) QT2 v2∗ (t− 2) .. .

⎞ ⎟ ⎟ ⎟ (24) ⎠

∗ QTK vK (t− K)

In practice, it is better to solve the Schur complement of equation (23) Aβ = −b, or in the notation of equation (22) ⎡ −R2 R1 R1T + I ⎢ −R2T R2 R2T + I −R3 ⎢ ⎢ .. .. ⎢ . . ⎢ T T ⎣ RK−1 RK−1 +I −RK−1 T −RK

A ≡ BB T

(25) ⎤

..

.

−RK−1 T R K RK +I

⎛ β1 ⎥ ⎥ ⎜ β2 ⎥⎜ ⎥ ⎜ .. ⎥⎝ . ⎦ βK

⎞

⎛

⎟ ⎜ ⎟ ⎜ ⎟=⎜ ⎠ ⎝

−QT1 v1∗ (t− 1) −QT2 v2∗ (t− 2) .. .

⎞ ⎟ ⎟ ⎟ ⎠

(26)

∗ (t− −QTK vK K)

The vector α can be computed by using the upper blocks of (23) α = BT β Since Ri is a p × p matrix, equation (25) has Kp unknowns. As stated previously p should be equal to the number of positive Lyapunov exponents for a given chaotic system, which has been found to be much smaller than the total number of DoFs, n, for a number of large chaotic simulations. For the channel flow simulations considered by Keefe et al., the number of positive exponents was found to be roughly 2% of the total DoFs of the simulation [38], A higher Reynolds number channel was estimated to have around 1500 positive exponents for a 1.23 million DoF simulation, so p is 0.1% n for this case [39]. Therefore, solving the KKT system for NILSS is considerably cheaper than solving the Kn × Kn KKT system for MSS. However, it remains to be seen if the total cost of computing the p extra tangents, K QR-decomposition, and solving equation (25) is actually less expensive than MSS in all cases. 7

3.1.2. Implementation NILSS uses much of the same implementation practices as MSS. Firstly, the effect of time dilation term η(t) is computed the same way as in MSS. Like η in equation (14), the scalars ηWip (t) and ηv∗ (t) ensure that Wip (t) and v ∗ (t) are orthogonal to f (u(t); s) for t0 ≤ t ≤ tK . In practice, these terms are not directly computed [15]. Instead, a projection operator Pt is used v(t) = Pt v  (t) ≡ v  (t) − f (u(t); s)

f (u(t); s), v  (t) f (u(t); s)22

(27)

ti−1 ≤ t ≤ ti

(28)

∗

where v  (t) = v  i (t) + Wi (t)αi and therefore ∗

v(t) = Pt (v  i (t) + Wi (t)αi ), The first component,

∗ v  i (t),

satisfies the inhomogeneous tangent equation  ∗ ∂f   ∗ ∂f dv  i = v + dt ∂u t i ∂s

and the pth column of Wi satisfies the following homogeneous tangent equation  p ∂f   p dW  i W = dt ∂u t i

(29)

(30)

It is important to note that equation (30) is identical to equation (4) and that equation (29) is equation (4) with ∂f ∂s set to zero. Therefore, one can use a conventional tangent solver to solve equations (28), (29), ∗ and (30) for v  i (t) and Wi (t), then apply the projection operator Pt to obtain v(t) without any explicit computation of ηWip (t) and ηvi∗ (t). The sensitivity of some time-averaged objective function J¯ to a parameter s can be computed from v  (t) using the MSS sensitivity equation [15, 18]   T K ti   f (u(ti ); s)T v  (ti ) ∂ J¯ 1  ∂J   dJ¯ ¯ v (t) dt + J − J(u(ti ); s) + = (31) 2  ds ΔT i=1 f (u(ti ); s)2 ∂s ti−1 ∂u t where 1 ∂ J¯ = ∂s ΔT



T1 T0

∂J dt ∂s

The second term on the right-hand side removes the contribution of the component of v  (t)  f (u(t); s) to the sensitivity, as shown in detail in the appendices of [15, 18]. ∗ As v  (t) = v  i (t) + Wi (t)αi on segment i, equation (31) can be written as follows dJ¯  T ∂ J¯ = gi αi + h + ds ∂s i=1 K

(32)

where T  1 ¯ ∂J  Wi (t) dt + J − J(u(ti ); s) Fti Wi (ti )  ΔT ti−1 ∂u t    K T t i   1  ∂J   ∗ ∗  ¯ h=  v i (t) dt + J − J(u(ti ); s) Fti v i (ti ) ΔT i=1 ti−1 ∂u t

giT ≡

Fti ≡

1 ΔT



ti

1 f (u(ti ); s)T f (u(ti ); s)22

(33) (34) (35)

8

3.1.3. Algorithm Outline Tangent NILSS is implemented as follows Tangent NILSS Solver Inputs: Initial condition for the governing equations u0 , Spin-up time t0 , Specified time horizon and checkpoints t0 , t1 , ..., tK ¯ Ouputs: Sensitivities dJ/ds 1. 2. 3. 4. 5. 6. 7. 8. 9.

Compute u(t) by solving equation (1) until t = t0 . Set v1∗ (t0 ) = 0 and W1 (t0 ) = Q0 , where Q0 is some n × p unitary, orthonormal matrix. Set i = 1. ∗ Solve for u(t), Wi (t), and v  i (t) from ti−1 to ti using equations (1), (29), and (30), respectively. Compute Qi Ri = Wi (ti ) = Pti Wi (ti ) with a QR-decomposition. Save Qi and Ri . ∗ Compute and save bi = QTi v ∗ (ti ) = QTi Pti v  (ti ). ∗ ∗ T Set Wi+1 (ti ) = Qi and vi+1 (ti ) = (I − Qi Qi )Pti v  i (ti ). Set i = i + 1 and repeat steps 3 to 6. If i = K, proceed to the next step Form the KKT Schur complement (25) and solve it. Compute α = B T β.

10. Compute sensitivities using equation (32). Step 3 of the algorithm requires a primal solver for u(t) and a tangent solver for Wi (t) and vi∗ (t). These three solvers could be run in parallel to limit the run-time of step 3 to that of a single primal or tangent solution. The tangent solver must be able to compute homogeneous or inhomogeneous solutions for a given initial condition and output the final tangent solution in a given segment. The time dilation operator Pti in steps 4 to 6 requires the vector f (u(t); s) at the end of a given segment. Finally, computing the gradient requires either computing giT and h on the fly or access to the tangent solution for t0 < tK and f (u(ti ); s) for 0 ≤ i ≤ K. It should be noted that this algorithm is very similar to the one proposed by Benettin et al. for computing Lyapunov exponents [21]: Lyapunov Exponent Computation Inputs: Initial condition for the governing equations u0 , Spin-up time t0 , Specified time horizon and checkpoints t0 , t1 , ..., tK ¯ Ouputs: Sensitivities dJ/ds 1. 2. 3. 4. 5. 6. 7.

Compute u(t) by solving equation (1) until t = t0 . Set W1 (t0 ) = Q0 , where Q0 is some n × p unitary, orthonormal matrix. Set i = 1. Solve for u(t) and Wi (t) from ti−1 to ti using equations (1) and (29), respectively. Compute Qi Ri = Wi (ti ) with a QR-decomposition. Save Ri . Set Wi+1 (ti ) = Qi Set i = i + 1 and repeat steps 3 to 5. If i = K, proceed to the next step Compute Lyapunov exponents Λ j with the following expression [21] Λj ≈

K 1  [Ri ]j,j ΔT i=1

where [Ri ]j,j is the jth term on the diagonal of Ri .

9

(36)

The above algorithm lacks the projection operator Pt and does not require the computation of vi∗ (t), but both NILSS and Benettin’s approach require computing p homogeneous tangents and a QR-decomposition. Like Benettin’s approach the cost of NILSS is highly dependent on the number of homogeneous tangents, p that are required. For NILSS, recall that p should be at least the number of positive Lyapunov exponents for a given system. Unfortunately there is no existing method to estimate p a priori in a computationally efficient manner to the author’s knowledge. Some studies have estimated the Lyapunov spectrum by extrapolating from a partial spectrum, but it is unlikely this yields a very accurate estimate of p in all cases [37, 39]. Benettin’s approach can be used, but it costs more than NILSS since it must be used to compute all positive and some negative exponents to confirm what p should be. More computations of Lyapunov exponents may yield some rough empirical rules, but precious few of these calculations have been done to date [ 37, 38, 39]. Until these rules are found, the total cost of using NILSS will be very high in part because of the need to estimate p a priori with Benettin’s approach. 3.2. Adjoint NILSS To derive a discrete adjoint NILSS, first rewrite equation (32) in the notation of equation (23) as    α ∂ J¯ dJ¯  T 0 +h+ = g . β ds ∂s   Next, add the product of the residual of equation (22) and the length 2Kp vector φT ψ T .      −I B T  α   0   T   α ∂ J¯ dJ¯  T T 0 − φ + +h+ = g (37) ψ β β b B ds ∂s This can be rearranged as follows dJ¯  T = α ds The sensitivity

dJ¯ ds

β

T



  −

−I

BT



B

φ ψ



 +

g 0

 − bT ψ + h +

∂ J¯ ∂s

is not a function of α or β if the following adjoint equation is satisfied   −I B T  φ   g  = ψ 0 B

(38)

(39)

As for tangent NILSS, the Schur complement can be formed and solved instead of (39) A = BB T

Aψ = Bg,

(40)

Now, equation (32) can be rewritten as ∂ J¯ dJ¯ = −bT ψ + h + (41) ds ∂s Recall from equations (24) and (34) that b and h both depend on v ∗ (t), the inhomogeneous tangent solution. Since v ∗ (t) depends on the parameter s, it must be removed from equation (41) to obtain an adjoint NILSS algorithm that costs less than tangent NILSS for multiple parameter sensitivities to a single ¯ objective function J. This can be done by finding the discrete adjoint of v ∗ (t), as shown in Appendix C. When the discrete adjoint of v ∗ (t) is used, equation (41) becomes the conventional adjoint sensitivity equation dJ¯ = ds



tK

t0

T ∂ J¯ ∂f  w(t) ˆ dt + ,  ∂s t ∂s

where the adjoint w(t) ˆ is governed by the following adjoint equations 10

(42)

T dw ˆ df  =  w ˆ+ dt du t T dw ˆ df  = − w ˆ+ dt du t −

 1 ∂J  , ΔT ∂u t  1 ∂J  , ΔT ∂u  t

w(t ˆ K ) = xK − PtK QK ψK ,

tK−1 < t < tK

w(t ˆ − ˆ + i ) = Pti Pi w(t i ) + xi − Pti Qi ψi ,

ti−1 < t < ti ,

(43) i = 1, 2, ..., K − 1,

where the matrix Pi and the vector xi defined as Pi = I − Qi QTi  f (u(ti ); s) 1 ¯ J − J(u(ti ); s) xi = ΔT f (u(ti ); s)22 3.2.1. Algorithm Outline Adjoint NILSS is implemented as follows Adjoint NILSS Solver Inputs: Initial condition for the governing equations u0 , Spin-up time t0 , Specified time horizon and checkpoints t0 , t1 , ..., tK ¯ Ouputs: Sensitivities dJ/ds 1. 2. 3. 4.

Compute u(t) by solving equation (1) until t = t0 . Set W1 (t0 ) = Q0 , where Q0 is some n × p unitary, orthonormal matrix. Set i = 1. Solve for u(t) and Wi (t) from ti−1 to ti using equations (1) and (29) respectively. Compute and save T  ti 1 ∂J  T Wi (t) dt gi,1 = ΔT ∂u  ti−1

5. 6. 7. 8.

t

Pti Wi (ti )

Compute Qi Ri = Wi (ti ) = with a QR-decomposition. Save Qi , Ri , and Wi (ti ).  Set Wi+1 (ti ) = Qi . Set i = i + 1 and repeat steps 3 to 6. If i = K, proceed to the next step For all checkpoints 1 to K, compute and save T giT = gi,1 +

 1 ¯ J − J(u(ti ); s) Fti Wi (ti ) ΔT

9. Form the adjoint KKT Schur complement (40) and solve it. 10. Solve the adjoint equation (43) from segment K back to the beginning of segment 1. 11. Compute sensitivities using equation (42).

Note that steps 1-7 of the adjoint NILSS algorithm are very similar to the tangent algorithm, the main difference being that no inhomogeneous tangent solution is computed for adjoint NILSS. Therefore, adjoint NILSS requires a primal solver for u(t), a homogeneous tangent solver that solves from a given initial condition to compute Wi (t), and the ability to compute f (u(ti ); s) for 0 ≤ i ≤ K. Additionally, adjoint NILSS requires an adjoint solver to compute w(t) ˆ in each segment from the terminal conditions given by equation (43). This adjoint solver should be the discrete adjoint of the tangent solver used, otherwise the Adjoint NILSS algorithm presented in this paper will not be the discretely consistent with the tangent NILSS algorithm and there is no guarantee that the algorithm will converge to the correct sensitivity. 11

Figure 3: LEFT: z¯ vs. ρ for two different time horizons ΔT with sensitivities computed using NILSS indicated by the slopes of the lines. RIGHT: d¯ z /dρ vs. ρ for the same case shown in the left plot. Both time horizons used 1.0 time unit segments. As the time horizon is increased, the sensitivities converge to the correct value of around 1.0.

4. Results The following section demonstrates adjoint NILSS on several chaotic dynamical systems. 4.1. Lorenz 63 System The first system considered is the Lorenz 63 system [40], a three equation model of Rayleigh-B´enard convection between a hot and cold plate: dx dy dz = σ(y − x), = xρ − xz − y, = xy − β z (44) dt dt dt where x, y, and z are state variables and the baseline parameter values are σ = 10.0, ρ = 28.0, and β = 8/3. The Lorenz system only has one positive Lyapunov exponent, so p = 1 for NILSS [41]. NILSS was also conducted for p = 2 and p = 3 and gave very similar results for p = 1, so they are omitted from this paper. The objective function for this case is time-averaged z z¯ =

1 ΔT



T1

z(t) dt.

(45)

T0

Sensitivities with respect to the parameter ρ were computed as in [13]. Equation (44) was solved from the initial condition (x, y, z) = (1.0, 1.0, 28.0) with forward Euler time integration and a time step of 0.001 units. The beginning of the time horizon, T0 , is 50.0 time units for all computations presented in this section. Figure 3 shows sensitivities of z¯ to ρ for two time horizons. The sensitivities are correctly computed as roughly 1.0, as in previous studies of Lorenz 63 [13, 18]. Additionally, the sensitivities computed with adjoint NILSS are identical those computed using tangent NILSS to machine precision for a sufficiently small time segment size, as seen in table 1. Table 1 also shows that using fewer time segments for a given time horizon can cause the difference in the sensitivities to grow and eventually become very inaccurate. 4.2. Kuramoto-Sivashinsky Equation The K-S equation is a 4th order chaotic PDE, that can be used to model a number of physical phenomena [42]. Numerical studies typically use the 1D version of the K-S equation: ∂u ∂ 2 u ∂ 4 u ∂u = −(u + c) − − ∂t ∂x ∂x2 ∂x4 12

(46)

ΔT 5.0 50.0 50.0 50.0 50.0

K 5 50 5 2 1

Tangent 0.972897207597036240 0.939491532496952719 0.928118558969352248 1.071439409240778184 -679.568942847139055630

Adjoint 0.972897207597036795 0.939491532496957937 0.928118558971240404 1.071437449373630502 253.434337897725782796

Difference 5.551e-16 5.218e-15 1.888e-12 1.960e-06 9.330e+02

d¯ z Table 1: Tangent and adjoint sensitivities dρ for the Lorenz 63 system with ρ = 28.0 for various time horizons ΔT and number of time segments K. The sensitivities diverge because a discrete Adjoint will only be consistent with the tangent to machine precision, which is a finite number of digits. Therefore, as the tangent and adjoint solutions become very large, they will match to fewer and fewer decimal places, as seen in the above table.

The c term is added to make the system ergodic [34]. The equation was solved on the domain 0 ≤ x ≤ 128, with the boundary conditions:   ∂u   = =0 u ∂x  x=0,128 x=0,128

These boundary conditions make the K-S equation ergodic, unlike the periodic boundary conditions typically used by other studies. Figure 4 shows a typical solution of the K-S equation from a randomized initial condition. The initial condition is randomized as follows: since equation (46) is discretized with a 2nd order central difference scheme, u(x, 0) at each node is set to a random number drawn from a uniform distribution between -1.0 and 1.0. More detailed discussion of the numerics can be found in [34]. The objective function J¯ is u averaged over both space and time: J¯ =



1 128(T1 − T0 )

T1



128

u(x, t) dx dt T0

(47)

0

all solutions presented in this section had a run-up time T0 = 1000.0 units. Figure 4 shows a NILSS adjoint field computed for c = 0.0 with K = 50, ΔT = 500.0, and p = 15, the number of non-negative Lyapunov exponents [34]. The adjoint is continuous across time segments and does not grow exponentially. Note that if p is less than the number of positive Lyapunov exponents, the NILSS adjoint norm grows exponentially backwards in time, as shown in figure 5. Otherwise, if p is sufficiently large the adjoint stays bounded and can even decrease backwards in time. In this case the adjoint can be used for a number of analyses. Figures 6 and 7 show sensitivities of J¯ computed with adjoint NILSS. Figure 6 shows that inaccurate sensitivities are obtained if the number of modes p is smaller than the number of positive Lyapunov exponents. Figure shows that if a sufficient number of modes is used, then the sensitivities computed by NILSS will converge as the time horizon ΔT used to compute the objective function is increased. 4.3. Minimal Flow Unit The final application of adjoint NILSS presented in this paper is the minimal flow unit. The minimal flow unit is a channel flow with the smallest possible domain that can sustain turbulence [ 19]. It allows for a less computationally expensive study of near-wall turbulence, as it only captures the near-wall flow structures essential for turbulence. The minimal flow unit has flow statistics that match those of larger channels very well up to around y + = 40 wall units normal to the channel walls. The case considered in this paper is Reynolds number Re = 3000 flow. As in [19], the Reynolds number is defined as Uh ν where ν is the kinematic viscosity, h is half the channel height, and U is the centerline velocity of a laminar parabolic profile with the same volume flux as the turbulent channel. The domain size is π × 2 × 0.34π in the streamwise, wall-normal, and spanwise directions, respectively. Re =

13

(a) Primal solution u(x, t)

(b) NILSS adjoint solution w(x, ˆ t) computed with p = 15 Figure 4: Space-time plots of the primal solution and corresponding NILSS adjoint solution for the K-S equation with c = 0.0 and a random initial condition u(x, 0). The objective function J¯ is space and time averaged u(x, t) as in equation (47).

14

Figure 5: L2 norm of the NILSS adjoint vs. p for the KS-equation with c = 0.0. There are 14 positive Lyapunov exponents for this case. The objective function J¯ is space and time averaged u(x, t) as in equation (47).

¯ Figure 6: Sensitivity magnitude |dJ/dc| versus number of modes p for the K-S equation with c = 0 and J¯ is defined as in equation (47). Adjoint NILSS was applied on a ΔT = 800 time horizon with K = 80 time segments. The solid line indicates the ¯ value of |dJ/dc| computed with LSS by Blonigan and Wang[34]. The sensitivity magnitudes computed by NILSS are inaccurate for p < 14 but match the LSS results well for p ≥ 14.

15

Figure 7: Sensitivities computed for the K-S equation by adjoint NILSS with p = 15 modes. LEFT: J¯ vs. c for two different ¯ time horizons ΔT with sensitivities computed using NILSS indicated by the slopes of the lines. RIGHT: dJ/dc vs. c for the same case shown in the left plot. Both time horizons used 10.0 time unit segments. The objective function J¯ is defined as in equation (47). As the time horizon is increased, the sensitivities converge to roughly -1.0, the same value found in previous studies of the K-S equation using LSS and MSS [34, 15].

To ensure that the flow can sustain turbulence and to obtain grid resolution similar to that in the original study, the friction Reynolds number is set to Reτ = 140.0, resulting a channel width of 140 wall units. The friction Reynolds number is Reτ =

uτ δ ν

(48)

 where uτ = τw /ρ is friction velocity and δ is the channel half-width, which is 1.0 in this study. The minimal flow unit is simulated with a space-time discontinuous Galerkin spectral-element (DGSEM) solver [43]. The domain is discretized with a 4 × 16 × 2 mesh with 8th order spatial elements, resulting in a 32×128×16 distribution of nodes, similar to the 32 ×129×16 grid used in the original study [19]. The choice of Reτ = 140.0 results in grid resolutions of Δx+ ≈ 14 and Δz + ≈ 10 per node. The wall-normal spacing corresponds to Δy + ≈ 0.8 for the nodes closest to the walls, which is well resolved. A nearly incompressible mean-flow Mach number of M = 0.3 was chosen for this study. The space-time elements used for the minimal flow unit are 4th order in time and the time slabs are Δt+ = 0.0023 eddy turnover time units, where uτ t (49) δ Figure 8 shows that this time slab size results in good agreement between the flow computed by the DGSEM solver and the original study. NILSS was implemented using the discrete tangent and adjoint solvers in the DGSEM solver [45]. Time segments of 50 time slabs, or 0.115 eddy turnover time units, were used to keep the tangent solutions from growing too large in magnitude. Before running adjoint NILSS, the Lyapunov exponents for the flow unit were computed using the approach of Benettin et al. [21]. As shown in figure 9, there are roughly 150-160 positive exponents for the minimal flow unit. The objective function considered in this study is the kinetic energy of the flow integrated over the entire volume t+ =

 J¯ =

T1

T0

 V

1 ρ(u2 + v 2 + w2 ) dx dy dz dt 2

(50)

where V is the simulation domain, ρ is fluid density, and u, v, w are the x, y, and z velocities, respectively. A run-up time of T0 = 735.0 eddy turnover time units was used for the computations presented in this section. 16

(a) Shear velocity u+ vs. wall units y +

(b) u u and v  v  vs. y +

Figure 8: Flow statistics for the minimal flow unit as computed in the current study with the DGSEM solver, the original study [19], and data from Wei and Willmarth’s experimental study of Re ≈ 3850 channel flow [44].

Figure 9: Lyapunov exponents for the mimimal flow unit computed using the method of Benettin et al. [ 21]. Computations used 500 time segments of 0.115 eddy turnover time units .

17

(a) Sensitivity of Volume integrated kinetic energy to Reτ vs. time horizon size. Minimum and maxium sensitivites after 10 eddy turnover units are indicated by the red and blue dotted line, respectively.

(b) Time-averaged, volume-intergrated kinetic energy vs. ¯ Reτ . The sensitivity dJ/dRe τ is indicated by the slope of the line. The red and blue dotted lines show the slopes corresponding to the red and blue lines in figure 10(a)

Figure 10: NILSS sensitivity computations for the minimal flow unit. Time averaged objectives in figure 10(b) were computed with time horizons of ΔT = 1600.0 eddy turnover time units. Error bars were computed using the batch mean approach with 8 means [46].

Adjoint NILSS was computed with p = 160 modes over K = 400 time segments. Table 2 shows computation times for the NILSS implementation used. Note that the cost is dominated by the primal, tangent, and adjoint solutions on each segment, which take 92% of the total wall time. The second largest contribution to the total wall time is the input and output of the data needed for NILSS, which took about 6% of the total time. The first four functions in table 2 correspond to steps 2-6 of Benettin et al.’s algorithm presented in section 3.1.3, so computing the Lyapunov exponent algorithm for 160 exponents over 400 time segments requires around 86% of the wall time required by adjoint NILSS. The sensitivity computed using NILSS shown in figure 10(a) has not converged for this time horizon due to the presence of a very long time scale in this flow [19]. However, the sensitivity is the correct order of magnitude, as shown in figure 10(b). As shown in figure 11(a), the adjoint remains bounded for the time horizon considered in this study, showing that p = 160 is a sufficient number of modes for this flow. Interestingly, the spikes in the adjoint occur prior to rapid increases in friction on one or both walls, for instance the spike at t+ ≈ 12 lines up with a rapid rise in friction on both walls shown in figure 11(b). The spike in the magnitude of the adjoint at t+ ≈ 12 corresponds to the flow field and adjoint shown in figures 12(a) and 12(b). The flow on the upper wall is nearly laminar when the adjoint magnitude spikes. The spike occurs just prior to the flow on the upper wall experiencing a bursting event and becoming turbulent again in figure 12(c). Evidence of the burst can also be seen from the rapid increase in wall shear stress on the top wall that occurs after t+ = 12 in figure 11(b). The increase in adjoint magnitude prior to bursts shows that the adjoint is revealing when the flow is most susceptible to instabilities. This is because flow instabilities ¯ the objective function the adjoint is computed will increase the turbulent kinetic energy and therefore J, for. This means that the wavy adjoint structure shown in figure 12(b) gives a pattern of perturbations to x-momentum to excite a flow instability at time t+ = 12. Similar structures in the adjoint field are observed at the same times as other adjoint spikes. Figure 13 shows a cut through several of these structures in the adjoint field superimposed on the x-momentum field, two on the top wall and one on the bottom, wall. These adjoint structures exist in conjunction with ejection events, where lower momentum fluid is transported away from the wall [ 47]. The adjoint also shows perturbations to other conserved quanities. Figure 14 shows a snapshot of the 18

Function Solve primal Solve tangent QR decomposition Write Qi ,Ri , and tangent solutions Wi Read data for KKT system Construct and solve adjoint KKT System Read data for adjoint terminal conditions Compute adjoint terminal condition Solve adjoint and compute sensitivities

Algorithm Step 3 3 5 5 8,9 9 10 10 10,11

Normalized wall time 1.0 3.5 0.069 0.27 0.069 0.46 0.069 0.0007 0.69

Number of calls 400 400 400 400 1 1 400 400 400

% total wall time 17.9 62.5 1.23 4.82 0.00308 0.0205 1.23 0.0125 12.3

Table 2: Wall times for the NILSS implementation used for the minimal flow unit. The wall time required by each part of the NILSS algorithm is listed in order of its execution along with the corresponding step numbers from the algorithm in section 3.2.1. All wall times were normalized with the average wall time of a primal solution on a single time segment. Wall time estimates were computed for the case where all p = 160 tangent solutions were computed simultaneously for each segment. Note that each unsteady tangent was computed with a quarter of the number of processors used for the primal and adjoint solutions, resulting in a longer wall time than the primal or adjoint. Finally, the contribution of each time segment to the ¯ τ was computed as part of the adjoint run. sensitivity J/

(a) L2 norm of the NILSS adjoint field vs. time.

(b) Normalized integrated wall friction τw /τw,ref for the upper and lower walls vs. time.

Figure 11: NILSS adjoint norm and wall friction for the minimal flow unit solution. The reference friction τw,ref is defined as τw,ref = u2τ where uτ is computed from the prescribed Reτ using equation (48).

19

(a) Primal at t+ = 12.50.

(b) Adjoint at t+ = 12.50.

(c) Primal at t+ = 13.18.

(d) Adjoint at t+ = 13.18.

Figure 12: Snapshots of the primal solution and the NILSS adjoint field. The primal plots show Q-criteria isocontours at |Q| = 0.02, colored by x-momentum magnitude, with high magnitudes in red, medium magnitudes in green, and low magnitude in blue. The adjoint plots show isosurfaces of adjoint x-momentum for values of ±2.0 where red is positive. The walls in all plots are colored by local wall friction, with high magnitudes in red, medium magnitudes in green, and low magnitude in blue. x,y, and z are the streamwise, wall normal, and spanwise directions, respectively.

20

Figure 13: Snapshot of primal x-momentum with adjoint x-momentum contour lines at t+ = 12.60. The two-dimensional cut was taken at z = 0.85 as shown on the right. High magnitudes of x-momentum are shown in red, medium magnitudes in green, and low magnitude in blue. The contour lines of adjoint x-momentum are for values of ±2.0 where red is positive.

z-momentum adjoint field. It shows that pertubations to z-momentum at both edges of a low velocity streak will have the most impact on the flow’s kinetic energy by causing or supressing the burst at t+ = 12 − 13. This is consistent with the results of Farano et al., who used a linearized analysis to show that the optimal perturbations to increase kinetic energy and cause bursting events in a turbulent channel flow lie on low velocity streaks [48]. 5. Conclusions This paper presents a discrete adjoint version of NILSS for conducting sensitivity analysis of chaotic dynamical systems. It is shown that the adjoint field computed with adjoint NILSS will remain bounded if the number of modes p is at least the number of positive exponents. Additionally, adjoint NILSS will compute similar sensitivities to previous studies using LSS and MSS for the Lorenz 63 system and the K-S equation. Finally, adjoint NILSS was used to study the minimal flow unit, producing the first shadowingbased adjoint field for a 3D turbulent flow. The adjoint was used to compute a fairly accurate sensitivity and was shown to provide some physical insight into the flow itself. The insights provided by the adjoint were consistent with the results of previous studies of the stability of wall-bounded flows. These results show the potential of the unsteady adjoint to be used as a tool to increase physical understanding of flow. Future work includes further study of the adjoint field of the minimal flow unit and its physical interpretation. The cost of the NILSS scales with the number of modes p, which is already O(100) for the small turbulent flow simulation presented in this paper. Although the minimization problem (39) took around a minute of wall time to solve on a single core, the tangent and adjoint solutions required many CPU hours, as indicated by the wall times in table 2. Therefore, an important direction of future research is on approaches to reduce the cost of NILSS and shadowing-based sensitivity analysis in general. The NILSS implementation used in this study could be made faster by computing the primal and tangent solutions in a staggered fashion: compute step i + 1 of the primal while step i of the tangent solution(s) are computed. Staggering reading and writing to hard disk would also improve performance. Parallel reading and writing along with a parallel QR decomposition could allow NILSS to scale to large problems than those presented in this study. Ultimately, most of the computational cost of NILSS comes from computing the unsteady tangent solutions. Therefore, large reductions in the cost of NILSS will most likely be achieved through large reductions in the cost and number of unsteady tangent solutions. Additionally, efficient methods to estimate the number of positive Lyapunov exponents a priori are required to drive down the overall cost of NILSS for large chaotic systems. Cost reductions of shadowing-based sensitivity analysis could also be found by exploring new implementations of MSS. Although MSS involves solving a larger system KKT system than NILSS, it does not 21

Figure 14: Snapshot of primal x-velocity with adjoint Z-momentum isocontours at t+ = 12.60. LEFT: same view as figures 12(a)-12(d), RIGHT: View from y = 0 in the positive y direction, the flow goes from top to bottom. The two-dimensional cut was taken at a distance of y + = 20 from the top wall. High magnitudes of x-velocity are shown in red, medium magnitudes in green, and low magnitude in blue. The contour lines of adjoint z-momentum are for values of ±15.0 where red is positive.

require any knowledge of the number of positive Lyapunov exponents and could be made less expensive with better preconditioning. Further reductions in the computational cost of NILSS or MSS would enable adjoint analyses of ever larger simulations, and potentially pave the way for efficient optimization, error estimation, and uncertainty quantification of these simulations. Appendix A. Tangent NILSS In terms of the decomposition in equation (17), the MSS minimization statement becomes 1 ∗ + + T ∗ [v (t+ ) + Wi+1 (t+ i )αi+1 ] [vi+1 (ti ) + Wi+1 (ti )αi+1 ] 2 i=0 i+1 i K

min αi

(A.1)

To find the corresponding KKT system, we consider the following Lagrangian 1 ∗ + + T ∗ T T ∗ − [v (t+ ) + W (t+ i )αi+1 ] [vi+1 (ti ) + Wi+1 (ti )αi+1 ] + βi (αi+1 − Ri αi − Qi vi (ti )) 2 i=0 i+1 i K

Λ = min αi

where β0 = 0. This can be written as follows

Λ = min αi

K  1

2 i=0

1 T + + T + + T + ∗ T ∗ ∗ (t+ [vi+1 i )] vi+1 (ti ) + [vi+1 (ti )] Wi+1 (ti )αi+1 + αi+1 Wi+1 (ti ) Wi+1 (ti )αi+1 2

+ βiT (αi+1 − Ri αi − QTi vi∗ (t− i )) + T ∗ Since Wi+1 (t+ i ) = Qi and Qi vi+1 (ti ) = 0

22

Λ = min αi

K  1 i=0

2

1 T + ∗ T ∗ T T ∗ − (t+ [vi+1 i )] vi+1 (ti ) + αi+1 αi+1 + βi (αi+1 − Ri αi − Qi vi (ti )) 2

∂Λ/∂βi = 0 is equation (20), and ∂Λ T = 0 = αiT + βi−1 − βiT Ri ∂αi

(A.2)

Appendix B. NILSS and MSS To show how NILSS is related to MSS, we need to consider the form of the tangent solution   t   t  ti ,t τ,t ∂f  ηi (τ ) dτ f (u(t); s) + φ v(t) = φ vi +  dτ , ti ≤ t < ti+1 , ∂s ti ti τ

(B.1)



where vi = v(ti ) and φt,t is the tangent propagator from t to t , defined as

φt,t = I for all t ,





φt ,τ · φt,t = φt,τ ,

 d t,τ ∂f  t,τ φ = φ , dτ ∂u τ

T d ∗ t,τ ∂f  φ = −  φ∗ t,τ , dt ∂u t

(B.2)

Note that φ∗ t,τ is the adjoint propagator from τ back to t. It is used to write the adjoint solution in Appendix C. The projection operator Pt can be used to eliminate the second term in equation (B.1)   t ∂f  φτ,t  dτ , ti ≤ t < ti+1 , (B.3) v(t) = Pt v  (t) = Pt φti ,t vi + Pt ∂s τ ti The solutions of Wi (t) and vi∗ (t) can also be written in terms of the tangent propagator and Pt . Wi (t) = Pt Wi (t) = Pt φti ,t Qi vi∗ (t)

=

∗ Pt v  i (t)

=

∗ Pt φti ,t v  i (t+ i )



t

+ Pt

φ ti

ti ≤ t < ti+1 ,    dτ , ∂s 

τ,t ∂f

(B.4)

ti ≤ t < ti+1 ,

(B.5)

τ

Consider the case where we choose p = n and set all Qi = I, instead of computing Qi with a QRdecomposition. In this case, equation (B.4) becomes ti ≤ t < ti+1 ,

Wi (t) = Pt φti ,t and therefore

Wi (ti ) = Ri = Pti φti−1 ,ti = Φi This means that the KKT system becomes ⎡ −Φ2 Φ1 ΦT1 + I T ⎢ −Φ Φ ΦT2 + I −Φ3 2 2 ⎢ ⎢ .. .. ⎣ . .

⎤⎛ ⎥⎜ ⎥⎜ ⎥⎜ ⎦⎝

..

. −ΦTK

ΦK ΦTK + I

β1 β2 .. . βK

⎞

⎛

⎟ ⎜ ⎟ ⎜ ⎟=⎜ ⎠ ⎝

−v1∗ (t− 1) −v2∗ (t− 2) .. .

⎞ ⎟ ⎟ ⎟ ⎠

(B.6)

∗ −vK (t− K)

∗ (t+ This is the MSS KKT Schur complement system, since when Qi = I equation (19) simplifies to vi+1 i )=0 and

23

vi∗ (t− i )

 = Pti

ti

φ ti−1

t,ti

 ∂f  dt, ∂s t

which is equal to the ith block of the MSS KKT Schur complement’s right hand side. Appendix C. Adjoint NILSS This appendix shows how equation (41), ∂ J¯ dJ¯ = −bT ψ + h + , ds ∂s can be made independent of v ∗ (t). Appendix C.1 shows how this is done for the first term, −bT ψ, Appendix C.2 shows how this is done for the second term, h, and Appendix C.3 combines the results of the previous two sections to obtain equation (43). All of these derivations use the tangent propagator notation introduced in Appendix B to describe the adjoint solution. For an adjoint equation on the ith segment T  df  1 ∂J  dw ˆ = w ˆ+ , − dt du t ΔT ∂u t

ˆ i, w(t ˆ i) = w

ti−1 < t < ti ,

(C.1)

 ∂J  dτ ∂u τ

ti−1 ≤ t ≤ ti

(C.2)

the solution is w(t) ˆ = φ∗

t,ti

ˆi + w

1 ΔT



ti

φ∗

t

t,τ

These identities can be used to eliminate dependence on v ∗ (t) by finding an equivalent expression for equation (41) using an adjoint w(t). ˆ Appendix C.1. Computing Direct Sensitivities The first term on the right of equation (41), −bT ψ, can also be written as −bT ψ = −

K 

T [vi∗ (t− i )] Qi ψi

i=1

The dependence on [v ∗ (t− i )] means that equation (41) still requires N solutions different parameters s, negating any benefit of using an adjoint approach. To ∗ T v ∗ (t) an expression for v  (t− i ) is needed. Defining Pi = I − Qi Qi , and using Appendix B,   ti  − ∗ − ti−1 ,ti ∗ t,ti ∂f  Pi−1 v i−1 (ti−1 ) + Pti φ (C.3) v i (ti ) = Pti φ  dt ∂s ti−1 t

Where Qi ψi is a n×1 vector. of v ∗ (t) for sensitivities to N eliminate the dependence on the tangent propagator from

Substituting equation (C.3) into itself ∗

∗

ti−1 ,ti Pi−1 Pti−1 φti−2 ,ti−1 Pi−2 v  (t− v  i (t− i ) =Pti φ i−2 )   ti−1 ∂f  φt,ti−1 dt + Pti φti−1 ,ti Pi−1 Pti−1 ∂s t ti−2   ti ∂f  φt,ti dt + Pti ∂s t ti−1

24

(C.4)

∗

Recursive substitutions into equation (C.4) can be made until a closed form expression for v  i (t− i ) if i > 1. ∗ Then, since v  0 (t− 0 ) = 0, ∗ v  i (t− i )

 = Pti

ti

φ

t,ti

ti−1

 j      i−1  ti−j   ∂f  ti−k ,ti−k+1 t,tj−1 ∂f  dt + Pti−k+1 φ Pi−k Pti−j φ  dt ∂s t ∂s ti−j−1 t j=1

(C.5)

k=1

∗

To find the transpose of v  (t− i ), care needs to be taken with the indicies used in the product operator in equation (C.5)

T [vi∗ (t− i )]



ti

= ti−1

⎛ ⎡ ⎤ ⎞ T T  tj i−1 i−1    ∗ t,t ∂f  ∗ t,ti ∂f ∗ t ,t  φ j−1 Pt ⎣ ⎝ φ Pti dt + Pk φ k k+1 Ptk+1 ⎦ dt⎠ j  ∂s t tj−1 ∂s t j=1

(C.6)

k=j

T Equations (C.6) and (C.2) can be used to rewrite −[v ∗ (t− i )] Qi ψi as    T T T  ti i−1  tj ti  ∂f   ∂f   ∂f   ∗ − T w ˆ (t) dt + w ˆ (t) dt = ˆ (t) dt −[vi (ti )] Qi ψi =    w ti−1 ∂s t tj−1 ∂s t t0 ∂s t j=1

where w ˆ  (t) = −φ∗ t,ti Pti Qi ψi , ti−1 < t < ti , ⎡ ⎤ i−1  Pk φ∗ tk ,tk+1 Ptk+1 ⎦ Qi ψi , w ˆ  (t) = −φ∗ t,tj−1 Ptj ⎣

tj−1 < t < tj

j = 1, 2, ..., i − 1,

k=j

or, written as differential equations T dw ˆ  df   =  w ˆ, dt du t T dw ˆ  df   − ˆ, =  w dt du t

−

w ˆ  (ti ) = −Pti Qi ψi ,

ti−1 < t < ti

w ˆ  (t− ˆ  (t+ j ) = Ptj Pj w j ),

tj−1 < t < tj

j = 1, 2, ..., i − 1,

∗

This is an adjoint equation for the sensitivity contribution of v  (t− i ). Since adjoint equations are linear, they can be superimposed. Therefore, the K adjoint solutions w ˆ  (t) can be added to obtain −

K 

T [vi∗ (t− i )] Qi ψi



tK

= t0

i=1

T ∂f  (0) w ˆ (t) dt ∂s t

(C.7)

where w ˆ (0) (t) is the solution of the following adjoint equation T dw ˆ (0) df  (0) − ˆ , =  w dt du t  T dw ˆ (0) df  (0) =  w − ˆ , dt du t

w ˆ (0) (tK ) = −PtK QK ψK ,

tK−1 < t < tK

w ˆ (0) (t− ˆ (0) (t+ i ) = Pti Pi w i ) − Pti Qi ψi ,

25

ti−1 < t < ti

i = 1, 2, ..., K − 1,

Appendix C.2. Computing Indirect Sensitivities Note that the second term in equation (41), h, also depends on v ∗ (t). Recall equation (34)   T K ti   1  ∂J   ∗ ∗ ¯ h=  v i (t) dt + J − J(u(ti ); s) Fti v i (ti ) ΔT i=1 ti−1 ∂u t To eliminate the dependence of h on v ∗ (t), start with an equation for v ∗ (t) using the tangent propagator   t  + ∗ ti ,t  ∗ τ,t ∂f  φ (C.8) vi+1 (t) = φ v i+1 (ti ) +  dτ ∂s ti τ Equation (C.8) can be substituted into the first term of equation (34) to obtain 1 ΔT



ti

ti−1

 T T   ti  T  t 1 ∂J  ti−1 ,t  ∗ + ∂J   ∗ ∂J  τ,t ∂f  v (t) dt = φ v i (ti−1 ) +  φ ∂s  dτ dt ∂u t i ΔT ti−1 ∂u t ti−1 ∂u t τ T T   ti  ti  ti 1 ∂J  ti−1 ,t  ∗ + ∂J  τ,t ∂f  1 φ v (t ) dt + φ dt dτ = i i−1 ΔT ti−1 ∂u t ΔT ti−1 τ ∂u t ∂s τ  T  ti   ti  ti   1 ∂f  ∗ T 1 ∗ ti−1 ,t ∂J  ∗ τ,t ∂J  )] φ dt + φ dt dτ = [v  i (t+ i−1   ΔT ti−1 ∂u t ΔT ti−1 ∂s τ τ ∂u t (C.9) ∗

∗

−  using equation (C.6) and the definition v  (t+ i ) = Pi v (ti ), equation (C.9) can be rewritten as



1 ΔT

ti

ti−1

  T T   ti ti−1  1 ∂J   ∗ ∂f  ∗ t,ti−1 ∗ ti−1 ,t ∂J  v (t) dt = φ P dt P φ dt t i−1 i i−1   ∂u t ΔT ∂u t ti−2 ∂s t ti−1 ⎛ ⎡ ⎤ ⎞⎤ ⎡ T   ti  i−2 i−2   1 ⎣ ⎝ tj ∂f  ∗ t,tj−1 ∗ tk ,tk+1 ∗ ti−1 ,t ∂J  ⎣ ⎦ ⎠ ⎦ Ptj Pk φ Ptk+1 dt φ + Pi−1  φ  dt ΔT j=1 ∂s ∂u tj−1 ti−1 t t k=j T  ti   ti 1 ∂f  ∂J  + φ∗ τ,t dt dτ  ΔT ti−1 ∂s τ τ ∂u t

using equation (C.2), this expression becomes

1 ΔT



ti

ti−1

   T T T T  ti−1 i−2  tj ti  ∂J   ∗ ∂f   ∂f   ∂f   v (t) dt = ˆ (t) dt + ˆ (t) dt + ˆ (t) dt  w  w  w ∂u t i ti−2 ∂s t tj−1 ∂s t ti−1 ∂s t j=1 T  ti ∂f   = ˆ (t) dt  w t0 ∂s t

where 1 w ˆ (t) = ΔT 



w ˆ (t) = φ



ti

φ t

∗ t,tj−1

∗ t,τ

⎡ Ptk ⎣

 ∂J  dτ, ∂u τ

i−2 

k=j

ti−1 < t < ti , ⎤

∗ tk ,tk+1

Pk φ

Ptk+1 ⎦ Pi−1



ti

φ ti−1

26

   dt, ∂u 

∗ ti−1 ,t ∂J

t

tj−1 < t < tj

j = 1, 2, ..., i − 1,

and w ˆ  (t) is the solution of the following adjoint equation T  1 ∂J  dw ˆ  df   ˆ + , w ˆ  (ti ) = 0, =  w − dt du t ΔT ∂u t T dw ˆ  df   − ˆ, w ˆ  (t− ˆ  (t+ =  w j ) = Ptj Pj w j ), dt du t

ti−1 < t < ti tj−1 < t < tj

j = 1, 2, ..., i − 1,

As in Appendix C.1, the i adjoint solutions w ˆ  (t) can be superimposed. Therefore, T T  tK K  ∂f  (1) 1  ti ∂J   ∗ v (t) dt = w ˆ (t) dt ΔT i=1 ti−1 ∂u t i ∂s t t0

(C.10)

where w ˆ (1) (t) is the solution of the following adjoint equation T dw ˆ (1) df  (1) ˆ + =  w − dt du t T dw ˆ (1) df  (1) − ˆ + =  w dt du t

 1 ∂J  , ΔT ∂u t  1 ∂J  , ΔT ∂u t

w ˆ (1) (tK ) = 0,

tK−1 < t < tK

w ˆ (1) (t− ˆ (1) (t+ i ) = Pti Pi w i ),

ti−1 < t < ti

i = 1, 2, ..., K − 1,

Next, consider the second term in the definition of h (equation (34))  1 ¯ ∗ ∗ J − J(u(ti ); s) Fti v  i (ti ) = xTi v  i (t− i ) ΔT ∗ Where xi is a n × 1 vector. As in section Appendix C.1, the term contains an inner product of v  i (t− i ) and some n × 1 vector. Therefore, the second term in h can be computed as follows T  tK K 1  T ∗ ∂f  (2) xi v (t) = w ˆ (t) dt ΔT i=1 ∂s t t0

(C.11)

where w ˆ (2) (t) is the solution of the following adjoint equation T dw ˆ (2) df  (2) ˆ , =  w dt du t T dw ˆ (2) df  (2) − ˆ , =  w dt du t −

w ˆ (2) (tK ) = xi ,

tK−1 < t < tK

w ˆ (2) (t− ˆ (2) (t+ i ) = Pti Pi w i ) + xi ,

ti−1 < t < ti

i = 1, 2, ..., K − 1,

Appendix C.3. Computing All Sensitivity components ¯ Finally, equations (C.7), (C.10), and (C.11) can be combined to form an expression for the sensitivity dJ/ds ∗ that is independent of v i (t) dJ¯ = ds



tK

t0

T ∂ J¯ ∂f  w(t) ˆ dt +  ∂s t ∂s

(C.12)

ˆ (1) + w ˆ (2) is the solution of the following adjoint equation where w ˆ=w ˆ (0) + w T dw ˆ df  − w ˆ+ = dt du t T dw ˆ df  ˆ+ =  w − dt du t

 1 ∂J  , ΔT ∂u t  1 ∂J  , ΔT ∂u t

w(t ˆ − K ) = xK − PtK QK ψK ,

tK−1 < t < tK

w(t ˆ − ˆ + i ) = Pti Pi w(t i ) + xi − Pti Qi ψi , 27

ti−1 < t < ti ,

(C.13) i = 1, 2, ..., K − 1,

where xi =

 f (u(ti ); s) 1 ¯ J − J(u(ti ); s) ΔT f (u(ti ); s)22

Acknowledgments This research was supported by an appointment to the NASA Postdoctoral Program at the NASA Ames Research Center, administered by Universities Space Research Association under contract with NASA. This research was sponsored by NASA‘s Transformational Tools and Technologies (TTT) Project of the Transformative Aeronautics Concepts Program under the Aeronautics Research Mission Directorate. I would like to thank Laslo Diosady, Corentin Carton de Wiart, Anirban Garai, Nicholas Burgess, and Scott Murman for their help with setting up the minimal flow unit case in the DGSEM solver. I would also like to thank Laslo Diosady and Marian Nemec for their helpful suggestions on this paper. Bibliography [1] A. Jameson, Aerodynamic design via control theory, Journal of Scientific Computing 3 (3) (1988) 233–260. [2] M. Giles, E. S¨ uli, Adjoint methods for pdes: a posteriori error analysis and postprocessing by duality, Acta Numerica 11 (2002) 145–236. [3] D. Venditti, D. Darmofal, Grid adaptation for functional outputs: Application to two-dimensional inviscid flow, Journal of Computational Physics 176 (2002) 40–69. [4] M. Gunzburger, Perspectives in Flow Control and Optimization, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2002. [5] Q. Wang, K. Duraisamy, J. Alonso, G. Iaccarino, Risk assessment of scramjet unstart using adjoint-based sampling methods, AIAA Journal 50 (3) (2012) 581–592. [6] D. Lea, M. Allen, T. Haine, Sensitivity analysis of the climate of a chaotic system, Tellus 52A (2000) 523–532. [7] G. Eyink, T. Haine, D. Lea, Ruelle’s linear response formula, ensemble adjoint schemes and L´ evy flights, Nonlinearity 17 (5) (2004) 1867–1889. [8] J. Thuburn, Climate sensitivities via a Fokker-Planck adjoint approach, Quarterly Journal of the Royal Meteorological Society 131 (605) (2005) 73–93. [9] P. Blonigan, Q. Wang, Probability density adjoint for sensitivity analysis of the mean of chaos, Journal of Computational Physics 270 (2014) 660–686. [10] C. Leith, Climate response and fluctuation dissipation, Journal of the Atmospheric Sciences 32 (10) (1975) 2022–2026. [11] R. Abramov, A. Majda, Blended response algorithms for linear fluctuation-dissipation for complex nonlinear dynamical systems, Nonlinearity 20 (12) (2007) 2793. URL http://stacks.iop.org/0951-7715/20/i=12/a=004 [12] R. Abramov, A. Majda, New approximations and tests of linear fluctuation-response for chaotic nonlinear forced-dissipative dynamical systems, Journal of Nonlinear Science 18 (2008) 303–341. [13] Q. Wang, R. Hui, P. Blonigan, Least squares shadowing sensitivity analysis of chaotic limit cycle oscillations, Journal of Computational Physics 267 (2014) 210–224. [14] Q. Wang, Convergence of the least squares shadowing method for computing derivative of ergodic averages, SIAM Journal of Numerical Analysis 52 (1) (2014) 156–170. [15] P. Blonigan, Q. Wang, Multiple shooting shadowing for sensitivity analysis of chaotic dynamical systems, journal of Computational Physics, (submitted). arXiv:1704.02047 (2017). [16] A. Ni, P. J. Blonigan, M. Chater, Q. Wang, Z. Zhang, Sensitivity analysis on chaotic dynamical system by non-intrusive least square shadowing (ni-lss), AIAA 2016-4399, 2016. [17] P. Blonigan, Q. Wang, E. Nielsen, B. Diskin, Least squares shadowing sensitivity analysis of chaotic flow around a twodimensional airfoil, AIAA 2017-1534, 2017. [18] A. Ni, Q. Wang, Sensitivity analysis on chaotic dynamical systems by non-intrusive least square shadowing (nilss), journal of Computational Physics, (submitted). arXiv:1611.00880v6 (2017). [19] J. Jim´ enez, P. Moin, The minimal flow unit in near-wall turbulence, Journal of Fluid Mechanics 225 (1991) 213–240. [20] J.-P. Eckmann, D. Ruelle, Ergodic theory of chaos and strange attractors, Reviews of Modern Physics 57 (3) (1985) 617–656. [21] G.Benettin, L. Galgani, A. Giorgilli, J. Strelcyn, Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. part 2: Numerical application, Meccanica 15 (1) (1980) 21–30. [22] F. Ginelli, P. Poggi, A. Turchi, H. Chat´e, R. Livi, A. Politi, Characterizing dynamics with covariant Lyapunov vectors, Physical Review Letters 99:130601. [23] S. Y. Pilyugin, Shadowing in dynamical systems, 1st Edition, Vol. 1706 of Lecture Notes in Mathematics, Springer-Verlag, New York, 1999.

28

[24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48]

D. Ruelle, Differentiation of SRB states, Communications in Mathematical Physics 187 (1997) 227–241. C. Bonatti, L. Diaz, M. Viana, Uniform Hyperbolicity: A Global Geometric and Probabilistic Perspective, Springer, 2005. G. Gallavotti, E. Cohen, Dynamical ensembles in stationary states, Journal of Statistical Physics 80 (1995) 931–970. G. Gallavotti, E. Cohen, Dynamical ensembles in nonequilibrium statistical mechanics, Physical Review Letters 74 (1995) 2694–2697. D. Albers, J. Sprott, Structural stability and hyperbolicity violation in high-dimensional dynamical systems, Nonlinearity 19 (2006) 1801–1847. E. J. Doedel, M. J. Friedman, Numerical computation of heteroclinic orbits, Journal of Computational and Applied Mathematics 26 (1-2) (1989) 155–170. J. Sanchez, M. Net, On the multiple shooting continuation of periodic orbits by Newton-Krylov methods, International Journal of Bifurcation and Chaos 20 (1) (2010) 43–61. J. D. Farmer, J. J. Sidorowich, Optimal shadowing and noise reduction, Physica D 47 (1991) 373–392. Q. Wang, Forward and adjoint sensitivity computation of chaotic dynamical systems, Journal of Computational Physics 235 (2013) 1–13. M. Chater, A. Ni, P. Blonigan, Q. Wang, Least squares shadowing method for sensitivity analysis of differential equations, arXiv:1509.02882v1 (2015). P. Blonigan, Q. Wang, Least squares shadowing sensitivity analysis of a modified Kuramoto–Sivashinsky equation, Chaos, Solitons, and Fractals 64 (2014) 16–25. S. Gomez, Parallel multigrid for large-scale least squares sensitivity, Master’s thesis, Massachusetts Institute of Technology, Cambridge, MA (2013). P. Blonigan, Q. Wang, E. Nielsen, B. Diskin, Least squares shadowing sensitivity analysis of chaotic flow around a twodimensional airfoil, AIAA 2016-0296, 2016. L. Sirovich, A. Deane, A computational study of Rayleigh-Benard convection. Part 2. Dimension considerations, Journal of Fluid Mechanics 222 (1991) 251–266. L. Keefe, P. Moin, J. Kim, The dimension of attractors underlying periodic turbulent poiseuille flow, Journal of Fluid Mechanics 242 (1992) 1–29. P. Blonigan, P. Fernandez, S. Murman, Q. Wang, G. Rigas, L. Magri, Toward a chaotic adjoint for les, in: Proceedings of the 2016 Summer Program, Center for Turbulence Research, Stanford, 2016, pp. 385–394. E. Lorenz, Deterministic nonperiodic flow, Journal of the Atmospheric Sciences 20 (1963) 130–141. S. Strogatz, Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering, Westview Press, Philadelphia, PA, 1994. J. M. Hyman, B. Nicolaenko, The Kuramoto-Sivashinsky equation: A bridge between PDE’s and dynamical systems, Physica D: Nonlinear Phenomena 18:1-3 (1986) 113–126. L. Diosady, S. Murman, DNS of flows over periodic hills using a discontinuous-Galerkin spectral-element method, AIAA paper # 2014-2784, 2014. T. Wei, W. Willmarth, Reynolds-number effects on the structure of a turbulent channel flow, Journal of Fluid Mechanics 204 (1989) 57–95. M. Ceze, L. Diosady, S. Murman, Development of a high-order space-time matrix-free adjoint solver, AIAA 2016-0833, 2016. G. S. Fishman, Grouping observations in digital simulation, Management Science 24 (5) (1978) 510–521. J. Jim´ enez, Cascades in wall-bounded turbulence, Annual Review of Fluid Mechanics 44 (2012) 27–45. M. Farano, S. Cherubini, J. C. Robinet, P. de Palma, Optimal bursts in turbulent channel flow, Journal of Fluid Mechanics 817 (2017) 35–60.

29