This is an old revision of the document!
====== 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?300}} 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.
