This shows you the differences between two versions of the page.
Both sides previous revision Previous revision Next revision | Previous revision | ||
docs:utils:findorbit [2009/02/02 10:39] gibson |
docs:utils:findorbit [2010/02/02 07:55] (current) |
||
---|---|---|---|
Line 1: | Line 1: | ||
====== findorbit ====== | ====== findorbit ====== | ||
- | ''findorbit'' calculates equilibria, traveling waves, and periodic orbits of | + | ''findorbit'': Newton-Krylov-hookstep search for (relative) periodic orbit of plane Couette flow. Algorithm by |
- | plane Couette flow using a Newton-GMRES-hookstep algorithm developed by | + | [[#references|D. Viswanath]]. |
- | Divakar Viswanath, U Michigan. | + | |
+ | 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 | ||
+ | -sigma --sigma <string> file containing symmetry of relative eqb/orb | ||
+ | -symms --symmetries <string> constrain u(t) to invariant subspace generated by symmetries in listed file, arg is filename | ||
+ | -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 finite-time phase shift variables ''ax'' and ''az'' to the search space | ||
+ | * -sigma <filename> sets the initial guess (if -rel) or fixed value (if not -rel) of the symmetry \sigma | ||
+ | * -symms <filename> constrains the search to an invariant subspace. The file is a list of generators for the isotropy group. | ||
+ | |||
+ | In all cases we search for a solution of | ||
+ | |||
+ | <latex> | ||
+ | \sigma f^T(u) - u = 0 | ||
+ | </latex> | ||
+ | |||
+ | where | ||
+ | |||
+ | <latex> | ||
+ | \sigma [u,v,w](x,y,z) \rightarrow s [s_x u, \, s_y v, \, s_z w](s_x x + a_x/L_x, \, s_y y, \, s_z z + a_z/L_z) | ||
+ | </latex> | ||
+ | |||
+ | (will continue later today.... /2009-06-08 jfg) | ||
+ | |||
+ | ====== 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 half a minute's thinking to determine the right -orb, -eqb, -rel, -ax, and -az options based on the type of search you want. | ||
====== Mathematics ====== | ====== Mathematics ====== | ||
Line 94: | Line 173: | ||
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 | 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(δ)| = δ |
+ | * dx(δ) aligns with the gradient of the residual |DG \, dx + G(x)|^2 as δ → 0. | ||
* dx(δ) equals the Newton step if δ ≥ |Newton step| | * dx(δ) equals the Newton step if δ ≥ |Newton step| | ||
* dx(δ) varies smoothly between those values as a function of δ | * dx(δ) varies smoothly between those values as a function of δ | ||
Line 102: | Line 182: | ||
Of course, you don't know what the best radius δ for the trust region. The hookstep algorithm uses a few | Of course, you don't know what the best radius δ for the trust region. The hookstep algorithm uses a few | ||
- | heuristics | + | 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. | ||
- | * start with the best δ from the previous Newton-GMRES-hookstep | + | 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. |
- | * 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. | + | In pseudo-code, the algorithm is this |
- | + | <code> | |
- | 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 | + | start with initial guess x |
- | Newton algorithm. | + | 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> | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
Line 132: | Line 235: | ||
Beyonce, "Single Ladies", http://www.youtube.com/watch?v=8mVEGfH4s5g | Beyonce, "Single Ladies", http://www.youtube.com/watch?v=8mVEGfH4s5g | ||
- |