#!/usr/bin/env python
# coding: utf-8
# #### Problem 1: the "lift-up" solution
#
#
# Kierran presented us the Rayleigh criterion for *modal* instability and from that we can deduce that unstratified Couette flow $U(z)=\Lambda z$ is modally stable.
#
# Despite that, we'll construct a solution to the linearized perturbation equations about the Couette flow that grows.
#
# #### Step 1
#
# Consider **three dimensional** perturbations about $\boldsymbol{u} = \Lambda z\,\widehat{\boldsymbol{x}} + \boldsymbol{u}'(x, y, z, t)$. Write the linearized equations the these perturbations obey. (Don't forget incompressibility.)
#
# #### Step 2
#
# Assume two-dimensional perturbations **but** not the ones we usually have been assuming... Consider that the perturbations depend on $y$ and $z$ but *not* on $x$. Then expand all perturbation fields as:
#
# $$ \phi'(y, z, t) = \widehat{\phi}(y, t)\,e^{\mathrm{i} m z} .$$
#
# (We don't write the perturbations' time-dependance as $e^{\sigma t}$ because we are not looking for modal instabilities. One need not bother with modal instabilities for this problem since Rayleigh's necessary criterion for instability is not satisfied.)
#
# #### Step 3
# Use the incompressibility equation and show that you can express $v'$ and $w'$ in terms of a streamfunction $\psi'(y, z, t)$. Thus our three-dimensional perturbation field is now described by two fields: $u'$ and $\psi'$.
#
# #### Step 4
# By manipulating the $v'$ and $w'$ equation deduce that one possible solution has both $v'$ and $w'$ not varying with time. (Form the equation for the $x$-component of the vorticity and show that a streamfunction $\psi'(y,z)$ satisfies it.)
#
# #### Step 5
# Finally, use the $u'$ equation to show that:
#
# \begin{align*}
# u'(y, z, t) & = u'_0 - w'_0\,\Lambda\, t,\\
# \upsilon'(y, z, t) &= \upsilon'_0,\\
# w'(y, z, t) &= w'_0.
# \end{align*}
#
# where $\boldsymbol{u}'_0 = \boldsymbol{u}'(y, z, t=0)$ is the initial (infinitesimal) flow perturbation.
#
# You have just discovered a flow field that *(i)* satifies the perturbation equations and *(ii)* grows linearly with time! Note that although it grows, it does not grows *exponentially* with time.
#
# This solution is often called the "lift-up" solution. Can you think why is called like that? Try to visualize what the perturbation field is doing to enable $u'$ to grow and you'll see why...
# ## Problem 2: the Kelvin-Orr shearing wave solution
#
# Consider again the Couette flow basic state and the linearized equations that govern a perturbation flow about it. Now take two-dimensional perturbations as we used to take up to now:
#
# $$\boldsymbol{u}(x, z, t) = \big[\Lambda z+ u'(x, z, t)\big]\widehat{\boldsymbol{x}} + w'(x, z, t)\,\widehat{\boldsymbol{z}}.$$ As usual, using the incompressibility we can express $u'$ and $w'$ through a streamfunction $\psi'$.
#
# #### Step 1
#
# Write down (or derive it from the $u'$ and $w'$ linearized momentum equations) the equation satisfied by $\psi'$. This is a simplified version Rayleigh's equation with $\partial_z^2U(z)=0$. You should get:
#
# $$(\partial_t + \Lambda z\partial_x)\,\zeta'=0, \tag{1}$$
#
# where $\zeta'$ is the vorticity $\zeta'=(\boldsymbol{\nabla}\times\boldsymbol{u}')\boldsymbol{\cdot}\widehat{\mathbf{y}} = \partial_z u'-\partial_x w'= -\nabla^2\psi'$. Solving (1) implies finding either $\zeta'$ or $\psi'$, since the two are related.
#
# (OK, in the end I also wrote down Rayleigh's equations because I need it to discuss further.)
#
#
# #### Step 2
#
# How can we find a solution $\psi'$ to the above? Let's start with something that is doomed to fail: search for solutions in the form
#
# $$\psi'= A(t)\,e^{\mathrm{i} (k x + m z)}.\tag{2}$$
#
# Show that this cannot be the case by plugging this solution into (1) and arriving in an inconsistency. (E.g., you will find that the only way to satisfy (1) using (2) is when the solutions for $\psi'$ has a form which is **inconsistent** with (2).)
#
#
# #### Step 3
#
# Let's try something else. First of all, because of homogeneity in $x$ we can *always* seek for a solution that has $e^{ikx}$ $x$-dependance, e.g.:
#
# $$\zeta'(x, z, t) = \widehat{\zeta}(z, t)\,e^{\mathrm{i}kx}.$$
#
# (See notes on exponential solutions if this is not obvious by now.)
#
# Can you find a solution of (1) in the above form? What equation should function $\widehat{\zeta}(z,t)$ satisfy? If initially the vorticity perturbation is a simple harmonic function
#
# $$
# \zeta'(x, z, t=0) = \zeta_0\,e^{\mathrm{i} (k_0 x + m_0 z)}
# $$
#
# then show that
#
# $$
# \zeta'(x, z, t) = \zeta_0\,e^{\mathrm{i} [k_0 x + (m_0-k_0 \Lambda t) z]}\tag{3}
# $$
#
# solves (1).
#
# We've done it. This solves (1). Of course don't forget that our *actual* solution is either the real or the imaginary part of (3). Note that this is a harmonic in which the $z$-wavenumber changes with time! Quite different from all the modal solutions we have considered so far...
#
# #### Step 3
#
# Let's look closer at solution (3). First of all, what's the streamfunction that corresponds to vorticity solution (3)?
# Show that:
#
# $$
# \psi'(x, z, t) = \frac{\zeta_0\,e^{\mathrm{i} [k_0 x + (m_0-k_0 \Lambda t) z]}}{k_0^2 + (m_0-k_0\Lambda t)^2}.\tag{4}
# $$
#
# This solution is called the "Kelvin-Orr shearing wave". Plot contour plots of snapshots of $\psi'$ (remember to take the real part) for various values of $t$. (Perhaps even better draw a movie.) Use any domain, e.g., $(x, z)\in [-2\pi, 2\pi]\times[-2\pi, 2\pi]$ and plot the solution for various initial wavenumbers $k_0$, $m_0$ and values of $\Lambda$.
#
# Note: The purpose here is **not** simply produce as many figures possible, but rather to understand what solution (3) or (4) looks like. Please don't fill in your reports with figures. Show me only 2-3 snapshots and describe what you think the solution is doing. (But I'm sure you'll have to plot more than 2-3 snapshots to understand what the solution is doing..)
#
# #### Step 4
#
# Now let's compute the perturbation energy, $E=\frac1{2}(u'^2+w'^2)$ when averaged over a wavelength in $x$. That is, compute:
#
# $$
# E(t) = \frac1{2\pi/k}\int_0^{2\pi/k} \frac1{2} \Big( \mathrm{Re}[u'(x,z,t)]^2 + \mathrm{Re}[w'(x,z,t)]^2 \Big)\,\mathrm{d}x
# $$
#
# Hint: The calculations are much easier if you first note that the average over a wavelength $2\pi/k$ for $\cos^2(kx + \phi)$, $\sin^2(kx + \phi)$ is $\frac1{2}$.
#
# Show that:
#
# $$
# \frac{E(t)}{E(0)} = \frac{k_0^2 + m_0^2}{k_0^2 + (m_0 - k_0\Lambda t)^2}. \tag{5}
# $$
#
# #### Step 5
#
# Let's stop and think for a minute. We obtain solution (3) or (4) that satisfies the linearized perturbation equations and whose kinetic energy evolves as (5).
#
# First of all we see that (5) implies that perturbation energy **can grow**. Remember these are perturbations about the Couette profile, which has no modal instabilities. Deduce from (4) the following:
# 1. For which initial conditions do we have energy growth? (i.e., relation between $k_0$ and $m_0$)
# 2. For what time is maximum perturbation energy achieved?
# 3. What's the maximum energy growth? Can you find initial perturbations that maximize the maximum energy growth?
# 4. How does the perturbation energy behave at $t\to\infty$? Is this compatible with modal instability of Couette flow?
# ## Problem 3: best way to excite a linear system
#
# Consider the linear $2\times2$ system:
#
# $$
# \frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{x} = \mathbb{A}\,\boldsymbol{x}, \quad \boldsymbol{x}(t=0)=\boldsymbol{x}_0, \tag{6}
# $$
#
# with $\|\boldsymbol{x}_0\|^2=1$. Here, $\|\boldsymbol{x}\|^2\equiv x_1^2+x_2^2$.
#
# Imagine that $\boldsymbol{x}$ describes our perturbations about a fixed basic state. We want to investigate what's the best initial condition to achieve maximum perturbation growth when there is an instability (i.e. $\mathbb{A}$ has an eigenvalue with positive real part).
#
# Take:
#
# \begin{equation*}
# \begin{pmatrix}0 & 1 \\ 1 & 0\end{pmatrix},
# \end{equation*}
#
# which has eigenvalues $\sigma_1=1$ and $\sigma_2 = -1$ and corresponding eigenvectors $\boldsymbol{u}_1 = \frac1{\sqrt{2}}\begin{pmatrix}1\\1\end{pmatrix}$ and $\boldsymbol{u}_2 = \frac1{\sqrt{2}}\begin{pmatrix}1\\-1\end{pmatrix}$. (Make sure you can do the eigenanalysis.)
#
# Plot these eigenvectors on the unit circle. What their relation? Perhpaps if you compute $\boldsymbol{u}_1^\mathrm{T}\boldsymbol{u}_2$ it wil become crystal clear?
#
#
# The solution of (6) is:
#
# $$
# \boldsymbol{x}(t) =\underbrace{\mathrm{exp}\big(\mathbb{A}t\big)}_{2\times2 \textrm{ matrix}}\,\boldsymbol{x}_0,
# $$
#
# where the exponential of a matrix is defined through a Taylor series:
#
# $$
# \mathrm{exp}\big(\mathbb{A}t\big) = \mathbb{I} + \frac{\mathbb{A}t}{1!} + \frac{(\mathbb{A}t)^2}{2!} + \dots
# $$
#
# This gives:
#
# $$
# \mathrm{exp}\Big[ \begin{pmatrix}0 & 1 \\ 1 & 0\end{pmatrix}t\Big] = \begin{pmatrix}\cosh t & \sinh t \\ \sinh t & \cosh t\end{pmatrix}.
# $$
#
# We want to determine what is the initial condition $\boldsymbol{x}_0$ that maximizes the energy growth $E(t)/E(0)$ at time $t$. Here, energy is simply $E(t)\equiv \|\boldsymbol{x}(t)\|^2$.
#
# #### Step 1
#
# Can you think of the best initial condition so that we have maximum energy growth at time $t$? Remember that this matrix has one stable and one unstable eigenvalue...
#
# #### Step 2
#
# Let's see if we can prove what you might have thought above. Write the initial condition as $\boldsymbol{x}_0=\begin{pmatrix}\cos\theta \\ \sin\theta\end{pmatrix}$; this way $\|\boldsymbol{x}_0\|^2=1$. Find the angle $\theta$ that gives you maximum energy growth at time $t$. Plot this optimal initial condition $\boldsymbol{x}_{0,\mathrm{opt}}$ onto the same diagram as the eigenvectors; how is $\boldsymbol{x}_{0,\mathrm{opt}}$ related to $\boldsymbol{u}_{1}$ and $\boldsymbol{u}_{2}$?
#
# #### Step 3
#
# Plot the energy evolution for different initial conditions. Here's a sample of how this energy evolution should be. (The initial conditions plotted below do not necessarily correspond to the one that optimally excites the perturbations.
#
#
# 
#
#
# #### Step 4
#
# Now repeat everything for matrix:
#
# \begin{equation*}
# \begin{pmatrix}1 & 10 \\ 0 & -1\end{pmatrix}
# \end{equation*}
#
# This matrix has **the same** eigenvalues but different eigenvectors.
#
# Hint:
#
# $$
# \mathrm{exp}\Big[ \begin{pmatrix}1 & R \\ 0 & -1\end{pmatrix}t\Big] = \begin{pmatrix} e^t & R\sinh t \\
# 0 & e^{-t}
# \end{pmatrix}.
# $$
#
# A sample of energy evolutions for different initial conditions is shown below.
#
#
# 
#
#
# #### Step 5
#
# What about if the matrices were stable?
#
# \begin{equation*}
# \begin{pmatrix}
# -\frac{3}{2} & -\frac{1}{2} \\
# -\frac{1}{2} & -\frac{3}{2} \\
# \end{pmatrix}
# \quad\textrm{and}\quad \begin{pmatrix}-2 & 10 \\ 0 & -1\end{pmatrix}
# \end{equation*}
#
# both of which have eigenvalues $\sigma_1=-2$ and $\sigma_2=-1$.
# In[ ]: