Flow

Here we consider flow past a weak barrier. The force on the barrier is

\[ F = -E'(x_0) = \int V'(x-x_0)n(x)\;d{x} \]

# %#matplotlib inline
import numpy as np, matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec

import mmfutils.performance.fft, mmfutils.plot
from pytimeode.evolvers import EvolverSplit

from gpe.bec import StateTwist_x, Units
from gpe import utils

mmfutils.performance.fft.set_num_threads(1)


class State(StateTwist_x):
    def __init__(self, experiment, **kw):
        self.experiment = experiment
        super().__init__(**kw)

    def get_xyz(self, gpu=False):
        if gpu:
            return self._xyz_
        else:
            return self.xyz

    def _get_Vext_(self):
        """Return the external potential."""
        return self.experiment.get_Vext(state=self, gpu=True)

    def get_force(self):
        dV = self.experiment.get_Vext(state=self, d=1)
        n = self.get_density()
        return self.integrate(dV * n)

    def get_momentum(self):
        psi = self.get_psi()
        return self.hbar * self.braket(self.basis.get_gradient(psi)[0]).imag

    def plot(self, log=False, label=None, comoving=False):  # pragma: nocover
        n = self.get_density()

        def scale(n):
            if log:
                return np.log10(n)
            return n

        n = scale(n)
        if comoving:
            x = self.xyz[0]
            xlabel = r"$x$ [$\xi$]"
        else:
            x = self.x_lab
            xlabel = r"$x_{\rm lab}$ [$\xi$]"
            inds = np.argsort(x)
            x, n = x[inds], n[inds]

        x_ = x / self.experiment.healing_length
        plt.plot(x_, n, label=label)
        ax = plt.gca()
        ax.set(xlabel=xlabel)
        E = self.get_energy()
        N = self.get_N()
        t = np.ravel(self.t)[0]
        t_ = t / self.experiment.t_unit
        t_name = self.experiment.t_name
        plt.suptitle(f"$t={t_:.4f}{t_name}$, ${N=:.4f}$, ${E=:.4f}$")


class Experiment:
    Nxyz = (256,)
    Lxyz = (10.0,)
    healing_length_dx = 2.0
    hbar = m = 1
    n0 = 1

    V0_mu = 0.1
    V_sigma_healing_length = 3.0
    V_t_ = 0.0  # How long to take to turn barrier on in units of t_unit
    V_v_c = 0.0  # Speed of potential in units of c_cound
    V_integrable_mu = 0.0  # How much to break integrability

    State = State
    phi0 = 0.01  # Size of random phases

    v_x_V_v = 1.0  # Speed of frame in units of V_v
    t_cool = 1.0

    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):
        dx = self.Lxyz[0] / self.Nxyz[0]
        self.healing_length = self.healing_length_dx * dx
        self.mu = self.hbar**2 / 2 / self.m / self.healing_length**2

        self.t_unit = self.hbar / self.mu
        self.t_name = r"\hbar/\mu"

        self.g = self.mu / self.n0
        self.V_sigma = self.V_sigma_healing_length * self.healing_length

        self.V0 = self.V0_mu * self.mu
        self.V0_t_ = utils.get_smooth_transition(
            [0, self.V0], durations=[0], transitions=[self.V_t_]
        )

        self.c_sound = np.sqrt(self.mu / self.m)
        self.V_v = self.V_v_c * self.c_sound
        self.v_x = self.v_x_V_v * self.V_v
        self.state_args = dict(
            experiment=self,
            Nxyz=self.Nxyz,
            Lxyz=self.Lxyz,
            g=self.g,
            hbar=self.hbar,
            m=self.m,
            v_x=self.v_x,
        )
        self.rng = np.random.default_rng(seed=2)

    def get_Vext(self, state, gpu=False, d=0):
        """Return the external potential or it's derivative."""
        t = state.t
        t_ = t / self.t_unit

        x_lab = state.get_xyz(gpu=gpu)[0] + state.v_x * t
        x0 = self.V_v * t

        Lx = self.Lxyz[0]
        x0 = (x0 - x_lab.min()) % Lx + x_lab.min()
        ks = np.arange(1, 5) * np.pi / Lx

        # Background to break integrability
        V = (
            self.V_integrable_mu
            * self.mu
            * sum(np.sin(k * x_lab) for a, k in zip(ks, [1.0, -2.0, 3.0, 1.0]))
        )
        V += self.V0_t_(t_) * sum(
            np.exp(-(((x_lab - x0) / self.V_sigma) ** 2) / 2)
            for x0 in [x0 - Lx, x0, x0 + Lx]
        )
        if d == 0:
            return V
        elif d == 1:
            dV = self.V0 * sum(
                -(x_lab - x0)
                / self.V_sigma**2
                * np.exp(-(((x_lab - x0) / self.V_sigma) ** 2) / 2)
                for x0 in [x0 - Lx, x0, x0 + Lx]
            )
            return dV
        else:
            raise NotImplementedError(f"{d=} not supported.")

    def get_state(self, cool=False):
        state = self.State(**self.state_args)
        V = self.get_Vext(state)
        n = np.maximum((self.mu - V) / self.g, 0)
        phi = 2 * self.phi0 * (self.rng.random(size=n.shape) - 0.5)
        state.set_psi(np.sqrt(n) * np.exp(1j * phi))
        if cool:
            state.cooling_phase = 1j
            dt = 0.1 * state.t_scale
            steps = max(int(np.ceil(self.t_cool / dt)), 2)
            ev = EvolverSplit(state, dt=self.t_cool / steps, normalize=True)
            ev.evolve(steps)
            psi = ev.y.get_psi()

            state = self.State(**self.state_args)
            state.set_psi(psi)
        return state


e = Experiment()
s = e.get_state()
s.plot()

KeyboardInterrupt

from IPython.display import clear_output, display
from mmfutils.contexts import FPS
from pytimeode.evolvers import EvolverABM

e = Experiment(V_v_c=0.0, phi0=0.2, t_cool=1)
s0 = e.get_state(cool=False)
s0.plot()
plt.figure()
s1 = e.get_state(cool=True)
s1.plot()

assert np.allclose(s0.get_N(), s1.get_N())
s0.get_energy(), s1.get_energy()

f = 2
e = Experiment(V_v_c=0.2, V_sigma_healing_length=10, V_integrable_mu=0.01,
               V0_mu=0.1, phi0=0.1*np.pi, t_cool=0.1, Nxyz=(256*f,), Lxyz=(10.*f,))
s = e.get_state(cool=True)
s.cooling_phase = 1
ev = EvolverABM(s, dt=0.1*s.t_scale)
steps = 200
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
Fs = [s.get_force()]
ps = [s.get_momentum()]
ts = [s.t]
skip = 10
for frame in FPS(frames=10000, timeout=60*5):
    ev.evolve(steps)
    ts.append(ev.y.t)
    Fs.append(ev.y.get_force())
    ps.append(ev.y.get_momentum())
    if frame % skip == 0:
        ax = axs[0]
        ax.cla()
        plt.sca(ax)
        ev.y.plot()

        ax = axs[1]
        ax.cla()
        plt.sca(ax)
        Fmax = np.abs(Fs).max()
        pmax = np.abs(ps).max()
        plt.plot(ts, np.divide(Fs, Fmax))
        plt.plot(ts, np.divide(ps, pmax))
        ax.set(title=f"{Fmax=:0.2g} ({np.mean(Fs):.4g}), {pmax=:0.2g} ({np.mean(ps):.4g})")
        #plt.axhline(np.mean(Fs))
        ax.set(ylim=(-1,1))
        display(fig)
        clear_output(wait=True)

s = ev.get_y()
plt.plot(s.x, np.gradient(s.get_Vext(), s.x))
plt.plot(s.x, s.experiment.get_Vext(s, d=1))

Things to Check (for Chance)

1. No matter what the state, if \(V(x)\) is independent of time, then the total energy is conserved. This is consistent with the fact that, even if there is a force, there is no displacement, so there is no work done. Derive this.

2. A time-independent \(V(x)\) does not imply that total momentum is conserved.

3. A moving potential \(V(x-vt)\) can do both work (energy changes) and can change the momentum.

4. Thus, a time-independent potential \(V(x)\) in a constantly moving frame must change both energy and momentum. Explain why the derivation in 1) fails.

f = 2
kw = dict(V_v_c=0.2, V0_mu=0.1, Nxyz=(256*f,), Lxyz=(10.*f,), phi0=0.01)
e0 = Experiment(v_x_V_v=0, **kw)
e1 = Experiment(v_x_V_v=1, **kw)
s0 = e0.get_state(cool=False)
s1 = e1.get_state(cool=False)
s0.get_energy(), s1.get_energy()

ev0 = EvolverABM(s0, dt=0.1 * s0.t_scale)
ev1 = EvolverABM(s1, dt=0.1 * s0.t_scale)
[ev.evolve(10000) for ev in [ev0, ev1]]
print([ev.y.get_energy() for ev in [ev0, ev1]])
[ev.y.plot(comoving=False) for ev in [ev0, ev1]]

Review Paper

Here we explore some of the features described in the paper https://arxiv.org/abs/2306.05048: “Stationary transport above the critical velocity in a one-dimensional superflow past an obstacle.” Instead of solving for the stationary solutions (which can be done with a nice classical mechanical analog as described in the paper), we solve the dynamic problem similar to Paris-Mandoki et al., “Superfluid flow above the critical velocity” Scientific Reports, 7 (2015) 9070 where we slowly turn on the potential.

def go(v=1.5, U0=2, sigma=10.0, f=8, **kw):
    kw.update(V_v_c=v, V0_mu=U0, V_t_=50.0, phi0=0.0, V_sigma_healing_length=sigma)
    kw.update(Nxyz=(256 * f,), Lxyz=(10.0 * f,))
    e = Experiment(**kw)
    s = e.get_state(cool=False)
    s.cooling_phase = 1

    T = e.Lxyz[0] / e.c_sound
    ev = EvolverABM(s, dt=0.1 * s.t_scale)
    steps = 200
    # fig, axs = plt.subplots(1, 2, figsize=(10, 5))
    fig, ax = plt.subplots()
    axs = [ax]
    skip = 10
    frames = int(np.ceil(T / ev.dt / steps))
    for frame in FPS(frames=frames, timeout=60 * 5):
        ev.evolve(steps)
        if frame % skip == 0:
            ax = axs[0]
            ax.cla()
            plt.sca(ax)
            ev.y.plot(comoving=True)
            display(fig)
            clear_output(wait=True)
    return ev.get_y()

# Some points from Fig. 1.
# b)
go(U0=2.0, v=1.0)

# a) Here g*n_min = mu - U0, v_c/c_sound = sqrt(1-U0/mu)

#go(U0=0.5, v=0.4*np.sqrt(0.5))
s = go(U0=2, v=0.1)

# d)
s = go(U0=-0.4, v=2.0, V_integrable_mu=0.001, f=8)

f = 2
e = Experiment(V_v_c=0.9, v_x_V_v=1, V0_mu=0.1, V_t_=100.0, phi0=0.0, t_cool=0.1,
               
               Nxyz=(256*f,), Lxyz=(10.*f,))
s = e.get_state(cool=True)
s.cooling_phase = 1
ev = EvolverABM(s, dt=0.1*s.t_scale)
steps = 200
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
Fs = [s.get_force()]
ps = [s.get_momentum()]
ts = [s.t]
skip = 10
for frame in FPS(frames=10000, timeout=60*5):
    ev.evolve(steps)
    ts.append(ev.y.t)
    Fs.append(ev.y.get_force())
    ps.append(ev.y.get_momentum())
    if frame % skip == 0:
        ax = axs[0]
        ax.cla()
        plt.sca(ax)
        ev.y.plot(comoving=True)

        ax = axs[1]
        ax.cla()
        plt.sca(ax)
        Fmax = np.abs(Fs).max()
        pmax = np.abs(ps).max()
        plt.plot(ts, np.divide(Fs, Fmax))
        plt.plot(ts, np.divide(ps, pmax))
        ax.set(title=f"{Fmax=:0.2g} ({np.mean(Fs):.4g}), {pmax=:0.2g} ({np.mean(ps):.4g})")
        #plt.axhline(np.mean(Fs))
        ax.set(ylim=(-1,1))
        display(fig)
        clear_output(wait=True)

from ipywidgets import interact
def g(n):
    return n

def G(n):
    return n**2/2

@interact(v=(0, 3.0), U0=(-1.0, 1.0))
def draw_W(v=0.1, U0=0.0):
    n = np.linspace(0.5, 2, 500)[1:]
    W = v**2/2*(n + 1/n) + n*g(1) - G(n)
    W_1 = v**2/2*(1 + 1/1) + 1*g(1) - G(1)
    W0 = W - U0*n
    W0_1 = W_1 - U0*1
    plt.plot(np.sqrt(n), W-W_1)   
    plt.plot(np.sqrt(n), W0-W0_1)