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

# Here we discus how material line elements evolve in a fluid and further we show that in an inviscid fluid vorticity lines move with the flow, that is vorticity lines move **as if they were** material lines.

# ## How material line elements move in a fluid?
# 
# <center><img src="assets/materialline.png" width=45%></center>
# 
# Initially at time $t$ our material line element is given as:
# 
# $$
# \delta\boldsymbol{\ell} = \boldsymbol{x}_B-\boldsymbol{x}_A. \tag{1}
# $$
# 
# 
# We are interested to find the rate of change of $\delta\boldsymbol{\ell}$ *as we follow it with the flow*. This is Lagrangian thinking.
# 
# In a Lagrangian framework we follow the fluid parcels as they move with the flow and compute in this manner the rate of change of any quantity we are interested *along with the flow*. In an Eulerian framework, stay put at a fixed point in space and we compute the rate of change of a quantity as fluid is being advected. We can always translate form one framework to the other.
# 
# For our calculation here, it is much easier to compute the rate of change of $\delta\boldsymbol{\ell}$ in the Lagrangian framework. We follow points $A$ and $B$ for time $\delta \tau$. At $t=t+\delta t$, the points $A$ and $B$ have moved to points $A'$ and $B'$ respectively. Thus, our material line at $t=t+\delta\tau$ is:
# 
# \begin{align*}
# \delta\boldsymbol{\ell}' & = \boldsymbol{x}_{B'}-\boldsymbol{x}_{A'} \\
# & = \big\{\boldsymbol{x}_B + \delta\tau\,\boldsymbol{u}(\boldsymbol{x}_B,t) + \mathcal{O}\big[(\delta\tau)^2\big]\big\} -\big\{\boldsymbol{x}_A + \delta\tau\,\boldsymbol{u}(\boldsymbol{x}_A,t) + \mathcal{O}\big[(\delta\tau)^2\big]\big\}\\
# & = \delta\boldsymbol{\ell} + \delta\tau\big[\boldsymbol{u}(\boldsymbol{x}_B,t) - \boldsymbol{u}(\boldsymbol{x}_A,t)\big] + \mathcal{O}\big[(\delta\tau)^2\big].  \tag{2}
# \end{align*}
# 
# But, since points $A$ and $B$ are very close to each other we can write:
# 
# $$
# \boldsymbol{u}(\boldsymbol{x}_B,t) = \boldsymbol{u}(\boldsymbol{x}_A,t) + (\underbrace{\boldsymbol{x}_B-\boldsymbol{x}_A}_{=\delta\boldsymbol{\ell}})\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}(\boldsymbol{x}_A,t) + \mathcal{O}\big[(\delta\boldsymbol{\ell})^2\big], \tag{3}
# $$
# 
# which implies:
# 
# \begin{align*}
# \delta\boldsymbol{\ell}' & = \delta\boldsymbol{\ell} + \delta\tau\,\delta\boldsymbol{\ell}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}(\boldsymbol{x}_A,t) + \mathcal{O}\big[(\delta\tau)^2, \delta \tau\,(\delta\boldsymbol{\ell})^2\big]. \tag{4}
# \end{align*}
# 
# From (4) we can compute the Lagrangian rate of change of the infinitesimal material element $\delta\boldsymbol{\ell}$. That is the rate of change $\delta\boldsymbol{\ell}$ as we follow it with the flow:
# 
# $$
# \left.\frac{\mathrm{d}}{\mathrm{d}t}\delta\boldsymbol{\ell}\right|_{\mathrm{L}} = \lim_{\delta\tau\to0} \frac{\delta\boldsymbol{\ell}' - \delta\boldsymbol{\ell}}{\delta\tau} = \delta\boldsymbol{\ell}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}(\boldsymbol{x}_A,t) + \mathcal{O}\big[(\delta\boldsymbol{\ell})^2\big]. \tag{5}
# $$
# 
# Subscript $\mathrm{L}$ in equation (5) denotes that the rate of change is computed in the Lagrangian framework.
# 
# 
# For infinitesimal line elements we can neglect terms of order $(\delta\boldsymbol{\ell})^2$ and thus (5) simplifies to:
# 
# $$
# \left.\frac{\mathrm{d}}{\mathrm{d}t}\delta\boldsymbol{\ell} \right|_{\mathrm{L}}= \delta\boldsymbol{\ell} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}, \tag{6}
# $$
# 
# with both $\delta\boldsymbol{\ell}$ and $\boldsymbol{u}$ evaluated at $\boldsymbol{x},t$.
# 
# How can we tranform the above to our usual beloved Eulerian frame? The time-derivative $\mathrm{d}/\mathrm{d}t\big|_{\mathrm{L}}$ becomes in the Eulerian framework the material derivative:
# 
# $$
# \big(\partial_t + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\big) \delta\boldsymbol{\ell}  = \delta\boldsymbol{\ell} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}. \tag{7}
# $$

# ## Vorticity equation for an inviscid fluid
# 
# Let's take the simplest case of a fluid with homogeneous density, $\rho=\rho_m=\textrm{const}$. Then the mass conservation equation $\partial_t\rho + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho\,\boldsymbol{u}) = 0$ simplifies to
# 
# $$
#  \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0. \tag{8}
# $$
# 
# The momentum equatio
# n for a fluid under the influence of gravity is then:
# 
# $$
# \partial_t \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} = -\frac{\boldsymbol{\nabla}p}{\rho_m} + \rho_m g\widehat{\boldsymbol{z}}.\tag{9}
# $$
# 
# Note that because *(i)* gravity is a conservative force (i.e., it derives from a potential) and *(ii)* $\rho_m=\textrm{const.}$, we can write the right-hand-side of the above as
# 
# $$
# -\boldsymbol{\nabla} \Big( \frac{p}{\rho_m} - \rho_m g z \Big).\tag{10}
# $$
# 
# The vorticity $\boldsymbol{\omega}$ of the fluid is defined as:
# 
# $$
# \boldsymbol{\omega} =\boldsymbol{\nabla} \times \boldsymbol{u}.\tag{11}
# $$
# 
# We would like to find how the vorticity evolves. To do so we take the curl of (9) and use (11) to get:
# 
# $$
# \partial_t \boldsymbol{\omega} + \boldsymbol{\nabla}\times\big(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}\big) = 0.\tag{12}
# $$
# 
# The term $\boldsymbol{\nabla}\times\big(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}\big)$ involves the flow field $\boldsymbol{u}$ as well as spatial derivatives of $\boldsymbol{u}$. After doing some (rather elaborate) vector calculus and use of equation (8) we can rewrite 
# 
# $$
#  \boldsymbol{\nabla}\times\big(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}\big) = \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{\omega} - \boldsymbol{\omega}\boldsymbol{\cdot}\boldsymbol{\nabla}\,\boldsymbol{u}.\tag{13}
# $$
# 
# (The steps to obtain (13) are described below. But if you find those difficult or confusing then just take equation (13) as an identity that holds whenever (8) is true.)
# 
# Therefore, the vorticity equation for a homogeneous inviscid fluid becomes:
# 
# $$
# \partial_t \boldsymbol{\omega} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{\omega} = \boldsymbol{\omega}\boldsymbol{\cdot}\boldsymbol{\nabla}\,\boldsymbol{u},\tag{14}
# $$
# 
# or more compactly written,
# 
# $$
# \frac{\mathrm{D}  \boldsymbol{\omega}}{\mathrm{D}t} = \boldsymbol{\omega}\boldsymbol{\cdot}\boldsymbol{\nabla}\,\boldsymbol{u}.\tag{15}
# $$
# 
# 
# **Exercise**: Assume the flow $\boldsymbol{u}$ is two-dimensional. What does that imply about the direction of $\boldsymbol{\omega}$? In this case, what does equation (15) imply?

# ## Vorticity lines move with the flow in an inviscid fluid
# 
# The resemblance of equations (7) and (14) demonstrates the above. Vorticity lines and material lines evolve in the same way in an inviscid fluid.
# 
# **This is an extremely powerful and useful result.**
# 
# 
# **Exercise**: Does this argument follow on in a viscous fluid? Consider a fluid with viscosity $\nu$. Does equation (7) change? Does equation (14) change?

# ***
# ### Steps to derive equation (13)
# 
# First let's bring this vector identity to our short-term memory:
# 
# $$
# \boldsymbol{a} \times (\boldsymbol{b} \times \boldsymbol{c}) = (\boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{c})\,\boldsymbol{b} - (\boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{b})\,\boldsymbol{c}.\tag{16}
# $$
# 
# Does the identity above make sense? Definitely the vector on the left-hand-side will be lying on the plane which is *(i)* normal to $\boldsymbol{a}$ and *(ii)*  normal to $(\boldsymbol{b}\boldsymbol{\times}\boldsymbol{c})$. But *(ii)* implies that $\boldsymbol{a} \times (\boldsymbol{b} \times \boldsymbol{c})$ will by lying on the plane that is spanned by vectors $\boldsymbol{b}$ and $\boldsymbol{c}$ and this is exactly what (16) means.
# 
# Identity (16) is often referred to as the "$a$-$b$-$c$ identity". We can derive (16) easily using index notation:
# 
# \begin{align*}
# [ \boldsymbol{a} \times (\boldsymbol{b} \times \boldsymbol{c}) ]_i &= \epsilon_{ijk} a_j [ \boldsymbol{b} \times \boldsymbol{c} ]_k \\
# &= \epsilon_{ijk} a_j \epsilon_{klm} b_l c_m \\
# &= \big(\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl}\big) a_j  b_l c_m \\
# &= a_j  b_i c_j  - a_j  b_j c_i \\
# &= (\boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{c}) b_i  - (\boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{b}) c_i
# \end{align*}
# 
# (We can also derive it using geometry if we consider two cases: first take $\boldsymbol{a}$ to be perpendicular to $\boldsymbol{b}$ and then $\boldsymbol{a}$ to be perpendicular to $\boldsymbol{c}$.)
# 
# So can we use the $a$-$b$-$c$ identity to compute $\boldsymbol{u} \times  \boldsymbol{\omega}$?
# 
# Yes, ***but*** we have to be careful on where the $\boldsymbol{\nabla}$ acts on since the differentiation does not commute with the vectors (i.e., $f\partial_x g \ne \partial_x (fg)$). For example, we have $\boldsymbol{a}\mapsto\boldsymbol{u}$, $\boldsymbol{b}\mapsto\boldsymbol{\nabla}$, and $\boldsymbol{c}\mapsto\boldsymbol{u}$. The term $(\boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{b})\,\boldsymbol{c}$ is rather trivial:
# 
# \begin{align*}
# (\boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{b})\,\boldsymbol{c} \mapsto (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\,\boldsymbol{u},\tag{17}
# \end{align*}
# 
# but the term $(\boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{c})\,\boldsymbol{b}$ is tricky. If we just replace everything we have
# 
# \begin{align*}
# (\boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{c})\,\boldsymbol{b} \mapsto (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{u})\,\boldsymbol{\nabla}.\tag{18}
# \end{align*}
# 
# But (18) makes no sense since $\boldsymbol{\nabla}$ is left acting on thin air that is on its right! What we should have written is 
# 
# \begin{align*}
# (\boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{c})\,\boldsymbol{b} \mapsto \boldsymbol{\nabla} (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{u}).\tag{19}
# \end{align*}
# 
# But wait again! Now $\boldsymbol{\nabla}$ acts on **both** vectors $\boldsymbol{u}$ (while initially it was acting only on one of them)! This is wrong... The correct way to do that is
# 
# \begin{align*}
# (\boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{c})\,\boldsymbol{b} \mapsto \tfrac1{2} \boldsymbol{\nabla} (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{u}),\tag{20}
# \end{align*}
# 
# which gives:
# 
# \begin{align*}
# \boldsymbol{u} \times  \boldsymbol{\omega} = \boldsymbol{\nabla}\big(\tfrac1{2}\left\|\boldsymbol{u}\right\|^2\big) - (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\,\boldsymbol{u}.\tag{21}
# \end{align*}
# 
# Rather mind boggling, right?!
# 
# I find it ***much*** easier to do these kind of calculations with index notation. Then everything is much clearer; we don't have to remember e.g. on which fields the differential operators were acting on. We just have to be careful with putting parenthesis at the right places. 
# 
# \begin{align*}
# [\boldsymbol{u} \times  \boldsymbol{\omega}]_i &= \big[ \boldsymbol{u} \times (\boldsymbol{\nabla} \times \boldsymbol{u}) \big]_i \\
# &= \epsilon_{ijk} u_j [ \boldsymbol{\nabla} \times \boldsymbol{u} ]_k \\
# &= \epsilon_{ijk} u_j \epsilon_{klm} \partial_{l}\, u_m \\
# &= \big(\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl}\big) u_j  \partial_{l}\, u_m \\
# &= u_j  \partial_i u_j  - u_j \partial_j u_i \\
# &= \partial_i \big(\tfrac1{2}\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{u}\big)   - (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) u_i.
# \end{align*}
# 
# Same as (21) in a straightforward manner! (Index notation and Einstein summation are great tools once you get the grips out of them.)
# 
# 
# The first step towards getting equation (13) is to use (21) to rewrite the advection term as:
# 
# $$
#  \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}  =  \boldsymbol{\nabla}\big(\tfrac1{2}\left\|\boldsymbol{u}\right\|^2\big) + \boldsymbol{\omega} \times  \boldsymbol{u}.\tag{22}
# $$
# 
# Next we want to compute the curl of (22). The curl of $\boldsymbol{\nabla}\big(\tfrac1{2}\left\|\boldsymbol{u}\right\|^2\big)$ vanishes; that's good. The real challenge is the term $\boldsymbol{\nabla}\times (\boldsymbol{\omega} \times \boldsymbol{u} )$. 
# 
# If we use again the $a$-$b$-$c$ identity without much thinking, we get:
# 
# \begin{align*}
# \boldsymbol{\nabla}\times (\boldsymbol{\omega} \times \boldsymbol{u} ) & = (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u})\,\boldsymbol{\omega} - (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\omega})\,\boldsymbol{u}.\tag{23}
# \end{align*}
# 
# **Wrong!** Why? Same reason as before; the $\boldsymbol{\nabla}$ operator on the left-hand-side now acts on ***both*** $\boldsymbol{\omega}$ and $\boldsymbol{u}$ and on the right-hand-side it only acts on one of them.
# 
# All these is very confusing (for me) so I'll do it with index notation to be safe!
# 
# \begin{align*}
# [\boldsymbol{\nabla}\times (\boldsymbol{\omega} \times \boldsymbol{u} )]_i &= \epsilon_{ijk} \partial_j [ \boldsymbol{\omega} \times \boldsymbol{u} ]_k \\
# &= \epsilon_{ijk} \partial_j \epsilon_{klm} \omega_l\, u_m \\
# &= \big(\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl}\big) \partial_j \big(\omega_l \, u_m \big) \\
# &= \partial_j (\omega_i \, u_j ) - \partial_j (\omega_j \, u_i ) \\
# &= u_j \partial_j \omega_i + \omega_i \partial_j u_j - u_i \partial_j \omega_j - \omega_j  \partial_j u_i  \\
# &= (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\omega_i + \omega_i (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}) - u_i (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\omega}) - (\boldsymbol{\omega}\boldsymbol{\cdot}\boldsymbol{\nabla}) u_i. \tag{24}
# \end{align*}
# 
# But, $\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\omega}=0$ (identically) and also for a homogeneous fluid we have $\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0$ (see equation (1)). Thus, (24) simplifies to equation (13):
# 
# $$
# \boldsymbol{\nabla}\times (\boldsymbol{\omega} \times \boldsymbol{u} ) = \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{\omega}  - \boldsymbol{\omega}\boldsymbol{\cdot}\boldsymbol{\nabla}\,\boldsymbol{u}.
# $$
# 
# ***
# 
# **Exercise**: Try generalizing the vorticity conservation for an inviscid fluid with varying density. Prove that for a *barotropic fluid*, that is, for a fluid for which $\boldsymbol{\nabla}\rho\times\boldsymbol{\nabla}p=0$ then:
# 
# $$
# \frac{\mathrm{D}}{\mathrm{D}t}\Big(\frac{\boldsymbol{\omega}}{\rho}\Big) = \Big(\frac{\boldsymbol{\omega}}{\rho}\Big)\boldsymbol{\cdot}\boldsymbol{\nabla}\,\boldsymbol{u}.\tag{25}
# $$
# 
# (Hint: You will have to revisit the mass conservation equation. Now, the mass conservation equation does *not* imply  $\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0$.)