User Tools

Site Tools


docs:utils:findorbit

Differences

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

Link to this comparison view

Both sides previous revision Previous revision
Next revision
Previous revision
docs:utils:findorbit [2009/02/02 11:04]
gibson
docs:utils:findorbit [2010/02/02 07:55] (current)
Line 2: Line 2:
  
 ''​findorbit'':​ Newton-Krylov-hookstep search for (relative) periodic orbit of plane Couette flow. Algorithm by  ''​findorbit'':​ Newton-Krylov-hookstep search for (relative) periodic orbit of plane Couette flow. Algorithm by 
-[[references|D. Viswanath]].+[[#references|D. Viswanath]]
 + 
 +Note: this algorithm can take hours or days to run to completion.
  
-===== usage ===== 
 <​code>​ <​code>​
 +usage:
   findorbit [-orb|-eqb] [-rel] [-T maptime] [-opts] guess.ff   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>​ </​code>​
  
-=== example 1: equilibrium ​ ===+[[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''​ Find an equilibrium from initial guess ''​guess.ff''​
 <​code>​ <​code>​
Line 15: Line 62:
 </​code>​ </​code>​
  
-=== example 2: traveling wave ===+===== traveling wave =====
 Find a traveling wave (relative equilibrium) from initial guess ''​guess.ff'',​ with initial guess for wavespeed 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)+specified in terms of finite shifts (ax/Lx, az/Lz) in finite time T. Equivalent to wavespeed ​(cx, cz) =  1/(ax/Lx, az/Lz).
 <​code>​ <​code>​
   findorbit -eqb -rel -ax 0.3  -az 0.4 -T 10 guess.ff   findorbit -eqb -rel -ax 0.3  -az 0.4 -T 10 guess.ff
 </​code>​ </​code>​
 +the ''​--rel''​ option adds the finite shift variables (ax, az) to the search space.
  
-=== example 3: periodic orbit ===+===== periodic orbit =====
 Find a periodic orbit with initial guess for period T=83 Find a periodic orbit with initial guess for period T=83
 <​code>​ <​code>​
Line 28: Line 76:
 </​code>​ </​code>​
  
-=== example 4: relative periodic orbit===+===== relative periodic orbit ====
 Find a relative periodic orbit with initial guess of a shift by Lz/4 every period T=83 Find a relative periodic orbit with initial guess of a shift by Lz/4 every period T=83
 <​code>​ <​code>​
Line 34: Line 82:
 </​code>​ </​code>​
  
-===== main options ===== +It takes half a minute'​s thinking to determine the right -orb-eqb, -rel, -axand -az options based on the type of search you want.
-<​code>​ +
-  -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>​+
  
 ====== Mathematics ====== ====== Mathematics ======
Line 145: 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 153: 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 smoothlybased 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 accurateincrease δ 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 solutionthe 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)|^2where 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 183: Line 235:
  
 Beyonce, "​Single Ladies",​ http://​www.youtube.com/​watch?​v=8mVEGfH4s5g Beyonce, "​Single Ladies",​ http://​www.youtube.com/​watch?​v=8mVEGfH4s5g
- 
docs/utils/findorbit.1233601454.txt.gz · Last modified: 2009/02/02 11:04 by gibson