User Tools

Site Tools


docs:utils:findorbit

This is an old revision of the document!


A PCRE internal error occured. This might be caused by a faulty plugin

====== findorbit ====== ''findorbit'': Newton-Krylov-hookstep search for (relative) periodic orbit of plane Couette flow. Algorithm by [[references|D. Viswanath]]. <code> usage: findorbit [-orb|-eqb] [-rel] [-T maptime] [-opts] guess.ff options: -eqb --equilibrium search for equilibrium or relative equilibrium (trav wave) -orb --periodicorbit search for periodic orbit or relative periodic orbit -rel --relative search for relative periodic orbit or relative equilibrium -T --maptime <real> default == 10 initial guess for orbit period or time of eqb/reqb map f^T(u) -Nn --Nnewton <int> default == 20 max number of Newton steps -Ng --Ngmres <int> default == 40 max number of GMRES iterations per restart -d --delta <real> default == 0.01 initial radius of trust region -sx --x-sign change u,x sign -sy --y-sign change v,y sign -sz --z-sign change w,z sign -ax --axshift <real> default == 0 translate by ax*Lx -az --azshift <real> default == 0 translate by az*Lz -R --Reynolds <real> default == 400 Reynolds number -vdt --variableDt adjust dt to keep CFLmin<=CFL<CFLmax -dt --timestep <real> default == 0.03125 (initial) integration timestep -ts --timestepping <string> default == SBDF3 timestep algorithm: See README. -nl --nonlinearity <string> default == Rotational nonlinearity method: See README. <flowfield> (trailing arg 1) initial guess for Newton search </code> === example 1: equilibrium === Find an equilibrium from initial guess ''guess.ff'' <code> findorbit -eqb guess.ff </code> === example 2: traveling wave === Find a traveling wave (relative equilibrium) from initial guess ''guess.ff'', with initial guess for wavespeed specified in terms of finite shifts (ax/Lx, az/Lz) in finite time T (giving wavespeeds are ax/(Lx T), az/(Lz T) <code> findorbit -eqb -rel -ax 0.3 -az 0.4 -T 10 guess.ff </code> === example 3: periodic orbit === Find a periodic orbit with initial guess for period T=83 <code> findorbit -orb -T 83 guess.ff </code> === example 4: relative periodic orbit=== Find a relative periodic orbit with initial guess of a shift by Lz/4 every period T=83 <code> findorbit -orb -rel -az 0.25 -T 83 guess.ff </code> ====== Mathematics ====== The following is an overview of the Viswanath Newton-GMRES-hookstep algorithm, meant as a reference for the options of the ''findorbit'' utility. It leaves out a number of gnarly details. And please excuse the mishmash of plain text, html, and dokuwiki-latex! Let u = discretized velocity field, as difference from laminar flow f^T(u) = time-T map of the Navier Stokes eqn, numerically integrated σ = a symmetry of plane Couette flow G(u,σ,T) = σ f^T(u) - u, a zero of G(u,σ,T) is an invariant solution of Navier-Stokes r(u,σ,T) = 1/2 ||G(u,σ,T)||^2 = residual of G(u,σ,T)=0 ||G|| = a norm of G, to be defined later. Equilibria, traveling waves, periodic orbits, and relative periodic orbits are all solutions of G(u,σ,T) = 0. * equilibria: σ is the identity and t is arbitrary. We set σ = 1, hold T fixed to a arbitrary value, and search over u * traveling waves: σ is a translation τ(cx T, cz T) and T is arbitrary. We hold T fixed and search over u and a finite phase shifts, (ax, az) = (cx T / Lx, cz T/ Lz). * periodic orbits: σ is fixed and we search over u and T. * relative periodic orbits: we search over u, σ, and T. Let x be the finite set of free variables being searched over, e.g. for a periodic orbit, 10^5 spectral coefficients of u plus the period T. We want to find a zero of G(x) or equivalently a minimum of r(x). ===== Newton algorithm ===== The Newton algorithm for solving G(x) = 0 is Let x* be the solution, i.e. G(x*) = 0. Suppose we start with an initial guess x near x*. Let x* = x + dx. Then <latex> G(x^*) = G(x) + DG(x) \, dx + O(|dx|^2) = 0 </latex> Drop higher order terms and solve <latex> DG \, dx = -G(x) </latex> Then start with new guess x' = x+dx and iterate. ===== GMRES ===== In the present case x is 10^5 dimensional or more, so it is not practical to compute the matrix DG in order to solve the Newton equation. So we use GMRES (Generalized Minimal Residuals), an iterative algorithm for solving linear systems A x = b by repeated calculations of the product A x for test values of x. Specifically, we look for solutions x within n-dimensional //Krylov subpaces// spanned by Ab, A^2b, A^3 b, ..., A^n b. If the eigenvalues of A are well-separated, you can often get a good solution x in a n-dimensional Krylov subspace with n << dim(x). In our case the product is DG dx, and it is calculated using numerical integration of Navier-Stokes. If x = (u,σ,T), then dx = (du,dσ,dT) and <latex> $ \begin{align*} DG \, dx &= G(u+du, \, \sigma+d\sigma, \, T+dT) \\ &= (\sigma+d\sigma) f^{T+dT}(u+du) - \sigma f^T(u) \end{align} $ </latex> which can be easily calculated by numerical integrations of f. For low Reynolds numbers (Re=400) and small cells [Lx,Ly,Lz] approx [pi, 2, pi], we can usually find an approximate solution to DG dx = -G(x) in forty or so iterations. That means forty or so integrations around the approximate periodic orbit. This accounts for the bulk of the CPU time consumed by ''findorbit''. For more on GMRES, see Trefethen and Bau in the [[#References]]. ===== Hookstep ===== Unfortunately, Newton iteration only works when you are very near the solution, i.e. within the radius of accurate linear approximation of dynamics around the solution. Otherwise the Newton step dx can be widely inacccurate and send you shooting off towards infinity (or "to infinity and beyond," as [[http://www.youtube.com/watch?v=8mVEGfH4s5g|Beyonce]] says). So, instead taking the Newton step directly, we limit the size of the step dx to a //trust region// <latex> \| dx \|^2 \leq \delta^2 </latex> where δ is an estimate of the radius of validity of the linearization about the current best solution x, and we look for the solution to the //constrained minimization problem// <latex> \text{minimize } \| DG \, dx + G(x) \|^2_{Krylov} \text{ subject to } \| dx \|^2 \leq \delta^2 </latex> The //Krylov// subscript indicates that the norm is the norm of the projection onto the low-d Krylov subspace. It turns out there's a unique, cheaply computable solution to this problem as a function of δ. Call it dx(δ). It has the following properties * dx(δ) aligns with the gradient of the residual |G(x)|^2 as δ → 0. * dx(δ) equals the Newton step if δ ≥ |Newton step| * dx(δ) varies smoothly between those values as a function of δ dx(δ) defines a 1d curve in the Krylov subspace. It is typically a smooth hook-shaped curve parallel to the gradient direction at one end and the Newton step direction at the other --hence the name //hookstep algorithm//. Of course, you don't know what the best radius δ for the trust region. The hookstep algorithm uses a few heuristics * start with the best δ from the previous Newton-GMRES-hookstep * evaluate the residual |G(u+dx(δ))| and compare to the reduction predicted Newton-GMRES-hookstep model * if the predicted reduction is very accurate, increase δ and evaluate the new residual |G(x+dx(δ))| * if the predicted reduction is inaccurate, decrease δ and evaluate the new residual |G(x+dx(δ))| * if the predicted reduction is moderately accurate, accept dx(δ), let x → x + dx(δ), and restart Newton-GMRES-hookstep iteration Each adjustment of δ takes two integrations around the periodic orbit. So, the upshot is this: When you are far from a solution, the hookstep algorithm creeps downhill like a gradient search. When you are near the solution, it switches to the rapidly converging Newton search. It adjusts between these algorithms smoothly, based on its estimate of the radius of validity of the linear approximations in the Newton algorithm. ====== References ====== Divakar Viswanath, "Recurrent motions within plane Couette turbulence", J. Fluid Mech. 580 (2007), http://arxiv.org/abs/physics/0604062 Lloyd N. Trefethen and David Bau, "Numerical Linear Algebra", SIAM, Philadelphia, 1997. Dennis and Schnabel, "Numerical Methods for Unconstrained Optimization and Nonlinear Equations", Prentice-Hall Series in Computational Mathematics, Englewood Cliffs, New Jersey, 1983. Beyonce, "Single Ladies", http://www.youtube.com/watch?v=8mVEGfH4s5g

docs/utils/findorbit.1233602188.txt.gz · Last modified: 2009/02/02 11:16 by gibson