Identifying invariant solutions of wall-bounded three-dimensional shear flows using robust adjoint-based variational techniques

Abstract Invariant solutions of the Navier–Stokes equations play an important role in the spatiotemporally chaotic dynamics of turbulent shear flows. Despite the significance of these solutions, their identification remains a computational challenge, rendering many solutions inaccessible and thus hindering progress towards a dynamical description of turbulence in terms of invariant solutions. We compute equilibria of three-dimensional wall-bounded shear flows using an adjoint-based matrix-free variational approach. To address the challenge of computing pressure in the presence of solid walls, we develop a formulation that circumvents the explicit construction of pressure and instead employs the influence matrix method. Together with a data-driven convergence acceleration technique based on dynamic mode decomposition, this yields a practically feasible alternative to state-of-the-art Newton methods for converging equilibrium solutions. We compute multiple equilibria of plane Couette flow starting from inaccurate guesses extracted from a turbulent time series. The variational method outperforms Newton(-hookstep) iterations in converging successfully from poor initial guesses, suggesting a larger convergence radius.


Introduction
Viewing fluid turbulence as a deterministic chaotic dynamical system has revealed new insights beyond what can be achieved through a purely statistical approach (see reviews by Kawahara et al. (2012) and Graham & Floryan (2021)).The idea for a dynamical description by envisioning turbulence as a chaotic trajectory in the infinite-dimensional state space of the Navier-Stokes equations dates back to the seminal work of Hopf (1948).A remarkable progress in bridging the gaps between ideas from dynamical systems theory and practically studying turbulence in this framework has been the numerical computation of invariant solutions -an advance that did not happen until the 1990's.Invariant solutions are nonchaotic solutions to the governing equations with simple dependence on time.This includes equilibria (Nagata (1990)), travelling waves (Faisst & Eckhardt (2003); Wedin & Kerswell (2004)), periodic and relative periodic orbits (Kawahara & Kida (2001); Chandler & Kerswell (2013); Budanur et al. (2017)) and invariant tori (Suri et al. (2019); Parker et al. (2023)).In the dynamical description, the chaotic trajectory of the turbulent dynamics transiently, yet recurringly, visits the neighbourhood of the unstable invariant solutions embedded in the state space of the evolution equations.In this picture, therefore, unstable invariant solutions serve as the building blocks supporting the turbulent dynamics, and extracting them is the key for studying turbulence in the dynamical systems framework.
Equilibria of plane Couette flow (PCF) numerically computed by Nagata (1990) were the first nontrivial invariant solutions discovered in a wall-bounded three-dimensional (3D) fluid flow.Despite their lack of temporal variation, equilibrium solutions can capture essential features of chaotic flows and play an important role in characterising their chaotic dynamics.In PCF for instance, Nagata (1990), Clever & Busse (1992), Waleffe (1998), Itano & Toh (2001), Wang et al. (2007) and others compute equilibrium solutions.These equilibria typically contain wavy streaks together with pairs of staggered counter-rotating streamwise vortices, and thus capture basic structures of near-wall turbulence.Gibson et al. (2008Gibson et al. ( , 2009) ) and Halcrow et al. (2009) demonstrate how the chaotic dynamics is organised by coexisting equilibrium solutions together with their stable and unstable manifolds; Schneider et al. (2010) and Gibson & Brand (2014) compute equilibria that capture localisation in the spanwise direction; Eckhardt & Zammert (2018) compute equilibria that capture localisation in the streamwise direction; Brand & Gibson (2014) compute equilibria that capture localisation in both streamwise and spanwise directions; and Reetz et al. (2019) identify an equilibrium solution underlying self-organised oblique turbulent-laminar stripes.While equilibrium solutions have been shown to capture features of the chaotic flow dynamics, their numerical identification in very high-dimensional fluid flow problems remains challenging.
One approach to computing equilibrium solutions is to consider a root finding problem.Irrespective of their dynamical stability, equilibria of the dynamical system / =  () are, by definition, roots of the nonlinear operator governing the time evolution,  () = 0.The root finding problem can be solved by Newton(-Raphson) iterations.Newton iterations are popular because of their locally quadratic convergence.However, employing Newton iterations for solving the root finding problem has two principal drawbacks: For a system described by  degrees of freedom, the update vector in each iteration is the solution to a linear system of equations whose coefficient matrix is the  ×  Jacobian.Solving this large system of equations and the associated quadratically scaling memory requirement are too costly for very high-dimensional, strongly coupled fluid flow problems.In addition to poor scaling, Newton iterations typically have a small radius of convergence, meaning that the algorithm needs to be initialised with an extremely accurate initial guess in order to converge successfully.Finding sufficiently accurate guesses is not simple even for weakly chaotic flows close to the onset of turbulence.Newton-GMRES-hookstep is the state-ofthe-art matrix-free variant of the Newton method commonly used for computing invariant solutions of fluid flows.This method defeats the  2 memory scaling drawback by employing the generalised minimal residual (GMRES) method and approximating the update vector in a Krylov subspace (Saad & Schultz (1986); Tuckerman et al. (2019)).In addition, the robustness of the convergence is improved via hook-step trust-region optimisation (Dennis & Schnabel (1996); Viswanath (2007Viswanath ( , 2009))).Newton-GMRES-hookstep thereby enlarges the basin of convergence of Newton iterations.Yet, requiring an accurate initial guess is still a bottleneck of this method, and identifying unstable equilibria remains challenging.
An alternative to the root finding setup is to view the problem of computing an equilibrium solution as an optimisation problem.Deviation of a flow field from being an equilibrium solution can be penalised by the norm of the to-be-zeroed right-hand side operator, ∥ () ∥.The absolute minima of this cost function, ∥ () ∥ = 0, correspond to equilibrium solutions of the system.Therefore, the problem of finding equilibria can be recast as the minimisation of the cost function.A matrix-free method is crucial for solving this minimisation problem in very high-dimensional fluid flows.Farazmand ( 2016) proposed an adjoint-based minimisa-tion technique to find equilibria and travelling waves of a 2D Kolmogorov flow.The adjoint calculations allow the gradient of the cost function to be constructed analytically as an explicit function of the current flow field.This results in a matrix-free gradient-descent algorithm whose memory requirement scales linearly with the size of the problem.The adjoint-based minimisation method is significantly more robust to inaccurate initial guesses in comparison to its alternatives based on solving a root finding problem using Newton iterations.This improvement, however, is obtained by sacrificing the quadratic convergence of the Newton iterations and exhibiting slow convergence.In the context of fluid mechanics, the variational approach has been successfully applied to the 2D Kolmogorov flows (see Farazmand (2016); Parker & Schneider (2022)).
Despite the robust convergence and favourable scaling properties of the adjoint-based minimisation method, it has not been applied to 3D wall-bounded flows.Beyond the highdimensionality of the 3D wall-bounded flows, the main challenge in the application of this method lies in handling the wall boundary conditions that cannot be imposed readily while evolving the adjoint-descent dynamics (see §3.3).This is in contrast to doubly periodic 2D (or triply periodic 3D) flows where the adjoint-descent dynamics is subject to periodic boundary conditions only, that can be imposed by representing variables in Fourier basis (Farazmand (2016); Parker & Schneider (2022)).To construct a suitable formulation for 3D flows in the presence of walls, we project the evolving velocity field onto the space of divergencefree fields and constrain pressure so that it satisfies the pressure Poisson equation instead of evolving independent of the velocity (see §3.4).However, solving the pressure Poisson equation with sufficient accuracy is not straightforward in wall-bounded flows.The challenge in computing the instantaneous pressure associated with a divergence-free velocity field stems from the absence of explicit physical boundary conditions on pressure at the walls (Rempfer (2006)).As a result, a successful implementation of the constrained dynamics hinges on resolving the challenge of accurately expressing the pressure in wall-bounded flows.
We propose an algorithm for computing equilibria of wall-bounded flows using adjointdescent minimisation in the space of divergence-free velocity fields.The proposed algorithm circumvents the explicit construction of pressure, thereby overcoming the inherent challenge of dealing with pressure in the application of the adjoint-descent method to wall-bounded flows.We construct equilibria of plane Couette flow, and discuss the application of the introduced method to other wall-bounded flows and other types of invariant solutions where the challenge of dealing with pressure exists analogously.To accelerate the convergence of the algorithm we propose a data-driven procedure which takes advantage of the almost linear behaviour of the adjoint-descent dynamics in the vicinity of an equilibrium solution.The acceleration technique approximates the linear dynamics using dynamic mode decomposition, and thereby approximates the asymptotic solution of the adjoint-descent dynamics.The large basin of convergence together with the improved convergence properties renders the adjoint-descent method a viable alternative to the state-of-the-art Newton method.
The remainder of the manuscript is structured as follows: The adjoint-based variational method for constructing equilibrium solutions is introduced in a general setting in §2.The adjoint-descent dynamics is derived for wall-bounded shear flows in §3, and an algorithm for numerically integrating the derived dynamics is presented in §4.The method is applied to plane Couette flow in §5 where the convergence of multiple equilibria is demonstrated.The data-driven procedure for accelerating the convergence is discussed in §6.Finally, the article is summarised and concluding remarks are provided in §7.

Adjoint-descent method for constructing equilibrium solutions
Consider a general autonomous dynamical system u  = r(u), (2.1) where u is an -dimensional real-valued field defined over a -dimensional spatial domain x ∈ Ω ⊆ R  and varying with time  ∈ R. Within the space of vector fields ℳ = {u : Ω → R  }, the evolution of u is governed by the smooth nonlinear operator r subject to time-independent boundary conditions at Ω, the boundary of Ω. Equilibrium solutions of this dynamical system are u * ∈ ℳ for which r(u * ) = 0. (2. 2) The residual of Equation (2.2) is not zero for non-equilibrium states u ≠ u * .We thus penalise non-equilibrium states by the non-negative cost function  2 defined as where ⟨•, •⟩ denotes an inner product defined on ℳ.The cost function takes zero value if and only if u = u * .We thereby recast the problem of finding equilibrium solutions u * as a minimisation problem over ℳ, and look for the global minima of  2 at which  2 = 0, following the arguments of Farazmand (2016).
In order to find minima of  2 , we construct another dynamical system in ℳ whose evolution monotonically decreases the cost function  2 .The objective is to define an evolution equation where the choice of the operator g guarantees Here,  is a fictitious time that parametrizes the evolution governed by the constructed dynamics.The rate of change of  2 along trajectories of the dynamical system (2.4) is where ℒ(u; g) is the directional derivative of r(u) along u/ = g: (2.7) We can rewrite Equation (2.6) as where ℒ † is the adjoint operator of the directional derivative ℒ, with the following definition: (2.9) To guarantee the monotonic decrease of  2 with  we choose g(u) = −ℒ † (u; r). (2.10) This choice results in monotonic decrease of  2 along solution trajectories of the adjoint

Focus on Fluids articles must not exceed this page length
Identifying invariant solutions of wall-bounded flows using variational techniques 5 (a) u/ = r(u) (b) u/ = g(u) Figure 1: Replacing the original dynamics with the gradient descent of the cost function  = ∥r(u) ∥ by the adjoint-descent method.Panel (a) schematically shows the trajectories and two equilibria of the original system parametrized by the physical time , while panel (b) shows contours of  and sample trajectories of its gradient flow parametrized by the fictitious time .Trajectories of the adjoint-descent dynamics converge to a stable fixed point, that is either an equilibrium of the original dynamics, where the global minimum value of  = 0 is achieved, or a state at which  takes a local minimum value.
dynamical system (2.4): (2.11) In summary, in order to find equilibria of u/ = r(u) the variational approach proposed by Farazmand (2016) constructs a globally contracting dynamical system u/ = g(u), that is essentially the gradient descent of the cost function  2 .Every trajectory of the constructed dynamical system eventually reaches a stable equilibrium corresponding to a minimum of the cost function.Equilibria of the original dynamics are equilibria of the adjoint dynamics at which the cost function takes its global minimum value  2 = 0.However, the adjoint dynamics might have other equilibria that correspond to a local minimum of the cost function with  2 > 0, and are not equilibria of the original dynamics.This is schematically illustrated in Figure 1.Finding equilibria of u/ = r(u) requires integrating the adjoint dynamics u/ = g(u) forward in the fictitious time .The solutions obtained at  → ∞ for which  2 = 0 are equilibria of the original system.Otherwise, when the trajectory gets stuck in a local minimum of the cost function, the search fails and the adjoint dynamics should be integrated from another initial condition.

Governing equations
We consider the flow in a three-dimensional rectangular domain Ω of non-dimensional size  ∈ [0,   ),  ∈ [−1, +1] and  ∈ [0,   ).The domain is bounded in  between two parallel plates, and is periodic in the lateral directions  and .Incompressible, isotherm flow of a Newtonian fluid is governed by the Navier-Stokes equations (NSE).The non-dimensional, perturbative form of the NSE reads Here,  is the Reynolds number and u  is the laminar base flow velocity field.The fields u and  are the deviations of the total velocity and pressure from the base flow velocity and pressure fields, respectively.For common driving mechanisms such as the motion of walls in the  plane, externally imposed pressure differences, or injection/suction through the walls, the laminar base flow satisfies the inhomogeneous boundary conditions (BCs), absorbs body forces, and is known a priori.Consequently, the perturbative Navier-Stokes equations (3.1) and (3.2) are subject to the boundary conditions This requires  to satisfy the Poisson equation with a velocity-dependent source term (Canuto et al. (2007); Rempfer (2006)).
We could not derive a variational dynamics based on expressing pressure explicitly as the solution to this Poisson equation.Therefore, instead of the state space ℳ of the NSE we define the search space such that 'velocity' and 'pressure' can evolve independently.Accordingly, we define the cost function such that residuals of both Equations (3.1) and (3.2) are included.Otherwise, the derivation follows §2.

The search space
We define the inner product space of general flow fields as where v and  are sufficiently smooth functions of space.Hereafter, the symbols u,  indicate physically admissible velocity and pressure, implying u ∈ ℳ and  satisfying the relevant Poisson equation.The space of general flow fields  is endowed with the real-valued inner Here • is the conventional Euclidean inner product in R 3 .Physically admissible velocity and pressure fields form the following subset of the general flow fields: where  0 is the subset of  whose vector-valued component satisfies the homogeneous Dirichlet BC at the walls: We aim to impose the zero-divergence constraint together with the defining property of an equilibrium solution via the variational minimisation discussed in §2.To that end, we consider an evolution in the space of general flow fields U = [v, ] ∈  0 in which the velocity and the pressure component are evolved independently.A flow field U ∈  0 neither necessarily satisfies the defining property of an equilibrium solution nor the zero-divergence constraint.Therefore, we define the residual field R ∈  associated with a general flow field as and the cost function  2 as At the global minima of the cost function,  2 = 0, the defining property of an equilibrium solution (3.12) and the incompressibility constraint (3.2) are both satisfied.The operator The operator G is derived following the adjoint-based method described in §2 to guarantee the monotonic decrease of the cost function along trajectories of the variational dynamics (3.15).

3.3.
Adjoint operator for the NSE The variational dynamics (3.15) must ensure that the flow field U remains within  0 , thus U is periodic in  and , and its velocity component v takes zero value at the walls for all .In order for these properties of U to be preserved under the variational dynamics, the operator G must be periodic in  and , and g 1 = v/ must take zero value at the walls, meaning that G ∈  0 .In addition, we choose the residual R to lie within  0 .The periodicity of R in  and  automatically results from the spatial periodicity of U in these two directions.However, at the walls we enforce the condition r 1 (v, )| =±1 = F (v, )| =±1 = 0.With the choice of U, R, G ∈  0 , the flow field remains within  0 as desired.Following this choice, all the boundary terms resulting from partial integrations in the derivation of the adjoint operator cancel out (see Appendix A), and the adjoint of the directional derivative of R(U) along G is obtained as Therefore, with G = −ℒ † (U; R) the variational dynamics takes the form Equation (3.18) is fourth order with respect to v and Equation (3.19) is second order with respect to .Therefore, four BCs for each component of v and two BCs for  are required in the inhomogeneous direction .The choice of U ∈  0 implies v = 0 and the choice of R ∈  0 implies r 1 = F (v, ) = 0 at each wall.Consequently, the adjoint-descent dynamics requires two additional wall BCs in order to be well-posed.As the additional BCs, we impose ∇ • v = 0 at each wall.Therefore, the adjoint-descent dynamics is subject to the following BCs: ) forward in the fictitious time is not straightforward.Consequently, instead of advancing the derived variational dynamics directly, we constrain the adjoint-descent dynamics to the subset of physical flow fields   .Within this subset, pressure does not evolve independently but satisfies the pressure Poisson equation.Thereby, we obtain an evolution equation for the velocity within the state space of the NSE.This allows us to employ the influence matrix method (Kleiser & Schumann (1980)) to integrate the constrained adjoint-descent dynamics.

Variational dynamics constrained to the subset of physical flow fields
To obtain a numerically tractable variational dynamics, we constrain the adjoint-descent dynamics (3.18)-(3.24) to the subset of physical flow fields   .Within   , the velocity component u is divergence-free over the entire domain.In addition, the pressure component  is no longer governed by an explicit evolution equation, but by Poisson equation with a velocity-dependent source term.Let  = P [u] denote the solution to the Poisson equation yielding pressure associated with an instantaneous divergence-free velocity u.In order for u to remain divergence-free, g 1 needs to be projected onto the space of divergence-free fields, yielding the evolution where P denotes the projection operator.The argument of the operator P is the right-hand side of Equation (3.18) with  2 = 0 and ∇ 2 = 0 that result from the zero divergence of u.According to the Helmholtz's theorem, a smooth 3D vector field can be decomposed into a divergence-free and a curl-free component.Thus, g 1 = u/ is decomposed as g 1 = f − ∇ where f = P {g 1 } is the divergence-free component and  is the scalar potential whose gradient gives the curl-free component.Therefore, the evolution of the divergence-free velocity is governed by where r 1 = r 1 (u, P [u]) and thus the BC (3.31) is automatically satisfied.It is necessary to verify that the constrained variational dynamics still guarantees a monotonic decrease of the cost function.For U ∈   , the scalar component of the steepest descent direction,  2 , vanishes (see Equation (3.19)).Therefore, according to the definition of the inner product on   (3.9), it is sufficient to verify that The Helmholtz decomposition is an orthogonal decomposition with respect to the  2 inner product defined on the state space of the NSE, ⟨f, ∇⟩ ℳ = 0. Therefore, ⟨f, ℳ ⩾ 0, and thus the evolution of u along f guarantees the monotonic decrease of the cost function, as desired.
The variational dynamics (3.26)-(3.31) is equivariant under continuous translations in the periodic directions  and .Furthermore, one can verify through simple calculations that this dynamics is also equivariant under the action of any reflection or rotation permitted by the laminar base velocity field u  .Consequently, the symmetry group generated by translations, reflections and rotations in the obtained variational dynamics is identical to that of the NSE (3.1)-(3.5).Therefore, to construct equilibria within a particular symmetry-invariant subspace of the NSE, one can use initial conditions from the same symmetry-invariant subspace to initialise the variational dynamics, and the variational dynamics preserves the symmetries of the initial condition.
In the variational dynamics, the scalar field  plays a role analogous to the pressure  in the incompressible NSE.The scalar fields  and  adjust themselves to the instantaneous physical velocity u such that ∇ • u = 0 and u( = ±1) = 0 are preserved under the evolution with the fictitious time  and the physical time , respectively.Similar to the pressure in the NSE,  satisfies a Poisson equation with a velocity-dependent source term.Solving the Poisson equation for  and  is a numerically challenging task in the present wall-bounded configuration (Rempfer (2006)).Therefore, instead of attempting to compute  and  and thereby advancing the variational dynamic (3.26), we formulate the numerical integration scheme based on the influence matrix method (Kleiser & Schumann (1980)) where the noslip BC and zero divergence are precisely satisfied while the explicit construction of  and  is circumvented.

Numerical implementation
To advance the variational dynamics (3.26)-(3.31)without explicitly computing  and , we take advantage of the structural similarity between the variational dynamics and the NSE.In order to evaluate the right-hand side of Equation (3.26), we consider the following PDE for the residual field r 1 : where N(r 1 ) = (∇r 1 ) (u  + u)−(∇(u  + u)) ⊤ r 1 with both u and u  being treated as constant fields.We use the dummy Equation (4.1) to evaluate the right-hand side of Equation (3.26) since the instantaneously evaluated right-hand side of these two systems are identically equal.For brevity, we are omitting the periodic BCs in  and  since spatial periodicity can be enforced via spectral representation in an appropriate basis, such as Fourier basis, that is periodic by construction.Equation (4.1) together with the BC (4.2) and the zero-divergence constraint (4.3) resembles the structure of the incompressible NSE: which is subject to The influence matrix (IM) algorithm has been developed to numerically advance this particular type of dynamical systems, which have a Laplacian linear term and gradient of a scalar on the right-hand side, and are subject to zerodivergence constraint and homogeneous Dirichlet BCs at the walls.This algorithm enforces zero divergence and the homogeneous Dirichlet BCs within the time-stepping process while the scalar field is handled implicitly and is not resolved as a separate variable (Kleiser & Schumann (1980)

Rapids articles must not exceed this page length
(ii) The residual field r 1 , which is by definition the right-hand side of the NSE (3.1), is approximated via finite differences Since both u and u Δ are divergence-free and satisfy homogeneous Dirichlet BCs at the walls, ∇ • r 1 = 0 and r 1 ( = ±1) = 0.
(v) Having approximated f, which is the descent direction at the current fictitious time , we advance the velocity for one step of size Δ using u Δ = u + Δ f. (4.9) Since both u and f are divergence-free and take zero value at the walls, the updated velocity satisfies ∇ • u Δ = 0 and u Δ ( = ±1) = 0.The finite differences (4.7) and (4.8) affect the accuracy of time-stepping the variational dynamics, but they do not interfere with imposing the boundary condition u( = ±1) = 0 and the constraint ∇ • u = 0 within machine precision.The low accuracy of the first-order finite differences does not affect the accuracy of the obtained equilibrium solution since both r 1 and f tend to zero when an equilibrium is approached.We are also not concerned about the low accuracy of the first-order forward Euler update rule (4.9) since the objective is to obtain the attracting equilibria of the adjoint-descent dynamics reached at  → ∞.Therefore, the introduced procedure is able to construct equilibrium solutions within machine precision.
We implement this procedure in Channelflow 2.0, an open-source software package for numerical analysis of the incompressible NSE in wall-bounded domains.In this software, an instantaneous divergence-free velocity field is represented by Chebyshev expansion in the wall-normal direction  and Fourier expansion in the periodic directions  and : where   () is the -th Chebyshev polynomial of the first kind,  is the imaginary unit, and indices 1 to 3 specify directions ,  and , respectively.Channelflow 2.0 employs the influence matrix algorithm for time-marching the NSE (4.4).With modification for the nonlinear term N(r 1 ), Equation (4.1) can also be advanced in time.

Application to plane Couette flow
We apply the introduced variational method to plane Couette flow (PCF), the flow between two parallel plates moving at equal and opposite velocities.PCF is governed by the general NSE (3.1)-(3.5)with the laminar base flow u  = [, 0, 0] ⊤ .Due to the periodicity in  and , PCF is equivariant under continuous translations in these directions: where ,  and  are the components of u in ,  and  directions, respectively.In addition, PCF is equivariant under two discrete symmetries as well, rotation around the line  =  = 0: (5. 3) The variational dynamics (3.26)-(3.31) is easily verified to be equivariant under the same continuous and discrete symmetry operators.Therefore, the variational dynamics preserves these symmetries, if present in the initial condition.In the following, we demonstrate the convergence of multiple equilibrium solutions from guesses both within a symmetry-invariant subspace and outside.

Results
We search for equilibria of PCF at  = 400 within a domain of dimensions   = 2/1.14and   = 2/2.5(see §3.1).This domain was first studied by Waleffe (2002).Several equilibrium solutions of PCF in this domain at  = 400 were computed by Gibson et al. (2008Gibson et al. ( , 2009)).These are available in the database on channelflow.org.Here, the flow field is discretised with   = 31 collocation points in the wall-normal direction and   =   = 32 points in the lateral directions.The adjoint-descent dynamics is numerically integrated by the forward Euler scheme (4.9) with Δ = 0.03, and r 1 and f are approximated via finite differences (4.7) and (4.8) with the step size Δ = 0.25 and Δ τ = 0.25, respectively (see §4).An accurate finite-difference approximation of r 1 and f suggests choosing Δ and Δ τ as small as possible.However, smaller values for these step sizes result in a less stable forward Euler integration scheme, requiring a smaller value of Δ to remain stable.Since for an equilibrium solution r 1 = f = 0, larger values of Δ and Δ τ do not diminish the accuracy of the obtained equilibrium solution.Consequently, we empirically choose values for Δ and Δ τ so that a reasonably large value for Δ can be used.
To verify the scheme and its implementation, we converge the so-called 'Nagata's lower branch' equilibrium solution (Nagata (1990)) at  = 400.As initial guess, we take an equilibrium solution on the same branch but at a significantly different .The Nagata's lower branch solution at  = 400 continued from Nagata's original domain dimensions to the ones considered here is available in the database on channelflow.org.We continue this equilibrium solution to  = 230, and use the resulting solution to initialise both the adjoint-descent variational method and the standard Newton iterations at  = 400.The standard Newton iterations, i.e. without optimisations such as hook steps, fail to converge.However, the adjoint-descent variational method successfully converges to the equilibrium solution at  = 400 on the same branch.
Along the trajectory of the adjoint-descent dynamics, the cost function initially drops rapidly and subsequently decreases with an exponential rate, as shown in Figure 2. The exponential decrease of the cost function is explained by the dynamical system picture of the adjoint descent: the adjoint-descent dynamics converges to a stable fixed point, hence the evolution is dominated by the slowest eigenmode of the linearised dynamics in the vicinity of that fixed point.The sharp initial drop and the following exponential decay of the cost function are reflected in fast and slow traversal, respectively, of the trajectory within the state space.Figure 3 presents a 2D projection of the trajectory, with markers indicating that the majority of the trajectory is traversed quickly in the beginning of the integration, and the majority of the integration time is spent on the remaining, much shorter portion of the trajectory.For instance, the portion of the trajectory traversed during the first 1.2 × 10 6 fictitious time units, that decreases the cost function from  = 5.9 × 10 −3 to  = 10 −5 , is against  1 = ℜ{ û0,3,0,1 }.The majority of the trajectory is traversed rapidly at the beginning, as indicated by a sharp drop of  in Figure 2, followed by a slow traversal of the remaining portion towards the asymptotic solution, reflected in Figure 2 as an exponential decay of the cost function.
considerably longer than the remaining portion which takes over 90 % of the integration time to be traversed.In Figure 3,  1 and  2 are the real parts of û0,3,0,1 and û0,5,0,1 , i.e. the coefficients of the third and the fifth Chebyshev polynomial in the expansion of the mean streamwise velocity in  (see Equation (4.10)).The visualisation of the trajectory in different projections of the state space yields a similar observation.Nagata's lower branch equilibrium solutions are symmetric under shift-and-rotate symmetry  1 = (  /2,   /2) 1 : (5.4) The  2 -norm of the velocity field against the physical time  in direct numerical simulation from a random initial condition.The snapshots corresponding to the local extrema of ∥u∥ are selected as guesses for an equilibrium solution.Table 1 summarises the result of the convergence from each guess using Newton-GMRES-hookstep and the adjoint-descent variational method.
and shift-and-reflect symmetry  2 = (  /2, 0) 2 : (5.5) Therefore, the initial guess in the present example, namely the Nagata's lower branch solution at  = 230, is symmetric under  1 and  2 that are preserved by the adjoint-descent dynamics.
The velocity field remains symmetric under  1 and  2 without explicitly enforcing them during the forward integration until the equilibrium solution on the same branch at  = 400 is converged.
To further investigate the robustness of the adjoint-descent variational method in successfully converging from inaccurate guesses, we initialise the method with guesses obtained from a direct numerical simulation.We construct a random divergence-free velocity field with  2 -norm ∥u∥ = 0.2, and time-march the NSE along a turbulent trajectory until the flow laminarises.The initial condition and therefore the entire trajectory are not symmetric under any of the symmetries allowed by the PCF.We extract the local extrema of ∥u∥ as a function of time , where  ∥u∥/ = 0, as guesses for potential equilibrium solutions.Figure 4 shows ∥u∥ plotted against  from which 26 guesses are extracted.The standard Newton iterations do not converge starting from any of the guesses.With hook-step optimisation, 5 of the searches converge within 50 Newton-GMRES-hookstep (NGh) iterations.The converged solutions include the trivial laminar solution u = 0 as well as two nontrivial solutions EQ1 and EQ3 (see Tables 1 and 2 for properties of the converged solutions).By integrating the adjointdescent dynamics, 11 of the guesses converge to an equilibrium solution.These solutions include the trivial solution as well as five nontrivial equilibria EQ1 to EQ5 (see Tables 1  and 2).Among these solutions, EQ1, EQ4 and EQ5 have been documented in the literature (Gibson et al. (2009)).Yet, to the best of our knowledge, the equilibria labelled EQ2 and EQ3 have not been previously reported.Snapshots that lead to a successful search via either NGh iterations or the adjoint-descent algorithm are marked on Figure 4.
The variational method succeeds in more than twice as many cases as the NGh method, and extracts three more non-trivial equilibria from a turbulent trajectory with a crude criterion Identifying invariant solutions of wall-bounded flows using variational techniques 15 snapshot NGh iterations NGh solution adjoint-descent solution Table 1: The list of the equilibrium solutions converged by Newton-GMRES-hookstep (NGh) and the adjoint-descent variational method from the guesses marked in Figure 4. See Table 2 2: Properties of the equilibrium solutions converged by Newton-GMRES-hookstep and the adjoint-descent variational method (see Table 1 and Figure 4).The second column contains the  2 norm of the solutions and the third column contains the total energy dissipation of the solutions normalised by that of the laminar base flow.
for selecting guesses.This suggests that the basin of attraction to converge an equilibrium solution is typically larger for the adjoint-descent variational method compared to the NGh method.However, the larger basin of attraction does not necessarily contain the smaller one.Notice, for instance, that the NGh iterations and the adjoint-descent algorithm converge to different equilibrium solutions when initialised with the snapshot 4, or the NGh iterations converge when initialised with the snapshot 5 while the adjoint-descent does not.Despite the advantage of the variational method in successfully converging from inaccurate guesses, this method exhibits a very slow rate of convergence.For instance, the convergence in our first example (Figures 2) takes near 650 hours of wall clock time on one core of a 2.60GHz Intel Xeon E5-2640 CPU.Besides the improvements on the computer programming side, such as parallel computations on multiple CPU cores, the convergence can be significantly accelerated by employing the inherent predictability of the variational dynamics, namely its almost linear behaviour when the trajectory reaches the vicinity of a solution.In the following, we introduce a data-driven technique for such an acceleration.

Accelerating the convergence
The variational dynamics evolves along the gradient descent of the cost function.As a result, this dynamics is globally contracting, and almost all its trajectories eventually converge to a stable fixed point where the cost function takes a minimum value.When the trajectory of the adjoint-descent dynamics has got sufficiently close to its destination fixed point, the cost function is well represented by a quadratic function and its gradient flow is almost linear.The approximately linear behaviour of the variational dynamics in the vicinity of an asymptotic fixed point inspires the idea of the following data-driven technique for accelerating the slow convergence of the variational method.
Our acceleration technique aims to approximate the expected linear dynamics and thereby approximate the equilibrium solution of the adjoint-descent dynamics.Since the destination fixed point is not known a priori, linearisation around the unknown fixed point is obviously not possible.Instead, we employ dynamic mode decomposition (DMD) to approximate the linear dynamics based on the available portion of the trajectory that has been traversed.DMD is a regression framework that constructs the best-fit linear model over a series of snapshots (Rowley et al. (2009); Schmid (2010); Kutz et al. (2016); Schmid (2022)).The equilibrium solution of the adjoint-descent dynamics is approximated by letting the fictitious time go to infinity in the approximated linear system.

Dynamic mode decomposition (DMD)
Suppose each instantaneous spatially resolved flow field u(x; ) is represented by an dimensional real-valued column vector (). snapshots   = (  );  = 1, . . .,  along a single trajectory can be related to the snapshots taken  later along the same trajectory,  ′  = (  + ), via the following linear relation: where   is the error in approximating  ′  by the linear map   ↦ → A  .DMD constructs the  ×  linear operator A which minimises the sum of squares of the elements of   over all  snapshot pairs: where . . ′  , and the superscript + denotes the Moore-Penrose pseudo-inverse.The dimensionality of the system can be prohibitively large for constructing A directly as defined in Equation (6.2), which is typically the case in a fluid dynamics problem.Therefore, we instead use a rank-reduced representation of this matrix.For this, the data matrix Ψ is factorised via singular value decomposition (SVD) as Ψ ≈ UΣV ⊤ with truncation rank .The  ×  projection of A on the POD modes U is The dynamic modes and their temporal behaviour are constructed from the eigendecomposition of Ã: Dynamic modes are   = Ψ ′ VΣ −1   with  = 1, . . ., , where   are eigenvectors of Ã; and the dynamic mode   evolves as     where   = ln(  )/ and   is the eigenvalue of Ã associated with   .Finally, the linear evolution of () is approximated as (6.4) where   are the amplitudes of the dynamic modes at a reference time, for instance at   .Based on this linear model we approximate the asymptotic equilibrium solution of the variational dynamics as follows.
6.2.Numerical implementation In order to approximate the linear dynamics using DMD, we collect snapshots   after the initial fast drop in the cost function, when it decays exponentially.The exponential decay implies that the dynamics is almost linear and dominated by only few of the slowest attracting eigenmodes.As a result, the snapshot matrices are of low column rank.In the vicinity of the yet to-be-found attracting fixed point, the Jacobian of the descent dynamics is the negative of the (discretised) second variation, or Hessian, of the cost function, and thus symmetric.Consequently, the eigenvalues of the linear dynamics approximated by the DMD are real.Suppose the dynamic modes are sorted in increasing order of |  |.Then, the exponent  1 is significantly closer to zero than the rest, and  2 , . . .,   are negative.By assuming  1 ≈ 0, the linear model (6.4) can be expressed as the superposition of the steady state   :=  1  1 , and the decaying terms     exp(  );  = 2, . . ., .The steady state   approximates the equilibrium solution of the almost linear adjoint-descent dynamics.The state vector   is mapped back to the corresponding flow field, from where the integration of the adjointdescent dynamics is restarted.Let  * denote the rank of the snapshot matrices.Then, the truncation rank  ⩽  * is chosen such that the cost function associated with the approximated equilibrium is the smallest.We consistently found the minimum for  =  * .In the following, we demonstrate the acceleration of the first test case presented in §5.
The snapshot vectors  are the (real-valued) state vectors containing the minimum number of independent variables required for describing a divergence-free velocity field in Fourier-Chebyshev-Fourier spectral representation (4.10).The vector  has  = 20 218 elements for the discretisation used in §5.Initially, we integrate the adjoint descent dynamics and let the cost function drop to log() = −4.5 before performing the first DMD extrapolation.The linear model is constructed using  = 100 snapshots uniformly spaced over an interval of 2 × 10 4 time units ( = 200).The next DMD extrapolations are performed using the same number of snapshots  and the same spacing  while the adjoint dynamics is integrated forward in time for 15 × 10 4 time units before starting to collect new snapshots.The acceleration technique allows to achieve the convergence criterion  = 10 −12 through  = 7.36 × 10 5 time units of total forward integration while without acceleration it takes  = 1.38 × 10 7 time units, that is almost 19 times longer (see Figure 5, compare with Figure 2).The time required for performing the extrapolation is negligible compared to the time required for the forward integration of the adjoint-descent dynamics.The first DMD extrapolation has resulted in a slight increase in the value of .The 2D projection of the state space, displayed in Figure 6, shows that the first extrapolated state is significantly closer to the destination fixed point, despite being located on a higher level of .By restarting the integration from the extrapolated state, the trajectory gets quickly attracted to the dominating eigendirection of the linearised dynamics resulting in a rapid drop in  (see Figures 5 and 6).
Exploiting the linear behaviour of the variational dynamics, the acceleration technique typically achieves an order of magnitude speed-up in converging equilibria of PCF.The linear behaviour in the vicinity of an equilibrium solution at sufficiently large  is a generic characteristic of the adjoint-descent variational method.Therefore, the introduced DMDbased acceleration technique is system-independent, and provided the snapshot vectors of the variational dynamics can be applied directly to any other problem.

Summary and concluding remarks
The unstable invariant solutions embedded within the chaotic attractor of the Navier-Stokes equations underpin the dynamics of a turbulent flow.Despite the significance of invariant solutions for a dynamical description of chaotic flows, the identification of these solutions remains a computational challenge, demanding robust algorithms.In this work, we have presented a matrix-free, adjoint-based variational method for computing equilibrium solutions of wall-bounded shear flows.We have applied the introduced method to plane Couette flow, and demonstrated the convergence of multiple equilibrium solutions.The variational method outperforms the state-of-the-art Newton iterations in successfully converging from inaccurate initial guesses, that suggests a larger basin of attraction.The present method employs the norm of the right-hand side of the evolution equation as a cost function to penalise the deviation of a flow field from the equilibrium state.Thereby, the problem of finding an equilibrium solution is recast as the minimisation of the cost function.To solve the minimisation problem, we adopted the variational approach of Farazmand (2016) where the gradient of the cost function is constructed analytically via adjoint calculations, and thereby a matrix-free gradient descent method is utilised.The cost function decreases monotonically along trajectories of the gradient descent dynamics until a minimum value is obtained.The global minima of the cost function, taking zero value, correspond to the equilibrium solutions of the flow.If a local minimum is obtained, the search for an equilibrium solution has failed.However, a local minimum of the cost function corresponds to the locally slowest state with respect to the chosen norm.This provides a means of characterising the so-called 'ghost' of a saddle-node bifurcation (Strogatz (2018)), which may influence the emerging spatiotemporal structures in chaotic flows (see, for example, Reetz et al. (2020), §3.1).
The present work describes two key contributions: First, we apply the adjoint-based variational method to 3D wall-bounded flows.Previously, the variational approach had only been successfully applied to a 2D Kolmogorov flow in a doubly periodic domain without walls (Farazmand (2016); Parker & Schneider (2022)).The primary challenge in extending the variational method for computing equilibria to wall-bounded flows lies in handling the nonlinear, nonlocal pressure in the presence of solid walls.To overcome this challenge, we have formulated the variational dynamics in a way that an explicit computation of pressure is avoided, allowing for application to 3D wall-bounded flows.We demonstrated the variational method for plane Couette flow.
The second contribution is addressing the slow convergence of the adjoint-based variational method, that poses a challenge in practically utilising this method for 3D Navier-Stokes equations.We propose a data-driven technique for accelerating the convergence by extrapolating the asymptotic fixed point of the variational dynamics based on the traversed portion of its trajectory.Since any trajectory of the variational dynamics converges to a stable fixed point, the dynamics behaves almost linearly when the trajectory has got close enough to the asymptotic solution.The extrapolation technique takes advantage of this predictability, and approximates the best-fit linear dynamics using dynamic mode decomposition (DMD).The asymptotic solution of the approximated linear system approximates the asymptotic solution of the variational dynamics.This results in an order-of-magnitude speed-up in the overall duration of the forward integration required to converge to a solution within machine accuracy.The proposed acceleration technique is based on the generic properties of gradient descent minimisation, and is therefore independent of the physical system of study.
The variational dynamics have been derived for the deviation of the velocity field from the laminar base flow.Consequently, an identical formulation and implementation directly translates to other wall-bounded flows such as plane Poiseuille flow (PPF) and asymptotic suction boundary layer (ASBL) flow as only the respective base velocity profiles in the variational dynamics (3.26)-(3.31)needs to be adapted.However, due to the mean advection in PPF and ASBL, steady solutions in these flows take the form of travelling waves or relative equilibria, that are equilibria in a moving frame of reference.The present method can be extended to compute travelling waves by expressing the NSE in a moving frame of reference.The speed of the Galilean transformation is an additional unknown in the cost function whose minimisation yields relative equilibrium solutions (Farazmand (2016)).The handling of the pressure and boundary conditions remains unchanged.
In the derivation of the adjoint-descent dynamics we assume the base velocity profile to be known and constant, which is the case when a fixed mean pressure gradient is imposed on the flow.For the alternative integral constraint of fixed mass flux, the base velocity profile, or more precisely its amplitude, needs to be determined together with the velocity and pressure perturbations.As for the generalisation for travelling waves, the additional unknown can be included in the variational formulation without modifying the handling of pressure and boundary conditions.
The advantages of the adjoint-based variational method have inspired its application in computing other invariant sets, such as periodic orbits (Azimi et al. (2022); Parker & Schneider (2022)) and connecting orbits (Ashtari & Schneider (2023)).These methods view the identification of a periodic or connecting orbit as a minimisation problem in the space of space-time fields with prescribed behaviour in the temporal direction.They then employ a similar adjoint-based technique to solve the minimisation problem.The robust convergence of these extensions has so far only been demonstrated in 2D flows in a doubly periodic domain and for 1D model systems.Like in computing equilibria, dealing with pressure is the key challenge in formulating the adjoint-based variational method for computing periodic or connecting orbits in 3D wall-bounded flows.In our ongoing research, the next step is to extend the introduced algorithm to the computation of more complex invariant solutions in wall-bounded flows via extensions of the adjoint-based variational method.

A.2. The adjoint operator
To derive the adjoint operator of the directional derivative of the residual, ℒ(U; G), we expand the inner product of ℒ(U; G) and the residual R as follows: therefore, the components of ℒ † (U; R) are obtained as ℒ † 2 = ∇ • r 1 .

JFigure 2 :Figure 3 :
Figure 2: Convergence of the adjoint-descent variational method for constructing an equilibrium solution of the plane Couette flow.The minimisation of the cost function  evolves the initial guess towards a true equilibrium solution at which  = 0.
Figure4: The  2 -norm of the velocity field against the physical time  in direct numerical simulation from a random initial condition.The snapshots corresponding to the local extrema of ∥u∥ are selected as guesses for an equilibrium solution.Table1summarises the result of the convergence from each guess using Newton-GMRES-hookstep and the adjoint-descent variational method.

Figure 5 :Figure 6 :
Figure5: Acceleration of the convergence of the adjoint-descent variational method by successive DMD-based extrapolations.The extrapolation employs DMD to construct a best-fit linear model for the dynamics in the vicinity of an equilibrium, and approximates the asymptotic solution of the adjoint-descent dynamics by the asymptotic solution of the linear model.The acceleration technique reduces the total duration of the forward integration by 95% in this example.The jumps in the state space associated with the first two extrapolations,  1 and  2 , are shown in Figure6.
for properties of the equilibria EQ0 to EQ5.