Differences

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

--- docs:tutorials:findsoln [2014/12/04 11:47]
gibson [periodic orbit of plane Couette flow]
+++ docs:tutorials:findsoln [2015/08/28 06:52] (current)
gibson
@@ Line 5: / Line 5: @@
 ===== periodic orbit of plane Couette flow =====
-. Generate a random velocity field with ''randomfield.''
+=== 1. Generate a random velocity field with ''randomfield.'' ===
 <code>
@@ Line 13: / Line 14: @@
 This will construct a random velocity field in file ''u0.h5'', in a  Lx x [-1, 1] x Lz periodic box with Lx = 2 lx pi = 1.75pi and Lz = 2 lz pi = 1.2 pi (the plane Couette minimal flow unit of Hamilton, Kim, Waleffe JFM 1995) with a 32 x 33 x 32 collocation grid (which is a lower resolution than I like for these calculations, but will result in a fast-running demo.) The velocity field u0 will be incompressible and have u=0 BCs at the wall, and will have an L2Norm of 0.50. You can check this by running ''fieldprops -n u0''. I chose this magnitude for the field after testing some smaller fields and seeing that they died to laminar pretty quickly.
-. Integrate the random field forward in time.
+=== 2. Integrate the random field forward in time. ===
 <code>
-couette -T1 1000 -dt 0.05 -R 400 -symms sxyz-sxytxz.asc u0.h5
+couette -T1 1100 -dt 0.05 -R 400 -symms sxyz-sxytxz.asc u0.h5
 </code>
-This will time-integrate the field u0 as a perturbation on top of a laminar base flow U(y) = y for 1000 time units, at Re=400, and with time integration step dt=0.05, and enforcing symmetries in the velocity field specified by the file ''sxyz-sxytxz.asc'' By default the velocity field will be saved into a ''data'' directory at intervals dT=1. See [[docs:tutorials:integration]] for more information on time integration, such as the base flow plus fluctuation decomposition, or changing from plane Couette to channel conditions.
+This will time-integrate the field u0 as a perturbation on top of a laminar base flow U(y) = y for 1100 time units, at Re=400, and with time integration step dt=0.05, and enforcing symmetries in the velocity field specified by the file ''sxyz-sxytxz.asc'' By default the velocity field will be saved into a ''data'' directory at intervals dT=1. See [[docs:tutorials:integration]] for more information on time integration, such as the base flow plus fluctuation decomposition, or changing from plane Couette to channel conditions.
 The symmetry file ''sxyz-sxytxz.asc'' has contents
@@ Line 28: / Line 30: @@
 </code>
 That specifies a symmetry group with two generators
 \begin{eqnarray*}
- &\\
   \sigma_{xyz} : [u,v,w] (x,y,z) &\rightarrow [-u,-v,-w] (-x,-y,-z) \\
-  \sigma_{xy} \tau_{xz} : [u,v,w] (x,y,z) &\rightarrow [-u,v,-w] (-x+L_x/2,y,-z+L_z/2)
+  \sigma_{xy} \tau_{xz} : [u,v,w] (x,y,z) &\rightarrow [-u,-v,w] (-x+L_x/2,-y,z+L_z/2)
 \end{eqnarray*}
 The pointwise inversion $\sigma_{xyz}$ of this group fixes the origin and prevents traveling waves
-and arbitrary relative periodic orbits. For more on symmetry groups of plane Couette flow see [[docs:math:symmetry]] and {{gibson:publications:gibson_jfm09b.pdf | Gibson, Halcrow, Cvitanovic JFM 2009}.
+and arbitrary relative periodic orbits. For more on symmetry groups of plane Couette flow see [[docs:math:symmetry]] and {{:gibson:publications:gibson_jfm09b.pdf | Gibson, Halcrow, Cvitanovic JFM 2009}}.
+=== 3. Select a good initial guess from a recurrence plot. ===
+The $\langle \sigma_{xyz}, \sigma_{xy} \tau_{xz} \rangle$ symmetry group allows periodic orbits of the form $\tau u - f^T(u) = 0$ for $\tau \in \{i, \tau_x, \tau_z, \tau_{xz} \}$, where $f^t$ is the forward-time evolution under Navier-Stokes, i.e. $u(t+T) = f^T(u(t))$, and where $\tau_x$ represents a half-box shift in the $x$ direction, etc.
+For this demo I'll just look for a solution of the form $u - f^T(u) = 0$. Thus local minima of $\| u(t) - u(t+T) \|$ should provide good initial guesses for the search.
+<code>
+seriesdist -T0 0 -T1 1000 -dT 1 -tmax 100 -da data/ -db data/ -as -kx 8 -kz 8
+</code>
+This reads in the saved velocity fields from the ''data'' directory, computes $\| u(t) - u(t+T) \|$ as a function of $t,T$, and saves this in the file ''dist.asc''. Here's the recurrence plot for the data in question.
+{{:docs:tutorials:2014-12-04-a.png?direct&600|}}
+You can see strikingly periodic behavior over the range $825 \leq t \leq 1000$ with nice horizontal blue streak at $T \approx 65$. This suggests the turbulent trajectory is shadowing a periodic orbit with period $T \approx 65$. The minimum of $\| u(t) - u(t+T) \|$ in this region occurs at $t=917$ and $T=63$. That's an unusually promising initial guess for a periodic orbit.
+=== 4. Find the periodic orbit with ''findsoln'' ===
+<code>
+mkdir findsoln-917-63
+cd findsoln-917-63
+findsoln -orb -T 63 -dt 0.05 -R 400 -symms ../sxyz_sxytxz.asc ../data/u917.h5
+</code>
+Since I tend to run many searches for many initial guesses, I like to do each search in a subdirectory
+named after the initial guess. The ''findsoln'' command runs a Newton-Krylov-hookstep search that finds a $u,T$ solution of the equation $u - f^T(u) = 0$ given the initial guess $T=63$ and $u = u(917)$. The search is restricted to the $\langle \sigma_{xyz}, \sigma_{xy} \tau_{xz} \rangle$ symmetry group. The symmetry restriction vastly reduces the search space and results in a faster and more robust search.
+ After about half an hour the search succeeds, finding a solution that satisfies the equation numerically to $O(10^{-13})$. The characteristics of the search routine are recorded in the file ''convergence.asc''
+<file>
+%-L2Norm(G)   r             delta         dT            L2Norm(du)    L2Norm(u)     L2Norm(dxN)   dxNalign      L2Norm(dxH)   GMRESerr      ftotal    fnewt     fhook     CFL
+.00888606    2.20611e-05   0.01          0             0             0.390225      0             0             0             0             1         0         1         0.609524
+.00663114    1.19704e-05   0.01          -0.0640999    0.0116139     0.38873       0.00902435    0             0.00902435    0.00076865    14        12        2         0.609524
+.000855143   2.14948e-07   0.01          0.190862      0.00429167    0.388192      0.0548914     -0.238276     0.01          0.000119926   29        25        4         0.608904
+.000597692   9.06641e-08   0.01          0.194443      0.00326027    0.387135      0.035679      0.998803      0.01          0.000613148   42        37        5         0.581667
+.000505451   6.2892e-08    0.01          0.18166       0.00364325    0.386258      0.00942818    0.997453      0.00942818    0.000289388   55        49        6         0.583459
+.77901e-05   8.04837e-11   0.01          -0.0438562    0.000824778   0.386507      0.00227387    -0.99388      0.00227387    0.000612772   67        60        7         0.604789
+.46682e-09   8.266e-18     0.01          -0.000571434  1.73752e-05   0.38651       3.13122e-05   0.981989      3.13122e-05   5.95861e-05   80        72        8         0.584728
+.51125e-13   1.58971e-25   0.01          -2.67003e-07  1.49819e-08   0.38651       1.70279e-08   0.940857      1.70279e-08   0.000134271   93        84        9         0.584723
+.72014e-14   2.358e-27     0.01          -1.08282e-10  2.94411e-12   0.38651       5.77394e-12   0.921371      5.77394e-12   0.000331508   106       96        10        0.584723
+</file>
+The solution $u$ is stored in the file ''ubest.h5'' and the value of $T$ in ''Tbest.asc.''
-. Compute

channelflow.org

User Tools

Site Tools

Differences

Page Tools