--- execution: timeout: 60 jupytext: cell_metadata_filter: -all formats: md:myst,ipynb main_language: python text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.16.7 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- ```{code-cell} :tags: [hide-cell] import mmf_setup; mmf_setup.nbinit() import os, sys from pathlib import Path FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/' os.makedirs(FIG_DIR, exist_ok=True) import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL) %matplotlib import numpy as np, matplotlib.pyplot as plt plt.rcParams['figure.constrained_layout.use'] = True try: from myst_nb import glue except: glue = None from matplotlib.animation import FuncAnimation from gpe.contexts import FPS ``` (sec:Viscosity5)= # Viscosity Notes 5 These are some supporting notes for the document {ref}`sec:EffectiveViscosity`. ## Microscopic Navier-Stokes :::{margin} As we shall see, the springs are not typical Hookean springs. This should be anticipated by the condition that the particles should not pass through each other. Thus $V(\delta q) \rightarrow \infty$ as $\delta q \rightarrow 0$. ::: To better understand the physics of the Navier-Stokes equations, consider representing the gas with a ball-and-"spring" model with $N+1$ positions $q_j$ and velocities $\dot{q}_j$ and inter-particle potential $V(\delta q)$. In the absence of dissipation we have the Lagrangian \begin{gather*} L[\vect{q}, \dot{\vect{q}}] = \sum_{i=0}^{N} \frac{m}{2}\dot{q}_i^2 - \sum_{i=0}^{N-1} V(q_{i+1} - q_{j}). \end{gather*} The corresponding Euler-Lagrange equations are \begin{gather*} m\ddot{q}_{i} = V'(q_{i+1} - q_{i}) - V'(q_{i} - q_{i-1}) = V'(\delta q_{i}) - V'(\delta q_{i-1}). \end{gather*} To interpret this as a fluid, we introduce the mass density. \begin{gather*} \rho_i = \frac{m}{q_{i+1} - q_{i}} = \frac{m}{\delta q_{i}}. \end{gather*} In the continuum limit the density (and velocity $u = \dot{q}_i$) must remain fixed so $m \sim \delta q_i \rightarrow 0$ in the same way. We can thus write \begin{gather*} \delta q_{i} = \frac{m}{\rho_i}, \qquad \ddot{q}_{i} = \frac{V'(\delta q_{i}) - V'(\delta q_{i-1})}{\rho_i \delta q_i}. \end{gather*} To take the continuum limit, we use the following notation \begin{gather*} f_{i} - f_{i-1} \equiv \delta f \rightarrow \diff{f}{i} = \diff{f}{q}\diff{q}{i} \equiv \delta q \diff{f}{q}. \end{gather*} Thus, the primary tool for effecting this limit are \begin{gather*} \delta \rightarrow \diff{}{i} = \delta q \diff{}{q} = \frac{1}{n}\diff{}{q}, \qquad \delta q \rightarrow \frac{1}{n}\diff{q}{q} = \frac{1}{n}. \end{gather*} Thus: \begin{gather*} \ddot{q}_{i} = \frac{1}{\rho_i \delta q_i}\delta V'(\delta q_i) \rightarrow \frac{1}{m} \frac{1}{n}\underbrace{\diff{V'(1/n)}{q}} _{V''(1/n)\frac{n_{,x}}{n^2}} = V''(1/n)\frac{n_{,x}}{mn^3} = -\frac{n_{,x}}{m}\mathcal{E}''(n), \end{gather*} Thus, in the continuum limit, we must have \begin{gather*} V''(1/n) \sim n^3 \mathcal{E}''(n). \end{gather*} E.g., for the GPE, we have $\mathcal{E}(n) = gn^2/2$, so $\mathcal{E}''(n) = g$. Here we set $r=\delta q = 1/n$ for notational convenience: :::{margin} I have not thought carefully about the meaning of the linear piece, but the overall constant can certainly be ignored. ::: \begin{gather*} V''(r) = \frac{g}{r^3},\qquad V'(r) = a - \frac{g}{2r^2},\qquad V(r) = \frac{g}{2r} + ar + c. \end{gather*} This means that \begin{gather*} \ddot{q} = \frac{-m\rho_{,x}}{\rho^3} (g n^3) = \frac{-m\rho_{,x}}{\rho^3} \frac{g \rho^3}{m^3} = \frac{-g}{m^2}\rho_{,x}. \end{gather*} Everything here must be finite in the continuum limit, so we must fixed $g/m^2$ as our interaction strength in order to compare with GPE simulations. :::{important} If you are used to working with the GPE, you might naturally think of $n$ as being finite, but this is because we usually work with a fixed particle mass $m$ that is already very small. When taking he continuum limit here, you must regard the mass density $\rho = mn$ as fixed to get the correct scaling. ::: :::{admonition} Strategy for taking the Continuum Limit :class: dropdown My strategy for taking the continuum limit is to express everything in terms of what will become physical quantities: \begin{align*} \rho(x_i) &\sim \rho_i, \\ \rho_{,x}(x_i) &\sim \frac{\rho_{i+1} - \rho_{i}}{q_{i+1} - q_{i}} = \frac{\rho_{i+1} - \rho_{i}}{\delta q_{i+1}}. \end{align*} It is helpful to think of our goal to match the Euler equation \begin{gather*} \ddot{q} = \frac{-\bigl(\mathcal{E}'(n)\bigr)_{,x}}{m} = \frac{-\bigl(\mathcal{E}'(\rho/m)\bigr)_{,x}}{m} = \frac{-\mathcal{E}''(\rho/m)\rho_{,x}}{m^2} = \frac{-\mathcal{E}''(n)n_{,x}}{m}. \end{gather*} Thus, we are trying to manipulate our equations into a form with a gradient: \begin{gather*} \nabla \sim \frac{1}{\delta q_i}. \end{gather*} We must also connect $\delta q_{i}$ and $\delta q_{i-1}$: \begin{gather*} \delta q_{i} = \frac{m}{\rho_{i}}, \qquad \rho_{i} \sim \rho_{i-1} + \delta q_{i}\rho_{,x}(x_{i-1}),\\ \delta q_{i-1} = \frac{m}{\rho_{i-1}} = \frac{m}{\rho_{i} - \delta q_{i}\rho_{,x}(x_{i-1})} = \frac{m}{\rho_{i}}\left( 1 + \delta q_{i}\frac{\rho_{,x}(x_{i-1})}{\rho_i} + O(\delta q)^2\right),\\ = \frac{m}{\rho_{i}}\left( 1 + \delta q_{i}\frac{\rho_{,x}}{\rho} + O(\delta q)^2\right). \end{gather*} In the last equation, we note that, since $\rho_{,x}$ will be finite in the continuum limit, $\rho_{,x}(x_{i-1}) = \rho_{,x}(x_{i}) + O(\delta q)$, as long as $\rho_{,xx}$ does not diverge. The key to these manipulations is to keep consistent orders of expansion. We use this to compute \begin{gather*} V'(\delta q_{i-1}) = V'\left( \delta q_i + \delta q_{i}^2\frac{\rho_{,x}}{\rho} + O(\delta q)^3\right),\\ = V'(\delta q_i) + \delta q_{i}^2\frac{\rho_{,x}}{\rho} V''(\delta q_i) + O(\delta q)^3, \end{gather*} from which we can compute the desired result above. ::: :::{admonition} Numerical Check :class: dropdown Here we numerically check this result. We start with $N$ particles separated such that the initial density is $n_0(x) \propto R^2-x^2$ where $R$ is the radius of the initial cloud (Thomas-Fermi radius): \begin{gather*} n(q) = \frac{R^2 - q^2}{\alpha m}, \\ \delta q = q'(i) = \frac{\alpha m}{R^2-q^2}, \\ \int_{-R}^{q} (R^2-q^2)\d{q} = \int_{0}^{i} \alpha m\d{i}, \\ \frac{3q_iR^2 - q_i^3 + 2R^3}{3} = \alpha m i\\ \alpha m = \frac{4R^3}{3N}, \qquad \frac{3\frac{q_i}{R} - \frac{q_i^3}{R^3} + 2}{4} = \frac{i}{N}. \end{gather*} This requires inverting a cubic, which can be done easily using Newton's method. One caveat: with these normalizations, we will end up with a fixed mass density $\rho$, and the number density $n = \rho/m$ will integrate to the number of particles $N$. This makes it a little hard to compare things since our physical equation of state $\mathcal{E}(n) = gn^2/2$ is expressed in terms of the number density. To allow us to compare states with different numbers of particles (i.e., to test the continuum limit), we consider the classical equations: \begin{gather*} \ddot{q} = -\frac{g}{m^2}\rho_{,x} + \cdots. \end{gather*} Thus, we should fix the ratio $g/m^2$, not $g$ like we are accustomed to doing with the GPE where $m$ is held fixed in experiments. ::: :::{admonition} Globally Convergent Newton's Method :class: dropdown To get our initial density profile, we need to solve \begin{gather*} 3q - q^3 + 2 = (q-2)(q+1)^2 = 4x, \quad q \in [-1, 1], \qquad x \in [0, 1]. \end{gather*} This can be done with Newton's method, but the current form has zero derivatives at the endpoints, which will cause problems. Using the symmetry of the problem, we restrict ourselves to $q<0$ and use something linear at the left endpoint: \begin{gather*} f(q) = (q+1) - \sqrt{\frac{4x}{2-q}},\\ f'(q) = 1 - \sqrt{\frac{x}{(2-q)^3}},\\ q_0 = (2x)^{0.6}-1. \end{gather*} This converges to machine precision in 3 steps for all $x$. ::: ```{code-cell} def get_q_R(x, n=3): """Return root of 3*q-q**3+2 = 4*x for 0 <= x < 0.5.""" def f(q, x, d=0): if d == 0: return q + 1 - np.sqrt(4*x/(2-q)) elif d == 1: return 1 - np.sqrt(x/(2-q)**3) def iter(q, x): return q - f(q, x) / f(q, x, d=1) x_ = np.where(x < 0.5, x, 1-x) q = (2*x_)**0.6 - 1 for n in range(n): q = iter(q, x_) return np.where(x < 0.5, q, -q) x = np.linspace(0, 1) q = get_q_R(x) assert abs(3*q - q**3 + 2 - 4*x).max() < 2.1*np.finfo(float).eps ``` As an additional check, consider the two particles at $q_{i\pm 1}$ neighboring $q_{i}$. The net force on the particle at $q_i$ is \begin{align*} F_{i} &= -V'(q_{i+1}-q_{i}) + V'(q_{i} - q_{i-1}) = -\mathcal{E}''(n)n_{,x},\\ F_{i} &= - V''(\delta q)\diff{\delta q}{i} = - V''(\delta q)\diff{\delta q}{q}\underbrace{\diff{q}{i}}_{\delta q}. \end{align*} Now replacing $\delta q = 1/n$ we have \begin{gather*} F = - V''(n^{-1})(n^{-1})_{,x}\frac{1}{n} = V''(n^{-1})\frac{n_{,x}}{n^3} = -\mathcal{E}''(n)n_{,x}\\ V''(n^{-1}) = -n^3\mathcal{E}''(n) \end{gather*} ## Viscosity We now model viscosity through a dissipation term in the "springs". The idea is that there is a resisting force $F_{\beta}(\delta q, \delta \dot{q})$ that depends on the extension of the spring $\delta q$ and the velocity $\delta \dot{q}$: :::{margin} Here we use \begin{gather*} \dot{q}_{,x} = n \delta \dot{q} = \frac{\delta \dot{q}}{\delta q}. \end{gather*} ::: \begin{gather*} F_i = \cdots + \bigl(f(\delta q_{i}, \delta \dot{q}_{i}) - f(\delta q_{i-1}, \delta \dot{q}_{i-1}) = \cdots + \diff{}{i}f(\delta q, \delta \dot{q}),\\ \diff{}{i}f(\delta q, \delta \dot{q}) = \left( f_{,\delta q}\diff{\delta q}{q} + f_{,\delta \dot{q}}\diff{\delta \dot{q}}{q}\right)\diff{q}{i} = \Bigl( f_{,\delta q}(n^{-1})_{,x} + f_{,\delta \dot{q}}(\dot{q}_{,x}n^{-1})_{,x} \Bigr)n^{-1},\\ = -f_{,\delta q}\frac{n_{,x}}{n^3} + f_{,\delta \dot{q}}\left( \frac{\dot{q}_{,xx}}{n^2} - \frac{\dot{q}_{,x}n_{,x}}{n^3} \right) = f_{,\delta \dot{q}}\frac{\dot{q}_{,xx}}{n^2} -(f_{,\delta q} + n\delta \dot{q} f_{,\delta \dot{q}})\frac{n_{,x}}{n^3} \end{gather*} :::{margin} \begin{align*} \frac{n}{m\nu}\diff{f}{i} &= \frac{\overbrace{n\delta q}^{1}}{m\nu}\diff{f}{q} \\ \frac{1}{m\nu}\diff{f}{q} &= \diff{(n^2\delta \dot{q})}{q},\\ \frac{f}{m\nu} &= n^2\delta \dot{q}. \end{align*} ::: E.g., we can generate a Lamé viscosity: \begin{gather*} \diff{f}{i} = m\nu \frac{(n\dot{q}_{,x})_{,x}}{n}, \qquad f(\delta q, \delta \dot{q}) = m \nu \frac{\delta \dot{q}}{(\delta q)^2}. \end{gather*} :::{admonition} Details :class: dropdown I was not sure how to first express this, so I went back to the idea \begin{gather*} \delta q_{i+1} - \delta q_{i} \sim \diff{\delta q}{i} = \diff{\delta q}{q}\diff{q}{i} = \diff{n^{-1} q}{q}\delta q = -\frac{n_{,x}}{n^2}\delta q = -\frac{\rho_{,x}}{\rho}\frac{\delta q}{n} = -\frac{\rho_{,x}}{\rho}\delta q^2. \end{gather*} In other words: \begin{gather*} \delta \sim \diff{}{i} = \diff{q}{i}\diff{}{q} = \delta q\diff{}{q}, \qquad \delta \dot{q} = \delta q \diff{\dot{q}}{q} = \frac{\dot{q}_{,x}}{n}, \end{gather*} etc. We consider the following two cases: \begin{align*} \frac{\beta f}{m\nu} &\approx \dot{q}_{,xx},\\ \frac{\beta f}{m\nu} &\approx \frac{(n \dot{q}_{,x})_{,x}}{n} = \dot{q}_{,xx} + \frac{n_{,x} \dot{q}_{,x}}{n}. \end{align*} The first is the simple viscosity, while the second is a Lamé viscosity with power $\alpha = 1$. In both cases: \begin{gather*} \frac{\beta}{m\nu} f_{,\delta \dot{q}} = n^2 = \frac{1}{(\delta q)^2}, \qquad \frac{\beta f}{m\nu} = \frac{\delta \dot{q}}{(\delta q)^2} + g(\delta q). \end{gather*} In the first case, we must have \begin{gather*} f_{,\delta q} + n\delta \dot{q} f_{,\delta\dot{q}} = -2\frac{\delta \dot{q}}{(\delta q)^3} + g'(\delta q) + \frac{\delta \dot{q}}{(\delta q)^3} = g'(\delta q) - \frac{\delta \dot{q}}{(\delta q)^3} = 0,\\ g'(\delta q) = \frac{\delta \dot{q}}{(\delta q)^3}. \end{gather*} This has no solutions since $g(\delta q)$ should be independent of $\delta \dot{q}$. In the second case \begin{gather*} \frac{\beta}{m\nu}(f_{,\delta q} + n\delta \dot{q} f_{,\delta\dot{q}}) = -n^2\dot{q}_{,x} = -n^3\delta \dot{q},\\ \frac{\beta}{m\nu}\left( g'(\delta q) - \frac{\delta \dot{q}}{(\delta q)^3} \right) = -\frac{\delta \dot{q}}{(\delta q)^3}, \beta = m \nu, \qquad g(\delta q) = 0. \end{gather*} This is as we stated above. ::: ```{code-cell} from scipy.integrate import solve_ivp from scipy.interpolate import InterpolatedUnivariateSpline class Particles: Ns = 32 # Number of springs g_m2 = 1.0 # Coupling constant over mass^2 R = 1.0 # Initial radius beta = 0.0 # Dissipation nu = 0.0 # Lame viscosity 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): x0 = np.arange(0, self.Ns+1)/self.Ns q0_R = get_q_R(x0) self.q0 = self.R * q0_R self.m = 4*self.R**3/3/self.Ns self.dq0 = 0*self.q0 self.x0 = np.linspace(-self.R, self.R, self.Ns+1) def pack(self, q, dq): return np.ravel([q, dq]) def unpack(self, y): return np.reshape(y, (2, y.shape[0]//2) + y.shape[1:]) def V(self, dq, d=0): g = self.g_m2 * self.m**2 if d == 0: return g/2/dq elif d == 1: return -g/2/dq**2 def compute_dy_dt(self, t, y): q, dq_dt = self.unpack(y) dq = np.concatenate([[0], np.diff(q), [0]]) dV = np.concatenate([[0], self.V(np.diff(q), d=1), [0]]) ddq_dt = np.concatenate([[0], np.diff(dq_dt), [0]]) f = self.m * self.nu * ddq_dt/dq**2 f[0] = f[-1] = 0 #ddq_dt2 = (dV[1:] - dV[:-1]) / self.m + self.beta * (dv[1:] - dv[:-1]) ddq_dt2 = (dV[1:] - dV[:-1] + (f[1:] - f[:-1])) / self.m return self.pack(dq_dt, ddq_dt2) p1 = Particles(Ns=32) p2 = Particles(Ns=256) ps = [p1, p2] y0s = [p.pack(p.q0, p.dq0) for p in ps] T = 2.0 Nt = 100 t_eval = np.linspace(0, T, Nt) ress = [solve_ivp(p.compute_dy_dt, y0=y0, t_span=(0,T), t_eval=t_eval, method="LSODA") for p, y0 in zip(ps, y0s)] for res, fmt in zip(ress, ['--C0', ':C1']): q, dq = p2.unpack(res.y) plt.plot(res.t, q.T, fmt); ``` ```{code-cell} Nt = len(t_eval) rhos = [] xs = [] xlim = [0,0] rho_max = 0 for p, res in zip(ps, ress): q, dq = p.unpack(res.y) rho = p.m/np.diff(q, axis=0) x = (q[1:] + q[:-1])/2 xlim = [min(xlim[0], x.min()), max(xlim[1], x.max())] rho_max = max(rho_max, rho.max()) rhos.append(rho) xs.append(x) args = dict(xlim=xlim, ylim=(0, rho_max)) fig, ax = plt.subplots() for i in FPS(Nt, fig=fig, embed=True, fps=10): ax.cla() for j, p in enumerate(ps): ax.plot(xs[j][:, i], rhos[j][:, i]) t = t_eval[i] ax.set(title=f"t={res.t[i]:.4f}", **args) ``` Here is the Lagrangian hydrodynamic description for comparison: ```{code-cell} import scipy.fft sp = scipy class Flow: Nx = 128 A = 0.5 # g/2/m^2 nu = 0.04 # Viscosity coefficients 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): self._th0 = (2*np.arange(self.Nx) + 1) * np.pi / 2 / self.Nx self._k = np.arange(self.Nx) self.x0 = np.cos(self._th0) # Initial density profile self.n0 = 0*self.x0 + 1.0 self.n0_x0 = 0*self.x0 self.n0 = (1-self.x0**2) self.n0_x0 = -2*self.x0 def get_y0(self): X0 = self.x0 Xt0 = 0*X0 return self.pack(X0, Xt0) def diff(self, f, d=1): ft = sp.fft.dct(f, type=2, norm='forward') th, k = self._th0, self._k s = np.sin(th) # Shift the coefficients and pad with zero. df_dth = sp.fft.dst(np.concatenate([-(ft*k)[1:], [0]]), type=3) if d == 1: df_dx = -df_dth / s return df_dx else: d2f_dth2 = sp.fft.idct(-ft*k**2, type=2, norm='forward') d2f_dx2 = (d2f_dth2 - df_dth / np.tan(th)) / s**2 if d == 2: return d2f_dx2 return self.diff(df_dx, d=d-2) def pack(self, X, X_t): return np.ravel([X, X_t]) def unpack(self, y): return np.reshape(y, (2, self.Nx)) def unpack_xnu(self, y): X, X_t = self.unpack(y) X_x0 = self.diff(X) n = self.n0 / X_x0 u = X_t x = X return x, n, u def get_a_ext(self): return 0 def compute_dy_dt(self, t, y): X, X_t = self.unpack(y) X_x0 = self.diff(X) X_x0x0 = self.diff(X, d=2) n = self.n0 / X_x0 n_x = self.n0_x0 / X_x0**2 - self.n0 * X_x0x0 / X_x0**3 X_tx0 = self.diff(X_t) u_x = X_tx0/X_x0 u_xx = (self.diff(X_t, d=2)) / X_x0**2 - self.diff(X_t) * X_x0x0 / X_x0**3 # Zero-laplacian boundary conditions u_xx[0] = u_xx[-1] = 0 a_ext = self.get_a_ext() X_tt = a_ext - 2*self.A * n_x + self.nu * (n_x * u_x / n + u_xx) return self.pack(X_t, X_tt) def solve(self, T=10.0, Nt=100): y0 = self.get_y0() t_eval = np.linspace(0, T, Nt) res = solve_ivp(self.compute_dy_dt, y0=y0, t_span=(0, T), t_eval=t_eval, method="BDF") if not res.success: print(res.message) self.res = res return res ``` ```{code-cell} %%time from scipy.integrate import solve_ivp p1 = Particles(Ns=32) p2 = Particles(Ns=256) ps = [p1, p2] y0s = [p.pack(p.q0, p.dq0) for p in ps] T = 2.0 Nt = 100 t_eval = np.linspace(0, T, Nt) ress = [solve_ivp(p.compute_dy_dt, y0=y0, t_span=(0,T), t_eval=t_eval, method="LSODA") for p, y0 in zip(ps, y0s)] rhos = [] xs = [] xlim = [0,0] rho_max = 0 for p, res in zip(ps, ress): q, dq = p.unpack(res.y) rho = p.m/np.diff(q, axis=0) x = (q[1:] + q[:-1])/2 xlim = [min(xlim[0], x.min()), max(xlim[1], x.max())] rho_max = max(rho_max, rho.max()) rhos.append(rho) xs.append(x) f = Flow(A=p1.g_m2/2, nu=0.0) res = f.solve(T=T, Nt=Nt) x, n, u = np.einsum('tix->itx', [f.unpack_xnu(y) for y in res.y.T]) args = dict(xlim=(min(xlim[0], x.min()), max(xlim[1], x.max())), ylim=(0, max(rho_max, n.max()))) Nt = len(res.t) fig, ax = plt.subplots() for i in FPS(Nt, fig=fig, embed=True, fps=10): t = res.t[i] ax.cla() ax.plot(x[i], n[i], '-') for j, p in enumerate(ps): ax.plot(xs[j][:, i], rhos[j][:, i], '--') t = t_eval[i] ax.set(title=f"t={res.t[i]:.4f}", **args) ``` ```{code-cell} %%time from scipy.integrate import solve_ivp nu = 1.0 #p1 = Particles(Ns=32, beta=0.5*32**2) #p2 = Particles(Ns=256, beta=0.5*256**2) p1 = Particles(Ns=32, nu=nu) p2 = Particles(Ns=256, nu=nu) ps = [p1, p2] y0s = [p.pack(p.q0, p.dq0) for p in ps] T = 2.0 Nt = 100 t_eval = np.linspace(0, T, Nt) ress = [solve_ivp(p.compute_dy_dt, y0=y0, t_span=(0,T), t_eval=t_eval, method="LSODA") for p, y0 in zip(ps, y0s)] rhos = [] xs = [] xlim = [0,0] rho_max = 0 for p, res in zip(ps, ress): q, dq = p.unpack(res.y) rho = p.m/np.diff(q, axis=0) x = (q[1:] + q[:-1])/2 xlim = [min(xlim[0], x.min()), max(xlim[1], x.max())] rho_max = max(rho_max, rho.max()) rhos.append(rho) xs.append(x) f = Flow(A=p1.g_m2/2, nu=0.0) res = f.solve(T=T, Nt=Nt) x, n, u = np.einsum('tix->itx', [f.unpack_xnu(y) for y in res.y.T]) f = Flow(A=p1.g_m2/2, nu=nu) res = f.solve(T=T, Nt=Nt) x, n, u = np.einsum('tix->itx', [f.unpack_xnu(y) for y in res.y.T]) args = dict(xlim=(min(xlim[0], x.min()), max(xlim[1], x.max())), ylim=(0, max(rho_max, n.max()))) Nt = len(res.t) fig, ax = plt.subplots() for i in FPS(Nt, fig=fig, embed=True, fps=10): t = res.t[i] ax.cla() ax.plot(x[i], n[i], '-') for j, p in enumerate(ps): ax.plot(xs[j][:, i], rhos[j][:, i], '--') t = t_eval[i] ax.set(title=f"t={res.t[i]:.4f}", **args) ``` [Rankine-Hugoniot conditions]: [Shallow water equations]: [conservative form of the shallow water equations]: [weir]: https://en.wikipedia.org/wiki/Weir [Shallow water equations]: [conservative form of the shallow water equations]: [Lagrangian fluid dynamics]: [van der Waals force]: [mole]: [smoothed-particle hydrodynamics]: