--- execution: timeout: 60 jupytext: cell_metadata_filter: -all formats: md:myst,ipynb main_language: python text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.16.7 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- ```{code-cell} :tags: [hide-cell] import mmf_setup; mmf_setup.nbinit() import os, sys from pathlib import Path FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/' os.makedirs(FIG_DIR, exist_ok=True) import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL) %matplotlib import numpy as np, matplotlib.pyplot as plt try: from myst_nb import glue except: glue = None from matplotlib.animation import FuncAnimation from gpe.contexts import FPS ``` (sec:Viscosity4)= # Viscosity Notes 4 These are some supporting notes for the document {ref}`sec:EffectiveViscosity`. ## Saint-Venant Equations: Vacuum Interface In {ref}`sec:Viscosity3` we explored the Saint-Venant equations in the case where the density $n > 0$ was finite everywhere. The case where $n = 0$ -- sometimes referred to as "dry bed" problems -- poses a challenge. To see this, we demonstrate the dam breaking problem: ```{code-cell} :tags: [hide-input] import viscosity f = viscosity.Flow() 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, embed=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', cla=False, label=f"$n_R=0.01$: nfev={nfevs1[-1]}") fig_axs[1][0].legend(loc='center left') print(f"Wet bed: nR=0.1 : nfev = {nfevs0[0]}") print(f"Dry bed: nR=0.01: nfev = {nfevs1[0]}") ``` Notice that the initial step for the dry bed $n_R=0.01$ takes an order of magnitude more function evaluations (`nfev`) than the wet bed $n_R=0.1$. As we make the bed dryer, this gets even worse. One approach to dealing with this is to reformulate the problem with [Lagrangian fluid dynamics][]. ## [Lagrangian Fluid Dynamics][] There are two alternative formulations of fluid dynamics. The Eulerian approach we have adopted so far describes the density and velocity as functions of fixed spatial position. The Lagrangian approach instead follows the particle motion. Thus, consider the position $x = X(x_0, t)$ of the fluid element that started at $x_0$ at time $t=0$: $X(x_0, 0) = x_0$. The relationship between the two is expressed in terms of the velocity \begin{gather*} \overbrace{u\bigl(X(x_0, t), t\bigr)}^{u} = \overbrace{X_{,t}(x_0, t)}^{X_{,t}}, \qquad \overbrace{n(x_0, 0)}^{n_0} = \overbrace{n\bigl(X(x_0, t), t\bigr)}^{n} \overbrace{X_{,x_0}(x_0, t)}^{X_{,x_0}}. \end{gather*} :::{admonition} Intuition The Lagrangian description is actually the more physically natural one. Taken to its extreme, it simply describes the motion of the fluid particles using Newton's law: $F = ma = m X_{,tt}$. One has particle conservation "for free" since the number of particles is explicitly conserved -- they just move around. What is tricky is computing the force. At a microscopic level, this is just the inter-particle interactions: I.e., the [van der Waals force][] between the fluid particles. Of course, we do not have the resources to describe the [moles][mole] of particles in a typical fluid, so we must instead break the fluid up into small chunks -- fluid elements -- that we assume stay together during the course of the evolution. The net force on such an element is obtained by averaging over the microscopic forces, and is expressed instead in terms of the equation of state and gradient of the external potential. For a recent attempt to derive these forces from an averaging procedure, see {cite}`Deng:2025`. In our Lagrangian approach, we characterize the locations $X(x_0, t)$ of the center-of-mass of the fluid element that started at position $x_0$, and track its motion over time, assuming that it "stays together", thus, the original density $n_0(x_0)$ is stretched by the factor $X_{,x_0}$. This is very much in the spirit of [smoothed-particle hydrodynamics][] (SPH). A more flexible -- but complicated -- approach would be to track the boundaries of each fluid element, allowing fluid to flow from one to the other. This approach would allow us to relax the mesh as a simulation evolved, but we will not discuss it here. ::: The following differential is useful: \begin{gather*} \d{X} = X_{,x_0}\d{x_0} + X_{,t}\d{t}. \end{gather*} This gives us the following differentiation rules \begin{align*} \pdiff{f(x, t)}{x} &= \pdiff{f(x_0, t)}{x_0}\frac{1}{X_{,x_0}}, & \pdiff{f(x, t)}{t} &= \pdiff{f(x_0, t)}{t} - \pdiff{f(x_0, t)}{x_0}\frac{X_{,t}}{X_{,x_0}},\\ \pdiff{f(x, t)}{x} &= \frac{f_{,x_0}}{X_{,x_0}}, & \pdiff{f(x, t)}{t} &= f_{,t} - f_{,x_0}\frac{X_{,t}}{X_{,x_0}}. \end{align*} *(In the second form, we suppress arguments -- everything on the right-hand side in the Lagrangian framework is a function of $x_0$ and $t$, while everything on the left-hand side is formulated in the Eulerian framework as a function of $x$ and $t$.)* Thus, we can convert from Eulerian to Lagrangian with the following formulae: \begin{align*} u &= X_{,t}, & n &= \frac{n_0}{X_{,x_0}},\\ u_{,x} &= \frac{X_{,tx_0}}{X_{,x_0}}, & n_{,x} &= \frac{[n_{0}]_{,x_0}}{X_{,x_0}^2} - \frac{n_0X_{,x_0 x_0}}{X_{,x_0}^3},\\ u_{,t} &= X_{,tt} - \underbrace{X_{,t}}_{u}\underbrace{\frac{X_{,tx_0}}{X_{,x_0}}}_{u_{,x}}, & n_{,t} &= - \underbrace{\frac{n_0}{X_{x_0}}}_{n}\underbrace{\frac{X_{,x_0 t}}{X_{x_0}}}_{u_{,x}} - \underbrace{X_{,t}}_{u}\underbrace{\left( \frac{[n_{0}]_{,x_0}}{X_{,x_0}^2} - \frac{n_0X_{,x_0 x_0}}{X_{,x_0}^3}\right)}_{n_{,x}}. \end{align*} As a check, note the physical content of the last two equations. The second is just the continuity equation $n_{,t} = -(nu)_{,x}$, while the first shows that the covariant (material) derivative is just the acceleration of the fluid elements: \begin{gather*} D_t(u) = u_{,t} + uu_{,x} = X_{,tt}(x_0, t). \end{gather*} For the viscosity term, we need one more: \begin{gather*} u_{,xx} = \frac{X_{,tx_0x_0}}{X_{,x_0}^2} - \frac{X_{,tx_0}X_{,x_0x_0}}{X_{,x_0}^3},\\ \frac{n_{,x}u_{,x}}{n} = \frac{X_{,tx_0}}{n_0}\left( \frac{[n_{0}]_{,x_0}}{X_{,x_0}^2} - \frac{n_0X_{,x_0 x_0}}{X_{,x_0}^3}\right). \end{gather*} ```{code-cell} :tags: [margin, hide-input] import sympy x, x0, t = sympy.var('x,x_0,t') X = x0**2/t n0 = sympy.sin(x0) ss = [(x0, sympy.solve(X-x, x0)[1])] u = X.diff(t).subs(ss) n = (n0/X.diff(x0)).subs(ss) X_t = X.diff(t) X_tt = X_t.diff(t) X_x0 = X.diff(x0) X_x0x0 = X_x0.diff(x0) X_x0t = X_x0.diff(t) assert (u.diff(x) - (X_x0t/X_x0)).subs(ss) == 0 assert (u.diff(t) - (X_tt - X_x0t*X_t/X_x0)).subs(ss) == 0 assert (n.diff(x) - (n0.diff(x0)/X_x0**2 - n0*X_x0x0/X_x0**3)).subs(ss).simplify() == 0 ``` Thus, the hydrodynamic formulation is \begin{gather*} X_{,tt} = -\frac{V'(X) + \mathcal{E}''(n)n_{,x}}{m} + \frac{\Bigl(\nu(n)u_{,x}\Bigr)_{,x}}{n}. \end{gather*} While one could write this out explicitly, we keep this form and compute the various parts independently. ## Vacuum Interface The advantage of the Lagrangian formulation is that we can now model the vacuum interface. A simple example is if we have a fluid that is initially confined to some compact region, then the boundary will follow the evolution of the outermost fluid elements. For this to work, we must be careful about the boundary conditions. :::{todo} Ed: Explore and make rigorous the boundary conditions. ::: ## Numerical Implementation :::{margin} The viscosity used here is relevant for a quantum gas, but not for shallow water flow which would slow due to shear viscosity from the stream bed. For some rigorous results in quantum systems, see {cite}`Wu:2014`. ::: Here we demonstrate an initial numerical simulation using a finite-difference approach. We model a harmonically trapped gas with a small bump in the middle, starting from the ground state but suddenly moving the potential. Without the bump, the solution will smoothly oscillate back and forth with a constant velocity profile, hence without any dissipation since our viscosity term contains $u_{,x} = 0$: the bump will break this uniformity, inducing drag, and slowing the motion. ```{code-cell} 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 # Free bc's (second derivative is zero). 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) ``` ```{code-cell} :tags: [hide-input] %%time from scipy.integrate import solve_ivp y0 = pack(X0, 0*X0) Tw = 2*np.pi / w # Trapping period T = 40*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) print(f"{res.nfev=}") 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$"); ``` Now we show the animation: ```{code-cell} :tags: [hide-input] %%time data = [] for i, (V0_m, res) in enumerate(ress): X, X_t = res.y.reshape((2, Nx, Nt)) 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 data.append((V0_m, x_, b_, X, b, h)) fig, axs = plt.subplots(1, 2, figsize=(6, 3), sharey=True, constrained_layout=True) for i in FPS(Nt, fig=fig, embed=True, fps=10): t = ts[i] for ax, (V0_m, x_, b_, X, b, h) in zip(axs, data): ax.cla() ax.plot(x_, b_, 'k-') ax.fill_between(X[:, i], b[:, i], h[:, i], color='cyan') ax.set(title=f"$V_0/m = {V0_m:.1f}$", xlabel="$x$", xlim=(X.min(), X.max()), ylim=(0, h.max())) fig.suptitle(rf"$t={t/Tw:.2f}T_{{\omega}}$") ``` ## Spectral Basis The finite difference approach used above works reasonably well, but we would like to see if we could apply a spectral method using a {ref}`sec:ChebyshevBasis`. We do this here, imposing zero Laplacian boundary conditions ```{code-cell} 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 # Initial density provile n0 = w**2/2/g_m * (1 - x0**2) n0_x0 = -w**2 / g_m * x0 # To do: Make this work with arrays f. def diff(f, d=1, k=k, th=th0): ft = sp.fft.dct(f, type=2, norm='forward') 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 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 unpack_xnu(y): X, X_t = unpack(y) X_x0 = np.array([diff(x) for x in 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 # Zero-laplacian boundary conditions u_xx[0] = u_xx[-1] = 0 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) ``` ```{code-cell} :tags: [hide-input] %%time from scipy.integrate import solve_ivp y0 = pack(X0, 0*X0) Tw = 2*np.pi / w # Trapping period T = 40*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) print(f"{res.nfev=}") 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$"); ``` Now we show the animation: ```{code-cell} :tags: [hide-input] %%time data = [] for i, (V0_m, res) in enumerate(ress): X, X_t = res.y.reshape((2, Nx, Nt)) X_x0 = np.transpose([diff(x) for x in X.T]) 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 data.append((V0_m, x_, b_, X, b, h)) fig, axs = plt.subplots(1, 2, figsize=(6, 3), sharey=True, constrained_layout=True) for i in FPS(Nt, fig=fig, embed=True, fps=10): t = ts[i] for ax, (V0_m, x_, b_, X, b, h) in zip(axs, data): ax.cla() ax.plot(x_, b_, 'k-') ax.fill_between(X[:, i], b[:, i], h[:, i], color='cyan') ax.set(title=f"$V_0/m = {V0_m:.1f}$", xlabel="$x$", xlim=(X.min(), X.max()), ylim=(0, h.max())) fig.suptitle(rf"$t={t/Tw:.2f}T_{{\omega}}$") ``` ## Comparison Let's now do a comparison between the codes using the code in `viscosity.py`. We start with flow over a barrier in a periodic box. ```{code-cell} import viscosity class Flow(viscosity.FlowBase): Nx = 128 xR = 5.0 xL = -5.0 Basis = viscosity.BasisFFT n0 = 1.0 v0 = 0.0 # Potential V0 = 0.1 V_sigma = 1.0 V_grad = 0.2 def init(self): self.basis = self.Basis(Nx=self.Nx, Lx=self.xR - self.xL, xL=self.xL) x = self.basis.x zero = np.zeros_like(x) self.y0 = self.pack(self.n0 + zero, self.v0 + zero) def get_Vext(self, x, t, d=0): if d == 0: return self.V0 * np.exp(-(x/self.V_sigma)**2/2) + x*self.V_grad elif d == 1: return - self.V0 * x / self.V_sigma**2 * np.exp(-(x/self.V_sigma)**2/2) + x else: raise NotImemplementedError f = Flow() ``` ## Restructuring the Grid The Lagragian approach seems to work quite well, but the resulting meshes can end up being quite distorted -- especially if shockwaves present -- leading to numerical instabilities. The problem is that the Lagrangian framework naturally puts particles where the density is high, but we need instead lots of grid points where the gradient or curvature of the density is high. Consider the state at time $t$ (which we will suppress now). The Lagrangian approach encodes the density in the functions $X(x_0)$ and $\dot{X}(x_0)$ of the initial coordinate $x_0$, along with the initial density profile $n_0(x_0)$: \begin{gather*} n\bigl(X(x_0)\bigr) = \frac{n_0(x_0)}{X'(x_0)}, \qquad u\bigl(X(x_0)\bigr) = \dot{X}(x_0). \end{gather*} The same state can be described by a new pair of functions $Y(y_0)$ and $\tilde{n}_0(y_0)$ if \begin{gather*} n\bigl(X(x_0)\bigr) = n\bigl(Y(y_0)\bigr)\\ \frac{n_0(x_0)}{X'(x_0)} = \frac{\tilde{n}_0(y_0)}{Y'(y_0)} \end{gather*} where \begin{gather*} X(x_0) = Y(y_0), \qquad x_0 = X^{-1}\bigl(Y(y_0)\bigr). \end{gather*} :::{admonition} Systematic Solution :class: dropdown It can be productive to develop intuitions where you can solve these problems quickly with shortcuts. This allows you to quickly explore different options and discard false starts. But when you are stuck, having a systematic approach is also useful. To proceed systematically, we note that $x_0$ is a dummy variable, thus the total variation of the physical density $\d{n} = 0$ should be zero: \begin{gather*} \d n\Bigl(X(x_0)\Bigr) = n'(X)X'(x_0)\d{x_0} = \left(\frac{{n}_0'(x_0)}{X_{,x_0}(x_0)} - \frac{{n}_0(x_0)X_{,x_0}'(x_0)}{X_{,x_0}^2(x_0)}\right)\d{x_0} = 0. \end{gather*} This requires that \begin{gather*} n_0'(x_0)X_{,x_0}(x_0) - n_0(x_0)X_{,x_0}'(x_0) = 0. \end{gather*} One can proceed systematically to solve this, but some insight might be gained by considering a similar type of equation: \begin{gather*} f'g + fg' = (fg)' = 0. \end{gather*} This is thus solved by constant $fg$. Here we have the wrong sign, but recall that signs come from the exponents: $(fg^a)' = f'g^a + afg^{a-1}g'$. Dividing by the ugly factor of $g^{a-1}$, we have \begin{gather*} \frac{(fg^{a})'}{g^{a-1}} = f'g + afg'. \end{gather*} Thus, this condition can be expressed as \begin{gather*} X_{,x_0}\left(\frac{n_0(x_0)}{X_{,x_0}(x_0)}\right)' = 0. \end{gather*} Of course, if we had been observant, we would have noticed this since we obtained our express by differentiating the original fraction. Such a systematic approach might help, however, in more complicated cases where you have multiple variables. ::: Practically speaking, we can use this freedom to adjust the spacing between particles, compensating through the function $n_0(x_0)$. From a numerical perspective, we would like keep the transformation $X^{-1} \circ Y$ smooth and close to the identity -- especially near the boundaries -- so we do not spoil the nice properties of the Chebyshev basis. To do this, we can construct a heuristic $h(x_0)$ which measures the sensitivity of the problem at position $X(x_0)$. We want to choose $Y(y_0)$ such that the spacing between the points $1/Y'(y_0) \sim h(y_0)$ while maintaining the properties of the basis. :::{margin} Since we are fixing $t$, $X_{,x_0}(x_0) \equiv X'(x_0)$ etc. We use both notations depending on which is clearer. ::: As simple example is scaling by a multiplicative factor: \begin{gather*} n_0(x_0) \rightarrow f(x_0)n_0(x_0), \qquad X_{,x_0}(x_0) \rightarrow f(x_0)X_{,x_0}(x_0),\\ X(x_0) \rightarrow \int^{x_0}f(x)X'(x)\d{x}. \end{gather*} To be careful, we want to ensure that the extent of the cloud is not modified, $X(x_{L,R}) \rightarrow X(x_{L,R})$, and that $X(x_0)$ remains monotonic. The latter is ensured if $f$ is positive. This can be ensured by expressing the transformation in terms of a new function $g(x_0)$ \begin{gather*} X(x_0) \rightarrow X(x_L) + \int_{x_L}^{x_0}f(x)X'(x)\d{x},\\ f(x_0) = 1 + g(x_0), \qquad \int_{x_L}^{x_R}g(x)X'(x)\d{x} = 0, \qquad \abs{g(x_0)} < 1. \end{gather*} One more constraint concerns the derivatives at the boundaries: \begin{align*} X'(x_0) &\rightarrow f(x_0)X'(x_0),\\ X''(x_0) &\rightarrow f'(x_0)X'(x_0) + f(x_0)X''(x_0),\\ &\vdots \end{align*} To keep these unchanged we need $f(x_0) = 1$, $f'(x_0) = 0$, etc. up to one order less than the desired invariance. For example, to keep up to second derivatives, we have the following constraints on $g(x_0)$ at the boundaries: \begin{gather*} g(x_{L,R}) = g'(x_{L,R}) = 0. \end{gather*} To complete this, we also need to store \begin{gather*} n'(x_0) \rightarrow f'(x_0)n(x_0) +f(x_0)n'(x_0). \end{gather*} Another simple example is adding a term: \begin{gather*} n_0(x_0) \rightarrow n(x_0) + f(x_0), \qquad X_{,x_0} \rightarrow \left(1 + \frac{f(x_0)}{n(x_0)}\right)X_{,x_0}(x_0),\\ X(x_0) \rightarrow \int^{x_0}\left(1 + \frac{f(x_0)}{n(x_0)}\right)X'(x_0)\d{x}. \end{gather*} To be careful, we want to ensure that the extent of the cloud is not modified, $X(x_{L,R}) \rightarrow X(x_{L,R})$, and that $X(x_0)$ remains monotonic. The latter is ensured if $f(x_0) + n(x_0) > 0$ is positive. \begin{gather*} X(x_0) \rightarrow X(x_L) + \int_{x_L}^{x_0}\left(1 + \frac{f(x_0)}{n(x_0)}\right)X'(x_0)\d{x},\\ X(x_0) \rightarrow X(x_0) + \int_{x_L}^{x_0}\frac{f(x_0)}{n(x_0)}X'(x_0)\d{x},\\ X(x_0) \rightarrow X(x_0) + \int_{x_L}^{x_0}\frac{n'(x_0)f(x_0) - n(x_0)f'(x_0)}{n^2(x_0)}X(x_0)\d{x},\\ \int_{x_L}^{x_R} \frac{f(x_0)}{n(x_0)}X'(x_0)\d{x} = 0,\qquad f(x_0) + n(x_0) > 0. \end{gather*} One more constraint concerns the derivatives at the boundaries: \begin{align*} X'(x_0) &\rightarrow X'(x_0) + \frac{f(x_0)}{n(x_0)}X'(x_0),\\ X''(x_0) &\rightarrow X''(x_0) + \frac{f'(x_0)}{n(x_0)}X'(x_0) + \frac{f(x_0)}{n(x_0)}\left( X''(x_0) - \frac{n'(x_0)}{n(x_0)}X'(x_0)\right),\\ &\vdots \end{align*} To keep these unchanged we need $f(x_0) = 0$, $f'(x_0) = 0$, etc. up to one order less than the desired invariance. For example, to keep up to second derivatives, we have the following constraints on $f(x_0)$ at the boundaries: \begin{gather*} f(x_{L,R}) = f'(x_{L,R}) = 0. \end{gather*} To complete this, we also need to store \begin{gather*} n'(x_0) \rightarrow n'(x_0) + f'(x_0). \end{gather*} For example, suppose we take \begin{gather*} f(x_0) = (1-x_0^2)^2P(x_0), \qquad f'(x_0) = (1-x_0^2)^2P'(x_0) - 4x_0(1-x_0^2)P(x_0). \end{gather*} :::{admonition} A Strategy A strategy is to define a heuristic like $n_{,x_0x_0}^2$ which is large in the regions where we need more grid points. We then form an appropriate $g(x_0)$ by scaling this and multiplying by a weighting function $w(x_0)$ which vanishes sufficiently fast at the boundaries to ensure an appropriate number of derivatives of $g$ are zero. For example, working with the Chebyshev basis with $x_0\in [-1, 1]$: \begin{gather*} w(x_0) = \cos^{2p}(\theta_0 - \tfrac{\pi}{2}) = (1-x_0^2)^{p} \end{gather*} will ensure that $p$ derivatives of $g$ vanish at the boundaries. This will give us a suitable function $g(x_0)$ tabulated at the current abscissa. We now need to consistently to get our transformations: \begin{align*} f(x_0) &= 1 + w(x_0)g(x_0),\\ n_0(x_0) &\rightarrow f(x_0)n_0(x_0),\\ X(x_0) &\rightarrow X(x_L) + \int_{x_L}^{x_0}f(x)X'(x)\d{x}\\ &= X(x_L) + f(x_0)X(x_0) - f(x_L)X(x_L) - \int_{x_L}^{x_0}f'(x)X(x)\d{x}\\ \end{align*} Once we have a suitable function $g(x_0)$, we can integrate it numerically to obtain the appropriate coordinate transform. Splines can be used to simplify this process. Assuming that $X(x_0)$ is reasonably well approximated in the basis, then ::: ## Hybrid Approach Can we develop a hybrid approach that resolves this? The idea here is to allow both the density and grid points to vary so that we may guide the grid points to a more appropriate distribution. We start with the Generally, we might consider evolving a density $N(x_0, t)$ and the grid $X(x_0, t)$ such that \begin{gather*} n(X, t) = f(X, N), \qquad u(X, t) = g(X, N). \end{gather*} One must then find suitable forms for $f$ and $g$ such that the continuity equation is satisfied. Perhaps there is a nice choice that simplifies things, but currently it seems like a bit of a mess. Possibly relevant: * [Dynamic-mesh finite element method for Lagrangian computational fluid dynamics](https://www.sciencedirect.com/science/article/pii/S0168874X02000884) [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]: