#!/usr/bin/env python
# coding: utf-8

# # The unstratified Kelvin-Helmholtz instability: analytical approach
# 
# 
# #### Navid C. Constantinou
# #### RSES, ANU, 2018

# The diagrams below show the basic "flavors" of Kelvin-Helmholtz instability: unstratified on the left and stratified on the right.
# <center><img src="assets/KH-profiles.png" width=45%></center>

# Here we will study the unstratified Kelvin-Helmholtz, that is we will study the stability of basic state:
# $$
# \boldsymbol{u} = U_0\,\mathrm{sign}{(z)}\,\widehat{\boldsymbol{x}}\quad,\quad\rho = \rho_0. \label{1a}\tag{1a}
# $$
# along with pressure
# $$
# p_0(z) = P_0-\rho_0 g z.\label{1b}\tag{1b}
# $$
# Above, $P_0$ is just a constant; our reference point for measuring the pressure.
# 
# For simplicity, we consider an infinite domain.
# 
# We have seen in the previous lecture that (1a)-(1b) satisfy the inviscid Boussinesq equations.
# 
# Let's compute the vorticity $\boldsymbol{\omega}$ of the basic state:
# 
# $$
#  \boldsymbol{\omega} = \boldsymbol{\nabla}\times\boldsymbol{u} = 2U_0\,\delta(z)\,\widehat{\boldsymbol{y}}.
# $$
# 
# That is, the vorticity of the basic state is everywhere zero except from $z=0$. There we have a, so-called, "vortex sheet" with homogeneously distributed vorticity density. 

# ### Physical mechanism of the instability
# 
# The above profile is unstable as we will show (spoiler alert). But before we start doing eigenanalysis of the basic state let's try to understand the basic physical mechanism of this instability.
# 
# There are two ways to understand why this flow will produce instabilities.
# 
# #### Vorticity dynamics (Batchelor (1967))
# 
# Rembember (e.g., from Physics of Fluid Flows) that *(i)* in an inviscid fluid (like here) vortex lines are material lines, that its they move with the flow. Also, remember that each vortex line induces a rotating flow with strength proportional to the circulation around a loop enclosing it. Our vortex sheet can be thought of continuum of vortex lines aligned in $\widehat{\boldsymbol{y}}$ direction.

# Now consider we slightly perturb the vortex sheet as shown below. The vortex lines in regions like point A will be swept away to the right (since they are in the $z>0$ region) and likewise, vorticity of regions like point C will be swept away to the left. This will result in vorticity being accumulated in points like B. Thus points like B will increase their vorticity and induces a rotation counterclockwise, while points like D will decrease their vorticity and induces a rotation clockwise. But these induced rotational motions of regions B and D will result in the initial perturbation of the vortex sheet to grow even more. Thus we expect instability.
# 
# <center><img src="assets/KHmechanismVorticity.png" width=65%></center>
# 
# 
# **Question**: Sure the vorticity dynamics above demonstrate a positive feedback mechanism. But does this result in the perturbation grown exponentially as is required to have a modal instability?

# #### Bernoulli's principle and lift
# 
# The other mechanistic view of Kelvin-Helmholtz instability involves Bernoulli's principle. Consider the upper region $z>0$. After the perturbation fluid above point A will increase its $u$-velocity by a bit to account for the slight reduction in the flow's cross-section. Likewise fluid above point C will recude its $u$-velocity slightly to account for the slight increase in the flow's cross-section. (The opposite occurs in the lower region.) But then, there is a upward lift in points like A (similar to the lift that an airplane's wing feel because of the different flow speeds on either side of the wing) and, for the same reason, there is a downward lift in points like C. These lifts will result in the initial perturbation of the vortex sheet to grow even more. Thus we expect instability.
# 
# <center><img src="assets/KHmechanismLift.png" width=65%></center>

# ### Study of the instability
# 
# Now that we have an idea of what the physical process is and we know what to expect, we move on to the mathematical study of the instability.
# 
# Assume small perturbations superimposed on the basic state (1a)-(1b). For simplicity we consider only two-dimensional perturbations perturbations in $x$, $z$:
# $$
# \boldsymbol{u}(x, z, t) = \big[U_0\,\mathrm{sign}{(z)} + u'(x, z, t)\big]\,\widehat{\boldsymbol{x}} + w'(x, z, t)\,\widehat{\boldsymbol{z}},\quad p = p_0 + p',\label{2}\tag{2}
# $$
# where $|u'|, |w'| \ll U_0$, and $|p'|\ll |p_0|$.
# Perturbations are denoted with primes.

# The procedure then goes as follows:
# 1. Think about the boundary conditions the perturbation fields should obey.
# 2. Insert the perturbed fields in the Boussinesq equations, use that $u_0$, $\rho_0$, $p_0$ are actual solutions to simplify the equations, and discard any terms that involve products of primed fields (in other words, linearize the equations). At this point you should end up with something in the form:
# $$
# \partial_t \phi = \mathcal{A} \phi,
# $$
# where $\phi$ denotes collectively all perturbation fields and $\mathcal{A}$ is a differential operator that depends on the basic state $u_0$, $\rho_0$, $p_0$.
# 5. Search for eigensolutions, i.e., solutions for which all perturbation flow fields have time-dependence $e^{\sigma t}$. In other words, assume $\phi = \widetilde{\phi}\,e^{\sigma t}$. This transforms $\partial_t \phi = \mathcal{A} \phi$ into an eigenvalue problem in which the allowed $\sigma$ are the eigenvalues. If there exist eigenvalues with $\mathrm{Re}(\sigma)>0$ we have instability.

# ### 1. Boundary conditions (part 1)
# 
# Since the perturbation fields are assumed infinitesimal and our is domain infinite, it is reasonable to assume that far away from $z=0$ the flow should not "feel" the perturbations. I.e.,
# 
# $$
# |u'|, |w'|, |p'| \to 0 \quad \textrm{for } |z|\to\infty.
# $$

# ### 2. Insert in EOM ($=$equations of motion) and linearize
# 
# Our basic state is discontinuous at $z=0$. Furthermore, at $z=0$ there exist initially our vortex sheet and this sheet behaves as a material surface. Thus, the regions in either side of the vortex sheet behave as two "different" immiscible fluid layers --- fluid cannot pennetrate from the upper region to the lower. Therefore, we consider separately the $z>0$ and $z<0$ regions. 
# 
# Since we are interested in the unstratified case, we neglect the density equation and any perturbations to the density field.
# 
# We label fields in region $z>0$ with a subscript 1 and fields in region $z<0$ with subscript 2. Then from the equations of motion we have:

# For $z>0$:
# \begin{gather}
# \partial_t u'_1 + U_0 \partial_x u'_1 = -\partial_x p'_1, \label{3}\tag{3}\\
# \partial_t w'_1 + U_0 \partial_x w'_1 = \underbrace{-\partial_z p_0 - g}_{=0} -\partial_z p'_1, \label{4}\tag{4}\\
# \partial_x u'_1 + \partial_z w'_1 = 0,  \label{6}\tag{6}\\
# \end{gather}
# while for $z<0$:
# \begin{gather}
# \partial_t u'_2 - U_0 \partial_x u'_2 = -\partial_x p'_2,  \label{7}\tag{7}\\
# \partial_t w'_2 - U_0 \partial_x w'_2 = \underbrace{-\partial_z p_0 - g}_{=0} -\partial_z p'_2,  \label{8}\tag{8}\\
# \partial_x u'_2 + \partial_z w'_2 = 0.  \label{10}\tag{10}\\
# \end{gather}

# ### 1. Boundary conditions (part 2)
# 
# Now that we split our perturbation field in two regions we have to revisit boundary conditions: there must be some boundary conditions for both upper and lower region fields at $z=\eta'(x, y, t)$.
# 
# The vortex sheet moves as a material surface. This means that at the vortex sheet $z=\eta'(x, y, t)$ the component of the flow which is *normal* to the vorticity sheet in the upper and the lower region should match the normal velocity component component of the vortex sheet itself.
# 
# This may read a bit cumbersome. Let's think about it a bit more...
# 
# 
# 
# 

# #### Kinematic conditions on material surfaces
# 
# <center><img src="assets/vortexsheet.png" width=65%></center>
# 
# The condition that the surface $\eta$ remains a material surface for *all times* implies a relationship among the velocity of each point of the material surface $\boldsymbol{V}(\boldsymbol{x},t)$ and the flow adjacent to either sides of the surface. Namely, the requirement is that the component of the velocity of the material surface that is normal to the surface at each point matches the corresponding component of the fluid at either side:
# 
# \begin{align*}
# \boldsymbol{u}_1(\boldsymbol{x},t)\boldsymbol{\cdot}\widehat{\boldsymbol{\kappa}} = \boldsymbol{V}(\boldsymbol{x},t)\boldsymbol{\cdot}\widehat{\boldsymbol{\kappa}}\quad\textrm{at }z=\eta(x, y, t),\\
# \boldsymbol{u}_2(\boldsymbol{x},t)\boldsymbol{\cdot}\widehat{\boldsymbol{\kappa}} = \boldsymbol{V}(\boldsymbol{x},t)\boldsymbol{\cdot}\widehat{\boldsymbol{\kappa}}\quad\textrm{at }z=\eta(x, y,t).
# \end{align*}

# First let's find $\widehat{\boldsymbol{\kappa}}$. The surface of the vortex sheet can be described as the zero-contour level of $H$:
# 
# $$
# H(\boldsymbol{x},t) \equiv z-\eta(x, y, t) = 0.
# $$
# 
# (*Is this always true?*)
# 
# The unit normal vector $\widehat{\boldsymbol{\kappa}}$ is parallel to $\boldsymbol{\nabla}H$, thus:
# 
# $$
# \widehat{\boldsymbol{\kappa}} = \frac{\boldsymbol{\nabla}H}{\|\boldsymbol{\nabla}H\|} = \frac{-\partial_x\eta\,\widehat{\boldsymbol{x}} -\partial_y\eta\,\widehat{\boldsymbol{y}} + \widehat{\boldsymbol{z}}}{\sqrt{1+(\partial_x\eta)^2}} .
# $$

# Now let's find the velocity of the points on the vortex sheet $\boldsymbol{V}$.
# 
# A trajectory of any point $P$ is:
# 
# $$
# \boldsymbol{x}_P(t) = x_P(t)\,\widehat{\boldsymbol{x}} + y_P(t)\,\widehat{\boldsymbol{y}} + z_P(t)\,\widehat{\boldsymbol{z}}.
# $$
# 
# If further point $P$ *lies on* the material surface then for any $t$ we have that $z_P = \eta$, i.e.,
# 
# $$
# \boldsymbol{x}_P(t) \Big|_{\textrm{on }\eta} = x_P(t)\,\widehat{\boldsymbol{x}} + y_P(t)\,\widehat{\boldsymbol{y}} + \eta\big(x_P(t), y_P(t), t\big)\,\widehat{\boldsymbol{z}}.
# $$
# 
# The velocity of the points $P$ that lie on surface $\eta$ is simply
# 
# $$
# \boldsymbol{V} \equiv \frac{\mathrm{d}}{\mathrm{d}t} \boldsymbol{x}_P(t) \Big|_{\textrm{on }\eta} =  \frac{\mathrm{d}x_P(t)}{\mathrm{d}t} \widehat{\boldsymbol{x}} + \frac{\mathrm{d}y_P(t)}{\mathrm{d}t} \widehat{\boldsymbol{y}} + \frac{\mathrm{d}}{\mathrm{d}t}\eta\big(x_P(t), y_P(t), t\big)\,\widehat{\boldsymbol{z}}
# $$

# Now, for points $P$ on the surface $\eta$ the value of $H$ is *always* $H=0$ (by definition). Thus, for points on the surface $\mathrm{d}H/\mathrm{d}t = 0$. But for points $P$ on the surface we have
# 
# $$
# H\big(x_P(t), y_P(t), \eta\big(x_P(t), y_P(t), t\big) \big)
# $$
# 
# from which we deduce that:
# 
# \begin{align*}
# 0 = \frac{\mathrm{d}H}{\mathrm{d}t} &= (\partial_x H) \frac{\mathrm{d}x_P}{\mathrm{d}t} + (\partial_y H) \frac{\mathrm{d}y_P}{\mathrm{d}t} + (\partial_z H) \frac{\mathrm{d}\eta(x_P, y_P, t)}{\mathrm{d}t} + \partial_t H\\
# &=\boldsymbol{\nabla}H\boldsymbol{\cdot}\boldsymbol{V} + \partial_t H
# \end{align*}
# 
# Using that $\boldsymbol{\nabla}H = \|\boldsymbol{\nabla}H\|\widehat{\boldsymbol{\kappa}}$ we get 
# 
# $$
# \boldsymbol{V}\boldsymbol{\cdot} \widehat{\boldsymbol{\kappa}} = \frac{-\partial_t H}{\|\boldsymbol{\nabla}H\|} = \frac{\partial_t\eta}{\sqrt{1+(\partial_x\eta)^2}}.
# $$

# At last we can combine the previous to write the kinematic boundary condition:
# 
# \begin{align*}
# -u_1\, \partial_x\eta -v_1\, \partial_y\eta+ w_1 = \partial_t \eta \quad\textrm{at }z=\eta(x, y, t),\\
# -u_2\, \partial_x\eta -v_2\, \partial_y\eta+ w_2 = \partial_t \eta \quad\textrm{at }z=\eta(x, y, t).
# \end{align*}
# 
# 
# 
# We are ready now to return to our problem. In our case the material surface in quastion is the vortex sheet.

# The kinematic boundary condition for our problem, in which $\eta = \eta'(x, t)$, are:
# 
# \begin{align*}
# -(U_0+u'_1)\,\partial_x\eta' + w_1' = \partial_t\eta' \quad\textrm{at }z=\eta'(x, t),\\
# -(-U_0+u'_2)\,\partial_x\eta' + w_2' = \partial_t\eta' \quad\textrm{at }z=\eta'(x, t)
# \end{align*}
# 
# These boundary conditions can be further simplified upon linearization. Because $\eta'$, $u'$, $w'$ are infinitesimal we can drop any terms that involve products of primed quantities. Further, all fields can be evaluated at $z=0$ rather than at $z=\eta'$ since, e.g.,
# 
# \begin{align*}
# w'_1(x, y,\eta'(x,t), t)&= w'_1(x, 0, t) + \underbrace{\eta'(x, t)\big[\partial_z w'_1(x, z, t)\big]_{z=0}+ \dots}_{\textrm{products of primed fields; drop them}} \\
# & = w'_1(x, 0, t). \label{11}
# \end{align*}
# 
# Therefore our *linearized* boundary conditions at the vortex sheet are:
# 
# \begin{align*}
# -U_0\partial_x\eta' + w_1' = \partial_t\eta' \quad\textrm{at }z=0,\tag{11}\\
# +U_0\partial_x\eta' + w_2' = \partial_t\eta' \quad\textrm{at }z=0.\tag{12}
# \end{align*}

# Additionally, we have the *total* pressure of the fluid must be continuous accros the vortex sheet:
# $$
# P_0 -\rho_0 g z + p'_1(x, \eta'(x, t), t) = P_0 -\rho_0 g z + p'_2(x, \eta'(x, t), t) 
# $$
# After linearization and some cancellations:
# $$
# p'_1(x, 0, t) = p'_2(x, 0, t) \label{13}\tag{13}
# $$

# ### 3. Dynamical equations Versus constraints
# 
# Our equations (3)-(10) comprise of dynamical equations (those that involve $\partial_t$) and also constraints like the incompressibility equations (6) and (10). Our constraints essentially imply that the $u'$ and $w'$ fields are not completely indepentent but rather they have to be such so that the constraints are statisfied.
# 
# (**Note**: The situation here is very similar case as solving for the motion of a mass on an inclided plane using Lagrange multipliers.)
# 
# Often we can write our flow fields in such a way that the constrait is always satisfied. Note, that in doing so we have to write our flow fields **in the most general way possible** otherwise we might loose some of the solutions we are looking for. For example, perturbations $u'=w'=0$ satisfy the incompressibility constraints (6) and (10), but they are obviously **not** the most general ones. The solution $u'=w'=0$ would argue that there is no instability but that's not the case.
# 
# A usual way to incorporate incompressibility in a two dimensional flow is to write our flow fields in terms of a streamfunction $\psi$. That is, in each region we consider a streamfunction $\psi_j'(x, z, t)$, $j=1,2$ and write the flow as:
# $$
# u'_j = -\partial_z\psi'_j,\quad w'_j = \partial_x\psi'_j\quad\textrm{for}\quad j=1,2.
# $$
# This way (6) and (10) are trivially satisfied so we can forget about them.

# ### Recap
# Let's see where do we stand for the moment. We have:
# 
# For $z>0$:
# \begin{gather}
# -\partial_t \partial_z\psi'_1 - U_0 \partial^2_{xz} \psi'_1 = -\partial_x p'_1, \label{14}\tag{14}\\
# \partial_t \partial_x\psi'_1 + U_0 \partial^2_x \psi'_1 = -\partial_z p'_1, \label{15}\tag{15}\\
# \end{gather}
# for $z<0$:
# \begin{gather}
# -\partial_t \partial_z\psi'_2 + U_0 \partial^2_{xz}\psi'_2 = -\partial_x p'_2, \label{16}\tag{16}\\
# \partial_t \partial_x\psi'_2 - U_0 \partial_x^2 \psi'_2 = -\partial_z p'_2. \label{17}\tag{17}\\
# \end{gather}
# Also,
# \begin{gather}
# (\partial_t + U_0\partial_x)\eta' = \partial_x\psi'_1\quad\textrm{at}\quad z=0 \label{18}\tag{18}\\
# (\partial_t - U_0\partial_x)\eta' = \partial_x\psi'_2\quad\textrm{at}\quad z=0 \label{19}\tag{19}\\
# p'_1 = p'_2\quad\textrm{at}\quad z=0  \label{20}\tag{20}
# \end{gather}
# and at $|z|\to\infty$
# \begin{gather}
# |\psi'_1|, |p'_1| \to 0\quad\textrm{for}\quad z\to +\infty, \label{21}\tag{21}\\
# |\psi'_2|, |p'_2| \to 0\quad\textrm{for}\quad z\to -\infty. \label{22}\tag{22}\\
# \end{gather}

# Equations (14)-(22) form a well-posed eigenvalue problem. ("well-posed" means that we have all equations and also all boundary conditions needed.)
# 
# Can we do any further simplifications?
# 
# With a closer look in, e.g., equation (14)-(15) we can see that we can derive a relationship that involves *only* $p_1'$. After $\partial_x(14) + \partial_z(15)$ we get:
# $$
# (\partial_x^2+\partial_z^2)p_1' = 0, \label{23}\tag{23}
# $$
# and similarly for the lower region:
# $$
# (\partial_x^2+\partial_z^2)p_2' = 0. \label{24}\tag{24}
# $$
# 
# (Note that $\partial_x(14) + \partial_z(15)$ is actually the incompressibility equation.)

# ### Homogeneity
# 
# Our eigenvalue problem (14)-(21) is, as we say, homogeneous in $x$. This means that all equations and boundary conditions remain unchanged under the transformation:
# $$
# x \mapsto x+\alpha \quad \textrm{for any } \alpha.
# $$
# 
# Similarly, the eigenvalue problem is homogeneous in $t$, i.e., equations don't change under $t\mapsto t+\tau$ for any $\tau$.
# 
# These symmetries allow us to search for solutions for which all our flow fields are in the form:
# 
# $$
#  \phi' = \widehat{\phi}\,e^{i k x} e^{\sigma t}. \label{25}\tag{25}
# $$
# 
# Above, $k$ is **real** and often called $x$-wavenumber or zonal wavenumber; $\sigma$ can be complex.
# 
# #### Parenthetical note: Homogeneity
# 
# What does "homogeneity in $x$ really means"? For a discussion on the matter refer to [Notes on exponential solution in linear systems](http://nbviewer.jupyter.org/format/slides/github/navidcy/Instabilities-in-Fluids/blob/master/lectures/exponentialsolutions.ipynb#/)

# ### Back to our eigenvalue problem
# 
# #### Pressure
# 
# Let's first investigate what expansion (25) implies or pressure fields. Now (23) and (24) read:
# 
# $$
# (\partial_z^2-k^2)\widehat{p}_1 = 0,\\
# (\partial_z^2-k^2)\widehat{p}_2 = 0.
# $$
# 
# Using boundary conditions at infinity (22) and (23) we have that $\widehat{p}_1 = A e^{-kz}$ and $\widehat{p}_2 = B e^{kz}$. Last, after taking boundary condition (20) into account we get
# 
# $$
# \widehat{p}_1 = A e^{-kz} \quad\text{and}\quad \widehat{p}_2 = A e^{kz}.
# $$
# 
# The constant $A$ can be anything. Actually, since $p$ has dimensions of $(\mathrm{lenght})^2 (\mathrm{time})^{-2}$ let us redefine $A$ to be non-dimensional and $p'$:
# 
# $$
# \widehat{p}_1 = A U_0^2 e^{-kz} \quad\text{and}\quad \widehat{p}_2 = A U_0^2 e^{kz}. \label{26}\tag{26}\\
# $$
# 
# (This last step is not needed. But having variables from which you can read off their dimensions easily can help you track down potential error or typos that occur along the way.)

# #### Streamfunction
# 
# From equations (15) and (26) we have
# 
# $$
# ik (\sigma + i k U_0) \widehat{\psi}_1 = k A U_0^2 e^{-k z} \quad\textrm{for}\quad z>0,
# $$
# 
# and by evaluating it at $z=0^+$ we have
# 
# \begin{gather}
# \widehat{\psi}_1(z=0^+) = \frac{-i A U_0^2}{ \sigma + i k U_0}. \label{27}\tag{27}
# \end{gather}
# 
# Similarly, from (17) and (26):
# 
# $$
# ik (\sigma - i k U_0) \widehat{\psi}_2 = -k A U_0^2 e^{k z} \quad\textrm{for}\quad z<0,
# $$
# 
# which implies
# 
# \begin{gather}
# \widehat{\psi}_2(z=0^-) = \frac{i A U_0^2}{ \sigma - i k U_0}. \label{28}\tag{28}\\
# \end{gather}
# 
# **Note**: If instead of (15) and (17) we use (14) and (16) we arrive at exactly the very same solutions for the streamfunction in each region.

# We are almost there! We haven't used equations (18) and (19) yet... 
# 
# #### Vortex sheet displacement
# 
# Taking into account (27) and (28), equations (18) and (19) read:
# 
# \begin{align}
# (\sigma + i k U_0) \widehat{\eta} = ik  \underbrace{\frac{-i AU_0^2}{ \sigma + i k U_0}}_{=\widehat{\psi}_1(z=0^+)} \quad\textrm{at}\quad z=0, \label{29}\tag{29}\\
# (\sigma - i k U_0) \widehat{\eta} = ik  \underbrace{\frac{i AU_0^2}{ \sigma - i k U_0}}_{=\widehat{\psi}_2(z=0^-)} \quad\textrm{at}\quad z=0. \label{30}\tag{30}
# \end{align}

# But wait! How can **both** equtions (29) and (30) **hold simultaneously**. That cannot possibly be true... At least it cannot be true *for any value of $\sigma$*. Only for some specital values of $\sigma$ both  (29) and (30) can hold simultaneously. These "special values" are the *eigenvalues $\sigma$* we were after.
# 
# After some fiddling around with (29) and (30) we get:
# 
# $$
# (\sigma + i k U_0)^2 = -(\sigma - i k U_0)^2
# $$
# 
# which can be true *only* if
# 
# $$
# \sigma = \pm k U_0.
# $$
# 
# These are the eigenvalues we were looking for.

# **Exercise**: At some point in the derivation we manipulated equations (14) and (15) in such a way so that $\psi_1$ dropped out. We could have instead did $\partial_z(14) - \partial_x(15)$ and this way the pressure $p_1'$ would drop out. Repeat the eigenvalue derivation doing the latter. You should, of course, get the same eigenvalues in the end.

# **Question**: But did we cheat by searching only for two-dimensional perturbations? Two-dimensional perturbations are only a subset of three-dimensional perturbations, thus, finding instability for 2D implies that 3D are unstable as well.
# 
# There is, unfortunately, a famous theorem by Squire (1933) that states that a flow will have its largest instability for two-dimensional perturbations. In that sense not only we did *not* cheat but actually we found the largest instability possible for this problem! 
# 
# Why "unfortunately"? In the words of Jacob Bedrossian (2015),
# > I say "unfortunately" because Squire's theorem is *horrifically* misleading as it suggests all interesting aspects of hydrodynamic stability can be found in 2D equations – this is extremely false.
# 
# More discussion on the issue when we talk about transient perturbation growth and non-normality.

# ### Final Recap
# 
# After we take a few deep breaths, let us now recap again.
# 
# We showed that our eigenvalue problem (14)-(21) has solutions proportional to:
# 
# $$
# e^{i k x}e^{\pm k U_0 t}
# $$
# 
# This means that *for every choice of $k$* there is always one eigenmode that grows (the one with $\sigma = kU_0$).
# 
# Through the procsess of determining the allowed eigenvalues we also determined relations between the various perturbation flow fields. All flow fields were determined up to a multiplying factor (the constant $A$ in (26)-(30)). This is because we have linearized the perturbation equations to begin with and thus any multiple of a solution is also a solution itself.

# For the unstable mode ($\sigma=k U_0$) we find that (setting $A=1$ for simplicity):
# 
# \begin{align*}
# p'_1(x, z, t)  = U_0^2e^{i k x} e^{-k z} e^{k U_0 t}\ &,\ \ \psi'_1 = \frac{e^{-3i\pi/4}}{\sqrt{2}} \frac{U_0}{k} e^{i k x} e^{-k z} e^{k U_0 t},\\
# p'_2(x, z, t) = U_0^2e^{i k x} e^{k z} e^{k U_0 t}\ &,\ \ \psi'_2 = \frac{e^{+3i\pi/4} }{\sqrt{2}} \frac{U_0}{k} e^{i k x} e^{-k z} e^{k U_0 t}.
# \end{align*}
# 
# From the above we can compute everything else. For example, $u'$, $w'$, or the vortex sheet elevation $\eta'$. For the unstable eigenvalue $\sigma=k U_0$ we have:
# 
# \begin{align*}
# u'_1 = \frac{e^{-3i\pi/4}}{\sqrt{2}}  U_0 e^{i k x} e^{-k z} e^{k U_0 t},&\ \ w'_1 = \frac{e^{-i\pi/4}}{\sqrt{2}} U_0 e^{i k x} e^{-k z} e^{k U_0 t},\\
# u'_2 = \frac{e^{-i\pi/4}}{\sqrt{2}}  U_0e^{i k x} e^{-k z} e^{k U_0 t},&\ \ w'_2 = \frac{e^{-3i\pi/4}}{\sqrt{2}} U_0 e^{i k x} e^{-k z} e^{k U_0 t},\\
# \eta' = \frac{e^{-i\pi/2}}{2 k} &e^{i k x} e^{k U_0 t}.
# \end{align*}

# ### All these is good; but how do the eigenfunctions look like?
# 
# Let us plot these fields to see how the look.
# 
# (Don't bother with the code below --- we will discuss numerics in details in the next lecture.)

# In[1]:


import numpy as np
import matplotlib.pyplot as plt
from matplotlib import gridspec

from matplotlib import rc
rc('font', **{'family': 'serif', 'serif': ['Computer Modern'], 'size':22})
rc('text', usetex=True)
rc('xtick', labelsize=16) 
rc('ytick', labelsize=16) 
rc('axes', labelsize=20)    # fontsize of the x and y labels

nz, Lz = 200, 6.0
nx, Lx = 200, 4*np.pi
x, z = np.linspace(0, Lx, nx), np.linspace(-Lz/2, Lz/2, nz)
X, Z = np.meshgrid(x, z)


# In[2]:


U0, k = 1.0, 1.0

psicomplex = U0/(np.sqrt(2)*k) * np.exp(-k*np.abs(Z)) * np.exp(1j*k*X)
psicomplex[Z>0] = psicomplex[Z>0]*np.exp(-1j*3*np.pi/4)
psicomplex[Z<0] = psicomplex[Z<0]*np.exp(1j*3*np.pi/4)

ucomplex = U0/(np.sqrt(2)) * np.exp(-k*np.abs(Z)) * np.exp(1j*k*X)
ucomplex[Z>0] = ucomplex[Z>0]*np.exp(-1j*3*np.pi/4)
ucomplex[Z<0] = ucomplex[Z<0]*np.exp(-1j*1*np.pi/4)

wcomplex = U0/(np.sqrt(2)) * np.exp(-k*np.abs(Z)) * np.exp(1j*k*X)
wcomplex[Z>0] = wcomplex[Z>0]*np.exp(-1j*1*np.pi/4)
wcomplex[Z<0] = wcomplex[Z<0]*np.exp(-1j*3*np.pi/4)

psi = np.real(psicomplex)
u = np.real(ucomplex)
w = np.real(wcomplex)
eta = np.real( 1/(2*k) * np.exp(-1j*np.pi/2) * np.exp(1j*k*x))


# In[6]:


f, axs = plt.subplots(1, 3, figsize=(18, 5))

gs = gridspec.GridSpec(1, 3) 
axs[0], axs[1], axs[2] = plt.subplot(gs[0, 0]), plt.subplot(gs[0, 1]), plt.subplot(gs[0, 2])
axs[0].pcolormesh(x/(2*np.pi/k), k*z, u)
axs[0].set_title("shading: $u$, contours: $\psi$", fontsize=22)
axs[0].set_ylabel("$k z$")
axs[1].pcolormesh(x/(2*np.pi/k), k*z, w)
axs[1].set_title("shading: $w$, contours: $\psi$", fontsize=22);
axs[2].pcolormesh(x/(2*np.pi/k), k*z, 0*w, vmin=-1, vmax=1)
axs[2].set_title("shading: vorticity, contours: $\psi$", fontsize=22);
for i in range(3):
    axs[i].contour(x/(2*np.pi/k), k*z, psi, 11, colors="k", linewidths=1.2)
    axs[i].plot(x/(2*np.pi/k), k*eta/s, 'r', linewidth=3, label="$\eta$")
    axs[i].set_xlabel("$x/(2\pi/k)$")
    axs[i].set_ylim(-Lz*k/4, Lz*k/4)
axs[1].legend();


# ### References
# 
# - Batchelor, G. K. An introduction to fluid dynamics. Cambridge University Press, 1967.
# - Bedrossian, J. "Hydrodynamic stability and mixing at high Reynolds number: MSRI Graduate Summer school 2015." Lecture notes available online at https://www.msri.org/ckeditor_assets/attachments/306/jacobcourseiii.pdf.
#