---
execution:
  timeout: 60
jupytext:
  cell_metadata_filter: -all
  formats: md:myst,ipynb
  main_language: python
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.7
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

```{code-cell}
:tags: [hide-cell]

import mmf_setup; mmf_setup.nbinit()
import os, sys
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)
%matplotlib
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None

from matplotlib.animation import FuncAnimation
from gpe.contexts import FPS
```

(sec:Viscosity3)=
# Viscosity Notes 3

These are some supporting notes for the document {ref}`sec:EffectiveViscosity`.

## Saint-Venant Equations

Here we consider the following 1D version of the Navier-Stokes equations called the
Saint-Venant equations:
\begin{gather*}
  \dot{n} + (nu)' = 0, \qquad
  D_{t}u = \dot{u} + uu_{,x} = -\frac{\Bigl(V(x) + \mathcal{E}'(n)\Bigr)_{,x}}{m} +
  \frac{\Bigl(\nu(n)u_{,x}\Bigr)_{,x}}{n}.
\end{gather*}

:::::{admonition} A Physical Analogy: Shallow Water
:class: dropdown

These equations follow from the flow of an incompressible fluid of mass density $\rho$
like water flowing along a channel of width $w$ whose bottom height $b(x)$ varies, but
not enough to excite appreciable vertical dynamics.  If the fluid height is $\eta(x)$,
then the internal energy density is given by the gravitational potential energy (here
$a_g$ is the acceleration due to gravity)
\begin{gather*}
  \int_{b(x)}^{b(x) + \eta(x)} a_g \rho w z\;\d{z}
              = a_g w\rho \frac{(b+\eta)^2 - b^2}{2}
              = a_g w\rho \frac{\eta^2}{2} + a_g w b \eta
\end{gather*}
This suggests identifying the following 1D number density and coupling constant to match the GPE.
\begin{gather*}
    n_{1}(x) \equiv \frac{\rho w \eta}{m}, \qquad
    g \equiv \frac{m^2 a_g}{w \rho^2}, \qquad
    V(x) \equiv \frac{m a_g b(x)}{\rho}.
\end{gather*}

To help with intuition, we will sometimes plot our flow $n(x, t)$ and $u(x, t)$ in terms
of this channel flow:
\begin{gather*}
  \eta(x) \equiv \frac{m}{\rho w}n(x), \qquad
  b(x) \equiv \frac{\rho}{ma_g} V(x) =
  \frac{m}{\rho w} \frac{V(x)}{g}
\end{gather*}
For the purposes of plotting, we will choose the width of the channel $w$ and fluid density
$\rho$ such that $m/\rho w = 1$.

:::{admonition} Comparison with [Wikipedia][Shallow water equations]
:class: dropdown

To compare with the [conservative form of the shallow water equations][], note that, in
2D, the following are equivalent, after application of the continuity equation:
\begin{gather*}
   \rho \eta D_t u = \pdiff{(\rho \eta u)}{t} 
   + \pdiff{\rho\eta u^2}{x} + \pdiff{\rho\eta u v}{y}.
\end{gather*}
Their form does not include the potential $V(x)$ (they consider a flat bed).  Explicitly
\begin{gather*}
  \dot{n} = - \sum_{i}\nabla_{i}(nv_{i}), \qquad
  D_t v_i = \dot{v}_{i} + \sum_{j} v_j \nabla_{i}v_i,\\
  \pdiff{n v_i}{t} + \sum_{j}\nabla_j(nv_iv_j) 
  = \underbrace{n \dot{v}_i + n\sum_{j}v_j\nabla_j(v_i)}_{nD_tv_i} 
  + \underbrace{\dot{n}v_i + v_i\sum_{j}\nabla_j(nv_j)}_{0}.
\end{gather*}
:::
:::::

## Numerical Implementation

We start with a simple numerical implementation of these equations on a periodic
lattice.  This should work well when everything is smooth, but will likely fail when the
density drops to zero which will produce a kink.  Another potential failure mode appears
if the local flow velocity exceeds the speed of sound.  Finite viscosity allows for mild
supersonic flow

:::{margin}
To develop this code, I started by working in a Notebook, tweaking things as I went
along.  This is about as far as I want to go before refactoring -- there are getting to
be too many knobs.
:::
```{code-cell}
:tags: [hide-cell]
from scipy.integrate import solve_ivp
from gpe.contexts import FPS

%pylab inline
from scipy.integrate import solve_ivp
from gpe.contexts import FPS

Lx = 10.0
Nx = 128
dx = Lx/Nx
x = np.arange(Nx) * dx - Lx/2
k = 2*np.pi * np.fft.fftfreq(Nx, dx)

# Fourier-space representation of derivatives
D = {1: None, 2: None}
D[1] = 1j*k
assert Nx % 2 == 0
D[1][Nx//2] = 0  # Kill highest unmatched frequency
D[2] = -k**2

def diff(f, d=1):
    """Compute the derivative of f."""
    global D
    return np.fft.ifft(D[d] * np.fft.fft(f)).real

def pack(n, u):
    return np.concatenate([n, u])
    
def unpack(y):
    return np.reshape(y, (2, Nx))

def get_V_m(t, x, d=0):
    """Return the dth derivative of V/m."""
    global V0_m, sigma, dV_m_dx
    V_m = V0_m * np.exp(-(x/sigma)**2/2)
    if d == 0:
        dV = dV_m_dx*x
        return V_m + dV - np.min(dV)  # Constant to make V_min = 0
    elif d == 1:
        return - x * V_m / sigma**2 + dV_m_dx

def get_E_m(n, d=0):
    """Return the dth derivative of E/m where E is the energy density."""
    global g_m
    if d == 0:
        return g_m * n**2 / 2 
    elif d == 1:
        return g_m * n
    elif d == 2:
        return g_m

def get_nu_n(n, d=0):
    """Return the dth derivative of the vicosity diff(nu(n), n, d)/n, divided by n."""
    if d == 0:
        return nu # * n / n
    elif d == 1:
        return nu / n
    else:
        return 0
        
def compute_dy_dt(t, y):
    global get_V_m, get_E_m, get_nu_n
    n, u = unpack(y)
    j = u*n    
    dn_dt = -diff(j)
    n_x = diff(n)
    u_x = diff(u)
    u_xx = diff(u, d=2)
    du_dt = -u*u_x
    du_dt -= get_V_m(t, x, d=1) + get_E_m(n, d=2)*n_x
    du_dt += get_nu_n(n, d=1)*n_x*u_x + get_nu_n(n) * u_xx / n
    
    # Correction needed to keep finite j = j_avg
    du = j_avg/n - u
    du_dt += j_corr_factor * du
    return pack(dn_dt, du_dt)
    
# Quick test of derivatives
k = 2*np.pi / Lx
s, c = np.sin(k*x), np.cos(k*x)
f = np.exp(s)
df = k*c*f
ddf = k**2*f*(c**2 - s)
assert np.allclose(diff(f), df)
assert np.allclose(diff(f, 2), ddf)
```
```{code-cell}
:tags: [hide-input]
# Here is a simulation of flow over a bump.
V0_m = 1.0
dV_m_dx = -0.05
sigma = 0.5
g_m = 1.0
n0 = 1.5
u0 = 0.3
j_avg = n0*u0
nu = 0.6
t = 0

b = get_V_m(t, x) / g_m
h = n0
eta = h - b
n = eta
u = j_avg / n
y0 = pack(n, u)

dt = dx / nu * 0.002
j_corr_factor = 0/dt/100

fig, axs = plt.subplots(2, 1, sharex=True)
ax = axs[0]
ax_stream = axs[1]
ax_j = ax_stream.twinx()
ax_c = ax.twinx()
y = y0.copy()
for n in FPS(120, fig=fig, embed=True, fps=20):
    for _ in range(500):
        dy_dt = compute_dy_dt(t, y)
        y += dy_dt * dt
        t += dt
    n, u = unpack(y)
    ax.cla()
    ax.plot(x, n)
    ax.set(ylim=(0, 2.0), title=f"{t=:.4g}, N={sum(n)*dx:.4g}",
           ylabel="$n$")
           
    # Speed of sound
    u_c = np.sqrt(n*get_E_m(n, d=2))
    ax_c.cla()
    ax_c.yaxis.set_label_position("right")
    ax_c.plot(x, u/u_c, '--C1')
    ax_c.axhline(1.0, color='k', ls=':')
    ax_c.set(ylim=(-0.1, 1.3), ylabel="Mach number $u/u_c$")
    
    ax_stream.cla()
    b = get_V_m(t, x) / g_m
    eta = n
    h = b + eta
    ax_stream.plot(x, b, '-k')
    ax_stream.fill_between(x, b, h, color='cyan')
    ax_stream.set(ylim=(0, 2.0), xlim=(x.min(), x.max()),
                  xlabel="$x$", ylabel="Shallow-water analog.")

    ax_j.cla()
    ax_j.plot(x, n*u)
    ax_j.plot(x, u)
    ax_j.yaxis.set_label_position("right")
    ax_j.set(ylim=(-0.2, 1.0), 
             ylabel="$j$ (blue) and $u$ (orange)")
```
:::{margin}
Originally I tried to enforce a finite (constant) current $j$ by adding the following
term:
\begin{gather*}
  \dot{u}(x) = \cdots + \lambda \left(\frac{j}{n(x)} - u(x)\right).
\end{gather*}
This is still in the code where $\lambda$ is called `j_corr_factor` but is not used in
either of these movies.
:::
In this example, we add a gradient to our potential, noting that only $V'(x)$ appears in
our equations of motion.  Thus, although the potential is not periodic, the force is.
This constant is needed to ensure that a finite flow remains.


[Rankine-Hugoniot conditions]: <https://en.wikipedia.org/wiki/Rankine%E2%80%93Hugoniot_conditions>
[Shallow water equations]: <https://en.wikipedia.org/wiki/Shallow_water_equations>
[conservative form of the shallow water equations]: <https://en.wikipedia.org/wiki/Shallow_water_equations#Conservative_form>
[weir]: https://en.wikipedia.org/wiki/Weir
[Shallow water equations]: <https://en.wikipedia.org/wiki/Shallow_water_equations>
[conservative form of the shallow water equations]: <https://en.wikipedia.org/wiki/Shallow_water_equations#Conservative_form>