Astrophysical Fluids Lec 7 | Introduction to Numerical Fluid Eq. Solvers

Peter Ma | Jan 5th 2025

Over the last few lectures we have shown how to model various kinds of astrophysical problems with fluid equations and finally getting out the famed Navier Stokes Equation. Sadly we have yet to tame such an equation analytically, and so until someone does that, we can instead solve it numerically and beat the fluid equations with compute alone. To get started let us introduce ourselves to how numerical simulatiors / integrators work, their draw backs and design. We start with the basics.

<< Home Page

Numerical errors and instabilities

Aside from errors in programming or discretizing the physical model, there are 2 types of errors in computations:

1. Rounding errors: errors in how the computer stores or manipulates numbers.
2. Approximation errors (sometimes called truncation errors): errors in approximations to various functions or methods.

Reason: computers are machines, with inherent limitations.

Machine error: round-off error.

Under what circumstances is the following possible? $$(x+y)+z \neq x + (y + z)$$

Let's try it in Python:

x = 1e16  # this will actually work
y = -1e16  # need to decrease
z = 1.

print("(x+y)+z = {}".format((x+y)+z))
print("x+(y+z) = {}".format(x+(y+z)))

(x+y)+z = 1.0
x+(y+z) = 0.0

Lmao, this is not right, what happened? Round-off error!

Algorithmic error: instability in Integration

• Consider this system representing phasor rotation in the complex plane: $$\dot Z = i\omega Z, \quad\text{given}\quad Z_0 = Z(t=0).$$

• Analytical solution is $Z(t) = Z_0 \exp(i\omega t).$

• How can we solve it numerically?

Taylor expansion: $$\dot Z(t) = \frac{Z(t+\Delta t)-Z(t)}{\Delta t} + H.O.T. = i\omega Z(t).$$

Suggests algorithm:

• Start with $Z(t=0)=Z_{old}$,
• $Z_{new} = Z_{old}+i\omega Z_{old}\Delta t$,
• repeat.

Let's code it up and see

# Unstable solution to dZ/dt
from numpy import pi

# initialize Z0, omega and dt below;
Z = complex(1., 0.)
nsteps = 200
omega = 1.
dt = 2*pi/(nsteps*omega)

# print initial value
print('t = {0}, Z = {1}, |Z| = {2}'.format(0., Z, abs(Z)))
# rotate by 2pi
for k in range(nsteps):
    Z = (1+1j*omega*dt) * Z
#print final value
print('t = {0}, Z = {1}, |Z| = {2}'.format(k*dt, Z, abs(Z)))
#print('t = {0:.3f}'.format(k*dt))
#print('Z = {0:.3f}'.format(Z))
print('|Z| = {0:.3f}'.format(abs(Z)))

t = 0.0, Z = (1+0j), |Z| = 1.0
t = 6.2517693806436885, Z = (1.1036746878104917-0.002280042726220316j), |Z| = 1.1036770429380236
|Z| = 1.104

What happened? What is the problem? Why did it happen?

Propagation of errors

Errors propagate statistically like they do in experimental physics.

These errors in differences become very important when taking numerical derivatives: $$\frac{df}{dt} \approx \frac{f_{i+1} - f_i}{\Delta t} = \text{danger zone}$$ as we will see later.

One important rule

Never, ever. Never ever ever. Do something like this:

if 7./3. - 4./3. - 1. == 0:
    print('7/3 - 4/3 - 1 == 0')
else:
    print('7/3 - 4/3 - 1 != 0')
    

7/3 - 4/3 - 1 != 0

Instead:

delta = 1e-15
if abs(7./3. - 4./3. - 1.) < delta:
    print('7/3 - 4/3 - 1 == 0 (or close enough anyway...)')
else:
    print('7/3 - 4/3 - 1 != 0')
    

7/3 - 4/3 - 1 == 0 (or close enough anyway...)

Ordinary Differential Equation Solvers

Before we get to PDE's let us crack the ODE problem first. With two methods the Euler Method and the Runge-Kutta Methods of integrations.

Euler Method

• Start at $t = t_0,~ x = x_0$
• Discretize into time-steps $t_i$ with constant spacing $h$
• At each time-step, find $x$, using $x$ at the previus time-step and $f$:

$\quad x_{i} = x_{i-1} + hf(x_i) $

Example Simple Harmonic Oscillator (Forward Euler)

Consider the simple problem $$\frac{d^2x}{dt^2} + w^2 x = 0$$ We can solve this analytically and obtain the solution set of $x(t) = C_0e^{iw t}$ To numerically integrate this we have the simple "forward" scheme. We define a velocity $v_i$ and displacement $x_i$ and write the second order equation as 2 first order equations coupled together. $$x_{i+1} = x_i + \Delta t v_i $$ $$v_{i+1} = v_i - \Delta t w^2 x_i$$ We call this the forward Euler method becuase of the order at which we updated the equations. First position then momenta.

import numpy as np
import matplotlib.pyplot as plt

# Parameters
omega = 1.0  # Angular frequency
dt = 0.01    # Time step
T = 10       # Total time

# Time array
time = np.arange(0, T, dt)

# Initial conditions
x0 = 1.0     # Initial displacement
v0 = 0.0     # Initial velocity

# Arrays to store results
x = np.zeros(len(time))
v = np.zeros(len(time))

# Initial values
x[0] = x0
v[0] = v0

# Forward Euler loop
for i in range(len(time) - 1):
    x[i + 1] = x[i] + dt * v[i]
    v[i + 1] = v[i] - dt * omega**2 * x[i]

# Plot results
plt.figure(figsize=(10, 5))
plt.plot(time, x, label="Displacement (x)")
plt.plot(time, v, label="Velocity (v)")
plt.title("Simple Harmonic Oscillator: Forward Euler Method")
plt.xlabel("Time (t)")
plt.ylabel("Amplitude")
plt.legend()
plt.grid()
plt.show()

Problems

The forward Euler method is NOT energy conserving. We can check the plot in the bottom for this case.

E_t = []
E_k_total = []
E_p_total = []
for i, t in enumerate(time):
    E_k = 0.5 * v[i]**2 
    E_p = 0.5 * x[i]**2
    E_k_total.append(E_k)
    E_p_total.append(E_p)
    E_t.append(E_k +E_p)

plt.plot(time, E_t , label='total energy')
plt.plot(time, E_k_total, label="kinetic energy")
plt.plot(time, E_p_total, label='spring potential')
plt.legend()
plt.grid()
plt.show()

We can further see that this is because the forward Euler method produces a phase space that is NOT closed. This is

plt.plot(x,v)
plt.xlabel("position")
plt.ylabel("velocity")
plt.grid()
plt.show()

Example Simple Harmonic Oscillator (Backward Euler)

We consider the Euler integration but back subsituting the values of $x_i$ after updating $v_i$. Let us first write out the forward case $$x_{i+1} = x_i + \Delta t v_{i+1}$$ $$v_{i+1} = v_i - w^2 \Delta t x_{i+1} $$ But let us say we only have access to the $v_i$ how do we update $v_{i+1}$ before we compute it? Well we can just rearrange to solve for the expressions we want (right hand side to not have any forward terms). Sub eq2 into eq1 $$x_{i+1} = x_i + \Delta t [v_i - w^2 \Delta t x_{i+1}]$$ $$x_{i+1} = x_i + \Delta t v_i - w^2 \Delta t^2 x_{i+1}$$ $$x_{i+1} + w^2 \Delta t^2 x_{i+1} = x_i + \Delta t v_i $$ $$x_{i+1} (1+ w^2 \Delta t^2 ) = x_i + \Delta t v_i $$ $$\boxed{x_{i+1} = \frac{x_i + \Delta t v_i}{1+ w^2 \Delta t^2 }} $$ $$\boxed{v_{i+1} = v_i - w^2 \Delta t x_{i+1}} $$ Now this is something we can write in a algorithm! The back subsitution makes it a "backward" euler.

import numpy as np
import matplotlib.pyplot as plt

# Parameters
omega = 1.0  # Angular frequency
dt = 0.01    # Time step
T = 10       # Total time

# Time array
time = np.arange(0, T, dt)

# Initial conditions
x0 = 1.0     # Initial displacement
v0 = 0.0     # Initial velocity

# Arrays to store results
x = np.zeros(len(time))
v = np.zeros(len(time))

# Initial values
x[0] = x0
v[0] = v0

# Backward Euler loop
for i in range(len(time) - 1):
    # Solve for v_{i+1} first
    v_next = (v[i] - dt * omega**2 * x[i]) / (1 + dt**2 * omega**2)
    # Solve for x_{i+1} using the updated v_{i+1}
    x_next = x[i] + dt * v_next
    
    # Update the arrays
    v[i + 1] = v_next
    x[i + 1] = x_next

# Plot results
plt.figure(figsize=(10, 5))
plt.plot(time, x, label="Displacement (x)")
plt.plot(time, v, label="Velocity (v)")
plt.title("Simple Harmonic Oscillator: Backward Euler Method")
plt.xlabel("Time (t)")
plt.ylabel("Amplitude")
plt.legend()
plt.grid()
plt.show()

E_t = []
E_k_total = []
E_p_total = []
for i, t in enumerate(time):
    E_k = 0.5 * v[i]**2 
    E_p = 0.5 * x[i]**2
    E_k_total.append(E_k)
    E_p_total.append(E_p)
    E_t.append(E_k +E_p)

plt.plot(time, E_t , label='total energy')
plt.plot(time, E_k_total, label="kinetic energy")
plt.plot(time, E_p_total, label='spring potential')
plt.legend()
plt.grid()
plt.show()

Energy is still not conserved !..... this is also wrong!

Example Simple Harmonic Oscillator (Symplectic Euler)

Okay so none of these really worked that well... We notice that one over estimates and the other underestimates the total energy. Let us combine them together through a central difference technique. Or better called the Symplectic Euler method. The key is to calculate "half step" corrections.

$$\boxed{v_{i+1/2} = v_i - w^2 \frac{\Delta t}{2} x_{i}} $$

Then using the half velocity correction we update the $x_i$ $$\boxed{x_{i+1} = x_i + \frac{\Delta t}{2} v_{i+1/2}} $$ Then we fully update the velocity $$\boxed{v_{i+1} = v_{i+1/2} - w^2 \frac{\Delta t}{2} x_{i+1}} $$ This is basically the combination of both the foreward and backward euler implementation! Here the the mention of symplectic means that this method preserves the symplectic structure of the system. This just means a phase space volume is preserved in a hamiltonian system. (i think)

import numpy as np
import matplotlib.pyplot as plt

# Parameters
omega = 1.0  # Angular frequency
dt = 0.01    # Time step
T = 10       # Total time

# Time array
time = np.arange(0, T, dt)

# Initial conditions
x0 = 1.0     # Initial displacement
v0 = 0.0     # Initial velocity

# Arrays to store results
x = np.zeros(len(time))
v = np.zeros(len(time))

# Initial values
x[0] = x0
v[0] = v0 - 0.5 * dt * omega**2 * x0  # Half-step velocity update

# Central difference loop
for i in range(len(time) - 1):
    # Update position
    x[i + 1] = x[i] + dt * v[i]
    # Update velocity at the next half-step
    v[i + 1] = v[i] - dt * omega**2 * x[i + 1]

# Plot results
plt.figure(figsize=(10, 5))
plt.plot(time, x, label="Displacement (x)")
plt.plot(time, v, label="Velocity (v)")
plt.title("Simple Harmonic Oscillator: Central Difference Method")
plt.xlabel("Time (t)")
plt.ylabel("Amplitude")
plt.legend()
plt.grid()
plt.show()

E_t = []
E_k_total = []
E_p_total = []
for i, t in enumerate(time):
    E_k = 0.5 * v[i]**2 
    E_p = 0.5 * x[i]**2
    E_k_total.append(E_k)
    E_p_total.append(E_p)
    E_t.append(E_k +E_p)

plt.plot(time, E_t , label='total energy')
plt.plot(time, E_k_total, label="kinetic energy")
plt.plot(time, E_p_total, label='spring potential')
plt.legend()
plt.grid()
plt.show()

Energy is conserved! (pretty much at least)

• The Euler method has error $O(h^2)$ at each step
• integrating across the whole interval: global error is $O(h)$: $$\text{Taylor expansion} \Rightarrow x(t+h) = x(t) + h\frac{dx}{dt} + \overbrace{ \frac{h^2}{2} \frac{d^2x}{dt^2} } ^{\epsilon} + O(h^3)$$ $$\sum\epsilon = \sum_{k=0}^{N-1}\frac{h^2}{2}\left. \frac{d^2x}{dt^2} \right|_{x_k, t_k} = \frac{h}{2}\sum_{k=0}^{N-1}h\left.\frac{df}{dt}\right|_{x_k, t_k}\\ \approx \frac{h}2\int_a^b\frac{df}{dt}d t = \frac{h}{2}\left[f_b - f_a\right]$$
• For some applications, this is good enough. But for others, we need to do better!

Runge-Kutta (RK) Methods (2nd order)

The RK method looks at using the mid points to evaluate Eulers method. Intuitvely the Euler method is just piecing together slopes together to approximate the general curve and trajectory given IC. Instead a better way is estimate the next timestep is the compute the slope midway between two points and use that slope to update you to the next step.

• Use the middle point $t+h/2$ and evaluate with Euler's method,

$\quad\displaystyle x\left(t + \frac{h}2\right) \approx x(t) + \frac{h}2 f[x(t), t]$

• Slope at $t + \frac{h}2 \approx f \! \left[ x(t) + \frac{h}2 f\left[x(t), t\right], t + \frac{h}2\right]$

$\displaystyle \Rightarrow \boxed{x(t+h) = x(t) + h f\! \left[x(t) + \frac{h}2 f[x(t), t], t+\frac{h}2\right]}$

Example Solving a 1st order ODE

Let us use a simple toy example $$\frac{dx}{dt} = -kx$$ And let us try to integrate out the solution

import numpy as np
import matplotlib.pyplot as plt

# Parameters
k = 1.0       # Decay constant
dt = 0.3      # Time step
T = 5.0       # Total time

# Time array
time = np.arange(0, T + dt, dt)

# Initial condition
y0 = 1.0  # Initial value of y

# Array to store results
y = np.zeros(len(time))
y[0] = y0

# RK2 loop
for i in range(len(time) - 1):
    t_i = time[i]
    y_i = y[i]
    
    # Compute k1 and k2
    k1 = -k * y_i
    k2 = -k * (y_i + k1 * dt)
    
    # Update y
    y[i + 1] = y_i + dt * (k1 + k2) / 2  # RK2 update formula

# Analytical solution for comparison
y_analytical = y0 * np.exp(-k * time)

# Plot results
plt.figure(figsize=(10, 5))
plt.plot(time, y, label="RK2 Numerical Solution")
plt.plot(time, y_analytical, label="Analytical Solution", linestyle="dashed")
plt.title("Runge-Kutta 2nd Order Method for ODE: dy/dt = -ky")
plt.xlabel("Time (t)")
plt.ylabel("y(t)")
plt.legend()
plt.grid()
plt.show()

Benefits:

• Improved Accuracy: More accurate than the Euler method for the same time step size.
• Intermediate Step: Takes the slope at an intermediate point for better estimation.
• Error: Local error is $\mathcal{O}(h^3)$ instead of $\mathcal{O}(h^2)$

For even higher accuracy, consider using Runge-Kutta 4th Order (RK4), which is a widely used numerical integration method.

RK (4th Order)

To get the 4th order expansion you can implement by hand but its annoying but the framework for doing so is the following: This is the common go-to integration scheme used!

• Various Taylor expansions at various points in the interval $\Rightarrow$ higher-order RK's.
• RK4 is reasonable to code yourself. Higher-order methods $\Rightarrow$ use canned routines.
• End result, after tedious algebra:
  1. $k_1 = hf(x, t)$,
  2. $k_2 = hf\left(x + \frac{k_1}{2}, t+\frac{h}2\right)$,
  3. $k_3 = hf\left(x + \frac{k_2}{2}, t+\frac{h}2\right)$,
  4. $k_4 = hf\left(x + k_3, t + h \right)$,
  5. $x(t+h) = x(t) + \frac{1}{6}(k_1 + 2 k_2 + 2k_3 + k_4)$.

Partial Differential Equation Solvers

Firstly, the classification of partial differential equations (PDEs) is critical for numerical solvers because different types of PDEs exhibit different stable conditions which is needed for accurate solutions. Misclassifying a PDE can lead to instability, inefficiency, or incorrect results. Some examples.

1. Elliptic PDEs (Laplace's equation: $\nabla^2u = 0$) Key Characteristics: No time dependence. Solutions are smooth if the boundary conditions are smooth. Boundary value problems. Methods: Finite difference methods, finite element methods, spectral methods.
2. Parabolic PDEs (Heat equation: $\frac{\partial u}{\partial t} = c\nabla^2 u$) Key Characteristics: Time-dependent, with a "smoothing" effect on solutions over time. Initial and boundary value problems. Methods: Explicit and implicit time-stepping methods (e.g., Forward Euler, Backward Euler, Crank-Nicolson).
3. Hyperbolic PDEs (Wave equation: $\frac{\partial^2 u}{\partial t^2} = c\nabla^2 u$) Key Characteristics: Time-dependent, with solutions that propagate information along characteristic lines.

Solutions may develop discontinuities (shocks) even if initial conditions are smooth. Initial value and boundary value problems. Methods: Explicit methods like finite difference, finite volume, or finite element methods.

We will see what makes the Navier Stokes Eq solvers to be so difficult is because it combines the challenges of multiple classes of PDEs. This requires hybrid solvers that adapt to each behavior.

We will not have time to explore each of these. But we will explore 2 example to understand the general framework of a PDE solver. Before we approach a full NSE problem.

Elliptical PDE Solver (2-D Laplace Equation)

Example problem: Modeling Electric potential. We want to model the electric potential for an empty 2D box, 10cm x 10cm in size, where the top wall is held at $V$ = 1.0V and the other walls at 0V.

$$0 = \nabla^2 \phi = \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2},$$

We have boundary conditions. We note that this is time independent. $$ \phi(y = 10) = 1.0 V$$ $$ \phi(y = 0) = \phi(x = 0) = \phi(x = 10) = 0$$

To solve this we will use a method called Jacobi Relaxation the update scheme is derived as follows:

Recall that finte element difference to computing derivative is

$$ \frac{\partial \phi(x,y)}{\partial x} \approx \frac{ \phi(x+\Delta x,y) - \phi(x,y)}{\Delta x}$$

We can compute the approximate second order version of this $$ \frac{\partial^2 \phi(x,y)}{\partial x^2} \approx \frac{\frac{\partial}{\partial x} \phi(x+\Delta x,y) - \frac{\partial}{\partial x}\phi(x,y)}{\Delta x}$$ $$ \frac{\partial^2 \phi(x,y)}{\partial x^2} \approx \frac{\frac{ \phi(x+2\Delta x,y) - \phi(x+\Delta x,y)}{\Delta x} - \frac{ \phi(x+\Delta x,y) - \phi(x,y)}{\Delta x}}{\Delta x}$$ $$ \frac{\partial^2 \phi(x,y)}{\partial x^2} \approx \frac{\phi(x+2\Delta x,y) - 2\phi(x+\Delta x,y) + \phi(x,y)}{\Delta x^2}$$ We can repeat for the y version $$ \frac{\partial^2 \phi(x,y)}{\partial y^2} \approx \frac{\phi(x ,y+2\Delta y) - 2\phi(x,y+\Delta y) + \phi(x,y)}{\Delta y^2}$$ Together we have the following according the Laplace equation $$ \frac{\phi(x+2\Delta x,y) - 2\phi(x+\Delta x,y) + \phi(x,y)}{\Delta x^2}+ \frac{\phi(x ,y+2\Delta y) - 2\phi(x,y+\Delta y) + \phi(x,y)}{\Delta y^2} = 0$$ We assume that the spacing is the same $\Delta x = \Delta y$ same size grids. $$ \frac{\phi(x+2\Delta x,y) - 2\phi(x+\Delta x,y) + \phi(x,y) +\phi(x ,y+2\Delta y) - 2\phi(x,y+\Delta y) + \phi(x,y)}{\Delta x^2} = 0$$ $$ \phi(x+2\Delta x,y) - 2\phi(x+\Delta x,y) + \phi(x,y) +\phi(x ,y+2\Delta y) - 2\phi(x,y+\Delta y) + \phi(x,y) = 0$$

$$\phi(x+a, y) + \phi(x-a, y) + \phi(x, y+a) + \phi(x, y-a) - 4\phi(x,y) = 0$$

• Iterate the rule

$\phi_{new}(x, y) = \frac14\left[\phi(x+a, y) + \phi(x-a, y) + \phi(x, y+a) + \phi(x, y-a)\right].$

• Much like the relaxation method for finding solutions of $f(x) = x$,
• For this problem it turns out that Jacobi Relaxation is always stable and so always gives a solution!

import numpy as np
from tqdm import tqdm
def initPhi(M):
    phi = np.zeros((M+1,M+1), np.float32)
    phi[0,:] = 1.0
    return phi

def jacobi(target, M, phi):
    delta = 1
    phinew = np.empty([M+1,M+1], np.float32)
    iterations = 0

    while delta > target:
        for i in range(M+1):
            for j in range(M+1):
                if (i==0 or i ==M or j==0 or j==M):
                    phinew[i,j] = phi[i,j]
                else:
                    phinew[i,j] = (phi[i+1,j] + phi[i-1,j] + phi[i,j+1]+ phi[i,j-1])/4
        delta = np.amax(abs(phi-phinew))
        phi, phinew = phinew, phi
        iterations += 1
    return iterations, phi

numIterations, phiJacobi = jacobi(1e-4, 100, initPhi(100))
print(numIterations)
print(phiJacobi)

1909
[[1.0000000e+00 1.0000000e+00 1.0000000e+00 ... 1.0000000e+00
  1.0000000e+00 1.0000000e+00]
 [0.0000000e+00 4.9966612e-01 6.9698507e-01 ... 6.9698513e-01
  4.9966609e-01 0.0000000e+00]
 [0.0000000e+00 3.0168000e-01 4.9866584e-01 ... 4.9866581e-01
  3.0168000e-01 0.0000000e+00]
 ...
 [0.0000000e+00 1.5559126e-05 3.1097243e-05 ... 3.1097246e-05
  1.5559128e-05 0.0000000e+00]
 [0.0000000e+00 7.7394761e-06 1.5468486e-05 ... 1.5468488e-05
  7.7394761e-06 0.0000000e+00]
 [0.0000000e+00 0.0000000e+00 0.0000000e+00 ... 0.0000000e+00
  0.0000000e+00 0.0000000e+00]]

import matplotlib.pyplot as plt
plt.imshow(phiJacobi, cmap='hot')

<matplotlib.image.AxesImage at 0x10f590a50>

Hyperbolic PDE Solver (Wave Equation)

Recall the 1D wave equation: $$\frac{\partial^2\phi}{\partial t^2}=v^2\frac{\partial^2\phi}{\partial x^2}$$

Consider a piano string of length $L$, initially at rest. At time $t=0$ the string is struck by the piano hammer a distance $d$ from the end of from the string. The string vibrates as a result of being struck, except at the ends, $x=0$, and $x=L$, where it is held fixed.

Consider the case $v=100\mathrm{ms^{-1}}$, with the initial condition that $\phi(x)=0$ everywhere but the velocity $\psi(x)$ is nonzero, with profile \begin{equation} \psi(x) = C\frac{x(L-x)}{L^2}\exp\left[-\frac{(x-d)^2}{2\sigma^2}\right], \end{equation} where $L=1$m, $d=10$cm, $C=1\mathrm{ms^{-1}}$, and $\sigma=0.3$m.

Approach using the FTCS method

import numpy as np
import matplotlib.pyplot as plt

# Parameters
L = 1.0  # Length of the domain
T = 0.11  # Total simulation time
c = 100

# Discretization
h = 1e-6
Nt = int(T/h)
Nx = 2**7
dx = L/Nx
x = np.linspace(0, L, Nx)
u = np.zeros((Nx, Nt + 1))
u_new = np.zeros(Nx)

# Initial condition
d = 0.1
C = 0.001 
sigma = 0.3

u0 =0
u1 = C*(x*(L-x))/(L**2) * np.exp( - ((x-d)**2/(2*sigma**2)))

u[:, 0] = u0
u[:, 1] = u0 + 0.5 * c * h * u1 

plt.figure(figsize=(8, 6))
for n in range(1, Nt):
    for i in range(1, Nx - 1):
        u_new[i] = 2 * u[i, n] - u[i, n - 1] + (c**2) * (h**2) / (dx**2) * (u[i + 1, n] - 2 * u[i, n] + u[i - 1, n])
        u_new[0] = 0
    u_new[-1] = 0
    u[:, n + 1] = u_new
    if n*h in [0.006, 0.004, 0.002, 0.012, 0.100002] :
        plt.plot(u_new)
plt.xlabel('Time')
plt.ylabel('Position')
plt.title('Solution of 1D Wave Equation using FTCS')

Text(0.5, 1.0, 'Solution of 1D Wave Equation using FTCS')

Toy Hydrodynamical Simulation (Lid Driven Cavity Example)

We recall that the full NSE looks like the following $$\frac{\partial \rho}{\partial t} + \nabla \rho v = 0 \tag{mass cons.}$$ $$\frac{\partial v}{\partial t} + v \cdot \nabla v = \frac{-1}{\rho}\nabla P - \nabla \Phi + \nu[\nabla^2 v + \frac{1}{3}\nabla(\nabla \cdot v)] \tag{mom. cons}$$

In our lid driven example, we assume a few things, namely the fluid is incompressible and so $\nabla \cdot u = 0 $ we also assume this is 2-D problem.

Lid Driven Cavity

In a cavity the top of the fluid has some velocity but everything inside is initially static. The top moving fluid drags the fluid internally in the cavity which creates swirly patterns (intuitively speaking). Initially velocity and pressure have zero initial condition. At the walls of the cavity we enforce a homogeneous dirichlet boundary conditions everywhere except for velocity at top. It is driven by an external flow.

Chorin's Projection

The primary challenge with these simulations is enforcing the incompressibility condition, which requires that the velocity field has zero divergence. Chorin’s Projection solves this by "projecting" the velocity field onto a divergence-free space. The method operates by splitting the computation of the velocity field into intermediate steps, bypassing the need to solve for pressure and velocity simultaneously loosely speaking.

Framework

1. Solve tentative velocity + velocity BC. Let the two velocity components be $u,v$.

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y}= \nu \frac{\partial^2 u}{\partial x^2} + \nu\frac{\partial^2 u}{\partial y^2} $$$$ \frac{\partial v}{\partial t} + v \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y}= \nu \frac{\partial^2 v}{\partial x^2} + \nu\frac{\partial^2 v}{\partial y^2} $$

2. Solve pressure poisson + pressure BC

$$ \frac{\partial^2 P}{\partial x^2} +\frac{\partial^2 P}{\partial y^2} = \frac{\rho}{\Delta t} (\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} )$$

3. Correct velocity + velocity BC project the velocity

$$u − \frac{\Delta t}{\rho}\frac{\partial P}{\partial x} \to u$$$$v −\frac{\Delta t}{\rho} \frac{\partial P}{\partial y} \to v$$

Implementation

Below is an implementation following code derived from video here. Credits to machine-learning-and-simulation youtube channel

N_POINTS = 41
DOMAIN_SIZE = 1.0
N_ITERATIONS = 500
TIME_STEP_LENGTH = 0.001
KINEMATIC_VISCOSITY = 0.1
DENSITY = 1.0
HORIZONTAL_VELOCITY_TOP = 1.0

N_PRESSURE_POISSON_ITERATIONS = 50
STABILITY_SAFETY_FACTOR = 0.5


element_length = DOMAIN_SIZE / (N_POINTS - 1)
x = np.linspace(0.0, DOMAIN_SIZE, N_POINTS)
y = np.linspace(0.0, DOMAIN_SIZE, N_POINTS)

X, Y = np.meshgrid(x, y)

u_prev = np.zeros_like(X)
v_prev = np.zeros_like(X)
p_prev = np.zeros_like(X)


def central_difference_x(f):
    diff = np.zeros_like(f)
    diff[1:-1, 1:-1] = (
        f[1:-1, 2:  ]
        -
        f[1:-1, 0:-2]
    ) / (
        2 * element_length
    )
    return diff

def central_difference_y(f):
    diff = np.zeros_like(f)
    diff[1:-1, 1:-1] = (
        f[2:  , 1:-1]
        -
        f[0:-2, 1:-1]
    ) / (
        2 * element_length
    )
    return diff

def laplace(f):
    diff = np.zeros_like(f)
    diff[1:-1, 1:-1] = (
        f[1:-1, 0:-2]
        +
        f[0:-2, 1:-1]
        -
        4
        *
        f[1:-1, 1:-1]
        +
        f[1:-1, 2:  ]
        +
        f[2:  , 1:-1]
    ) / (
        element_length**2
    )
    return diff


maximum_possible_time_step_length = (
    0.5 * element_length**2 / KINEMATIC_VISCOSITY
)
if TIME_STEP_LENGTH > STABILITY_SAFETY_FACTOR * maximum_possible_time_step_length:
    raise RuntimeError("Stability is not guarenteed")


for _ in tqdm(range(N_ITERATIONS)):
    d_u_prev__d_x = central_difference_x(u_prev)
    d_u_prev__d_y = central_difference_y(u_prev)
    d_v_prev__d_x = central_difference_x(v_prev)
    d_v_prev__d_y = central_difference_y(v_prev)
    laplace__u_prev = laplace(u_prev)
    laplace__v_prev = laplace(v_prev)

    # Perform a tentative step by solving the momentum equation without the
    # pressure gradient
    u_tent = (
        u_prev
        +
        TIME_STEP_LENGTH * (
            -
            (
                u_prev * d_u_prev__d_x
                +
                v_prev * d_u_prev__d_y
            )
            +
            KINEMATIC_VISCOSITY * laplace__u_prev
        )
    )
    v_tent = (
        v_prev
        +
        TIME_STEP_LENGTH * (
            -
            (
                u_prev * d_v_prev__d_x
                +
                v_prev * d_v_prev__d_y
            )
            +
            KINEMATIC_VISCOSITY * laplace__v_prev
        )
    )

    # Velocity Boundary Conditions: Homogeneous Dirichlet BC everywhere
    # except for the horizontal velocity at the top, which is prescribed
    u_tent[0, :] = 0.0
    u_tent[:, 0] = 0.0
    u_tent[:, -1] = 0.0
    u_tent[-1, :] = HORIZONTAL_VELOCITY_TOP
    v_tent[0, :] = 0.0
    v_tent[:, 0] = 0.0
    v_tent[:, -1] = 0.0
    v_tent[-1, :] = 0.0


    d_u_tent__d_x = central_difference_x(u_tent)
    d_v_tent__d_y = central_difference_y(v_tent)

    # Compute a pressure correction by solving the pressure-poisson equation
    rhs = (
        DENSITY / TIME_STEP_LENGTH
        *
        (
            d_u_tent__d_x
            +
            d_v_tent__d_y
        )
    )

    for _ in range(N_PRESSURE_POISSON_ITERATIONS):
        p_next = np.zeros_like(p_prev)
        p_next[1:-1, 1:-1] = 1/4 * (
            +
            p_prev[1:-1, 0:-2]
            +
            p_prev[0:-2, 1:-1]
            +
            p_prev[1:-1, 2:  ]
            +
            p_prev[2:  , 1:-1]
            -
            element_length**2
            *
            rhs[1:-1, 1:-1]
        )

        # Pressure Boundary Conditions: Homogeneous Neumann Boundary
        # Conditions everywhere except for the top, where it is a
        # homogeneous Dirichlet BC
        p_next[:, -1] = p_next[:, -2]
        p_next[0,  :] = p_next[1,  :]
        p_next[:,  0] = p_next[:,  1]
        p_next[-1, :] = 0.0

        p_prev = p_next
    

    d_p_next__d_x = central_difference_x(p_next)
    d_p_next__d_y = central_difference_y(p_next)

    # Correct the velocities such that the fluid stays incompressible
    u_next = (
        u_tent
        -
        TIME_STEP_LENGTH / DENSITY
        *
        d_p_next__d_x
    )
    v_next = (
        v_tent
        -
        TIME_STEP_LENGTH / DENSITY
        *
        d_p_next__d_y
    )

    # Velocity Boundary Conditions: Homogeneous Dirichlet BC everywhere
    # except for the horizontal velocity at the top, which is prescribed
    u_next[0, :] = 0.0
    u_next[:, 0] = 0.0
    u_next[:, -1] = 0.0
    u_next[-1, :] = HORIZONTAL_VELOCITY_TOP
    v_next[0, :] = 0.0
    v_next[:, 0] = 0.0
    v_next[:, -1] = 0.0
    v_next[-1, :] = 0.0


    # Advance in time
    u_prev = u_next
    v_prev = v_next
    p_prev = p_next


# The [::2, ::2] selects only every second entry (less cluttering plot)
plt.figure()
plt.contourf(X[::2, ::2], Y[::2, ::2], p_next[::2, ::2], cmap="coolwarm")
plt.colorbar()

plt.quiver(X[::2, ::2], Y[::2, ::2], u_next[::2, ::2], v_next[::2, ::2])
# plt.streamplot(X[::2, ::2], Y[::2, ::2], u_next[::2, ::2], v_next[::2, ::2], color="black")
plt.xlim((0, 1))
plt.ylim((0, 1))
plt.show()

100%|███████████████████████████████████████| 500/500 [00:00<00:00, 1392.49it/s]