User Tools

Site Tools


gtspring2009:howto:poincare

How to produce a Poincare section of plane Couette flow

around the Nagata EQ2 upper branch. This is procedure is too special-case and kludgy to put in channelflow documentation. It is a horrific mixture of general-purpose channnnelflow utilities, specialized channelflow programs, Unix utilities, and bash shell programming.

Integrate perturbations

I'll assume you have the Nagata upper-branch eqb EQ2.ff and the eigenfunctions of the complex instability, ef2.ff and ef3.ff. These fields are actually real and imaginary parts of the complex eigenfunctions, and the evolution of perturbations goes like


u(t) = \text{EQ}_2 + e^{\text{Re} \lambda t} [\text{ef}_2 \cos (\text{Im } \lambda t) - \text{ef}_3 \sin (\text{Im } \lambda t)]

The eigenfunctions are non-orthogonal so the first thing to do is to make an orthogonal basis from them. Throw in the next leading S-symmetric eigenfunctions for good measure.

 makebasis ef2 ef3 ef11 ef12

The produces e0, e1, e2, e3. The first two will span ef2, ef3.

Construct perturbations of the form


u(0) = \text{EQ}_2 + \epsilon \; \Lambda^{n/N} e_0

where Λ = exp(Re λ * 2 π / Im λ) = 6.7549 is the expansion multiplier for one period of oscillation. This will produce trajectories uniformly distributed under the iterated unstable oscillation.

I set ε = 1e-05 and started with N=16, and named my initial condition fields after the digits in Λ^(n/N) (using digits rather than integer labels will scale if I later increase N to 32 or 64). E.g Λ^(n/N) for n/N = 1/16 is 1.1268, and for 2/16 is 12696.

 addfields 1 EQ2.ff 1.1267e-05 e0.ff eq2_11268e0.ff
 addfields 1 EQ2.ff 1.2696e-05 e0.ff eq2_12696e0.ff
 ...

Integrate these 16 fields for a few hundred time units and save.

 couette -T0 0 -T1 400 -o data-11268 eq2_11268e0.ff
 couette -T0 0 -T1 400 -o data-12696 eq2_12696e0.ff
 ...

Instead of typing each of these out, you can use a bash for-loop,

 for i in eq2_*e0.ff ; do tag=${i#eq2_} ; couette -T0 0 -T1 400 -o data-${tag%.e0.ff} $i ; done

The ${...} stuff is bash string manipulation syntax to extract the numerical part of the input file names.
This produces data directories data-12696, etc. I made symbolic links to these directories with simpler 
labels to make some future processing simpler.

   ln -s data-11268 data-a
   ln -s data-12696 data-b

===== Compute the Poincare crossings =====

{{gtspring2009:howto:poincare:eq2poincare.cpp}} is a special program I wrote to compute crossings of a Poincare section defined by

  (u(t) - EQ2, e(θ)) == 0

where e(θ) = e0 cos θ + e1 sin θ, and where u(t) is always mapped into a canonical 1st quadrant defined by
(u(t), etx) ≥ 0 and (u(t), etz) ≥ 0. Here etx and etz are τx and τz antisymmetric basis vectors. In my calculations
I chose these to be the τx and τz antisymmetric basis vectors of the EQ2 translation basis described in our 2008 
JFM paper. 

To compile eq2poincare.cpp, use this {{gtspring2009:howto:poincare:makefile.txt}}. Edit the Makefile to so that CHANNELDIR
is set to you channelflow installation

  CHANNELDIR = /home/gibson/channelflow-1.3.5

then run "make eq2poincare.x". Then execute 

  eq2poincare.x -d data-a -o section-pi4 --theta 0.7854 -tag a -T0 0 -T1 400 etx etz e0 e1 EQ2.ff

That will produce files uM0a.ff, uP0a.ff, uM1a.ff, uP1a.ff in directory section-pi4/. The P/M indicates
whether (u(t) - EQ2, e(θ)) is increasing (P) or decreasing (M) at the crossing. The following integer 
indicates the number of the crossing (incremented once per M,P pair), and the 'a' is a label indicating
which trajectory the crossing came from. 

The point of this crazy labeling scheme is that it gives the right lexical ordering to the filenames of 
multiple crossings of multiple trajectories. For example, after computing the crossings of all trajectories
using the bash-for loop

   for i in data-[a-p] ; do tag=${i#data-} ; eq2poincare.x -d data-$tag -o section-pi4 --theta 0.7854 -tag ${tag} -T0 0 -T1 400 etx etz e0 e1 EQ2.ff; done

The files in directory section-pi4 will be ordered in increasing distance along the unstable manifold.

 gibson@tansen$ ls section-foo-pi4/
 crossingtimes-a.asc  uM0d.ff  uM1g.ff  uM3a.ff  uM4d.ff  uP0h.ff  uP2b.ff  uP3e.ff
 crossingtimes-b.asc  uM0e.ff  uM1h.ff  uM3b.ff  uM4e.ff  uP0i.ff  uP2c.ff  uP3f.ff
 crossingtimes-c.asc  uM0f.ff  uM1i.ff  uM3c.ff  uM4f.ff  uP1a.ff  uP2d.ff  uP3g.ff
 crossingtimes-d.asc  uM0g.ff  uM2a.ff  uM3d.ff  uM4g.ff  uP1b.ff  uP2e.ff  uP3h.ff
 crossingtimes-e.asc  uM0h.ff  uM2b.ff  uM3e.ff  uM4h.ff  uP1c.ff  uP2f.ff  uP4a.ff
 crossingtimes-f.asc  uM0i.ff  uM2c.ff  uM3f.ff  uP0a.ff  uP1d.ff  uP2g.ff  uP4b.ff
 crossingtimes-g.asc  uM1a.ff  uM2d.ff  uM3g.ff  uP0b.ff  uP1e.ff  uP2h.ff  uP4c.ff
 crossingtimes-h.asc  uM1b.ff  uM2e.ff  uM3h.ff  uP0c.ff  uP1f.ff  uP2i.ff  uP4d.ff
 crossingtimes-i.asc  uM1c.ff  uM2f.ff  uM3i.ff  uP0d.ff  uP1g.ff  uP3a.ff  uP4e.ff
 uM0a.ff              uM1d.ff  uM2g.ff  uM4a.ff  uP0e.ff  uP1h.ff  uP3b.ff  uP4f.ff
 uM0b.ff              uM1e.ff  uM2h.ff  uM4b.ff  uP0f.ff  uP1i.ff  uP3c.ff  uP4g.ff
 uM0c.ff              uM1f.ff  uM2i.ff  uM4c.ff  uP0g.ff  uP2a.ff  uP3d.ff  uP4h.ff

etc.

Project the crossing points

Next we project the crossing points of the Poincare section using the channelflow projectfields utility.

projectfields -b basis-EQ2ef/ -Nb 4 -or EQ2.ff -o portrait-pi4 section-pi4/*.ff

The above command assumes you have put the basis vectors e0,e1,e2,e3 into a directory named basis-EQ2ef/; if not you might use '-b .' instead.

This creates a new directory 'portrait-pi4' containing ASCII files uM0a.asc, uM0b.asc, …, each of which contains four numbers corresponding to the projection of the fields in uM0a.ff, etc onto the four basis vectors e0,e1,e2,e3. E.g.

gibson@tansen$ cat portrait-pi4/uM0a.asc
% e0 e1 e2 e3
-9.6755005579877134e-06 9.6754650183359171e-06 3.8296888337351566e-10 -2.5429703558485441e-08

It's easier to plot this data if it's in one file and in the right order. To do that, we strip out the comments with grep and save the results to two files.

gibson@tansen$ grep -v -h ^% portrait-foo-pi4/uM*.asc > uM.asc
gibson@tansen$ grep -v -h ^% portrait-foo-pi4/uP*.asc > uP.asc

That produces two files of incoming/outgoing crossings that can be plotted with matlab or whatever.

gtspring2009/howto/poincare.txt · Last modified: 2010/02/02 07:55 (external edit)