channelflow.org

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== Continuation ======

===== 2009-2-13 =====

Continuation is finding a solution to an equation by using a solution to a nearby equation
as an initial guess. 

E.g. **poor man's continuation**: Let ''u'' be a solution to the equation ''f(u,p) = 0'' with 
parameters ''p''.  To find a solution of ''f(u',p') = 0'' with parameters ''p' = p + dp,'' 
take ''u'' as an initial guess and solve by Newton search. 

Even better is **quadratic continuation**:  Take few solutions ''u'' for a few values of ''p'', 
make a local quadratic approximation to ''u(p)'', and then extrapolate to the parameter value 
you want to get a good initial guess ''u(p')''. Of course you can do linear, cubic, etc., too.

Jonathon and I did our continuation mostly the poor man's way, with ''p'' usually being the
Reynolds number. Sometimes I did linear continuation using the ''addfields'' utility.
Jonathan wrote a quadratic continuation utility ''branchswitch.cpp'' that tracked a solution 
around a saddle-node bifurcation with Reynolds number as the continuation parameter.  We also 
continued EQBs and POs in cell size ''(Lx, Lz)''  
or equivalently (α,γ) = (2π/Lx, 2π/Lz). This is more troublesome because you have to satisfy 
the divergence constraint. Typically I would use ''changegrid'' to change the cell size and 
fix the divergence of the stretched field. The combination of procedures we used was somewhat
cumbersome and didn't help us continue around saddle-node bifurcations in spatial parameters.

Because of these limitations and because I anticipate people in the study group wanting to
continue solutions, I wrote a better algorithm that does quadratic continuation in an arbitrary
parameter and handles geometric continuation seamlessly. This allowed me to continue EQ1 and
EQ2 (Nagata LB and UB, blue in figure) around the saddle-nodes at γ=1.7ish and EQ5 (green) 
around the saddle-node at γ=2.9ish.

{{n00bsie_gamma1.png?300}} {{n00bsie_gamma2.png?340}} 

**Left** as of last fall, **Right** same fig with new continuations. The figure shows energy
dissipation D of plane Couette equilibria as a function of γ=2π/Lz, with α=2π/Lx=1.14 fixed and 
Re=400. EQ1,2 blue; EQ3,4 red; EQ5,6 green, EQ7,8 black; EQ10,11 magenta.

I have yet to add the ''continuefields'' utility to the channelflow distribution and document it.
===== 2009-03-09 =====

We really need continuation studies of our periodic orbits. It was the prominent question
at the Newton Institute meeting, and it came up again at my talks at Harvard and McGill
and in discussion with Lennaert van Veen. I can do continuation now, with current utilities, 
but it requires a lot of manual oversight. I would really like to automate the continuation
procedure or help someone else automate it.

Here's an example of the problem. We have a periodic orbit solution P19p02 at a given Reynolds 
number (Re=400) and given cell size (HKW cell). Let ''u(Re)'' be a parameterization of the solution 
as a function of Re. We can use ''u(400)'' as an initial guess for the solution at ''Re=390''
and ''Re=410'' and then initiate searches at thsoe two Reynolds numbers. E.g.

  jg356@pacemaker$ ls  Re400/
  couette.args  data/  orbstats/  orbstats.args  P19p02.ff  taubest.asc
  jg356@pacemaker$ grep T Re400/taubest.asc
  19.018545664858518  %T
  jg356@pacemaker$ mkdir Re410
  jg356@pacemaker$ cd Re410
  jg356@pacemaker$ qsubmit P19p02-Re410 findorbit -orb -T 19.018545664858518 -R 410 ../Re400/P19p02.ff

and similarly for Re400. ([[[[gtspring2009:howto:pace#submit_a_job_to_the_queue|qsubmit]] is a shell function I wrote to simplify 
submitting jobs to the PACE job queue). Assuming those searches converge, we get three points on the function ''u(Re)'' and
we can start continuing''u(Re)'' via quadratic extrapolation for initial guesses and then solving the nonlinear equations.
E.g.

  jg356@pacemaker$ mkdir Re420
  jg356@pacemaker$ continuefields 390 Re390/ubest.ff 400 Re400/P19p02.ff 410 Re410/ubest.ff 420 Re420/uguess.ff
  jg356@pacemaker$ grep T Re390/taubest.asc Re4[01]0/taubest.asc
  Re390/taubest.asc: 1.8882243932716367e+01 %T
  Re400/taubest.asc: 1.9018545664858518e+01 %T
  Re410/taubest.asc: 1.9154708055204622e+01 %T
  jg356@pacemaker$ continueparam.x 390 1.8882243932716367e+01 400 1.9018545664858518e+01 410 1.9154708055204622e+01 420
  extrapolated p(mu) == 19.290731103754673
  jg356@pacemaker$ cd Re420
  jg356@pacemaker$ qsubmit P19p02-Re420 findorbit -T 19.290731103754673 -R 420 uguess.ff

This is great, and you can even automate the procedure with a shell script, if you keep a fixed steps
in Re. But the trouble is that solution curve starts to change too rapidly and you need to take smaller 
steps in Re, or you reach a bifurcation point and have to follow the solution through a turn and go
backwards in the continuation parameter. To do this you have to switch the continuation parameter to 
the quantity on the y axis (here average dissipation ''D'') and treat both ''u'' and Re as functions 
of ''D'' be estimated by quadratic extrapolation. All this is time consuming, especially when you have
many solutions to track. So it would be nice to automate it. The [[http://indy.cs.concordia.ca/auto/|AUTO]]
software system does continuation nicely for low to moderate dimensional systems but it doesn't handle
CFD and Krylov methods. 

(a) {{p19p02continue.png?200}} (b) {{p19p02continueA.png?200}} %%(c)%% {{p19p02continueB.png?200}}

(a) continuation of P19p02 in Reynolds number (b) close-up of saddle-node bifurcation near Re=330
%%(c)%% close-up of possible bifurcation point near Re=440
===== 2009-03-11 =====

{{p19p02continueRe.png?400}}

Managed to complete the continuation of P19p02 in Re via the ''continuefields'', ''continueparams'', and ''findorbit'' utilities as shown above.  It took a lot of work (running three or four utilities on PACE for each data point), way too much to do for the whole set of periodic orbits in HKW. Last night I wrote a Matlab arclength continuation code that can branch-switch through pitchfork and saddle-node bifurcations. It'll take a little time to reorganize ''findorbit'' as a function to be called from a continuation routine, but it'll be worth it.


===== 2009-04-06 HKW orbit continuation =====

{{hkworbitsred1.png?400}} {{hkworbitsred2.png?400}}

Continuation in Re of HKW orbits, via automated ''continuesoln.cpp''. Continuation showed that 

   * P19p02 and P19p06 are upper/lower branches of same solution
   * P31p81, P41p36, and P46p23 have upper (lower) branch partners at Re=400
   * P76p85 branches off of P41p36 at about Re=275 in a symmetry-breaking, period-doubling bifurcation.

The dissipation D shown is the mean over the orbit. It would be nice to plot the mean of turbulent simulations
as well but this harder to define and compute, due to occasional decay to laminar. We need an escape rate.

PACE is still cooking on these calculations. I'm burning up all the CNS CPUs and quite a few of the other
nodes. (2 x 20ish solutions is 40ish computations). Continuation still requires manual intervention now and 
then. Anyone is interested in using this utility please let me know and I will describe its ins and outs.
channelflow.org

User Tools

Site Tools

Page Tools
