import mmf_setup; mmf_setup.nbinit()
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)
%matplotlib
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None
This cell adds /home/docs/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:
- Choose "Trust Notebook" from the "File" menu.
- Re-execute this cell.
- Reload the notebook.
Using matplotlib backend: module://matplotlib_inline.backend_inline
%load_ext autoreload
%autoreload 2
from importlib import reload
from gpe.contexts import FPS
import time
import viscosity
f = viscosity.Flow(n_min=0)
n0 = np.where(f.basis.x<0, 1.0, 0.1)
n1 = np.where(f.basis.x<0, 1.0, 0.01)
y0 = f.pack(n0, 0*n0)
y1 = f.pack(n1, 0*n1)
dt = 0.01
t = 0
nfevs0 = []
nfevs1 = []
fig, axs = fig_axs = f.plot(y0)
for n in FPS(50, fig=fig, display=True):
y0 = f.evolve(y0, t=t, dt=dt)
nfevs0.append(f._res.nfev)
y1 = f.evolve(y1, t=t, dt=dt)
nfevs1.append(f._res.nfev)
t += dt
fig_axs = f.plot(y=y0, t=t, fig_axs=fig_axs, c='C0', label=f"$n_R=0.1$: nfev={nfevs0[-1]}")
fig, axs = f.plot(y=y1, t=t, fig_axs=fig_axs, c='C1', label=f"$n_R=0.01$: nfev={nfevs1[-1]}", cla=False)
axs[0].legend(loc='center left')
import viscosity; reload(viscosity)
f = viscosity.Flow(n_min=0)
n = np.where(f.basis.x<0, 1.0, 0.01)
y = f.pack(n, 0*n)
dt = 0.01
t = 0
nfevs = []
for n in FPS(50, fig=fig, display=True):
#y, t = f.step_euler(y, t=t, dt=dt), t+dt
y, t = f.evolve(y, t=t, dt=dt), t+dt
nfevs.append(f._res.nfev)
fig_axs = f.plot(y=y, t=t, fig_axs=fig_axs)
Viscosity 4: Lagrangian#
Nx = 128
x0 = np.linspace(-1.0, 1.0, Nx)
dx0 = np.diff(x0).mean()
g_m = 1.0 # Coupling constant g/m
nu = 0.2 # Viscosity coefficient
w = np.sqrt(2) # Trap frequency
V0_m = 0.5 # Bump height
sigma = 0.1 # Bump width
# Initial density provile
X0_shift = 0.5 # Initial shift of the cloud
X0 = x0 + X0_shift
n0 = w**2/2/g_m * (1 - x0**2)
n0_x0 = -w**2 / g_m * x0
# Finite difference approximations for the derivatives
o = np.ones(Nx)
Df = (np.diag(o[1:], 1) - np.diag(o, 0))/dx0
Db = (np.diag(o, 0) - np.diag(o[:-1], -1))/dx0
Dc = (Df + Db)/2
D2 = Db @ Df
# Neumann bc's
Dc[0, :2] = [-1/dx0, 1/dx0]
Dc[-1, -2:] = [-1/dx0, 1/dx0]
D2[0, :] = 0
D2[-1, :] = 0
def pack(X, X_t):
return np.ravel([X, X_t])
def unpack(y):
return np.reshape(y, (2, Nx))
def unpack_xnu(y):
X, X_t = unpack(y)
X_x0 = Dc @ X
n = n0 / X_x0
u = X_t
x = X
return x, n, u
def compute_dy_dt(t, y):
X, X_t = unpack(y)
X_x0 = Dc @ X
X_x0x0 = D2 @ X
n_x = n0_x0 / X_x0**2 - n0 * X_x0x0 / X_x0**3
u_xx = (D2 @ X_t) / X_x0**2 - Dc @ X_t * X_x0x0 / X_x0**3
a_ext = - w**2 * X
a_ext += V0_m * np.exp(-X**2/2/sigma**2) * X/sigma**2
X_tt = a_ext - g_m * n_x + nu * u_xx
return pack(X_t, X_tt)
y0 = pack(X0, 0*X0)
X, X_t = unpack(compute_dy_dt(0, y0))
plt.plot(x0, X_t)
%%time
from scipy.integrate import solve_ivp
Tw = 2*np.pi / w # Trapping period
T = 40*Tw
Nt = 500
t_eval = np.linspace(0, T, Nt)
res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, T),
t_eval=t_eval, method="BDF")
if not res.success:
print(res.message)
ts = res.t
Nt = len(ts)
X, X_t = res.y.reshape((2, Nx, Nt))
fig, ax = plt.subplots()
ax.plot(ts / Tw, X[::10].T)
ax.set(xlabel="$t/T_w$", ylabel="X");
from scipy.integrate import solve_ivp
T = 4.0
Nt = 500
t_eval = np.linspace(0, T, Nt)
res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, T),
t_eval=t_eval, method="BDF")
if not res.success:
print(res.message)
ts = res.t
Nt = len(ts)
X, X_t = res.y.reshape((2, Nx, Nt))
plt.plot(ts, X[::10].T);
X_x0 = Dc @ X
n = n0[:, None] / X_x0
x_ = np.linspace(X.min(), X.max(), 200)
V_ext = w**2 * X**2 / 2
V_ext += V0_m * np.exp(-X**2/2/sigma**2)
b = V_ext / g_m
b_ = w**2 * x_**2 / 2
b_ += V0_m * np.exp(-x_**2/2/sigma**2)
b_ /= g_m
h = b + n
fig, ax = plt.subplots()
for i in FPS(np.arange(Nt)[::1], fig=fig, display=True):
t = ts[i]
ax.cla()
#ax.plot(X[:, i], n[:, i])
ax.plot(x_, b_, 'k-')
ax.fill_between(X[:, i], b[:, i], h[:, i], color='cyan')
ax.set(title=f"{t=:.4f}",
xlim=(X.min(), X.max()),
ylim=(0, h.max()))
# Spectral version
def s():
import scipy.fft
sp = scipy
Nx = 128
th0 = (2*np.arange(Nx) + 1) * np.pi / 2 / Nx
k = np.arange(Nx)
x0 = np.cos(th0)
X0_shift = 0.5 # Initial shift of the cloud
X0 = x0 + X0_shift
g_m = 1.0 # Coupling constant g/m
nu = 0.2 # Viscosity coefficient
w = np.sqrt(2) # Trap frequency
V0_m = 0.5 # Bump height
sigma = 0.1 # Bump width
# Initial density provile
n0 = w**2/2/g_m * (1 - x0**2)
n0_x0 = -w**2 / g_m * x0
def diff(f, d=1, _k=k, _th=th0):
ft = sp.fft.dct(f, type=2, norm='forward')
# Shift the coefficients and pad with zero.
df_dth = sp.fft.dst(np.concatenate([-(ft*_k)[1:], [0]]), type=3)
s = np.sin(_th)
df_dx = -df_dth / s
if d == 1:
return df_dx
elif d == 2:
d2f_dth2 = sp.fft.idct(-ft*_k**2, type=2, norm='forward')
d2f_dx2 = (d2f_dth2 - df_dth / np.tan(_th)) / s**2
return d2f_dx2
return diff(df_dx, d=d-1)
def pack(X, X_t):
return np.ravel([X, X_t])
def unpack(y):
return np.reshape(y, (2, Nx))
def unpack_xnu(y):
X, X_t = unpack(y)
X_x0 = diff(X)
n = n0 / X_x0
u = X_t
x = X
return x, n, u
def compute_dy_dt(t, y):
X, X_t = unpack(y)
X_x0 = diff(X)
X_x0x0 = diff(X, d=2)
n_x = n0_x0 / X_x0**2 - n0 * X_x0x0 / X_x0**3
u_xx = (diff(X_t, d=2)) / X_x0**2 - diff(X_t) * X_x0x0 / X_x0**3
a_ext = - w**2 * X
a_ext += V0_m * np.exp(-X**2/2/sigma**2) * X/sigma**2
X_tt = a_ext - g_m * n_x + nu * u_xx
return pack(X_t, X_tt)
y0 = pack(X0, 0*X0)
return x0, y0, compute_dy_dt, n0_x0
# Spectral version
import scipy.fft
sp = scipy
Nx = 128
th0 = (2*np.arange(Nx) + 1) * np.pi / 2 / Nx
k = np.arange(Nx)
x0 = np.cos(th0)
X0_shift = 0.5 # Initial shift of the cloud
X0 = x0 + X0_shift
g_m = 1.0 # Coupling constant g/m
nu = 0.2 # Viscosity coefficient
w = np.sqrt(2) # Trap frequency
V0_m = 0.5 # Bump height
sigma = 0.1 # Bump width
# Initial density provile
n0 = w**2/2/g_m * (1 - x0**2)
n0_x0 = -w**2 / g_m * x0
def diff(f, d=1, _k=k, _th=th0):
ft = sp.fft.dct(f, type=2, norm='forward')
# Shift the coefficients and pad with zero.
df_dth = sp.fft.dst(np.concatenate([-(ft*_k)[1:], [0]]), type=3)
s = np.sin(_th)
df_dx = -df_dth / s
if d == 1:
return df_dx
elif d == 2:
d2f_dth2 = sp.fft.idct(-ft*_k**2, type=2, norm='forward')
d2f_dx2 = (d2f_dth2 - df_dth / np.tan(_th)) / s**2
return d2f_dx2
return diff(df_dx, d=d-1)
def pack(X, X_t):
return np.ravel([X, X_t])
def unpack(y):
return np.reshape(y, (2, Nx))
def unpack_xnu(y):
X, X_t = unpack(y)
X_x0 = diff(X)
n = n0 / X_x0
u = X_t
x = X
return x, n, u
def compute_dy_dt(t, y):
X, X_t = unpack(y)
X_x0 = diff(X)
X_x0x0 = diff(X, d=2)
n_x = n0_x0 / X_x0**2 - n0 * X_x0x0 / X_x0**3
u_xx = (diff(X_t, d=2)) / X_x0**2 - diff(X_t) * X_x0x0 / X_x0**3
a_ext = - w**2 * X
a_ext += V0_m * np.exp(-X**2/2/sigma**2) * X/sigma**2
X_tt = a_ext - g_m * n_x + nu * u_xx
return pack(X_t, X_tt)
y0 = pack(X0, 0*X0)
x0_, y0_, compute_dy_dt_, n0_x0_ = s()
dX, dX_t = unpack(compute_dy_dt(0, y0))
X, X_t = unpack(y0)
plt.plot(x0, X)
plt.plot(x0, dX_t)
dX_, dX_t_ = unpack(compute_dy_dt_(0, y0_))
X_, X_t_ = unpack(y0_)
plt.plot(x0_, X_)
plt.plot(x0_, dX_t_)
res_ = solve_ivp(compute_dy_dt_, y0=y0_, t_span=(0, 4), method="BDF")
res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, 4), method="BDF")
plt.plot(x0, unpack(res.y[:, -1])[0])
plt.plot(x0_, unpack(res_.y[:, -1])[0])
%%time
from scipy.integrate import solve_ivp
Tw = 2*np.pi / w # Trapping period
T = Tw
Nt = 500
t_eval = np.linspace(0, T, Nt)
ress = []
for V0_m in [0.0, 0.5]:
res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, T),
t_eval=t_eval, method="BDF")
if not res.success:
print(res.message)
ress.append((V0_m, res))
ts = res.t
Nt = len(ts)
fig, axs = plt.subplots(2, 1, figsize=(6, 3), sharex=True, gridspec_kw=dict(hspace=0))
for ax, (V0_m, res) in zip(axs, ress):
X, X_t = res.y.reshape((2, Nx, Nt))
ax.plot(ts / Tw, X[::10].T)
ax.set(ylabel=f"$X$ ($V_0/m = {V0_m:.1f})$")
axs[1].set(xlabel="$t/T_w$");
y = f.get_ground_state(x0=0.0)
fig_axs = fig, axs = f.plot(y)
dt = 0.01
t = 0
for n in FPS(200, fig=fig, display=True):
#y, t = f.step_euler(y, t=t, dt=dt), t+dt
y, t = f.evolve(y, t=t, dt=dt), t+dt
fig_axs = f.plot(y=y, t=t, fig_axs=fig_axs)
import viscosity; reload(viscosity)
f = viscosity.FlowTrap(n_min=0.1)
y = f.get_ground_state(x0=0.0)
fig_axs = fig, axs = f.plot(y)
dt = 0.01
t = 0
for n in FPS(200, fig=fig, display=True):
#y, t = f.step_euler(y, t=t, dt=dt), t+dt
y, t = f.evolve(y, t=t, dt=dt), t+dt
fig_axs = f.plot(y=y, t=t, fig_axs=fig_axs)
%load_ext autoreload
%autoreload 2
from gpe.Examples.hydro_1d import Experiment, u
res0 = []
for n in [6, 7, 8, 9, 10]:
e = Experiment(V0p=0.5, Nx=2**n)
s = e.get_state()
#s.plot()
hist = s.evolve(steps=2000, t_final=10*u.ms, hist=True, show_plot=True);
res0.append((np.array([s.t for s in hist]),
np.array([(s.get_density()*s.x).mean() / s.get_density().mean() for s in hist])))
plt.subplots(figsize=(20,3))
[plt.plot(t/s.experiment.t_unit, x, label=str(n)) for n, (t, x) in enumerate(res0)];plt.legend()
plt.plot(ts0/s.experiment.t_unit, x_cms0)
plt.plot(ts1/s.experiment.t_unit, x_cms1)
from scipy.integrate import solve_ivp
import scipy as sp
class Flow:
m = 1.2
g = 1.3
nu = 0.1 # Viscosity
V0 = 0.2 # Potential barrier height
L = 0.2 # Length of barrier
n0 = 1.3 # Steady-state density on the left.
v0 = 0.14 # Steady-state velocity on the left.
dV = 0.15 # Potential drop (magnitude)
potential = 'gaussian'
potential_k = 1
# Boundary conditions on the left
def __init__(self, **kw):
for key in kw:
if not hasattr(self, key):
raise ValueError(f"Unknown {key=}")
setattr(self, key, kw[key])
self.init()
def init(self):
if self.potential == 'gaussian':
f = 10.0
s = self.L/2
self._V = lambda x: self.V0*np.exp(-f*(x/s-1)**2)
self._dV = lambda x: -2*f*(x/s-1)/s*self.V0*np.exp(-f*(x/s-1)**2)
else:
xp, fp = [0, self.L/2, self.L], [self.dV, self.V0, 0]
if self.potential_k > 1:
xp = self.L*np.linspace(0, 1, 7)
fp = [self.dV]*3 + [self.V0] + [0]*3
xp = self.L*np.array([0, 0.01, 0.5, 0.99, 1])
fp = [self.dV]*2 + [self.V0] + [0]*2
self._V = sp.interpolate.InterpolatedUnivariateSpline(
xp, fp, k=self.potential_k, ext='const')
self._dV = self._V.derivative()
# Compute the nu=0 bounds assuming the GPE
j = self.v0*self.n0
m = self.m
nc = self.nc = (m * j**2/self.g)**(1/3) # Relies on GPE EoS.
self.mu0 = self.get_Eint(n=self.n0, d=1)
self.E0 = m*self.v0**2/2 + self.mu0
self.Vmax = self.E0 - m*(j/nc)**2/2 - self.get_Eint(nc, d=1)
def get_Eint(self, n, d=0):
if d == 0:
return self.g * n**2 /2
elif d == 1:
return self.g * n
elif d == 2:
return self.g
else:
return 0
def get_P(self, n):
E, mu = [self.get_Eint(n, d=d) for d in [0, 1]]
P = n*mu - E
return P
def compute_dy_dx(self, x, y):
n, u, du, F = y
dn = -n*du/u
ddu = n * (u*du + (self._dV(x) + self.get_Eint(n, d=2)*dn)/self.m) / self.nu
dF = n * self._dV(x)
return (dn, du, ddu, dF)
def solve(self, dv0=0.0, x0_L=-0.1, x_L=2.0,
method='BDF', max_step=0.1, x_eval=None, **kw):
x_span = (x0_L*self.L, x_L*self.L)
y0 = (self.n0, self.v0, dv0, 0.0)
res = solve_ivp(self.compute_dy_dx, y0=y0, t_span=x_span, t_eval=x_eval,
method=method, max_step=max_step, **kw)
if not res.success:
raise Exception(res.message)
res.x = res.t
res.n, res.u, res.du, res.F = res.y
res.dn, du, res.ddu, res.dF = self.compute_dy_dx(res.x, res.y)
self.res = res
return res
f = Flow()
res = f.solve()
class Flow1(Flow):
def plot(self, x0_L=0.0, x_L=1.2, ax=None):
"""Plot the density and velocity of flow over a barrier."""
res = self.solve(x0_L=x0_L, x_L=x_L)
x, n = res.x, res.n
Vx = self._V(x)
n_TF = (self.mu0 - Vx)/self.g
if ax is None:
fig, ax = plt.subplots(1, 1, figsize=(3, 1.5))
ax.plot(x, n, label="$n$")
ax.plot(x, n_TF, '--', label='$n_{TF}$')
ax.axhline(self.nc, color='y', ls=':')
ax.set(ylim=(0, n.max()*1.05),
xlabel="$x$", ylabel="$n$")
ax1 = ax.twinx()
ax1.plot(x, Vx/self.V0*self.nc, 'C2-.', label='$V(x)$')
ax.legend(loc='upper right')
plt.sca(ax)
res.ax = ax
res.ax1 = ax1
return res
Flow1(nu=0.0001, n0=2.0, V0=1.0, v0=0.3).plot()
Flow1(nu=0.0001, n0=2.0, V0=1.0, v0=0.425).plot()
Flow1(nu=0.003, n0=2.0, V0=1.0, v0=0.455).plot();
f = Flow1(nu=0.0001, n0=2.0, V0=1.0, v0=0.3)
nc = 1.0
jc = nc*np.sqrt(nc*f.get_Eint(nc, d=2)/f.m)
nus = np.linspace(0.001, 0.2)
fs = []
Fs = []
for nu in nus:
f = Flow1(n0=2.0, v0=1.0, V0=0.19, nu=nu)
res = f.solve(atol=1e-12, rtol=1e-12)
fs.append(f)
Fs.append(res.F[-1]-res.F[0])
plt.plot(nus, Fs)
nu = 0.05
V0s = np.linspace(0, 0.19, )
fs = []
Fs = []
for V0 in V0s:
f = Flow1(n0=2.0, v0=1.0, V0=V0, nu=nu)
res = f.solve(x_L=5.0, atol=1e-12, rtol=1e-12)
fs.append(f)
Fs.append(res.F[-1]-res.F[0])
fig, ax = plt.subplots(1, 1)
ax.plot(V0s, Fs)
ax.set(xlabel=r"$V_0$", ylabel="$F$")
ind = 25
fig, axs = plt.subplots(1, 3, figsize=(12, 3))
for f, ax in zip([fs[0], fs[ind], fs[-1]], axs):
f.plot(ax=ax)
def plot_flow(x0_L=0.0, x_L=2.0, **kw):
f = Flow(**kw)
res = f.solve(x0_L=x0_L, x_L=x_L, max_step=0.001, rtol=1e-12, atol=1e-12)
m, nu = f.m, f.nu
x, u, du, ddu, n, dn = res.x, res.u, res.du, res.ddu, res.n, res.dn
e, de, dde = [f.get_Eint(n, d=d) for d in [0, 1, 2]]
V, dV = f._V(x), f._dV(x)
j = f.v0*f.n0
n0 = f.n0
K0 = m*f.v0**2/2
mu0 = f.get_Eint(f.n0, d=1)
E0 = K0 + mu0
nc = (f.m * j**2 / f.g)**(1/3)
Vmax = E0 - 3/2*(f.m*f.g**2*j**2)**(1/3)
print(f.nc, f.Vmax, f.v0, Vmax)
nc0 = 2*(E0 - f.V0)/3/f.g
n_TF = (mu0 - V)/f.g
n_min = (mu0 - f.V0)/f.g
ns = [n_TF]
for _i in range(10):
_n = ns[-1]
dn = -m*j**2/2/f.g * (1/_n**2 - 1/n0**2)
ns.append(n_TF+dn)
α_min = (mu0 - f.V0)/(m*j**2/2*(1/n_min**2-1/n0**2))
print(n_min, n.min(), nc, nc0, n_min**2*(E0-f.V0)/(m*j**2/2))
#(de/(m*u**2/2)).min()
print(f"{α_min=:.2g}")
print(E0, f.V0)
fig, ax = plt.subplots(1, 1, figsize=(3, 2))
ax.plot(x, n)
for _n in ns:
ax.plot(x, _n, ':')
ax.set(ylabel="n")
ax1 = ax.twinx()
ax1.plot(x, f._V(x), 'C8--')
res.E0, res.j = E0, j
return f, res
f, res = plot_flow(dV=0.0, V0=1.0, nu=0.0001, n0=10.0, v0=1.795, #v0=1.795,
g=1.0, x_L=0.6)
Pn = [f.g, f.V0-res.E0, 0, f.m/2*res.j**2]
n = np.linspace(0, f.n0)
plt.plot(n, np.polyval(Pn, n))
Pu = [f.m/2, 0, f.V0-res.E0, f.g*res.j]
u = np.linspace(f.v0, 2*f.v0)
plt.plot(u, np.polyval(Pu, u))
Flow Explorer#
from ipywidgets import interact
kw = dict(L=0.1, x0_L=-100.0, x_L=100.2)
#plot_flow(dV=0.0, V0=8.46, nu=0.1, n0=10.0, v0=0.1, g=1.0);
plot_flow(dV=-0.1, V0=0.001, nu=0.1, n0=1.0, v0=1.0, g=1.0, potential='step', **kw);
plot_flow(dV=0.0, V0=1.0, nu=0.1, n0=10.0, v0=0.1, g=1.0)
plt.title("Slow flow over a barrier");
plot_flow(dV=0.0, V0=1.0, nu=0.1, n0=10.0, v0=1.9, g=1.0)
plt.title("Fast flow over a barrier");
plot_flow(dV=0.0, V0=1.0, nu=0.01, n0=10.0, v0=10.9, g=1.0, x_L=0.13)
plt.title("Fast flow over a barrier");
fig, ax = plt.subplots()
f = Flow(dV=0.05, V0=0.1, nu=0.1, n0=10.0, v0=2.743, g=1.0)
res = f.solve(x_L=2.0, max_step=0.001, rtol=1e-12, atol=1e-12)
m, nu = f.m, f.nu
x, u, du, ddu, n, dn = res.x, res.u, res.du, res.ddu, res.n, res.dn
E, dE, ddE = [f.get_Eint(n, d=d) for d in [0, 1, 2]]
V, dV = f._V(x), f._dV(x)
#ax.plot(x, n)
if False:
ax.plot(x, u*du)
ax.plot(x, -(dV + ddE*dn)/m)
ax.plot(x, nu*ddu/n)
else:
ax.plot(x, n)
#ax.set(ylim=(-0.1, 0.1))
ax1 = ax.twinx()
ax1.plot(x, f._V(x), 'C8--')
fig, ax = plt.subplots()
f = Flow(dV=0, V0=0.1, nu=0.1, n0=1.0, v0=0.15, g=0.0)
res = f.solve(x_L=2.0, max_step=0.001, rtol=1e-12, atol=1e-12)
m, nu = f.m, f.nu
x, u, du, ddu, n, dn = res.x, res.u, res.du, res.ddu, res.n, res.dn
E, dE, ddE = [f.get_Eint(n, d=d) for d in [0, 1, 2]]
V, dV = f._V(x), f._dV(x)
#ax.plot(x, n)
if False:
ax.plot(x, u*du)
ax.plot(x, -(dV + ddE*dn)/m)
ax.plot(x, nu*ddu/n)
else:
ax.plot(x, n)
#ax.set(ylim=(-0.1, 0.1))
ax1 = ax.twinx()
ax1.plot(x, f._V(x), 'C8--')
ec = (res.n*res.u*(
f.m*res.u**2/2 + f._V(res.x) + f.get_Eint(res.n, d=1))
-m*f.nu*res.u*res.du
)
plt.plot(x, ec)
from scipy.integrate import solve_ivp
from scipy.interpolate import InterpolatedUnivariateSpline as Spline
v0 = f.v0
L = f.L
T = L/v0
print(v0**2/2, f.V0)
def V_m(x, f=10.0, L=L, V0=f.V0, d=0):
s = L/2
if d == 0:
return V0 * np.exp(-f*(x/s-1)**2)
elif d == 1:
return -2*f*(x/s-1)/s*V0*np.exp(-f*(x/s-1)**2)
def compute_dy_dt(t, y):
x, v = y
dx = v
dv = -V_m(x, d=1)
return (dx, dv)
y0 = [0, v0]
res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, T),
max_step=0.001,
t_eval=np.linspace(0, T, 300))
t = res.t
x, v = res.y
n = f.n0/np.diff(x)*(x[1]-x[0])
x, v = [(v[1:]+v[:-1])/2 for v in [x, v]]
plt.plot(x, n)
from scipy.integrate import solve_ivp
from scipy.interpolate import InterpolatedUnivariateSpline as Spline
def V_m(x, f=10.0, L=1.0, V0=0.2, d=0):
s = L/2
if d == 0:
return V0 * np.exp(-f*(x/s-1)**2)
elif d == 1:
return -2*f*(x/s-1)/s*V0*np.exp(-f*(x/s-1)**2)
def compute_dy_dt(t, y):
x, v = y
dx = v
dv = -V_m(x, d=1)
return (dx, dv)
y0 = [0, 1]
res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, 1),
max_step=0.001,
t_eval=np.linspace(0, 1, 300))
t = res.t
x, v = res.y
n = f.n0/np.diff(x)*(x[1]-x[0])
x, v = [(v[1:]+v[:-1])/2 for v in [x, v]]
plt.plot(x, n)
f.V0
from scipy.integrate import solve_ivp
from scipy.interpolate import InterpolatedUnivariateSpline as Spline
def V_m(x, f=10.0, L=1.0, V0=0.2, d=0):
s = L/2
if d == 0:
return V0 * np.exp(-f*(x/s-1)**2)
elif d == 1:
return -2*f*(x/s-1)/s*V0*np.exp(-f*(x/s-1)**2)
def compute_dy_dt(t, y):
x, v = y
dx = v
dv = -V_m(x, d=1)
return (dx, dv)
y0 = [0, 1]
res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, 1),
max_step=0.001,
t_eval=np.linspace(0, 1, 300))
t = res.t
x, v = res.y
n = f.n0/np.diff(x)*(x[1]-x[0])
x, v = [(v[1:]+v[:-1])/2 for v in [x, v]]
plt.plot(x, n)
m = 1
v0 = 1
t = np.linspace(0, 2)
x = v0*t - 0.1*np.exp(-20*(t-1)**2)
v = v0 + 0.1*np.exp(-20*(t-1)**2)*40*(t-1)
V = m*(v0**2 - v**2)/2
a = 0.1*np.exp(-20*(t-1)**2)*40*(1 - 40*(t-1)**2)
plt.plot(x, a)
fig, ax = plt.subplots()
axn = ax.twinx()
js = np.linspace(0.001, 10)
nus = [0.025, 0.05, 0.1, 0.2]
nL = 1.0
g = 1.0
F_unit = g*nL**2/2
for nu in nus:
Fs = []
nRs = []
for j in js:
f = Flow(n0=nL, g=g, v0=j/nL, dV=0.0, nu=nu)
res = f.solve(method='LSODA', max_step=0.1, rtol=1e-8)
Fs.append(res.F[-1])
nRs.append(res.n[-1])
Fs, nRs = map(np.asarray, [Fs, nRs])
l, = ax.plot(js/nu, Fs/F_unit, label=f"{nu=}")
axn.plot(js/nu, nRs/nL, ls='--', c=l.get_c())
ax.set(xlabel=r'$j/\nu$', ylabel='$F$ [$gn_L^2/2$] (solid)',
xlim=(0, 100))
axn.set(ylabel='$n_R/n_L$ (dashed)')
ax.legend(loc='center left')
ax.grid('on');
fig, ax = plt.subplots()
axn = ax.twinx()
js = np.linspace(0.001, 3)
nus = [0.025, 0.05, 0.1, 0.2]
nL = 1.0
g = 1.0
F_unit = g*nL**2/2
for nu in nus:
Fs = []
nRs = []
for j in js:
f = Flow(n0=nL, g=g, v0=j/nL, dV=0.0, nu=nu)
res = f.solve(method='LSODA', max_step=0.1, rtol=1e-8)
Fs.append(res.F[-1])
nRs.append(res.n[-1])
Fs, nRs = map(np.asarray, [Fs, nRs])
l, = ax.plot(js, Fs/F_unit, label=f"{nu=}")
axn.plot(js, nRs/nL, ls='--', c=l.get_c())
ax.set(xlabel='$j$', ylabel='$F$ [$gn_L^2/2$] (solid)')
axn.set(ylabel='$n_R/n_L$ (dashed)')
ax.legend(loc='center left')
ax.grid('on');
In the previous plot, we look at static flow across a barrier as a function of the flow-rate (current) \(j\). We plot the total force on the barrier (solid), and the density drop (dashed).
We start with the two conservation laws for density and momentum far to the left and right of the potential, assuming that \(u_{,x}=0\) in these regions, and setting \(F(x_L) = 0\) and \(F(x_R) = F\):
\[\begin{gather*}
n_Lu_L = n_Ru_R = j = \text{const.}\\
mn_Lu_L^2 + n_L\mathcal{E}'(n_L) - \mathcal{E}(n_L) =
mn_Ru_R^2 + F + n_R\mathcal{E}'(n_R) - \mathcal{E}(n_R)\\
\frac{m j^2}{n_L} + n_L\mathcal{E}'(n_L) - \mathcal{E}(n_L) =
\frac{m j^2}{n_R} + F + n_R\mathcal{E}'(n_R) - \mathcal{E}(n_R)\\
\end{gather*}\]
\[\begin{gather*}
nu^2 + \frac{P}{m} - ν u_{,x}
= nu^2 + \frac{F(x) + n\mathcal{E}' - \mathcal{E}}{m} - ν
u_{,x}
= \text{const.},\\
P = F(x) + n\mathcal{E}'(n) - \mathcal{E}(n), \qquad
F(x) = \int_{x_0}^{x} n(x)V'(x)\d{x}.
\end{gather*}\]
Assuming that \(u_{,x} = 0\) far from the potential (static flow), and setiing \(V(x_R) = 0\), we have the following after integrating the energy current:
\[\begin{gather*}
\left.
u\left(
\frac{m}{2}nu^2 + nV + n\mathcal{E}'(n) - m\nu u_{,x}
\right)
\right|^{x_R}_{x_L}
= - \int_{x_L}^{x_R} m\nu u_{,x}^2 \d{x} \leq 0,\\
\tfrac{m}{2}(u_{R}^2 - u_{L}^2) + V_L + \mathcal{E}'(n_R) -
\mathcal{E}'(n_L) =
-\frac{1}{j}\int_{x_L}^{x_R} m\nu u_{,x}^2 \d{x}
\end{gather*}\]
Correction Needed Below#
The following analysis is incorrect. It is missing the factor of \(1/n\) discussed above, and neglects to include a difference in the potential from one side to the other (DC offset) required to maintain a persistent flow in the presence of dissipation.
Steady State (Incorrect)#
Now consider steady-state flow \(n(x, t) = n(x)\), \(u(x, t) = u(x)\). These become:
\[\begin{gather*}
(nu)_{,x} = 0, \qquad
muu_{,x} + V_{,x} + \mu'(n)n_{,x} = mνu_{,xx}.
\end{gather*}\]
The first equation establishes that the flux \(Φ = nu\) is the same at every point, allowing us to exchange velocity \(u\) for density \(n\):
\[\begin{gather*}
n = \frac{Φ}{u}, \qquad
n_{,x} = -\frac{Φ}{u^2}u_{,x}.
\end{gather*}\]
Let’s assume that \(V(x)\) has compact support over \([0, L]\), then we can give boundary conditions \(u(0) = u_0\), \(u_{,x}(0) = 0\). The second equation and be integrated once to give the non-linear first-order equation
\[\begin{gather*}
mνu_{,x} = \frac{mu^2}{2} + V(x) + \overbrace{\mu\left(\frac{Φ}{u}\right)}^{\mu(n)}.
\end{gather*}\]
For general \(V(x)\), this has no analytic solution, but we can compute a perturbative solution if we let \(V(x) \rightarrow \lambda V(x)\) where \(\lambda \ll 1\). We note that \(u(x) = u_0\) is a solution for \(\lambda = 0\) and write $\(u(x) = u_0 + \lambda u_1(x) + \lambda^2 u_2(x) + O(\lambda^3).\)\(_Note: In what follows, we revert back to the notation_ \)u_1’ = [u_1]{,x}$…
I originally wrote it this way, and have not changed it yet._
As shown below, we originally tried expanding to first order in \(\lambda\), but found this was insufficient, so we expand to second order in \(\lambda\) obtaining the following linear equations for \(u_1\) and \(u_2\):
\[\begin{align*}
mνu_1' + \overbrace{\left(\frac{n_0}{u_0}\mu'(n_0) - mu_0\right)}^{-mνc}u_1
&= V(x),\\
mνu_2' + \left(\frac{n_0}{u_0}\mu'(n_0) - mu_0\right)u_2
&= \underbrace{\left(\frac{m}{2}
+ \frac{2n_0\mu'(n_0)+n_0^2\mu''(n_0)}{2u_0^2}
\right)}_{mνb}u_1^2,\\
u_1' - c u_1
&= V(x)/mν,\\
u_2' - cu_2
&= bu_1^2,
\end{align*}\]
where we have introduced the constants
\[\begin{gather*}
c = \frac{u_0}{ν}-\frac{n_0\mathcal{E}''(n_0)}{mνu_0}
=
\frac{u_0}{ν}-\frac{n_0\mu'(n_0)}{mνu_0},\\
b = \frac{1}{2ν} + \frac{2n_0\mu'(n_0)+n_0^2\mu''(n_0)}{2mνu_0^2}
= \frac{1}{2ν} + \frac{\bigl(nP'(n)\bigr)'|_{n=n_0}}{2mνu_0^2}.
\end{gather*}\]
Do it!#
To solve this problem, one needs to find a solution to the following inhomogeneous differential equation:
\[\begin{gather*}
au_1'(x) - u_1(x) = f(x),
\end{gather*}\]
where \(a>0\). We are particularly interested in the limit where \(a\rightarrow 0\), which will correspond to zero viscosity.
You should find that
\[\begin{gather*}
u_1(x) = e^{x/a}\int_{x_0}^{x} d{y}\; f(y)e^{-y/a}
\end{gather*}\]
Notice that this has a non-analytic dependence \(e^{x/a}\) on the parameter \(a \propto ν\). Show that, when expanding from the appropriate limit, all coefficients of the Taylor series for \(e^{x/a}\) in powers of \(a\) are zero. This makes it somewhat tricky to expand for small viscosities since all terms with such a factor will vanish.
Argue, however, that there may be an analytic particular solution \(u_1(x, a)\) whose form can be found by expanding in powers of \(a\). We can then write the general solution as
\[\begin{gather*}
u_1(x) = u_1(x, a) + Ce^{x/a}.
\end{gather*}\]
The non-analytic piece, thus describes transients (exponentially decaying in the appropriate direction) induced by the initial (boundary) conditions.
Show that this particular solution is:
\[\begin{gather*}
u_1(x, a) = - \sum_{n=0}^{\infty} a^{n} f^{(n)}(x).
\end{gather*}\]
For example, this allows us to easily solve for \(f(x) = \sin(kx)\):
\[\begin{gather*}
-u_1(x, a) = \sin(x) + ak \cos(x) - (ak)^2 \sin(x) - (ak)^3 \cos(x) + \cdots\\
= \Bigl(1 - (ak)^2 + (ak)^4 - \cdots\Bigr)\Bigl(\sin(kx) + ak \cos(kx)\Bigr)\\
= \frac{\sin(kx) + ak \cos(kx)}{1 + (ak)^2},\\
u_1(x) = Ce^{x/a} - \frac{\sin(kx) + ak \cos(kx)}{1+(ak)^2}.
\end{gather*}\]
We can easily check by differentiating that this satisfies the original equation.
Force#
The (buoyant) force acting on the potential \(V(x)\) is
\[\begin{gather*}
F = -\int_{0}^{L}d{x}\;nV'
= -Φ\int_{0}^{L}d{x}\;\frac{V'}{u}
\end{gather*}\]
Returning to our original equation, we have
\[\begin{gather*}
-nV'(x) = m\overbrace{nu}^{Φ}u_{,x} + \overbrace{nμ_{,x}}^{[P(n)]_{,x}} - mνnu_{,xx},\\
F = \Bigl.(mΦu + P)\Bigr|_{0}^{L} - mν\int_0^{L}d{x}\;nu_{,xx}
\end{gather*}\]
A Numerical Example#
To make sure what we are doing makes sense, we give a concrete example. To minimize the chance for error, we implement the original equations in the form \(y=(n, u, u', F)\), allowing for a density-dependent viscosity:
\[\begin{gather*}
(nu)' = 0, \qquad
muu' + \Bigl(V + \mathcal{E}'(n)\Bigr)' = \frac{(mν(n)u')'}{n}
= m\frac{ν'(n)n'u'+ν(n)u''}{n}\\
\begin{pmatrix}
n\\
u\\
u'\\
F
\end{pmatrix}'
=
\begin{pmatrix}
-\frac{nu'}{u}\\
u'\\
\frac{n\Bigl(uu' + \frac{V'(x) + \mathcal{E}''(n)n'}{m}\Bigr)-ν'n'u'}{ν}\\
nV'(x)
\end{pmatrix}, \qquad
y(0) = \begin{pmatrix}
n_0\\
u_0\\
0\\
0
\end{pmatrix}.
\end{gather*}\]
assert False
%matplotlib inline
import matplotlib
from functools import partial
import numpy as np, matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
matplotlib.rcParams['figure.figsize'] = (4,3)
m = 1.8
nu0 = 0.032
L = 1.3
n0 = 0.85
u0 = 0.773/2
Φ = n0*u0
g_m = 1.5
lam = 0.02
def mu(n, d=0):
"""Return the GPE chemical potential derivatives."""
if d == 0:
return m*g_m*n
elif d == 1:
return m*g_m+0*n
elif d == 2:
return 0*n
def mu(n, d=0):
"""Return the UFG chemical potential derivatives."""
if d == 0:
return m*g_m*n**(5/3)
elif d == 1:
return 5/3*m*g_m*n**(2/3)
elif d == 2:
return 10/9*m*g_m*n**(-1/3)
def mu_(n, d=0):
"""Return the GPE chemical potential derivatives."""
if d == 0:
return m*g_m*n**3
elif d == 1:
return 3*m*g_m*n**2
elif d == 2:
return 6*m*g_m*n
def V(x, d=0):
"""Return the external potential or derivative."""
if d == 0:
return np.where(abs(2*x/L-1) < 1, m*x*(x-L), 0)
elif d == 1:
return np.where(abs(2*x/L-1) < 1, m*(2*x-L), 0)
elif d == 2:
return np.where(abs(2*x/L-1) < 1, m*2, 0)
def V(x, d=0):
"""Return the external potential or derivative."""
if d == 0:
return np.where(abs(2*x/L-1) < 1,
m*np.sin(np.pi*(x-L)/L)**2, 0)
elif d == 1:
return np.where(abs(2*x/L-1) < 1,
np.pi/L*m*np.sin(2*np.pi*(x-L)/L), 0)
elif d == 2:
return np.where(abs(2*x/L-1) < 1,
2*np.pi**2/L**2*m*np.cos(2*np.pi*(x-L)/L), 0)
elif d == 3:
return np.where(abs(2*x/L-1) < 1,
-4*np.pi**3/L**3*m*np.sin(2*np.pi*(x-L)/L), 0)
def V(x, d=0):
"""Return the external potential or derivative."""
s = np.sin(np.pi*(x-L)/L)
ds = np.cos(np.pi*(x-L)/L)*np.pi/L
ds2 = -np.sin(np.pi*(x-L)/L)*(np.pi/L)**2
ds3 = -np.cos(np.pi*(x-L)/L)*(np.pi/L)**3
if d == 0:
return np.where(abs(2*x/L-1) < 1,
m*s**4, 0)
elif d == 1:
return np.where(abs(2*x/L-1) < 1,
m*4*s**3*ds, 0)
elif d == 2:
return np.where(abs(2*x/L-1) < 1,
m*(12*s**2*ds**2 + 4*s**3*ds2), 0)
elif d == 3:
return np.where(abs(2*x/L-1) < 1,
m*(24*s*ds**3 + 36*s**2*ds*ds2 + 4*s**3*ds3), 0)
def nu(n, nu0=1.0, gamma=0.0, d=0):
"""Return the density-dependent viscosity."""
if d == 0:
return nu0*n**gamma
elif d == 1:
return nu0*gamma*n**(gamma-1)
def compute_dy_dt(x, y, lam=lam, nu0=nu0,
u0=u0, n0=n0, m=m):
n, u, du, F = y
dn = -n*du/u
dV = lam*V(x, d=1)
#ddu = (m*u*du + dV + mu(n, d=1)*dn)/m/nu
nu_, dnu_ = [nu(n, nu0=nu0, d=d) for d in [0, 1]]
ddu = (n*(m*u*du + dV + mu(n, d=1)*dn)/m - du*dnu_*dn)/nu_
dF = n*dV
return (dn, du, ddu, dF)
y0 = (n0, u0, 0, 0)
res = solve_ivp(
partial(compute_dy_dt, nu0=0.2, lam=0.5),
y0=y0,
t_span=(-L, 5*L), max_step=0.001)
x = res.t
n, u, du, F = res.y
assert np.allclose(n*u, n0*u0)
fig, ax = plt.subplots()
ax.plot(x/L, n/n0, label="$n$")
ax.plot(x/L, u/u0, label="$u$")
ax.axvspan(0, 1, color='y', alpha=0.5)
ax.set(xlabel="$x/L$", ylabel="$n/n_0$", title=f"{n[0]=}, {n[-1]=}", );
plt.legend(loc="upper left")
ax = plt.twinx()
ax.plot(x/L, F, '--C1', label="$F$")
ax.legend(loc="lower right")
ax.set(ylabel="$F$");
Notice that the potential causes the density to change, then the viscosity allows this to relax outside of the shaded region where the potential is active. We now check some relationships.
# Check that the potential and mu derivatives are correct
fig, axs = plt.subplots(1, 2)
ax = axs[0]
ax.plot(x/L, V(x), 'C0')
ax.plot(x/L, V(x, d=1), 'C1-')
ax.plot(x/L, np.gradient(V(x), x), 'C1:')
ax.plot(x/L, V(x, d=2), 'C2-')
ax.plot(x/L, np.gradient(V(x, d=1), x), 'C2:')
ax.plot(x/L, V(x, d=3), 'C3-')
ylim = ax.get_ylim()
ax.plot(x/L, np.gradient(V(x, d=2), x), 'C3:')
ax.set(ylim=ylim)
ax = axs[1]
ns = np.linspace(n0/2, 2*n0)
ax.plot(ns/n0, mu(ns), 'C0')
ax.plot(ns/n0, mu(ns, d=1), 'C1-')
ax.plot(ns/n0, np.gradient(mu(ns), ns), 'C1:')
ax.plot(ns/n0, mu(ns, d=2), 'C2-')
ax.plot(ns/n0, np.gradient(mu(ns, d=1), ns), 'C2:')
from functools import partial
def get_F(lam=lam, nu=nu, u0=u0, n0=n0, m=m, L=L, tol=1e-7):
y0 = (n0, u0, 0, 0)
dy_dt = partial(compute_dy_dt, lam=lam, nu=nu,
u0=u0, n0=n0, m=m)
res = solve_ivp(dy_dt,
y0=y0, t_span=(0, L),
method='BDF',
rtol=tol, atol=tol)
x = res.t
n, u, du, F = res.y
return F[-1]
lams = np.linspace(0, 0.5)[1:-1]
nus = np.linspace(0, 0.5)[1:-1]
Fs_lam = [get_F(lam=lam) for lam in lams]
Fs_nu = [get_F(nu=nu) for nu in nus]
fig, axs = plt.subplots(1, 2, figsize=(5, 2))
axs[0].plot(lams**2, Fs_lam)
axs[0].set(xlabel=r"$\lambda^2$", ylabel="$F$")
axs[1].plot(nus, Fs_nu)
axs[1].set(xlabel=r"$\nu$", ylabel="$F$");
Now let’s check the small \(λ\) approximations:
from functools import partial
from scipy.integrate import cumtrapz
c = u0/nu - n0*mu(n0, d=1)/m/nu/u0
b = (1 + (2*n0*mu(n0, d=1) + n0**2*mu(n0, d=2))/m/u0**2)/2/nu
a = c*nu
print(f"{a=}, {nu*b}")
def f(x, us):
u1, u2 = us
du1 = V(x)/m/nu + c*u1
du2 = b*u1**2 + c*u2
return (du1, du2)
res1 = solve_ivp(f, y0=[0, 0], t_span=(-L, 2*L), t_eval=x,
method='BDF', atol=1e-8, rtol=1e-8)
u1, u2 = res1.y
n1 = Φ/(u0+lam*u1)
F1 = cumtrapz(lam*n1*V(x, d=1), x, initial=0)
# Small nu approximations
u10 = -V(x)/m/a
u11 = -V(x, d=1)/m/a**2
u12 = -V(x, d=2)/m/a**3
fig, axs = plt.subplots(2, 2, figsize=(10, 5))
ax = axs[0,0]
ax.plot(x/L, u1/u0, label='$u_1/u_0$')
ax.plot(x/L, u10/u0, label=r"small $\nu$ approx")
ax.plot(x/L, (u-u0)/u0/lam, label="full")
ax.set(xlabel=r"$x/L$", ylabel="$u_1/u_0$")
ax.legend()
ax = axs[0, 1]
ax.plot(x/L, F1, label='$F_1$')
ax.plot(x/L, F, label="full")
ax.legend()
ax = axs[1, 0]
ax.plot(x/L, u1/u0, 'C0-')
ax.plot(x/L, u10/u0, 'C0:')
ax.plot(x/L, (u1-u10)/nu/u0, 'C1-')
ax.plot(x/L, u11/u0, 'C1:')
ax.plot(x/L, u12/u0, 'C2:')
ylim = ax.get_ylim()
ax.plot(x/L, (u1-u10-nu*u11)/nu**2/u0, 'C2-', scaley=False)
ax.set(ylim=ylim)
ax = axs[1, 1]
ax.plot(x/L, (u-u0)/lam/u0, 'C0-')
ax.plot(x/L, u1/u0, 'C0:')
ax.plot(x/L, (u-u0-lam*u1)/lam**2/u0, 'C1-')
ax.plot(x/L, u2/u0, 'C1:')
The following code demonstrates that the order \(λ\) solution \(u_1(x)\) gives the correct force when integrated. (There is no need to compute \(u_2(x)\).)
def get_F1(lam=lam, nu=nu, u0=u0, n0=n0, m=m, L=L, tol=1e-7):
a = u0 - n0*mu(n0, d=1)/m/u0
c = a/nu
def f(x, y):
u1, F1 = y
du1 = V(x)/m/nu + c*u1
dF = -lam**2 * n0 / u0 * u1 * V(x, d=1)
return (du1, dF)
res1 = solve_ivp(f,
y0=[0, 0], t_span=(0, L),
method='BDF', atol=tol, rtol=tol)
u1, F1 = res1.y
return F1[-1]
lams = np.linspace(0.01, 0.3, 20)
qs = [get_F1(lam=lam, tol=lam**2 * 1e-8) /
get_F(lam=lam, tol=lam**2 * 1e-8) for lam in lams]
plt.plot(lams, qs, '-+')
plt.plot([0], [1], 'o')
plt.gca().set(xlabel=r'$\lambda$', ylabel='$F_1/F$');
Thus, the complete solution to order \(λ\) is to solve the following problem:
\[\begin{gather*}
u_1'(x) + cu_1(x) = V(x)/mν,\qquad
c = \frac{u_0 - P'_0/mu_0}{ν}\\
F %= λ\int_{0}^{L}d(x)\;\frac{Φ V'(x)}{u_0+λu_1}
= \frac{-λ^2n_0}{u_0} \int_{0}^{L}d{x}\;u_1(x) V'(x)
= \frac{λ^2n_0}{u_0} \int_{0}^{L}d{x}\;u_1'(x) V(x)
\end{gather*}\]
Note that the coefficient \(c\) here is generally less than zero. We can see this by noting that the speed of sound in a material is \(c_s\) is:
\[\begin{gather*}
c_s = \sqrt{\frac{n_0\mu'(n_0)}{m}} = \sqrt{\frac{P'(n_0)}{m}},\\
c = \frac{u_0-c_s^2/u_0}{\nu} = \frac{u_0}{\nu}\left(1 - \frac{c_s^2}{u_0^2}\right).
\end{gather*}\]
Thus, for subsonic flow \(u_0<c_s\), \(c<0\). We keep this sign convention to avoid many other signs below. The homogeneous solution is \(u_1(x) = Ae^{cx}\). To find the general solution, there are several options:
One is the method of variation of parameters can be used, promoting \(A → A(x)\) to find
\[\begin{gather*} u_1(x) = A(x)e^{cx}, \qquad A'(x) e^{cx} = V(x)/{mν}, \qquad A(x) = \frac{1}{mν}\int_0^xd{y}\; e^{-cy}V(y), \\ u_1(x) = \frac{1}{mν}\int_0^xd{y}\; e^{c(x-y)}V(y),\\ F = \frac{-λ^2n_0}{mu_0}\int_{0}^{L}d{x}\;\int_0^xd{y}\; e^{c(x-y)}V(y)V'(x). \end{gather*}\]The integral is over the lower triangle \(y<x\) of the box \((x, y) ∈[0, L]×[0,L]\). To better express the dimensionality of the solution, we write the potential as \(V(x) = V_0U(x/L)\) so that the shape is described by the dimensionless function \(U(\xi)\), and we introduce the dimensionless quantity \(\kappa \propto \nu\):
\[\begin{gather*} λV(x) = λV_0U(x/L), \qquad \kappa = \frac{1}{-cL} = \frac{\nu}{-aL},\\ \begin{aligned} F &= \frac{-λ^2n_0V_0^2L}{mu_0} \underbrace{\int_{0}^{1}d{x}\;\int_0^xd{y}\; e^{(y-x)/\kappa}U(y)U'(x)}_{-\kappa^2 f(\kappa, U)},\\ &= \frac{Lu_0\kappa^2}{\nu}\frac{λ^2V_0^2n_0}{mu_0^2}f(\kappa, U). \end{aligned} \end{gather*}\]While formally correct, this solution obscures the viscosity dependence. We will see below that the normalization chosen here for \(f(\kappa, U)\) makes sense in the limit of small viscosity:
\[\begin{gather*} f(\kappa, U) = -\int_{0}^{1}d{x}\;\int_0^xd{y}\; \frac{e^{(y-x)/\kappa}}{\kappa^2}U(y)U'(x) \end{gather*}\]Another is formal manipulation. We first introduce the parameter \(a=cν\) which is independent of the viscosity, then solve algebraically treating the derivative operator as a matrix:
\[\begin{align*} νu_1'(x) &= au_1(x) + V(x)/m, \\ u_1(x) &= \frac{1}{\left(1-\frac{ν}{a}\frac{d}{d{x}}\right)}\frac{-V(x)}{am},\\ &= -\frac{V(x)}{am} - \frac{νV'(x)}{a^2m} - \frac{ν^2V''(x)}{a^3m} + O(ν^3). \end{align*}\]This is nice because it immediately allows us to expand in powers of \(ν\) by expanding \(1/(1-x) =1+x+x^2+x^3+\cdots\):
\[\begin{gather*} F = \frac{λ^2n_0 V_0^2}{mau_0} \int_{0}^{1}d{x}\;\left( U(x)U'(x) - \kappa[U'(x)]^2 + \kappa^2U''(x)U'(x) + O(ν^3) \right) \end{gather*}\]Since \(U(0)=U(1)=0\), the first term, which is a total derivative \((U^2)'/2\), vanishes. This gives:
\[\begin{align*} F &= \frac{mλ^2n_0}{au_0}\frac{V_0^2}{m^2}(-\kappa)\overbrace{ \Biggl( \int_{0}^{1}d{x}\;[U'(x)]^2 -\kappa\left.\frac{[U'(x)]^2}{2}\right|_{0}^{1} +\kappa^2\int_{0}^{1}d{x}\;U'(x)U'''(x) + O(ν^3) \Biggr)}^{f(\kappa, U)}\\ &= \frac{mλ^2n_0}{au_0}\frac{ν}{a} \frac{V_0^2}{m^2L} f(\kappa, U). \end{align*}\]This justifies the form of \(f(\kappa, U)\) chosen above. If we choose a symmetric potential so that \(|U'(0)| = |U'(1)|\), then the second term vanishes, allowing us to precisely define what a small viscosity means:
\[\begin{gather*} \frac{ν^2}{a^2} \frac{\int_{0}^{L}d{x}\;V'(x)V'''(x)} {\int_{0}^{L}d{x}\;[V'(x)]^2} = \kappa^2 \frac{\int_{0}^{1}d{x}\;U'(x)U'''(x)} {\int_{0}^{1}d{x}\;[U'(x)]^2} \ll 1. \end{gather*}\]
u10 = -V(x)/a/m
u11 = -1/a * V(x, d=1)/a/m
u12 = -1/a**2 * V(x, d=2)/a/m
u13 = -1/a**3 * V(x, d=3)/a/m
plt.plot(x/L, (u1 - u10)/nu)
plt.plot(x/L, (u1 - u10 - nu*u11)/nu**2)
plt.plot(x/L, (u1 - u10 - nu*u11 - nu**2*u12)/nu**3)
coeff = m*lam**2*n0/a/u0*nu/a
t1 = np.trapz(V(x, d=1)**2, x)/m**2
t2 = np.trapz(V(x, d=1)*V(x, d=2), x)/m**2/a
t3 = np.trapz(V(x, d=1)*V(x, d=3), x)/m**2/a**2
print(f"{t1=}, {t2=}, {t3=}")
F_ = get_F(tol=1e-13)/coeff
F1_ = get_F1(tol=1e-13)/coeff
print(F_, (F_-F1_)/lam)
print(F1_, t1)
print((F1_ - t1)/nu,
(F1_ - t1 - nu*t2)/nu**2,
(F1_ - t1 - nu*t2 - nu**2*t3)/nu**3,
)
We now replace the parameter \(a\) with the speed of sound \(c_s\). The final result, to leading order in \(λ\) is:
\[\begin{gather*}
F = \overbrace{\frac{ν}{u_0L}}^{1/\mathrm{Re}}\frac{1}{\left(\frac{c_s^2}{u_0^2}-1\right)^2}\frac{λ^2V_0^2n_0}{mu_0^2}
f(\kappa, U), \qquad
\mathrm{Re} = \frac{u_0L}{ν}, \qquad
\kappa = \overbrace{\frac{\nu}{u_0L}}^{1/\mathrm{Re}}\frac{1}{\left(\frac{c_s^2}{u_0^2}-1\right)}, \qquad
c_s = \sqrt{\frac{P'(n_0)}{m}},\\
F = \frac{1}{\left(\frac{c_s^2}{u_0^2}-1\right)}\frac{λ^2V_0^2n_0}{mu_0^2}
\Bigl(\kappa f(\kappa, U)\Bigr)
\end{gather*}\]
Here we have expressed the results in terms of the dimensionless Reynold’s number \(\mathrm{Re}\) and the ratio of the speed of sound \(c_s\) to the flow \(u_0\). Note that the force diverges as the flow speed approaches the speed of sound. The dependence on the shape appears with the dimensionless parameter \(\kappa = ν/(-aL)\):
\[\begin{align*}
f(\kappa, U) &= \frac{-1}{\kappa^2}\int_0^{L} dx\int_0^{x}d{y}\;e^{(y-x)/\kappa}
U(y)U'(x),\\
&\approx \sum_{n=1}^{\infty}(-\kappa)^{n-1}
\int_0^1d{x}\;U^{(n)}(x)U'(x).
\end{align*}\]
Caution: The expansion in \(\kappa \propto \nu\) here justifies our choice for the form of \(f(\kappa, U)\) and the corresponding separation of dimensionful parameters, but it turns out that this series does not converge to the solution. We now demonstrate this. Consider for example, the dimensionless potential
\[\begin{align*}
U(x) &= \sin(\pi(x-1)), &
f(\kappa, U) &= \frac{\pi^2}{2}
\frac{1+\pi^2\kappa^2 - 2\kappa(1+e^{-1/\kappa})}
{(1+\pi^2\kappa^2)^2},\\
U(x) &= \sin^2(\pi(x-1)), &
f(\kappa, U) &= \frac{\pi^2}{2}
\frac{\Bigl(
1 + 4\pi^2\kappa^2
+ 8\pi^2\kappa^3(1 - e^{-1/\kappa})\Big)}
{(1+ 4\pi^2\kappa^2)^2}
\end{align*}\]
This solution contains the non-analytic term \(e^{-1/\kappa}\) whose Taylor series for positive \(\kappa\) is zero. Thus, the series expansion converges to the incorrect form of:
\[\begin{gather*}
\sum_{n=1}^{\infty}(-\kappa)^{n-1}
\int_0^1d{x}\;U^{(n)}(x)U'(x) =
\frac{\pi^2}{2}\left(1 - (\pi\kappa)^2 + (\pi\kappa)^4 - \cdots\right)\\
= \frac{\pi^2}{2}\frac{1}{1+\pi^2\kappa^2}
\end{gather*}\]
fails to converge.
def V(x, d=0):
"""Return the external potential or derivative."""
if d == 0:
return np.where(abs(2*x/L-1) < 1,
np.sin(np.pi*(x-L)/L), 0)
elif d == 1:
return np.where(abs(2*x/L-1) < 1,
np.pi/L*np.cos(np.pi*(x-L)/L), 0)
elif d == 2:
return np.where(abs(2*x/L-1) < 1,
-np.pi**2/L**2*np.sin(np.pi*(x-L)/L), 0)
elif d == 3:
return np.where(abs(2*x/L-1) < 1,
-np.pi**3/L**3*np.cos(np.pi*(x-L)/L), 0)
def get_f(nu, u0=u0, n0=n0, V0=1, L=L, m=m, lam=lam):
cs2 = n0*mu(n0, d=1)/m
coeff = nu/u0/L / (cs2/u0**2-1)**2 * lam**2 * V0**2 *n0/m/u0**2
F1 = get_F1(nu=nu, u0=u0, n0=n0, L=L, m=m, lam=lam, tol=1e-8)
kappa = nu/u0/L/(cs2/u0**2-1)
return kappa, F1/coeff
nus = np.linspace(0.1, 10.0)
kappas, fs = np.transpose([get_f(nu=nu) for nu in nus])
xs = (np.pi*kappas)**2
plt.plot(xs, fs, '.')
plt.plot(xs, np.pi**2/2/(1+xs), label='Series')
plt.plot(xs, np.pi**2/2
* (1 + xs - 2*kappas*(1+np.exp(-1/kappas)))
/ (1+xs)**2, label='Full integral')
plt.semilogy(xs, np.pi**2/2
* (1 + xs - 2*kappas*(1+0*np.exp(-1/kappas)))
/ (1+xs)**2, label='Integral w/o non-analytic piece')
plt.legend()
plt.gca().set(xlabel=r"$\pi^2\kappa^2$", ylabel="$f(\kappa, U)$");
This figure shows that the integral gives the correct formula, and that for large \(\kappa\) the non-analytic piece is correct. It also shows that we have made a mistake in our series approximation. Please fix!
Here is \(U(x) = \sin^2\bigl(\pi(x-1)\bigr)\):
def V(x, d=0):
"""Return the external potential or derivative."""
if d == 0:
return np.where(abs(2*x/L-1) < 1,
np.sin(np.pi*(x-L)/L)**2, 0)
elif d == 1:
return np.where(abs(2*x/L-1) < 1,
np.pi/L*np.sin(2*np.pi*(x-L)/L), 0)
elif d == 2:
return np.where(abs(2*x/L-1) < 1,
2*np.pi**2/L**2*np.cos(2*np.pi*(x-L)/L), 0)
elif d == 3:
return np.where(abs(2*x/L-1) < 1,
-4*np.pi**3/L**3*np.sin(2*np.pi*(x-L)/L), 0)
nus = np.linspace(0.1, 10.0)
kappas, fs = np.transpose([get_f(nu=nu) for nu in nus])
xs = (2*np.pi*kappas)**2
plt.plot(xs, fs, '.')
plt.plot(xs, np.pi**2/2/(1+xs), label='Series')
plt.plot(xs, np.pi**2/2
* (1 + xs + 2*xs*kappas*(1-np.exp(-1/kappas)))
/ (1+xs)**2, label='Full integral')
plt.semilogy(xs, np.pi**2/2
* (1 + xs + 2*xs*kappas*(1 - 0*np.exp(-1/kappas)))
/ (1+xs)**2, label='Integral w/o non-analytic piece')
plt.legend()
plt.gca().set(xlabel=r"$\pi^2\kappa^2$", ylabel="$f(\kappa, U)$");
import sympy
x, y, k = sympy.var('x,y,kappa', positive=True)
def U(x):
return sympy.sin(sympy.pi*(x-1))
integrand = -1/k**2 * sympy.exp((y-x)/k)*U(x).diff(x)*U(y)
f = integrand.integrate((y, 0, x)).integrate((x, 0, 1)).simplify()
display(f.factor())
display(sympy.series((1+(sympy.pi*k)**2 - 2*k)/(1+(sympy.pi*k)**2), k, 0, 10))
f_series = sum(sympy.simplify(
(U(x).diff(*(x,)*n)*U(x).diff(x))
.integrate((x, 0, 1))*(-k)**(n-1))
for n in range(1, 9)[::2])
display(f_series)
If we choose \(U(\xi) = U(1-\xi)\) to be symmetric, then the even terms vanish, and we can write:
\[\begin{gather*}
f(\kappa, U) = \sum_{n=0}^{\infty}c_{n}\kappa^{2n}, \qquad
c_n = \int_0^1d{\xi}\;U^{(2n+1)}(\xi)U'(\xi).
\end{gather*}\]
Here are a couple of special cases. Note: For numerical work, it might be better to use a potential that is sufficiently flat at the boundaries, hence the higher power. This is what I have used in the code. (Please check these.)
\[\begin{gather*}
f(\kappa, U) = \pi^2g(\pi^2\kappa^2),\qquad
\kappa = \frac{\nu/u_0L}{\Bigl(1-\frac{c_s^2}{u_0^2}\Bigr)},\\
\begin{aligned}
U(\xi) = \sin\bigl(\pi(\xi-1)\bigr)&:&
c_n &= \frac{\pi^2}{2}(-\pi^2)^{n}, &
g(x) &= \frac{1}{2(1+x)},\\
U(\xi) = \sin^4\bigl(\pi(\xi-1)\bigr)&:&
c_n &= \frac{\pi^2}{2}(4^{2n-1} + 4^{n})(-\pi^2)^{n}, &
g(x) &= \frac{68x + 5}{8(1+4x)(1+16x)}.
\end{aligned}
\end{gather*}\]
Thus, for a sinusoidal potential (please check):
\[\begin{gather*}
\begin{aligned}
\frac{F}{L} &=
\frac{n_0 λ^2V_0^2}{mu_0^2}\frac{1}{\Bigl(1 - \frac{c_s^2}{u_0^2}\Bigr)}\frac{\pi^2\kappa}{2(1+\pi^2\kappa^2)}, &
U(\xi) &= \sin\bigl(\pi(\xi-1)\bigr),\\
\frac{F}{L} &=
\frac{n_0 λ^2V_0^2}{mu_0^2}\frac{1}{\Bigl(1 - \frac{c_s^2}{u_0^2}\Bigr)}\frac{\pi^2\kappa(68\pi^2\kappa^2+5)}{8(1+4\pi^2\kappa^2)(1+16\pi^2\kappa^2)}, &
U(\xi) &= \sin^4\bigl(\pi(\xi-1)\bigr),\\
\end{aligned}
\end{gather*}\]
These can also be obtained by :
\[\begin{gather*}
%-integrate(integrate(e^{(x-y)/kappa}*sin(pi*(y-1))*cos(pi*(x-1))*pi/kappa^2, y=0..x), x=0..1)
f(\kappa, U) = \frac{-1}{\kappa}\int_0^L d{x}\int_0^x d{y}\; e^{(x-y)/\kappa} U'(x)U(y),\\
\frac{\pi^2}{2}\frac{(1 + 2 (1 + e^{1/\kappa}) κ + π^2 κ^2)}{(1 + \pi^2 \kappa^2)^2}
\end{gather*}\]
import sympy
x, y, k = sympy.var('x,y,kappa', positive=True)
def U(x):
return sympy.sin(sympy.pi*(x-1))**2
integrand = -1/k*sympy.exp((x-y)/k)*U(x).diff(x)*U(y)
display(integrand)
f = integrand.integrate((y, 0, x)).integrate((x, 0, 1))
f.expand().factor()
sympy.series(sympy.exp(1/k), k, 0, 5)
(I found these from the expansion https://oeis.org/A092812):
\[\begin{gather*}
f(\kappa, U) = \pi^2 \sum_{n=0}^{\infty}
\left(\frac{16^{n}}{8} + \frac{4^{n}}{2}\right)(-\pi^2\kappa^2)^{n}
= \pi^2 g(\pi^2\kappa^2)
, \\
g(x) %= \frac{1 + 16x + 24x^2}{(1+4x)(1+16x)} - \frac{3}{8}
= \frac{68x + 5}{8(1+4x)(1+16x)}
\end{gather*}\]
import numpy as np
from scipy.integrate import quad
import sympy
x, k = sympy.var('x, kappa', real=True)
V = sympy.sin(sympy.pi*(x-1))
def get_c(n):
return k**(2*n)*(V.diff(*(x,)*(2*n+1))*V.diff(x)).integrate((x, 0, 1))
[get_c(n) for n in [0, 1, 2, 3, 4, 5]]
import numpy as np
from scipy.integrate import quad
import sympy
x, k = sympy.var('x, kappa', real=True)
V = sympy.sin(sympy.pi*(x-1))**4
def get_c(n):
f = sympy.lambdify([x], V.diff(*(x,)*(n+1))*V.diff(x))
return int(np.round(quad(f, 0, 1)[0]) / np.pi**(n+2)) * k**(n)*sympy.pi**(n+2)
get_c(4)
4**(0-1)/2 + 2**(0-1), 5/8
sympy.simplify((1+16*x+24*x**2)/(1+4*x)/(1+16*x) - sympy.S(3)/8).factor()
from scipy.integrate import quad
def f(x):
return V(x, d=1)**2
cs = np.sqrt(n0*mu(n0, d=1)/m)
Re = u0*L/nu
tmp = lam**2 * quad(f, 0, L)[0]
print(a/u0, (1 - (cs/u0)**2))
print(
(1/Re) * (1/(1-(cs/u0)**2)**2) * (n0*L/m/u0**2) * tmp,
get_F())
from scipy.integrate import quad
def f(x):
return V(x)*V(x,d=1)**2
tmp = quad(f, 0, L)[0]
lam**2*n0*nu/m/a**2/u0 * tmp, get_F()
>>>>>>> variant B
---
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.13.7
kernelspec:
display_name: Python 3 (ipykernel)
language: python
name: python3
---
```{code-cell}
import mmf_setup; mmf_setup.nbinit()
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)
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None
Effective Viscosity (Old)#
Conservation Laws#
Particle Number#
The first hydrodynamic equation demonstrates the conservation of particle number
\[\begin{gather*}
N = \int\d{x}\; n(x),\\
\dot{N} = \int\d{x}\; \dot{n}(x)
= \int\d{x}\; (-nu)_{,x}
= -\Bigl.nu\Bigr|_{x_L}^{x_R}
= \Phi_{n}(x_L) - \Phi_n(x_R).
\end{gather*}\]
This says that that total rate of change of particles in a region between \(x_L\) and \(x_R\) is equal to the flux of particles \(\Phi_{n}(x_L)\) entering the region from the left at \(x_L\) minus the flux of particles \(\Phi_{n}(x_R)\) leaving the region at \(x_R\) on the right. Notice that a positive flux corresponds to a flow to the right, consistent with the signs in this equation.
To Do. What is the meaning of \(D_t n = -nu_{,x}\)?
What is the physical meaning of the convective derivative of the density?
\[\begin{gather*}
D_t n = \dot{n} + u n_{,x} = -nu_{,x}.
\end{gather*}\]
Why is this not zero?
Solution
The convective derivative characterizes how a quantity changes as you follow the particle flow. If you are following a particle, how does the density change? Well, this depends on how close the neighbouring particles are. The average distance to the neighbours changes iff there is a gradient in the velocity field. To see this, consider two particles at \(x_0 < x_1\) so that
\[\begin{gather*}
n \propto \frac{1}{x_1-x_0} = \frac{1}{\delta},\\
\dot{n} \propto -\frac{\dot{x}_1-\dot{x}_0}{(x_1-x_0)^2}
= -\frac{\dot{x}_1-\dot{x}_0}{\delta^2} \approx -\frac{u_{,x}}{\delta} = -n u_{,x}.
\end{gather*}\]
I.e. the density only changes if the velocity of the relative particles is different – this is exactly what a gradient in the velocity indicates.
Momentum Density#
Multiply the continuity equation by \(u\) and add it to \(n\) times the Bernoulli equation:
\[\begin{gather*}
u\dot{n} + u(nu)_{,x} +
n\dot{u} + nuu_{,x} + \frac{n}{m}\Bigl(
V(x) + \mathcal{E}'(n)\Bigr)_{,x}
- (ν u_{,x})_{,x} = 0\\
(nu)_{,t} + \left(nu^2 + \frac{P}{m} - ν u_{,x}\right)_{,x} = 0,\\
\end{gather*}\]
This means that, if we have a static solution, then the following is constant everywhere:
\[\begin{gather*}
nu^2 + \frac{P}{m} - ν u_{,x}
= nu^2 + \frac{F(x) + n\mathcal{E}' - \mathcal{E}}{m} - ν
u_{,x}
= \text{const.},\\
P = F(x) + n\mathcal{E}'(n) - \mathcal{E}(n), \qquad
F(x) = \int_{x_0}^{x} n(x)V'(x)\d{x}.
\end{gather*}\]
Unfortunately, as mentioned above, the force depends on the solution through the density, so we can use this to check and interpret results, but it does not help us solve anything.
Energy#
We can similarly consider the energy flux. The energy density is the sum of the kinetic, potential, and internal energies.
\[\begin{gather*}
\varepsilon = \frac{m}{2}nu^2 + nV + \mathcal{E}(n),\\
\dot{\varepsilon} + \left[
u\left(\frac{m}{2}nu^2 + nV + n\mathcal{E}'(n)\right)
\right]_{,x} = n \dot{V} + mu(\nu u_{,x})_{,x},\\
\dot{\varepsilon} + \left[
u\left(\frac{m}{2}nu^2 + nV + n\mathcal{E}'(n)
- m\nu u_{,x}
\right)
\right]_{,x} = n \dot{V} - m\nu u_{,x}^2.
\end{gather*}\]
Here we see explicitly that energy conservation is violated by time-dependent potentials, and by a viscosity if the velocity field has a gradient \(u_{,x} \neq 0\). Thus, there will be no dissipation from this viscosity if we have simple sloshing in a harmonic trap where the velocity is constant. The second form gives a current density whose gradient is strictly decreasing if viscosity is active and \(\dot{V}=0\).
An Analytic Example#
Consider static flow: \(\dot{n} = \dot{u} = 0\). The density and momentum continuity equations give the following in terms of the constants \(j\) and \(p\):
\[\begin{gather*}
un = j, \qquad
mnu^2 + F(x) + n \mathcal{E}'(n) - \mathcal{E}(n) - m\nu u_{,x} = mp
= m \frac{j^2}{n_L} + P(n_L).
\end{gather*}\]
Consider a step-like potential with small \(\epsilon\rightarrow 0\):
\[\begin{gather*}
V(x) = \begin{cases}
V_L & x < 0,\\
0 & \epsilon < x,
\end{cases} \qquad
V_{,x}(x) \approx -V_L\delta(x).
\end{gather*}\]
Working backwards, let \(u(x)\) be linear in the region \([0, \epsilon]\)
\[\begin{gather*}
u'(x) = \begin{cases}
\frac{u_R - u_L}{\epsilon} & 0 < x < \epsilon\\
0 & \text{otherwise}
\end{cases}.
\end{gather*}\]
In this interval, we have
\[\begin{gather*}
u(x) = u_L + x\frac{u_R-u_L}{\epsilon}, \qquad
n(x) = \frac{j}{u(x)} = \frac{u_L n_L}{u_L + x\frac{u_R-u_L}{\epsilon}},\\
\end{gather*}\]
For \(x \gg \epsilon\) we have
\[\begin{gather*}
n(x \gg \epsilon) = n_L + \delta n, \qquad
u(x \gg \epsilon) = u_L + \delta u = \frac{j}{n_L + \delta n}, \quad
F(x \gg \epsilon) = \delta F,
\end{gather*}\]
where
\[\begin{gather*}
m\frac{j^2}{n_L + \delta n} + \delta F + P(n_L + \delta n) = m \frac{j^2}{n_L} + P(n_L),\\
\delta F = \int_{0}^{\epsilon} nV'(x)\d{x}
= -V_L(n_L + \alpha\delta n)
\in [-n_LV_L, -(n_L + \delta n)V_L],\\
0 \leq \alpha \leq 1.
\end{gather*}\]
For small \(\delta n\) we can write
\[\begin{gather*}
P(n_L + \delta n) - P(n_L) \approx \delta n P'(n_L) + \frac{(\delta n)^2}{2}P''(n_L),
\end{gather*}\]
which is exact for the GPE \(P(n) = gn^2/2\). This gives:
\[\begin{gather*}
\delta n = \frac{V_L n_L}{
P'(n_L) + \frac{\delta n}{2} P''(n_L)
- \frac{mj^2}{(n_L + \delta n)n_L}
- V_L\alpha}.
\end{gather*}\]
This can be iterated a couple of times to converge. If \(\delta n\) is small, we can expand:
\[\begin{gather*}
\frac{\delta n}{n_L} \approx \frac{1}{
\frac{n_L\mathcal{E}''(n_L) - mu_L^2}{V_L} - \alpha}
\end{gather*}\]
Let’s start with a perfect fluid without viscosity \(\nu=0\) flowing over a step The momentum continuity equation gives
\[\begin{gather*}
mnu^2 + F(x) + n\mathcal{E}'(n) - \mathcal{E}(n) - m\nu u_{,x} = \text{const.}
\end{gather*}\]
Starting with \(n(x<0) = n_L\), \(u(x<0) = u_L\), and \(F(x<0) = 0\) with \(j=u_Ln_L\) constant, we have
\[\begin{gather*}
m\frac{j^2}{n} + F(x) + n\mathcal{E}'(n) - \mathcal{E}(n) - m\nu u_{,x} =
m\frac{j^2}{n_L} + n_L\mathcal{E}'(n_L) - \mathcal{E}(n_L).
\end{gather*}\]
Integrating the momentum conservation equation across the potential we have:
\[\begin{gather*}
mj_0^2 n_{,x} = n^3 (V+\mathcal{E}'')_{,x}, \qquad
mj_0^2\delta n = n_0^2 (-V_L + \mathcal{E}''(n)\delta
\end{gather*}\]
A Numerical Example: Steady-State Flow#
We start with steady-state flow over a barrier. To make sure what we are doing makes sense, we give a concrete example. To minimize the chance for error, we implement the original equations in the form \(y=(n, u, u', F)\), looking for stationary solutions \(\dot{n} = \dot{u} = 0\):
\[\begin{gather*}
(nu)_{,x} = 0, \qquad
muu_{,x} = -(V+\mathcal{E}')_{,x} + m(\nu u_{,x})_{,x}/n\\
n_{,x} = -\frac{nu_{,x}}{u}, \qquad
u_{,xx} = \frac{mnuu_{,x} + n(V+\mathcal{E}')_{,x}}{m\nu}\\
\begin{pmatrix}
n\\
u\\
u_{,x}\\
F
\end{pmatrix}_{,x}
=
\begin{pmatrix}
-\frac{nu_{,x}}{u}\\
u_{,x}\\
n\frac{uu_{,x} + \frac{(V+\mathcal{E}')_{,x}}{m}}{\nu}\\
nV_{,x}
\end{pmatrix}, \qquad
y(0) = \begin{pmatrix}
n_0\\
u_0\\
0\\
0
\end{pmatrix}.
\end{gather*}\]
Removing the Velocity
Alternatively, we can use the conservation equation to remove the velocity from our equations since \(nu = n_0u_0 = j_0 = \text{const.}\) everywhere:
\[\begin{gather*}
u = \frac{j_0}{n}, \quad
u_{,x} = \frac{-j_0 n_{,x}}{n^2}, \quad
u_{,xx} = \frac{-j_0 n_{,xx}}{n^2}+\frac{2j_0 (n_{,x})^2}{n^3},\\
n_{,xx} = \frac{2(n_{,x})^2}{n} - \frac{-j_0 n_{,x} + n^3\frac{(V+\mathcal{E}')_{,x}}{mj_0}}{\nu}.
\end{gather*}\]
The momentum current is
\[\begin{gather*}
j_{p} = \frac{j_0^2}{n} + \frac{P}{m} + \nu\frac{j_0 n_{,x}}{n^2}
\end{gather*}\]
We will use the GPE equation of state:
\[\begin{gather*}
\mathcal{E} = \frac{gn^2}{2}, \qquad
\mathcal{E}' = gn, \qquad
\mathcal{E}'' = g,
\end{gather*}\]
and a piecewise linear potential.
import scipy as sp
x = np.linspace(-0.1, 1.1)
V0 = 1.2
dV = 0.2
xp, fp = [0, 0.5, 1.0], [dV, V0, 0]
V = sp.interpolate.InterpolatedUnivariateSpline(xp, fp, k=1, ext='const')
plt.plot(x, V(x), label='$V(x)$')
plt.plot(x, V.derivative()(x), label="$V'(x)$")
plt.legend()
from scipy.integrate import solve_ivp
import scipy as sp
class Flow:
m = 1.0
g = 1.0
nu = 0.1 # Viscosity
V0 = 0.2 # Potential barrier height
L = 0.2 # Length of barrier
n0 = 1.0 # Steady-state density on the left.
v0 = 0.1 # Steady-state velocity on the left.
dV = 0.1 # Potential drop (magnitude)
# Boundary conditions on the left
def __init__(self, **kw):
for key in kw:
if not hasattr(self, key):
raise ValueError(f"Unknown {key=}")
setattr(self, key, kw[key])
self.init()
def init(self):
xp, fp = [0, self.L/2, self.L], [self.dV, self.V0, 0]
self._V = sp.interpolate.InterpolatedUnivariateSpline(xp, fp, k=1, ext='const')
self._dV = self._V.derivative()
def get_Eint(self, n, d=0):
if d == 0:
return self.g * n**2 /2
elif d == 1:
return self.g * n
elif d == 2:
return self.g
else:
return 0
def compute_dy_dx(self, x, y):
n, u, du, F = y
dn = -n*du/u
ddu = n * (u*du + (self._dV(x) + self.get_Eint(n, d=2)*dn)/self.m) / self.nu
dF = n * self._dV(x)
return (dn, du, ddu, dF)
def solve(self, dv0=0.0, x0_L=-0.1, x_L=2.0):
x_span = (x0_L*self.L, x_L*self.L)
y0 = (self.n0, self.v0, dv0, 0.0)
res = solve_ivp(self.compute_dy_dx, y0=y0, t_span=x_span, max_step=0.001)
assert(res.success)
res.x = res.t
res.n, res.v, res.dv, res.F = res.y
res.dn, dv, res.ddv, res.dF = self.compute_dy_dx(res.x, res.y)
self.res = res
return res
f = Flow()
res = f.solve()
Analysis#
We start with the two conservation laws for density and momentum far to the left and right of the potential, assuming that \(u_{,x}=0\) in these regions, and setting \(F(x_L) = 0\) and \(F(x_R) = F\):
\[\begin{gather*}
n_Lu_L = n_Ru_R = j = \text{const.}\\
mn_Lu_L^2 + n_L\mathcal{E}'(n_L) - \mathcal{E}(n_L) =
mn_Ru_R^2 + F + n_R\mathcal{E}'(n_R) - \mathcal{E}(n_R)\\
\frac{m j^2}{n_L} + n_L\mathcal{E}'(n_L) - \mathcal{E}(n_L) =
\frac{m j^2}{n_R} + F + n_R\mathcal{E}'(n_R) - \mathcal{E}(n_R)\\
\end{gather*}\]
\[\begin{gather*}
nu^2 + \frac{P}{m} - ν u_{,x}
= nu^2 + \frac{F(x) + n\mathcal{E}' - \mathcal{E}}{m} - ν
u_{,x}
= \text{const.},\\
P = F(x) + n\mathcal{E}'(n) - \mathcal{E}(n), \qquad
F(x) = \int_{x_0}^{x} n(x)V'(x)\d{x}.
\end{gather*}\]
Assuming that \(u_{,x} = 0\) far from the potential (static flow), and setiing \(V(x_R) = 0\), we have the following after integrating the energy current:
\[\begin{gather*}
\left.
u\left(
\frac{m}{2}nu^2 + nV + n\mathcal{E}'(n) - m\nu u_{,x}
\right)
\right|^{x_R}_{x_L}
= - \int_{x_L}^{x_R} m\nu u_{,x}^2 \d{x} \leq 0,\\
\tfrac{m}{2}(u_{R}^2 - u_{L}^2) + V_L + \mathcal{E}'(n_R) -
\mathcal{E}'(n_L) =
-\frac{1}{j}\int_{x_L}^{x_R} m\nu u_{,x}^2 \d{x}
\end{gather*}\]
Correction Needed Below#
The following analysis is incorrect. It is missing the factor of \(1/n\) discussed above, and neglects to include a difference in the potential from one side to the other (DC offset) required to maintain a persistent flow in the presence of dissipation.
Steady State (Incorrect)#
Now consider steady-state flow \(n(x, t) = n(x)\), \(u(x, t) = u(x)\). These become:
\[\begin{gather*}
(nu)_{,x} = 0, \qquad
muu_{,x} + V_{,x} + \mu'(n)n_{,x} = mνu_{,xx}.
\end{gather*}\]
The first equation establishes that the flux \(Φ = nu\) is the same at every point, allowing us to exchange velocity \(u\) for density \(n\):
\[\begin{gather*}
n = \frac{Φ}{u}, \qquad
n_{,x} = -\frac{Φ}{u^2}u_{,x}.
\end{gather*}\]
Let’s assume that \(V(x)\) has compact support over \([0, L]\), then we can give boundary conditions \(u(0) = u_0\), \(u_{,x}(0) = 0\). The second equation and be integrated once to give the non-linear first-order equation
\[\begin{gather*}
mνu_{,x} = \frac{mu^2}{2} + V(x) + \overbrace{\mu\left(\frac{Φ}{u}\right)}^{\mu(n)}.
\end{gather*}\]
For general \(V(x)\), this has no analytic solution, but we can compute a perturbative solution if we let \(V(x) \rightarrow \lambda V(x)\) where \(\lambda \ll 1\). We note that \(u(x) = u_0\) is a solution for \(\lambda = 0\) and write $\(u(x) = u_0 + \lambda u_1(x) + \lambda^2 u_2(x) + O(\lambda^3).\)\(_Note: In what follows, we revert back to the notation_ \)u_1’ = [u_1]{,x}$…
I originally wrote it this way, and have not changed it yet._
As shown below, we originally tried expanding to first order in \(\lambda\), but found this was insufficient, so we expand to second order in \(\lambda\) obtaining the following linear equations for \(u_1\) and \(u_2\):
\[\begin{align*}
mνu_1' + \overbrace{\left(\frac{n_0}{u_0}\mu'(n_0) - mu_0\right)}^{-mνc}u_1
&= V(x),\\
mνu_2' + \left(\frac{n_0}{u_0}\mu'(n_0) - mu_0\right)u_2
&= \underbrace{\left(\frac{m}{2}
+ \frac{2n_0\mu'(n_0)+n_0^2\mu''(n_0)}{2u_0^2}
\right)}_{mνb}u_1^2,\\
u_1' - c u_1
&= V(x)/mν,\\
u_2' - cu_2
&= bu_1^2,
\end{align*}\]
where we have introduced the constants
\[\begin{gather*}
c = \frac{u_0}{ν}-\frac{n_0\mathcal{E}''(n_0)}{mνu_0}
=
\frac{u_0}{ν}-\frac{n_0\mu'(n_0)}{mνu_0},\\
b = \frac{1}{2ν} + \frac{2n_0\mu'(n_0)+n_0^2\mu''(n_0)}{2mνu_0^2}
= \frac{1}{2ν} + \frac{\bigl(nP'(n)\bigr)'|_{n=n_0}}{2mνu_0^2}.
\end{gather*}\]
Do it!#
To solve this problem, one needs to find a solution to the following inhomogeneous differential equation:
\[\begin{gather*}
au_1'(x) - u_1(x) = f(x),
\end{gather*}\]
where \(a>0\). We are particularly interested in the limit where \(a\rightarrow 0\), which will correspond to zero viscosity.
You should find that
\[\begin{gather*}
u_1(x) = e^{x/a}\int_{x_0}^{x} d{y}\; f(y)e^{-y/a}
\end{gather*}\]
Notice that this has a non-analytic dependence \(e^{x/a}\) on the parameter \(a \propto ν\). Show that, when expanding from the appropriate limit, all coefficients of the Taylor series for \(e^{x/a}\) in powers of \(a\) are zero. This makes it somewhat tricky to expand for small viscosities since all terms with such a factor will vanish.
Argue, however, that there may be an analytic particular solution \(u_1(x, a)\) whose form can be found by expanding in powers of \(a\). We can then write the general solution as
\[\begin{gather*}
u_1(x) = u_1(x, a) + Ce^{x/a}.
\end{gather*}\]
The non-analytic piece, thus describes transients (exponentially decaying in the appropriate direction) induced by the initial (boundary) conditions.
Show that this particular solution is:
\[\begin{gather*}
u_1(x, a) = - \sum_{n=0}^{\infty} a^{n} f^{(n)}(x).
\end{gather*}\]
For example, this allows us to easily solve for \(f(x) = \sin(kx)\):
\[\begin{gather*}
-u_1(x, a) = \sin(x) + ak \cos(x) - (ak)^2 \sin(x) - (ak)^3 \cos(x) + \cdots\\
= \Bigl(1 - (ak)^2 + (ak)^4 - \cdots\Bigr)\Bigl(\sin(kx) + ak \cos(kx)\Bigr)\\
= \frac{\sin(kx) + ak \cos(kx)}{1 + (ak)^2},\\
u_1(x) = Ce^{x/a} - \frac{\sin(kx) + ak \cos(kx)}{1+(ak)^2}.
\end{gather*}\]
We can easily check by differentiating that this satisfies the original equation.
Force#
The (buoyant) force acting on the potential \(V(x)\) is
\[\begin{gather*}
F = -\int_{0}^{L}d{x}\;nV'
= -Φ\int_{0}^{L}d{x}\;\frac{V'}{u}
\end{gather*}\]
Returning to our original equation, we have
\[\begin{gather*}
-nV'(x) = m\overbrace{nu}^{Φ}u_{,x} + \overbrace{nμ_{,x}}^{[P(n)]_{,x}} - mνnu_{,xx},\\
F = \Bigl.(mΦu + P)\Bigr|_{0}^{L} - mν\int_0^{L}d{x}\;nu_{,xx}
\end{gather*}\]
A Numerical Example#
To make sure what we are doing makes sense, we give a concrete example. To minimize the chance for error, we implement the original equations in the form \(y=(n, u, u', F)\), allowing for a density-dependent viscosity:
\[\begin{gather*}
(nu)' = 0, \qquad
muu' + \Bigl(V + \mathcal{E}'(n)\Bigr)' = \frac{(mν(n)u')'}{n}
= m\frac{ν'(n)n'u'+ν(n)u''}{n}\\
\begin{pmatrix}
n\\
u\\
u'\\
F
\end{pmatrix}'
=
\begin{pmatrix}
-\frac{nu'}{u}\\
u'\\
\frac{n\Bigl(uu' + \frac{V'(x) + \mathcal{E}''(n)n'}{m}\Bigr)-ν'n'u'}{ν}\\
nV'(x)
\end{pmatrix}, \qquad
y(0) = \begin{pmatrix}
n_0\\
u_0\\
0\\
0
\end{pmatrix}.
\end{gather*}\]
%matplotlib inline
import matplotlib
from functools import partial
import numpy as np, matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
matplotlib.rcParams['figure.figsize'] = (4,3)
m = 1.8
nu0 = 0.032
L = 1.3
n0 = 0.85
u0 = 0.773/2
Φ = n0*u0
g_m = 1.5
lam = 0.02
def mu(n, d=0):
"""Return the GPE chemical potential derivatives."""
if d == 0:
return m*g_m*n
elif d == 1:
return m*g_m+0*n
elif d == 2:
return 0*n
def mu(n, d=0):
"""Return the UFG chemical potential derivatives."""
if d == 0:
return m*g_m*n**(5/3)
elif d == 1:
return 5/3*m*g_m*n**(2/3)
elif d == 2:
return 10/9*m*g_m*n**(-1/3)
def mu_(n, d=0):
"""Return the GPE chemical potential derivatives."""
if d == 0:
return m*g_m*n**3
elif d == 1:
return 3*m*g_m*n**2
elif d == 2:
return 6*m*g_m*n
def V(x, d=0):
"""Return the external potential or derivative."""
if d == 0:
return np.where(abs(2*x/L-1) < 1, m*x*(x-L), 0)
elif d == 1:
return np.where(abs(2*x/L-1) < 1, m*(2*x-L), 0)
elif d == 2:
return np.where(abs(2*x/L-1) < 1, m*2, 0)
def V(x, d=0):
"""Return the external potential or derivative."""
if d == 0:
return np.where(abs(2*x/L-1) < 1,
m*np.sin(np.pi*(x-L)/L)**2, 0)
elif d == 1:
return np.where(abs(2*x/L-1) < 1,
np.pi/L*m*np.sin(2*np.pi*(x-L)/L), 0)
elif d == 2:
return np.where(abs(2*x/L-1) < 1,
2*np.pi**2/L**2*m*np.cos(2*np.pi*(x-L)/L), 0)
elif d == 3:
return np.where(abs(2*x/L-1) < 1,
-4*np.pi**3/L**3*m*np.sin(2*np.pi*(x-L)/L), 0)
def V(x, d=0):
"""Return the external potential or derivative."""
s = np.sin(np.pi*(x-L)/L)
ds = np.cos(np.pi*(x-L)/L)*np.pi/L
ds2 = -np.sin(np.pi*(x-L)/L)*(np.pi/L)**2
ds3 = -np.cos(np.pi*(x-L)/L)*(np.pi/L)**3
if d == 0:
return np.where(abs(2*x/L-1) < 1,
m*s**4, 0)
elif d == 1:
return np.where(abs(2*x/L-1) < 1,
m*4*s**3*ds, 0)
elif d == 2:
return np.where(abs(2*x/L-1) < 1,
m*(12*s**2*ds**2 + 4*s**3*ds2), 0)
elif d == 3:
return np.where(abs(2*x/L-1) < 1,
m*(24*s*ds**3 + 36*s**2*ds*ds2 + 4*s**3*ds3), 0)
def nu(n, nu0=1.0, gamma=0.0, d=0):
"""Return the density-dependent viscosity."""
if d == 0:
return nu0*n**gamma
elif d == 1:
return nu0*gamma*n**(gamma-1)
def compute_dy_dt(x, y, lam=lam, nu0=nu0,
u0=u0, n0=n0, m=m):
n, u, du, F = y
dn = -n*du/u
dV = lam*V(x, d=1)
#ddu = (m*u*du + dV + mu(n, d=1)*dn)/m/nu
nu_, dnu_ = [nu(n, nu0=nu0, d=d) for d in [0, 1]]
ddu = (n*(m*u*du + dV + mu(n, d=1)*dn)/m - du*dnu_*dn)/nu_
dF = n*dV
return (dn, du, ddu, dF)
y0 = (n0, u0, 0, 0)
res = solve_ivp(
partial(compute_dy_dt, nu0=0.2, lam=0.5),
y0=y0,
t_span=(-L, 5*L), max_step=0.001)
x = res.t
n, u, du, F = res.y
assert np.allclose(n*u, n0*u0)
fig, ax = plt.subplots()
ax.plot(x/L, n/n0, label="$n$")
ax.plot(x/L, u/u0, label="$u$")
ax.axvspan(0, 1, color='y', alpha=0.5)
ax.set(xlabel="$x/L$", ylabel="$n/n_0$", title=f"{n[0]=}, {n[-1]=}", );
plt.legend(loc="upper left")
ax = plt.twinx()
ax.plot(x/L, F, '--C1', label="$F$")
ax.legend(loc="lower right")
ax.set(ylabel="$F$");
Notice that the potential causes the density to change, then the viscosity allows this to relax outside of the shaded region where the potential is active. We now check some relationships.
# Check that the potential and mu derivatives are correct
fig, axs = plt.subplots(1, 2)
ax = axs[0]
ax.plot(x/L, V(x), 'C0')
ax.plot(x/L, V(x, d=1), 'C1-')
ax.plot(x/L, np.gradient(V(x), x), 'C1:')
ax.plot(x/L, V(x, d=2), 'C2-')
ax.plot(x/L, np.gradient(V(x, d=1), x), 'C2:')
ax.plot(x/L, V(x, d=3), 'C3-')
ylim = ax.get_ylim()
ax.plot(x/L, np.gradient(V(x, d=2), x), 'C3:')
ax.set(ylim=ylim)
ax = axs[1]
ns = np.linspace(n0/2, 2*n0)
ax.plot(ns/n0, mu(ns), 'C0')
ax.plot(ns/n0, mu(ns, d=1), 'C1-')
ax.plot(ns/n0, np.gradient(mu(ns), ns), 'C1:')
ax.plot(ns/n0, mu(ns, d=2), 'C2-')
ax.plot(ns/n0, np.gradient(mu(ns, d=1), ns), 'C2:')
from functools import partial
def get_F(lam=lam, nu=nu, u0=u0, n0=n0, m=m, L=L, tol=1e-7):
y0 = (n0, u0, 0, 0)
dy_dt = partial(compute_dy_dt, lam=lam, nu=nu,
u0=u0, n0=n0, m=m)
res = solve_ivp(dy_dt,
y0=y0, t_span=(0, L),
method='BDF',
rtol=tol, atol=tol)
x = res.t
n, u, du, F = res.y
return F[-1]
lams = np.linspace(0, 0.5)[1:-1]
nus = np.linspace(0, 0.5)[1:-1]
Fs_lam = [get_F(lam=lam) for lam in lams]
Fs_nu = [get_F(nu=nu) for nu in nus]
fig, axs = plt.subplots(1, 2, figsize=(5, 2))
axs[0].plot(lams**2, Fs_lam)
axs[0].set(xlabel=r"$\lambda^2$", ylabel="$F$")
axs[1].plot(nus, Fs_nu)
axs[1].set(xlabel=r"$\nu$", ylabel="$F$");
Now let’s check the small \(λ\) approximations:
from functools import partial
from scipy.integrate import cumtrapz
c = u0/nu - n0*mu(n0, d=1)/m/nu/u0
b = (1 + (2*n0*mu(n0, d=1) + n0**2*mu(n0, d=2))/m/u0**2)/2/nu
a = c*nu
print(f"{a=}, {nu*b}")
def f(x, us):
u1, u2 = us
du1 = V(x)/m/nu + c*u1
du2 = b*u1**2 + c*u2
return (du1, du2)
res1 = solve_ivp(f, y0=[0, 0], t_span=(-L, 2*L), t_eval=x,
method='BDF', atol=1e-8, rtol=1e-8)
u1, u2 = res1.y
n1 = Φ/(u0+lam*u1)
F1 = cumtrapz(lam*n1*V(x, d=1), x, initial=0)
# Small nu approximations
u10 = -V(x)/m/a
u11 = -V(x, d=1)/m/a**2
u12 = -V(x, d=2)/m/a**3
fig, axs = plt.subplots(2, 2, figsize=(10, 5))
ax = axs[0,0]
ax.plot(x/L, u1/u0, label='$u_1/u_0$')
ax.plot(x/L, u10/u0, label=r"small $\nu$ approx")
ax.plot(x/L, (u-u0)/u0/lam, label="full")
ax.set(xlabel=r"$x/L$", ylabel="$u_1/u_0$")
ax.legend()
ax = axs[0, 1]
ax.plot(x/L, F1, label='$F_1$')
ax.plot(x/L, F, label="full")
ax.legend()
ax = axs[1, 0]
ax.plot(x/L, u1/u0, 'C0-')
ax.plot(x/L, u10/u0, 'C0:')
ax.plot(x/L, (u1-u10)/nu/u0, 'C1-')
ax.plot(x/L, u11/u0, 'C1:')
ax.plot(x/L, u12/u0, 'C2:')
ylim = ax.get_ylim()
ax.plot(x/L, (u1-u10-nu*u11)/nu**2/u0, 'C2-', scaley=False)
ax.set(ylim=ylim)
ax = axs[1, 1]
ax.plot(x/L, (u-u0)/lam/u0, 'C0-')
ax.plot(x/L, u1/u0, 'C0:')
ax.plot(x/L, (u-u0-lam*u1)/lam**2/u0, 'C1-')
ax.plot(x/L, u2/u0, 'C1:')
The following code demonstrates that the order \(λ\) solution \(u_1(x)\) gives the correct force when integrated. (There is no need to compute \(u_2(x)\).)
def get_F1(lam=lam, nu=nu, u0=u0, n0=n0, m=m, L=L, tol=1e-7):
a = u0 - n0*mu(n0, d=1)/m/u0
c = a/nu
def f(x, y):
u1, F1 = y
du1 = V(x)/m/nu + c*u1
dF = -lam**2 * n0 / u0 * u1 * V(x, d=1)
return (du1, dF)
res1 = solve_ivp(f,
y0=[0, 0], t_span=(0, L),
method='BDF', atol=tol, rtol=tol)
u1, F1 = res1.y
return F1[-1]
lams = np.linspace(0.01, 0.3, 20)
qs = [get_F1(lam=lam, tol=lam**2 * 1e-8) /
get_F(lam=lam, tol=lam**2 * 1e-8) for lam in lams]
plt.plot(lams, qs, '-+')
plt.plot([0], [1], 'o')
plt.gca().set(xlabel=r'$\lambda$', ylabel='$F_1/F$');
Thus, the complete solution to order \(λ\) is to solve the following problem:
\[\begin{gather*}
u_1'(x) + cu_1(x) = V(x)/mν,\qquad
c = \frac{u_0 - P'_0/mu_0}{ν}\\
F %= λ\int_{0}^{L}d(x)\;\frac{Φ V'(x)}{u_0+λu_1}
= \frac{-λ^2n_0}{u_0} \int_{0}^{L}d{x}\;u_1(x) V'(x)
= \frac{λ^2n_0}{u_0} \int_{0}^{L}d{x}\;u_1'(x) V(x)
\end{gather*}\]
Note that the coefficient \(c\) here is generally less than zero. We can see this by noting that the speed of sound in a material is \(c_s\) is:
\[\begin{gather*}
c_s = \sqrt{\frac{n_0\mu'(n_0)}{m}} = \sqrt{\frac{P'(n_0)}{m}},\\
c = \frac{u_0-c_s^2/u_0}{\nu} = \frac{u_0}{\nu}\left(1 - \frac{c_s^2}{u_0^2}\right).
\end{gather*}\]
Thus, for subsonic flow \(u_0<c_s\), \(c<0\). We keep this sign convention to avoid many other signs below. The homogeneous solution is \(u_1(x) = Ae^{cx}\). To find the general solution, there are several options:
One is the method of variation of parameters can be used, promoting \(A → A(x)\) to find
\[\begin{gather*} u_1(x) = A(x)e^{cx}, \qquad A'(x) e^{cx} = V(x)/{mν}, \qquad A(x) = \frac{1}{mν}\int_0^xd{y}\; e^{-cy}V(y), \\ u_1(x) = \frac{1}{mν}\int_0^xd{y}\; e^{c(x-y)}V(y),\\ F = \frac{-λ^2n_0}{mu_0}\int_{0}^{L}d{x}\;\int_0^xd{y}\; e^{c(x-y)}V(y)V'(x). \end{gather*}\]The integral is over the lower triangle \(y<x\) of the box \((x, y) ∈[0, L]×[0,L]\). To better express the dimensionality of the solution, we write the potential as \(V(x) = V_0U(x/L)\) so that the shape is described by the dimensionless function \(U(\xi)\), and we introduce the dimensionless quantity \(\kappa \propto \nu\):
\[\begin{gather*} λV(x) = λV_0U(x/L), \qquad \kappa = \frac{1}{-cL} = \frac{\nu}{-aL},\\ \begin{aligned} F &= \frac{-λ^2n_0V_0^2L}{mu_0} \underbrace{\int_{0}^{1}d{x}\;\int_0^xd{y}\; e^{(y-x)/\kappa}U(y)U'(x)}_{-\kappa^2 f(\kappa, U)},\\ &= \frac{Lu_0\kappa^2}{\nu}\frac{λ^2V_0^2n_0}{mu_0^2}f(\kappa, U). \end{aligned} \end{gather*}\]While formally correct, this solution obscures the viscosity dependence. We will see below that the normalization chosen here for \(f(\kappa, U)\) makes sense in the limit of small viscosity:
\[\begin{gather*} f(\kappa, U) = -\int_{0}^{1}d{x}\;\int_0^xd{y}\; \frac{e^{(y-x)/\kappa}}{\kappa^2}U(y)U'(x) \end{gather*}\]Another is formal manipulation. We first introduce the parameter \(a=cν\) which is independent of the viscosity, then solve algebraically treating the derivative operator as a matrix:
\[\begin{align*} νu_1'(x) &= au_1(x) + V(x)/m, \\ u_1(x) &= \frac{1}{\left(1-\frac{ν}{a}\frac{d}{d{x}}\right)}\frac{-V(x)}{am},\\ &= -\frac{V(x)}{am} - \frac{νV'(x)}{a^2m} - \frac{ν^2V''(x)}{a^3m} + O(ν^3). \end{align*}\]This is nice because it immediately allows us to expand in powers of \(ν\) by expanding \(1/(1-x) =1+x+x^2+x^3+\cdots\):
\[\begin{gather*} F = \frac{λ^2n_0 V_0^2}{mau_0} \int_{0}^{1}d{x}\;\left( U(x)U'(x) - \kappa[U'(x)]^2 + \kappa^2U''(x)U'(x) + O(ν^3) \right) \end{gather*}\]Since \(U(0)=U(1)=0\), the first term, which is a total derivative \((U^2)'/2\), vanishes. This gives:
\[\begin{align*} F &= \frac{mλ^2n_0}{au_0}\frac{V_0^2}{m^2}(-\kappa)\overbrace{ \Biggl( \int_{0}^{1}d{x}\;[U'(x)]^2 -\kappa\left.\frac{[U'(x)]^2}{2}\right|_{0}^{1} +\kappa^2\int_{0}^{1}d{x}\;U'(x)U'''(x) + O(ν^3) \Biggr)}^{f(\kappa, U)}\\ &= \frac{mλ^2n_0}{au_0}\frac{ν}{a} \frac{V_0^2}{m^2L} f(\kappa, U). \end{align*}\]This justifies the form of \(f(\kappa, U)\) chosen above. If we choose a symmetric potential so that \(|U'(0)| = |U'(1)|\), then the second term vanishes, allowing us to precisely define what a small viscosity means:
\[\begin{gather*} \frac{ν^2}{a^2} \frac{\int_{0}^{L}d{x}\;V'(x)V'''(x)} {\int_{0}^{L}d{x}\;[V'(x)]^2} = \kappa^2 \frac{\int_{0}^{1}d{x}\;U'(x)U'''(x)} {\int_{0}^{1}d{x}\;[U'(x)]^2} \ll 1. \end{gather*}\]
u10 = -V(x)/a/m
u11 = -1/a * V(x, d=1)/a/m
u12 = -1/a**2 * V(x, d=2)/a/m
u13 = -1/a**3 * V(x, d=3)/a/m
plt.plot(x/L, (u1 - u10)/nu)
plt.plot(x/L, (u1 - u10 - nu*u11)/nu**2)
plt.plot(x/L, (u1 - u10 - nu*u11 - nu**2*u12)/nu**3)
coeff = m*lam**2*n0/a/u0*nu/a
t1 = np.trapz(V(x, d=1)**2, x)/m**2
t2 = np.trapz(V(x, d=1)*V(x, d=2), x)/m**2/a
t3 = np.trapz(V(x, d=1)*V(x, d=3), x)/m**2/a**2
print(f"{t1=}, {t2=}, {t3=}")
F_ = get_F(tol=1e-13)/coeff
F1_ = get_F1(tol=1e-13)/coeff
print(F_, (F_-F1_)/lam)
print(F1_, t1)
print((F1_ - t1)/nu,
(F1_ - t1 - nu*t2)/nu**2,
(F1_ - t1 - nu*t2 - nu**2*t3)/nu**3,
)
We now replace the parameter \(a\) with the speed of sound \(c_s\). The final result, to leading order in \(λ\) is:
\[\begin{gather*}
F = \overbrace{\frac{ν}{u_0L}}^{1/\mathrm{Re}}\frac{1}{\left(\frac{c_s^2}{u_0^2}-1\right)^2}\frac{λ^2V_0^2n_0}{mu_0^2}
f(\kappa, U), \qquad
\mathrm{Re} = \frac{u_0L}{ν}, \qquad
\kappa = \overbrace{\frac{\nu}{u_0L}}^{1/\mathrm{Re}}\frac{1}{\left(\frac{c_s^2}{u_0^2}-1\right)}, \qquad
c_s = \sqrt{\frac{P'(n_0)}{m}},\\
F = \frac{1}{\left(\frac{c_s^2}{u_0^2}-1\right)}\frac{λ^2V_0^2n_0}{mu_0^2}
\Bigl(\kappa f(\kappa, U)\Bigr)
\end{gather*}\]
Here we have expressed the results in terms of the dimensionless Reynold’s number \(\mathrm{Re}\) and the ratio of the speed of sound \(c_s\) to the flow \(u_0\). Note that the force diverges as the flow speed approaches the speed of sound. The dependence on the shape appears with the dimensionless parameter \(\kappa = ν/(-aL)\):
\[\begin{align*}
f(\kappa, U) &= \frac{-1}{\kappa^2}\int_0^{L} dx\int_0^{x}d{y}\;e^{(y-x)/\kappa}
U(y)U'(x),\\
&\approx \sum_{n=1}^{\infty}(-\kappa)^{n-1}
\int_0^1d{x}\;U^{(n)}(x)U'(x).
\end{align*}\]
Caution: The expansion in \(\kappa \propto \nu\) here justifies our choice for the form of \(f(\kappa, U)\) and the corresponding separation of dimensionful parameters, but it turns out that this series does not converge to the solution. We now demonstrate this. Consider for example, the dimensionless potential
\[\begin{align*}
U(x) &= \sin(\pi(x-1)), &
f(\kappa, U) &= \frac{\pi^2}{2}
\frac{1+\pi^2\kappa^2 - 2\kappa(1+e^{-1/\kappa})}
{(1+\pi^2\kappa^2)^2},\\
U(x) &= \sin^2(\pi(x-1)), &
f(\kappa, U) &= \frac{\pi^2}{2}
\frac{\Bigl(
1 + 4\pi^2\kappa^2
+ 8\pi^2\kappa^3(1 - e^{-1/\kappa})\Big)}
{(1+ 4\pi^2\kappa^2)^2}
\end{align*}\]
This solution contains the non-analytic term \(e^{-1/\kappa}\) whose Taylor series for positive \(\kappa\) is zero. Thus, the series expansion converges to the incorrect form of:
\[\begin{gather*}
\sum_{n=1}^{\infty}(-\kappa)^{n-1}
\int_0^1d{x}\;U^{(n)}(x)U'(x) =
\frac{\pi^2}{2}\left(1 - (\pi\kappa)^2 + (\pi\kappa)^4 - \cdots\right)\\
= \frac{\pi^2}{2}\frac{1}{1+\pi^2\kappa^2}
\end{gather*}\]
fails to converge.
def V(x, d=0):
"""Return the external potential or derivative."""
if d == 0:
return np.where(abs(2*x/L-1) < 1,
np.sin(np.pi*(x-L)/L), 0)
elif d == 1:
return np.where(abs(2*x/L-1) < 1,
np.pi/L*np.cos(np.pi*(x-L)/L), 0)
elif d == 2:
return np.where(abs(2*x/L-1) < 1,
-np.pi**2/L**2*np.sin(np.pi*(x-L)/L), 0)
elif d == 3:
return np.where(abs(2*x/L-1) < 1,
-np.pi**3/L**3*np.cos(np.pi*(x-L)/L), 0)
def get_f(nu, u0=u0, n0=n0, V0=1, L=L, m=m, lam=lam):
cs2 = n0*mu(n0, d=1)/m
coeff = nu/u0/L / (cs2/u0**2-1)**2 * lam**2 * V0**2 *n0/m/u0**2
F1 = get_F1(nu=nu, u0=u0, n0=n0, L=L, m=m, lam=lam, tol=1e-8)
kappa = nu/u0/L/(cs2/u0**2-1)
return kappa, F1/coeff
nus = np.linspace(0.1, 10.0)
kappas, fs = np.transpose([get_f(nu=nu) for nu in nus])
xs = (np.pi*kappas)**2
plt.plot(xs, fs, '.')
plt.plot(xs, np.pi**2/2/(1+xs), label='Series')
plt.plot(xs, np.pi**2/2
* (1 + xs - 2*kappas*(1+np.exp(-1/kappas)))
/ (1+xs)**2, label='Full integral')
plt.semilogy(xs, np.pi**2/2
* (1 + xs - 2*kappas*(1+0*np.exp(-1/kappas)))
/ (1+xs)**2, label='Integral w/o non-analytic piece')
plt.legend()
plt.gca().set(xlabel=r"$\pi^2\kappa^2$", ylabel="$f(\kappa, U)$");
This figure shows that the integral gives the correct formula, and that for large \(\kappa\) the non-analytic piece is correct. It also shows that we have made a mistake in our series approximation. Please fix!
Here is \(U(x) = \sin^2\bigl(\pi(x-1)\bigr)\):
def V(x, d=0):
"""Return the external potential or derivative."""
if d == 0:
return np.where(abs(2*x/L-1) < 1,
np.sin(np.pi*(x-L)/L)**2, 0)
elif d == 1:
return np.where(abs(2*x/L-1) < 1,
np.pi/L*np.sin(2*np.pi*(x-L)/L), 0)
elif d == 2:
return np.where(abs(2*x/L-1) < 1,
2*np.pi**2/L**2*np.cos(2*np.pi*(x-L)/L), 0)
elif d == 3:
return np.where(abs(2*x/L-1) < 1,
-4*np.pi**3/L**3*np.sin(2*np.pi*(x-L)/L), 0)
nus = np.linspace(0.1, 10.0)
kappas, fs = np.transpose([get_f(nu=nu) for nu in nus])
xs = (2*np.pi*kappas)**2
plt.plot(xs, fs, '.')
plt.plot(xs, np.pi**2/2/(1+xs), label='Series')
plt.plot(xs, np.pi**2/2
* (1 + xs + 2*xs*kappas*(1-np.exp(-1/kappas)))
/ (1+xs)**2, label='Full integral')
plt.semilogy(xs, np.pi**2/2
* (1 + xs + 2*xs*kappas*(1 - 0*np.exp(-1/kappas)))
/ (1+xs)**2, label='Integral w/o non-analytic piece')
plt.legend()
plt.gca().set(xlabel=r"$\pi^2\kappa^2$", ylabel="$f(\kappa, U)$");
import sympy
x, y, k = sympy.var('x,y,kappa', positive=True)
def U(x):
return sympy.sin(sympy.pi*(x-1))
integrand = -1/k**2 * sympy.exp((y-x)/k)*U(x).diff(x)*U(y)
f = integrand.integrate((y, 0, x)).integrate((x, 0, 1)).simplify()
display(f.factor())
display(sympy.series((1+(sympy.pi*k)**2 - 2*k)/(1+(sympy.pi*k)**2), k, 0, 10))
f_series = sum(sympy.simplify(
(U(x).diff(*(x,)*n)*U(x).diff(x))
.integrate((x, 0, 1))*(-k)**(n-1))
for n in range(1, 9)[::2])
display(f_series)
If we choose \(U(\xi) = U(1-\xi)\) to be symmetric, then the even terms vanish, and we can write:
\[\begin{gather*}
f(\kappa, U) = \sum_{n=0}^{\infty}c_{n}\kappa^{2n}, \qquad
c_n = \int_0^1d{\xi}\;U^{(2n+1)}(\xi)U'(\xi).
\end{gather*}\]
Here are a couple of special cases. Note: For numerical work, it might be better to use a potential that is sufficiently flat at the boundaries, hence the higher power. This is what I have used in the code. (Please check these.)
\[\begin{gather*}
f(\kappa, U) = \pi^2g(\pi^2\kappa^2),\qquad
\kappa = \frac{\nu/u_0L}{\Bigl(1-\frac{c_s^2}{u_0^2}\Bigr)},\\
\begin{aligned}
U(\xi) = \sin\bigl(\pi(\xi-1)\bigr)&:&
c_n &= \frac{\pi^2}{2}(-\pi^2)^{n}, &
g(x) &= \frac{1}{2(1+x)},\\
U(\xi) = \sin^4\bigl(\pi(\xi-1)\bigr)&:&
c_n &= \frac{\pi^2}{2}(4^{2n-1} + 4^{n})(-\pi^2)^{n}, &
g(x) &= \frac{68x + 5}{8(1+4x)(1+16x)}.
\end{aligned}
\end{gather*}\]
Thus, for a sinusoidal potential (please check):
\[\begin{gather*}
\begin{aligned}
\frac{F}{L} &=
\frac{n_0 λ^2V_0^2}{mu_0^2}\frac{1}{\Bigl(1 - \frac{c_s^2}{u_0^2}\Bigr)}\frac{\pi^2\kappa}{2(1+\pi^2\kappa^2)}, &
U(\xi) &= \sin\bigl(\pi(\xi-1)\bigr),\\
\frac{F}{L} &=
\frac{n_0 λ^2V_0^2}{mu_0^2}\frac{1}{\Bigl(1 - \frac{c_s^2}{u_0^2}\Bigr)}\frac{\pi^2\kappa(68\pi^2\kappa^2+5)}{8(1+4\pi^2\kappa^2)(1+16\pi^2\kappa^2)}, &
U(\xi) &= \sin^4\bigl(\pi(\xi-1)\bigr),\\
\end{aligned}
\end{gather*}\]
These can also be obtained by :
\[\begin{gather*}
%-integrate(integrate(e^{(x-y)/kappa}*sin(pi*(y-1))*cos(pi*(x-1))*pi/kappa^2, y=0..x), x=0..1)
f(\kappa, U) = \frac{-1}{\kappa}\int_0^L d{x}\int_0^x d{y}\; e^{(x-y)/\kappa} U'(x)U(y),\\
\frac{\pi^2}{2}\frac{(1 + 2 (1 + e^{1/\kappa}) κ + π^2 κ^2)}{(1 + \pi^2 \kappa^2)^2}
\end{gather*}\]
import sympy
x, y, k = sympy.var('x,y,kappa', positive=True)
def U(x):
return sympy.sin(sympy.pi*(x-1))**2
integrand = -1/k*sympy.exp((x-y)/k)*U(x).diff(x)*U(y)
display(integrand)
f = integrand.integrate((y, 0, x)).integrate((x, 0, 1))
f.expand().factor()
sympy.series(sympy.exp(1/k), k, 0, 5)
(I found these from the expansion https://oeis.org/A092812):
\[\begin{gather*}
f(\kappa, U) = \pi^2 \sum_{n=0}^{\infty}
\left(\frac{16^{n}}{8} + \frac{4^{n}}{2}\right)(-\pi^2\kappa^2)^{n}
= \pi^2 g(\pi^2\kappa^2)
, \\
g(x) %= \frac{1 + 16x + 24x^2}{(1+4x)(1+16x)} - \frac{3}{8}
= \frac{68x + 5}{8(1+4x)(1+16x)}
\end{gather*}\]
import numpy as np
from scipy.integrate import quad
import sympy
x, k = sympy.var('x, kappa', real=True)
V = sympy.sin(sympy.pi*(x-1))
def get_c(n):
return k**(2*n)*(V.diff(*(x,)*(2*n+1))*V.diff(x)).integrate((x, 0, 1))
[get_c(n) for n in [0, 1, 2, 3, 4, 5]]
import numpy as np
from scipy.integrate import quad
import sympy
x, k = sympy.var('x, kappa', real=True)
V = sympy.sin(sympy.pi*(x-1))**4
def get_c(n):
f = sympy.lambdify([x], V.diff(*(x,)*(n+1))*V.diff(x))
return int(np.round(quad(f, 0, 1)[0]) / np.pi**(n+2)) * k**(n)*sympy.pi**(n+2)
get_c(4)
4**(0-1)/2 + 2**(0-1), 5/8
sympy.simplify((1+16*x+24*x**2)/(1+4*x)/(1+16*x) - sympy.S(3)/8).factor()
from scipy.integrate import quad
def f(x):
return V(x, d=1)**2
cs = np.sqrt(n0*mu(n0, d=1)/m)
Re = u0*L/nu
tmp = lam**2 * quad(f, 0, L)[0]
print(a/u0, (1 - (cs/u0)**2))
print(
(1/Re) * (1/(1-(cs/u0)**2)**2) * (n0*L/m/u0**2) * tmp,
get_F())
from scipy.integrate import quad
def f(x):
return V(x)*V(x,d=1)**2
tmp = quad(f, 0, L)[0]
lam**2*n0*nu/m/a**2/u0 * tmp, get_F()
PDE#
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
%load_ext autoreload
from gpe.contexts import FPS
from IPython.display import clear_output
%autoreload 2
import viscosity
import itertools
for (periodic, fft) in itertools.product(*([[False, True]]*2)):
f = viscosity.FlowStep(V_dV=0.2, V_sigma=1.0, V0=1.0, n0=1.0, v0=0.0,
nu=1.0, x0=0, x1=2*np.pi, Nx=200,
periodic=periodic, fft=fft,
ghost_points=1)
y0 = f.get_y0()
x = f.x
y = np.sin(x)
dy = np.cos(x)
ddy = -np.sin(x)
#assert np.allclose(f.get_derivative(y)[1:-1], dy[1:-1], rtol=1.7e-4)
#assert np.allclose(f.get_derivative(y, d=2)[1:-1], ddy[1:-1], rtol=9e-5)
#plt.plot((f.get_derivative(y, d=1) / dy - 1)[1:-1])
#plt.plot((f.get_derivative(y, d=2) / ddy - 1)[1:-1])
#plt.plot(f.x, f.get_derivative(y, d=1))
#plt.plot(f.x, dy)
#plt.ylim(-1, 1)
Dam Breaking#
from scipy.integrate import solve_ivp
f = viscosity.Flow(x0=-1.0, x1=1.0, Nx=100)
assert np.allclose(f._D2x, -f._Dx.T @ f._Dx)
x = f.x
n0 = np.where(x < 0, 1.0, 0.00001)
v0 = 0*n0
y0 = f.pack(n0, v0)
t_span = (0, 0.1)
res = solve_ivp(f.compute_dy_dt, y0=y0, t_span=t_span, method="BDF", rtol=1e-12, atol=1e-12)
if not res.success:
print(res.message)
print(res.t.max())
nF, vF = f.unpack(res.y[:, -1])
plt.plot(x, nF)
Burger’s Equation#
from importlib import reload
reload(viscosity)
f = viscosity.FlowBurgers(x0=0, x1=2*np.pi, Nx=256, nu=0.07)#, fft=True)
n0, u0 = f.get_exact_sol(f.x, t=0)
y0 = f.pack(n0*0, u0)
dx = 2*np.pi/99
dt = f.nu * dx
T = dt * 100
print(T)
t_span = (0, T)
res = solve_ivp(f.compute_dy_dt, y0=y0, t_span=t_span, method="BDF")
if not res.success:
print(res.message)
print(res.t.max())
nF, vF = f.unpack(res.y[:, -1])
plt.plot(f.x, vF)
plt.plot(f.x, nF)
n_, v_ = f.get_exact_sol(f.x, t=T)
plt.plot(f.x, v_)
fig, ax = plt.subplots()
ax1 = ax.twinx()
axs = [ax, ax1]
n_range = [0, 0]
u_range = [0, 0]
for t, y in zip(res.t, res.y.T):
n, u = f.unpack(y)
n_range = [min(n.min(), n_range[0]), max(n.max(), n_range[1])]
u_range = [min(u.min(), u_range[0]), max(u.max(), u_range[1])]
skip = 1
fps = FPS(zip(res.t[::skip], res.y.T[::skip]), timeout=1000, display=fig, embed=True, filename="tmp.m4v")
for t, y in fps:
n, u = f.unpack(y)
[ax.cla() for ax in axs]
ax.plot(f.x, n)
ax1.plot(f.x, u, '--C1')
ax.set(title=f"{t=:.4f}", ylim=n_range)
ax1.set(ylim=u_range)
def step0(u, dt=dx * f.nu):
nx = len(u)
un = u
u = u.copy()
for i in range(1, nx-1):
u[i] = (un[i] + dt * ((-un[i] * (un[i] - un[i-1]) / dx
+ f.nu * (un[i+1] - 2 * un[i] + un[i-1]) / dx**2)))
u[0] = un[0] - un[0] * dt / dx * (un[0] - un[-2]) + f.nu * dt / dx**2 *\
(un[1] - 2 * un[0] + un[-2])
u[-1] = u[0]
return u
def step0(u, dt=dx * f.nu):
nx = len(u)
un = u
u = u.copy()
du = u.copy()
ddu = u.copy()
for i in range(1, nx-1):
du[i] = (un[i+1] - un[i]) / dx
ddu[i] = (un[i+1] - 2 * un[i] + un[i-1]) / dx**2
for i in range(1, nx-1):
du[i] = (un[i] - un[i-1]) / dx
#ddu[i] = (un[i+1] - 2 * un[i] + un[i-1]) / dx**2
u[i] = (un[i] + dt * (-un[i] * du[i] + f.nu * ddu[i]))
u[0] = un[0] - un[0] * dt / dx * (un[0] - un[-2]) + f.nu * dt / dx**2 *\
(un[1] - 2 * un[0] + un[-2])
u[-1] = u[0]
return u
def step1(u, dt=dx * f.nu):
du = f.get_derivative(u)
ddu = f.get_derivative(u, d=2)
du_dt = (-u * du + f.nu * ddu)
return u + du_dt * dt
def step(u, dt=dx * f.nu):
du_dt = f.unpack(f.compute_dy_dt(0, f.pack(u,u)))[1]
return u + du_dt*dt
us = [u0]
t = 0
for n in range(100):
t += dt
us.append(step0(us[-1]))
plt.plot(f.x, us[-1])
plt.plot(f.x, f.get_exact_sol(f.x, t=t)[1])
from scipy.integrate import solve_ivp
f = viscosity.FlowBurgers(x0=-1.0, x1=1.0, Nx=100)
x = f.x
n0 = np.where(x < 0, 1.0, 0.5)
v0 = 0*n0
y0 = f.pack(n0, v0)
t_span = (0, 1)
res = solve_ivp(f.compute_dy_dt, y0=y0, t_span=t_span, method="BDF", rtol=1e-12, atol=1e-12)
if not res.success:
print(res.message)
print(res.t.max())
nF, vF = f.unpack(res.y[:, -1])
plt.plot(x, nF)
fig, ax = plt.subplots()
ax1 = ax.twinx()
axs = [ax, ax1]
n_range = [0, 0]
u_range = [0, 0]
for t, y in zip(res.t, res.y.T):
n, u = f.unpack(y)
n_range = [min(n.min(), n_range[0]), max(n.max(), n_range[1])]
u_range = [min(u.min(), u_range[0]), max(u.max(), u_range[1])]
skip = 1
fps = FPS(zip(res.t[::skip], res.y.T[::skip]), timeout=1000, display=fig, embed=True, filename="tmp.m4v")
for t, y in fps:
n, u = f.unpack(y)
[ax.cla() for ax in axs]
ax.plot(f.x, n)
ax1.plot(f.x, u, '--C1')
ax.set(title=f"{t=:.4f}", ylim=n_range)
ax1.set(ylim=u_range)
fig, axs = plt.subplots(2, 1, height_ratios=(3, 1), gridspec_kw=dict(hspace=0))
ax = axs[0]
axV = axs[1]
ax1 = ax.twinx()
axs = [ax, axV, ax1]
n_range = [0, 0]
u_range = [0, 0]
for t, y in zip(res.t, res.y.T):
n, u = f.unpack(y)
n_range = [min(n.min(), n_range[0]), max(n.max(), n_range[1])]
u_range = [min(u.min(), u_range[0]), max(u.max(), u_range[1])]
skip = 1
for t, y in FPS(zip(res.t[::skip], res.y.T[::skip]), timeout=1000, display=fig):
n, u = f.unpack(y)
[ax.cla() for ax in axs]
ax.plot(f.x, n)
axV.plot(f.x, f.get_Vext(f.x, t))
ax1.plot(f.x, u, '--C1')
ax.set(title=f"{t=:.4f}", ylim=n_range)
ax1.set(ylim=u_range)
import viscosity
f = viscosity.Flow(V_dV=-0.2, V_sigma=1.0, V0=1.0, n0=2.0, v0=0.0)
f = viscosity.Flow(V_dV=0, V_sigma=0.1, V0=0, n0=2.0, v0=0.0, V_w=1, nu=0.1, nu_alpha=1)
y0 = f.solve_ivp(x0=-3, x1=3, Nx=200, ghost=1)
x = f.x
n0 = np.maximum(0, 1 - f.get_Vext(x, t=0))/f.g + 1e-5
v0 = 0*n0
y0 = f.pack(n0, v0)
from scipy.integrate import solve_ivp
t0 = 0
t1 = 0.5
res = solve_ivp(f.compute_dy_dt, y0=y0, t_span=(t0, t1), method='BDF')
if not res.success:
print(res.message)
print(res.t.max())
\[\begin{split}
f(x \pm h) = f \pm h f' + \frac{h^2}{2} f'' \pm \frac{h^3}{6} f'''\\
8\frac{f(x + h) - f(x-h)}{12h} - \frac{f(x + 2h) - f(x-2h)}{12h} = f'\\
\end{split}\]
%pylab inline
from scipy.integrate import solve_ivp
from gpe.contexts import FPS
Lx = 10.0
Nx = 256//2
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 for constant j = j_avg
du = j_avg/n - u
du_dt += j_corr_factor * du
return pack(dn_dt, du_dt)
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)
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.8
t = 0
b = get_V_m(t, x) / g_m
h = n0
eta = h - b
n = eta
u = 0*x + u0
y0 = pack(n, u)
dt = dx / nu * 0.002
j_corr_factor = 0
fig, axs = plt.subplots(2, 1, sharex=True)
ax = axs[0]
ax_c = ax.twinx()
ax_stream = axs[1]
ax_j = ax_stream.twinx()
y = y0.copy()
for n in FPS(80, fig=fig, display=True, 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, 1.0),
ylabel="$j$ (blue) and $u$ (orange)")
dt = dx / nu * 0.002
j_corr_factor = 0
fig, axs = plt.subplots(2, 1, sharex=True)
ax = axs[0]
ax_c = ax.twinx()
ax_stream = axs[1]
ax_j = ax_stream.twinx()
for n in FPS(80, fig=fig, display=True, 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, 1.0),
ylabel="$j$ (blue) and $u$ (orange)")
import sympy
x = sympy.symbols('x', real=True)
u_exact = sympy.exp((x**4-1)/4)
def deq(u):
return (u.diff(x, x) - (x**6 + 3*x**2)*u)
assert deq(u_exact).simplify() == 0
assert u_exact.subs([(x, -1)]) == 1
assert u_exact.subs([(x, 1)]) == 1
N = 2
an = sympy.symbols(f'a:{N}', real=True)
kn = tuple(sympy.pi * (2*sympy.S(n) + 1)/2 for n in range(N))
u = 1 + sum(a_ * sympy.cos(k_ * x) for a_, k_ in zip(an, kn))
import numpy as np
from scipy.optimize import curve_fit
def model(x, *a):
n_ = np.arange(len(a), float)[None, :]
k_ = np.pi * (2*n_ - 1)/2
a_ = np.array(a)[None, :]
x_ = np.asarray(x)[:, None]
return 1 + (a_ * np.cos(k_*x_)).sum(axis=-1)
xdata = np.linspace(-1, 1, 200)
ydata = sympy.lambdify(x, u_exact, "numpy")(xdata)
curve_fit(model, xdata, ydata, p0=[1,1])
import numpy as np
from scipy.optimize import curve_fit
from uncertainties import correlated_values, unumpy as unp
N = 2
def model(x, *a):
n_ = np.arange(len(a), dtype=float)[None, :]
k_ = np.pi * (2*n_ + 1)/2
a_ = np.array(a)[None, :]
x_ = np.asarray(x)[:, None]
return 1 + (a_ * np.cos(k_*x_)).sum(axis=-1)
xdata = np.linspace(-1, 1, 200)
ydata = sympy.lambdify(x, u_exact, "numpy")(xdata)
p0, pcov = curve_fit(model, xdata, ydata, p0=np.ones(N))
p = correlated_values(p0, pcov)
print(p)
y_ = model(xdata, *p)
y, dy = unp.nominal_values(y_), unp.std_devs(y_)
f = 3
plt.plot(xdata*f, ydata, '-')
plt.plot(xdata*f, y, ':')
plt.fill_between(xdata*f, y-dy, y+dy, color='C1', alpha=0.5)
plt.grid('on')
plt.xticks([0, 1, 2, 3, 4, 5])
import numpy as np
from scipy.optimize import curve_fit
def model(x, *a):
n_ = np.arange(len(a), dtype=float)[None, :]
k_ = np.pi * (2*n_ + 1)/2
a_ = np.array(a)[None, :]
x_ = np.asarray(x)[:, None]
return 1 + (a_ * np.cos(k_*x_)).sum(axis=-1)
exact = sympy.lambdify(x, u_exact, "numpy")
fig, ax = plt.subplots(figsize=(4,3), constrained_layout=True)
xdata = np.linspace(-1, 1, 200)
ydata = exact(xdata)
ax.plot(xdata, ydata, '-k')
for N in range(2, 5):
p, pcov = curve_fit(model, xdata, ydata, p0=np.ones(N))
ax.plot(xdata, model(xdata, *p), '--', label=f"{N=}")
fig.legend(loc='upper center');
plt.plot(xdata, model(xdata, *get_a(np.array([0, 0.999])).ravel()))
plt.plot(xdata, ydata, '-k')
import sympy
def get_spectral_approx(N, xi=None, show=False):
an = sympy.symbols(f'a:{N}', real=True)
kn = tuple(sympy.pi * (2*n + 1)/2 for n in sympy.S(np.arange(N)))
if xi is None:
xi = tuple((2*n+1)/(2*N+1) for n in sympy.S(np.arange(N)))
xi = tuple((2*n+1)/(2*N+1) for n in sympy.S(np.arange(N)))
#print(xi)
#xi = tuple(sympy.cos(sympy.pi * (N-n)/2/N) for n in sympy.S(np.arange(N)))
x = sympy.symbols('x')
u = 1 + sum(a_ * sympy.cos(k_ * x) for a_, k_ in zip(an, kn))
Eqs = [(u.diff(x,x) - (x**6 + 3*x**2)*u).subs([(x, x_)]).evalf() for x_ in xi]
sol = sympy.solve(Eqs, an)
u_x = u.subs(sol)
if show:
display(sol)
display(u.evalf())
u_ = sympy.lambdify(x, u_x, "numpy")
return u_
fig, ax = plt.subplots()
x = np.linspace(-1, 1)
u_exact = np.exp((x**4-1)/4)
ax.plot(x, u_exact, label="Exact: $u(x) = e^{(x^4-1)/4}$")
for N in [2, 3, 4, 5, 6]:
u = get_spectral_approx(N=N)
ax.plot(x, u(x), label=f"Spectral: ${N=}$")
u = get_spectral_approx(N=2, xi=[0.1,0.2])
ax.plot(x, u(x), label=f"Spectral: ${N=}$")
fig.legend()
from scipy.optimize import minimize
res = minimize(get_err, [0.3, 0.75])
res = minimize(get_err2, [0.3, 0.75])
u = get_spectral_approx(N=2, xi=res.x)
X = np.linspace(-1, 1, 200)
plt.plot(X, u(X))
plt.plot(X, exact(X))
from tqdm.auto import tqdm
from itertools import product
def get_a(x, exact=exact):
n_ = np.arange(len(x), dtype=float)[None, :]
k_ = np.pi * (2*n_ + 1)/2
x_ = np.asarray(x)[:, None]
A = np.cos(k_*x_)
a = np.linalg.solve(A, exact(x) - 1)
return a
def get_err(x, weight=False):
x = np.asarray(x)
a = get_a(x)
X = np.linspace(-1, 1, 201)[1:-1]
w = 1/np.sqrt(1-X**2) if weight else 1
return (w*(model(X, *a) - exact(X))**2).mean()
def get_err2(x, weight=False):
u = get_spectral_approx(N=2, xi=x)
X = np.linspace(-1, 1, 201)[1:-1]
w = 1/np.sqrt(1-X**2) if weight else 1
return (w*(u(X) - exact(X))**2).mean()
fig, axs = plt.subplots(1, 2, figsize=(9, 3))
x0, x1 = np.meshgrid(np.linspace(0.001, 0.99), np.linspace(0.002, 0.999), indexing='ij')
err = [get_err([x0, x1]) for (x0, x1) in tqdm(list(zip(x0.ravel(), x1.ravel())))]
err = np.reshape(err, x0.shape)
mesh = axs[0].pcolormesh(x0, x1, 1/err.T)
fig.colorbar(mesh, ax=axs[0])
x0, x1 = np.meshgrid(np.linspace(0.001, 0.99, 30), np.linspace(0.002, 0.999, 30), indexing='ij')
err = [get_err2([x0, x1]) for (x0, x1) in tqdm(list(zip(x0.ravel(), x1.ravel())))]
err = np.reshape(err, x0.shape)
mesh = axs[1].pcolormesh(x0, x1, 1/err.T)
fig.colorbar(mesh, ax=axs[1])
def model(x, *a):
n_ = np.arange(1, len(a)+1, dtype=float)[None, :]
a_ = np.array(a)[None, :]
x_ = np.asarray(x)[:, None]
return 1 + (a_ * (np.cos(2*n_*np.arccos(x_))-1)).sum(axis=-1)
def get_xi(N):
"""Return the N collocation points."""
i = np.arange(N)
xi = np.cos((i + 0.5) * np.pi / 2 / N)
return xi
def get_a(N, exact=exact):
x_ = get_xi(N)[:, None]
n_ = np.arange(1, N+1, dtype=float)[None, :]
k_ = 2*n_
A = np.cos(k_*np.arccos(x_)) - 1
a = np.linalg.solve(A, exact(x_.ravel()) - 1)
return a
xi = get_xi(2)
a = get_a(2)
plt.plot(xdata, model(xdata, *a))
plt.plot(xdata, exact(xdata))
list(map(plt.axvline, xi))
x_ = get_xi(N)[:, None]
n_ = np.arange(1, N+1, dtype=float)[None, :]
k_ = 2*n_
A = np.cos(k_*np.arccos(x_)) - 1
a = np.linalg.solve(A, exact(x_.ravel()) - 1)
A @ a, exact(x_.ravel()) - 1
def model(x, *a):
n_ = np.arange(1, len(a)+1, dtype=float)[None, :]
a_ = np.array(a)[None, :]
x_ = np.asarray(x)[:, None]
return 1 + (a_ * (np.cos(2*n_*np.arccos(x_))-1)).sum(axis=-1)
model(x_, *a)
model(x_, *a)
n_ = np.arange(1, len(a)+1, dtype=float)[None, :]
a_ = np.array(a)[None, :]
(a_ * (np.cos(k_*np.arccos(x_))-1)).sum(axis=-1), exact(x_) - 1
A @ np.ravel(a_), (a_ * (np.cos(k_*np.arccos(x_))-1)).sum(axis=-1), exact(x_) - 1
def phi(x, n, d=0):
"""Return the dth derivative of phi_n(x)"""
th = np.arccos(x)
cn = np.cos(2*n*th)
if d == 0:
return cn - 1
sn = np.sin(2*n*th)
s = np.sin(th)
if d == 1:
return (2*n) * sn / s
c = np.cos(th)
if d == 2:
return -((2*n)**2 * cn * s - (2*n) * sn * c) / s**3
x0 = 0.3
h = 1e-8
n = 3
assert np.allclose((phi(x0+h, n) - phi(x0-h, n))/2/h, phi(x0, n, d=1))
assert np.allclose((phi(x0+h, n, d=1) - phi(x0-h, n, d=1))/2/h, phi(x0, n, d=2))
def get_a(N):
i = np.arange(N)
xi = np.cos((i + 0.5) * np.pi / 2 / N)
phis = np.array([phi(xi, n=n) for n in 1+np.arange(N)])
d2phis = np.array([phi(xi, n=n, d=2) for n in 1+np.arange(N)])
rhs = xi**6 + 3*xi**2
A = np.transpose(d2phis - rhs[None, :]*phis)
print(A.shape)
a = np.linalg.solve(A, rhs)
return a
plt.plot(xdata, ydata, ':')
plt.plot(xdata, model(xdata, *get_a(6)))
Stiff Example#
https://www.mathworks.com/company/technical-articles/stiff-differential-equations.html
%pylab inline
from ipywidgets import interact
from scipy.special import lambertw
from scipy.integrate import solve_ivp
def F(t, y):
return y**2*(1-y)
def exact(t, delta):
a = 1/delta - 1
return 1/(lambertw(a*np.exp(a-t)) + 1)
@interact(method=['RK23', 'RK45', 'DOP853', 'Radau', 'BDF', 'LSODA'],
log_delta=(-6, -1))
def go(log_delta=-2.0, method="RK45"):
delta = 10**(log_delta)
res = %time solve_ivp(F, y0=[delta], t_span=[0, 2/delta], atol=1e-12, rtol=1.e-12, method=method)
plt.plot(res.t*delta, exact(res.t, delta=delta))
plt.plot(res.t*delta, res.y.T)
plt.title(f"{res.nfev=}")
display(plt.gcf())
plt.close('all')
solve_ivp??