====== 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.