---
jupytext:
  cell_metadata_filter: -all
  formats: md:myst
  text_representation:
    extension: .md
    format_name: myst
    format_version: '0.13'
    jupytext_version: 1.16.6
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

```{code-cell}
:tags: [hide-input]
import mmf_setup; mmf_setup.nbinit(quiet=True)
import os
from pathlib import Path
FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/'
os.makedirs(FIG_DIR, exist_ok=True)
import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL)
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None
```

(sec:Hydrodynamics)=
# Classical Hydrodynamics
This notebook describes the feature of the code to solve simple classical hydrodynamics
problems.  We use these to match to coarse-grained superfluid hydrodynamics simulations,
for example, to extract the effective viscosity.

## Example

As a sample problem, we consider here a harmonically trapped gas oscillating back and
forth in a harmonic trap with a small barrier potential inserted.  The idea is that the
barrier will produce some drag, slowing the oscillations as energy is dissipated.  For
comparison, a superfluid will not lose energy, but will generate phonons (wave
turbulence) resulting in a slow decay of the center of mass oscillations.


```{code-cell}
from gpe.Examples.hydro_1d import Experiment, u

expt = Experiment()
state = expt.get_state()
state.plot()
state.evolve(steps=200, t__final=10*u.ms);
```

### Linearization of Equation

Linearization of the above equation is little difficult in the current form. So let's do some redefinitions: 
$$
 h=-\frac{V}{g} \qquad \eta = n-h \qquad g' = \frac{g}{m}
$$
Then the shallow water equation in 1D takes the form
$$
\dot{\eta} = -\frac{\partial (\eta+h)u}{\partial x} \\
\dot{u} = -\frac{\partial (\frac{u^2}{2}+g'\eta)}{\partial x} + \nu \frac{\partial^2 u}{\partial x^2}
$$

The jacobian matrix for transforming fluxes for $X = (n, nu)$is given by:
\begin{align}
\mathbf{J} = \begin{bmatrix}
                u & h\\
                g' & u\\
                \end{bmatrix}
\end{align}
The eigenvalues are $u \pm \sqrt{g'h}$

Now, we can linearize this equation to get
$$
\dot{\eta} = -\frac{\partial hu}{\partial x} \\
\dot{u} = -\frac{\partial (g\eta)}{\partial x} + \nu \frac{\partial^2 u}{\partial x^2}
$$
I am little unclear on whether to include the double derivative in linear version. For simpler version of stability analysis, we drop the double derivative
$$
\dot{\eta} = -\frac{\partial hu}{\partial x} \\
\dot{u} = -\frac{\partial (g\eta)}{\partial x}
$$
i.e.
$$
\dot{Y} = -\mathbf{C} \frac{\partial Y}{\partial x}
$$
where 
\begin{align}
\mathbf{C} = \begin{bmatrix}
                0 & h\\
                g & 0\\
                \end{bmatrix}
\end{align}

\begin{align}
\mathbf{Y} = \begin{bmatrix}
                \eta\\
                u\\
                \end{bmatrix}
\end{align}


### Von-Neumann Stability Analysis: Convection Equation

Consider the Finite Difference discretization
$$
\frac{Y_{i}^{r+1} - Y_{i}^{r}}{\delta t} = - C \frac{Y_{i+1}^{r} - Y_{i}^{r}}{\delta x}
$$
Putting, $Y_{i}^{r} = A^r (1, 1) e^{ikx}$
$$
\frac{A-1}{\delta t} = -C \frac{e^{ik\delta x} - 1}{\delta x} \\
A = 1 - \frac{\delta t}{\delta x}C\left(e^{ik\delta x} - 1\right)
$$
$$
|A|^2 = (1-\frac{\delta t}{\delta x}C (\cos(k \delta x)-1))^2 + (\frac{\delta t}{\delta x}C \sin(k \delta x))^2
$$

$$
|A|^2 = (1 + \frac{|C| \delta t}{\delta x})^2 - 2 \frac{\delta t}{\delta x}|C| \cos(k \delta x) (1 +\frac{\delta t}{\delta x}|C|) + (\frac{\delta t}{\delta x}C)^2
$$
Now, $|A|^2 \leq 1$, for stability
$$
(1+\frac{\delta t}{\delta x}|C|)(1-\cos(k \delta x)) \leq 0
$$
Hence it is unconditionally unstable


### The Formulation (Finite Volume Method):

\begin{gather}
\dot{\rho} = -\frac{\partial F}{\partial x} + \mathcal{S} \\
\frac{d \rho}{dt} =  - \frac{1}{\bigtriangleup x} (F_{xR} - F_{xL}) + \mathcal{S}_{x} \\
F_{xR} = \frac{1}{2} \left((F_{x+dx} + F_{x}) - |\alpha_{R}|(\rho_{x+dx} - \rho_{x}) \right) \\
F_{xL} = \frac{1}{2} ((F_{x} + F_{x-dx}) - |\alpha_{L}|(\rho_{x} - \rho_{x-dx})) \\
\alpha_I = \frac{\bigtriangleup F}{\bigtriangleup \rho} \quad \bigtriangleup \rho > 1e-6 \\
 = \frac{\partial F}{\partial \rho} \quad otherwise
\end{gather}

\begin{gather}
\frac{\partial j}{\partial t} = -\frac{1}{2*\bigtriangleup x} \left(F_{x+dx} - F_{x-dx} - (\alpha_R j_{x+dx} - (\alpha_R + \alpha_L)j_{x} + \alpha_L j_{x-dx}) \right) + \mathcal{S_x}
\end{gather}

```python

```