docs:utils:findorbit

This shows you the differences between two versions of the page.

docs:utils:findorbit [2009/02/11 08:38] gibson |
docs:utils:findorbit [2010/02/02 07:55] |
||
---|---|---|---|

Line 1: | Line 1: | ||

- | ====== findorbit ====== | ||

- | ''findorbit'': Newton-Krylov-hookstep search for (relative) periodic orbit of plane Couette flow. Algorithm by | ||

- | [[#references|D. Viswanath]]. | ||

- | |||

- | Note: this algorithm can take hours or days to run to completion. | ||

- | |||

- | <code> | ||

- | usage: | ||

- | findorbit [-orb|-eqb] [-rel] [-T maptime] [-opts] guess.ff | ||

- | |||

- | main 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: (docs to be written) | ||

- | -nl --nonlinearity <string> default == Rotational nonlinearity method: (docs to be written) | ||

- | <flowfield> (trailing arg 1) initial guess for Newton search | ||

- | </code> | ||

- | |||

- | [[docs:utils:findorbit:alloptions|complete option list]] | ||

- | |||

- | The most important options are -orb, -eqb, -rel, and -ax and -az. You use these to specify the | ||

- | kind of solution you are searching for. | ||

- | |||

- | * -orb and -eqb are mutually exclusive and you must pick one (search for periodic orbit or equilibrium) | ||

- | * -eqb leaves the integration time fixed | ||

- | * -orb adds the integration time T to the search space | ||

- | * -rel modifies -orb / -eqb, makes the search for a relative periodic orbit or relative equilibrium (traveling wave), by adding the finite-time phase shift variables ''ax'' and ''az'' to the search space | ||

- | * -ax and -az specify a fixed or initial guess finite-time phase shifts for the solution (depending on whether -rel is specified) | ||

- | |||

- | In all cases we search for a solution of | ||

- | |||

- | <latex> | ||

- | \tau_{a_x, a_z} f^T(u) - u = 0 | ||

- | </latex> | ||

- | |||

- | where | ||

- | |||

- | <latex> | ||

- | \tau_{a_x, a_z} [u,v,w](x,y,z) \rightarrow [u,v,w](x+a_x/L_x, \, y, \, z+a_z/L_z) | ||

- | </latex> | ||

- | |||

- | Further symmetry operations are available but not yet documented. | ||

- | ====== Usage examples ====== | ||

- | |||

- | ==== equilibrium ==== | ||

- | Find an equilibrium from initial guess ''guess.ff'' | ||

- | <code> | ||

- | findorbit -eqb guess.ff | ||

- | </code> | ||

- | |||

- | ==== 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. Equivalent to wavespeed (cx, cz) = 1/T (ax/Lx, az/Lz). | ||

- | <code> | ||

- | findorbit -eqb -rel -ax 0.3 -az 0.4 -T 10 guess.ff | ||

- | </code> | ||

- | the ''--rel'' option adds the finite shift variables (ax, az) to the search space. | ||

- | |||

- | ==== periodic orbit ==== | ||

- | Find a periodic orbit with initial guess for period T=83 | ||

- | <code> | ||

- | findorbit -orb -T 83 guess.ff | ||

- | </code> | ||

- | |||

- | ==== 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> | ||

- | |||

- | It takes a small amount of thinking to determine the right -orb, -eqb, -rel, -ax, and -az options based on the type of | ||

- | search you want and | ||

- | |||

- | ====== 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(δ)| = δ | ||

- | * dx(δ) aligns with the gradient of the residual |DG \, dx + 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. It starts with the best δ from the previous Newton-GMRES-hookstep and increases or decreases | ||

- | δ based on the accuracy of the local quadratic model of the residual, |G(x) + DG dx(δ)|^2, compared to the actual residual |G(x + dx(δ))|^2. Each adjustment of δ takes two evaluations of G (i.e. two integrations around the periodic orbit). So the cost of the hookstep adjustment is small compared to the computation of the Krylov subspace and the Newton step. | ||

- | |||

- | 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. | ||

- | |||

- | In pseudo-code, the algorithm is this | ||

- | <code> | ||

- | start with initial guess x | ||

- | calculate residual r(x) = 1/2 |G(x)|^2, where G(x) = sigma f^t(u) - u | ||

- | |||

- | while residual r(x) is too big and # newton steps < max newton steps { | ||

- | | ||

- | // Compute Newton-Krylov step dxNewton via GMRES | ||

- | while GMRES residual > max GMRES residual { | ||

- | integrate Navier-Stokes to compute DG^n G(x) | ||

- | find least-squares solution to DG dxNewton = G(x) in n-d Krylov subspace | ||

- | GMRES residual = |DG dxNewton + G(x)| / |G(x)| | ||

- | } | ||

- | |||

- | while delta (the region of validity of linearization) is uncertain { | ||

- | compute dxHookstep as a function of delta (cheap) | ||

- | integrate Navier-Stokes to compute quadratic model of residual, |DG dxHookstep + G(x)|^2 | ||

- | integrate Navier-Stokes to compute actual residual, |G(x+dxHookstep)|^2 | ||

- | if quadratic model is | ||

- | very accurate: increase delta and try again | ||

- | inaccurate: decrease delta and try again | ||

- | alright: break | ||

- | } | ||

- | | ||

- | set x = x + dxHookstep | ||

- | } | ||

- | </code> | ||

- | | ||

- | |||

- | | ||

- | | ||

- | | ||

- | |||

- | |||

- | |||

- | ====== 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.txt · Last modified: 2010/02/02 07:55 (external edit)