Show code cell source
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('svg')
import seaborn as sns
sns.set_context("paper")
sns.set_style("ticks");
DMD Example#
We are going to replicate the DMD example shown in Fig. 7.3 of the book. The example is flow past a circular cylinder at Reynolds number 100. The flow is simulated using the incompressible Navier-Stokes equations. The data are available to download from DMD Book.
!curl -O http://dmdbook.com/DATA.zip
Show code cell output
% Total % Received % Xferd Average Speed Time Time Time Current
Dload Upload Total Spent Left Speed
100 417M 100 417M 0 0 3703k 0 0:01:54 0:01:54 --:--:-- 1603k59k 0 0:01:03 0:00:12 0:00:51 4013k 0 0:01:10 0:00:16 0:00:54 3733k4509k 0 0:01:34 0:00:34 0:01:00 2502k:51 0:00:46 5255k00:44 4963k1:38 0:00:56 0:00:42 3471k02 0:00:39 3396k:00:35 3989k 0 0:01:47 0:01:24 0:00:23 2813k 0 0 3851k 0 0:01:50 0:01:37 0:00:13 3061k 0 0:01:55 0:01:55 --:--:-- 1720k
Unzip:
!unzip -o DATA.zip
Show code cell output
Archive: DATA.zip
creating: DATA/
inflating: DATA/.DS_Store
creating: __MACOSX/
creating: __MACOSX/DATA/
inflating: __MACOSX/DATA/._.DS_Store
creating: DATA/FLUIDS/
inflating: DATA/FLUIDS/.DS_Store
creating: __MACOSX/DATA/FLUIDS/
inflating: __MACOSX/DATA/FLUIDS/._.DS_Store
inflating: DATA/FLUIDS/CYLINDER_ALL.mat
inflating: __MACOSX/DATA/FLUIDS/._CYLINDER_ALL.mat
inflating: DATA/FLUIDS/CYLINDER_basis.mat
inflating: __MACOSX/DATA/FLUIDS/._CYLINDER_basis.mat
creating: DATA/FLUIDS/ibpm_scripts/
inflating: DATA/FLUIDS/ibpm_scripts/.DS_Store
creating: __MACOSX/DATA/FLUIDS/ibpm_scripts/
inflating: __MACOSX/DATA/FLUIDS/ibpm_scripts/._.DS_Store
inflating: DATA/FLUIDS/ibpm_scripts/cylinder.geom
inflating: __MACOSX/DATA/FLUIDS/ibpm_scripts/._cylinder.geom
inflating: DATA/FLUIDS/ibpm_scripts/cylinder01.sh
inflating: __MACOSX/DATA/FLUIDS/ibpm_scripts/._cylinder01.sh
inflating: DATA/FLUIDS/ibpm_scripts/cylinder02.sh
inflating: __MACOSX/DATA/FLUIDS/ibpm_scripts/._cylinder02.sh
inflating: DATA/FLUIDS/ibpm_scripts/cylinder03.sh
inflating: __MACOSX/DATA/FLUIDS/ibpm_scripts/._cylinder03.sh
inflating: DATA/FLUIDS/ibpm_scripts/cylinder04.sh
inflating: __MACOSX/DATA/FLUIDS/ibpm_scripts/._cylinder04.sh
inflating: DATA/FLUIDS/ibpm_scripts/cylinder05.sh
inflating: __MACOSX/DATA/FLUIDS/ibpm_scripts/._cylinder05.sh
inflating: __MACOSX/DATA/FLUIDS/._ibpm_scripts
inflating: DATA/FLUIDS/ibpmUNSTEADY.plt
inflating: __MACOSX/DATA/FLUIDS/._ibpmUNSTEADY.plt
inflating: DATA/FLUIDS/README.txt
inflating: __MACOSX/DATA/FLUIDS/._README.txt
creating: DATA/NEURO/
inflating: DATA/NEURO/.DS_Store
creating: __MACOSX/DATA/NEURO/
inflating: __MACOSX/DATA/NEURO/._.DS_Store
inflating: DATA/NEURO/demo_ecog.mat
inflating: __MACOSX/DATA/NEURO/._demo_ecog.mat
inflating: DATA/NEURO/ecog_window.mat
inflating: __MACOSX/DATA/NEURO/._ecog_window.mat
inflating: __MACOSX/._DATA
Let’s load the data:
import scipy
data = scipy.io.loadmat('DATA/FLUIDS/CYLINDER_ALL.mat')
Let’s see what is inside:
data.keys()
dict_keys(['__header__', '__version__', '__globals__', 'UALL', 'UEXTRA', 'VALL', 'VEXTRA', 'VORTALL', 'VORTEXTRA', 'm', 'n', 'nx', 'ny'])
m = data['m'][0][0]
n = data['n'][0][0]
nx = data['nx'][0][0]
ny = data['ny'][0][0]
u_all = data['UALL']
u_extra = data['UEXTRA']
v_all = data['VALL']
v_extra = data['VEXTRA']
vort_all = data['VORTALL']
vort_extra = data['VORTEXTRA']
print(f"m\t\t= {m}")
print(f"n\t\t= {n}")
print(f"nx\t\t= {nx}")
print(f"ny\t\t= {ny}")
print(f"u_all\t\t= {u_all.shape}")
print(f"u_extra\t\t= {u_extra.shape}")
print(f"v_all\t\t= {v_all.shape}")
print(f"v_extra\t\t= {v_extra.shape}")
print(f"vort_all\t= {vort_all.shape}")
print(f"vort_extra\t= {vort_extra.shape}")
m = 199
n = 449
nx = 199
ny = 449
u_all = (89351, 151)
u_extra = (89351, 1)
v_all = (89351, 151)
v_extra = (89351, 1)
vort_all = (89351, 151)
vort_extra = (89351, 1)
We have 151 snapshots of the flow field. Each snapshot is \(449 \times 199\). We have \(x\) and \(y\) velocities as well as the voriticity field. We are going to focus on the vorticity field only. In 2D the vorticity field is defined as:
So, the vorticity is just a scalar.
where \(u\) and \(v\) are the \(x\) and \(y\) velocities, respectively. Working with the vorticity field we automatically satisfy the incompressibility condition \(\nabla \cdot \mathbf{u} = 0\).
Let’s visualize the vorticity at a fiew snapshots. Because the vorticity on the cylinder doesn’t make sense, we are going to mask it out by setting it to zero.
vorticity_max = 10.0
vort_all = np.clip(vort_all, -vorticity_max, vorticity_max)
for i in [0, 50, 100]:
fig, ax = plt.subplots(figsize=(5, 1.8))
cylinder = plt.Circle((50, 99), 27, color='gray')
c = ax.contourf(vort_all[:,i].reshape((nx, ny), order='F'), cmap='jet',
levels=np.linspace(-5, 5, 100))
#c = ax.contourf(u_all[:,i].reshape((nx, ny), order='F'), cmap='jet',
# levels=np.linspace(-2, 2, 100))
# plot a gray disk to represent the cylinder
ax.add_patch(cylinder)
plt.colorbar(c)
plt.title(f"Vorticity at time step {i}")
plt.show()
Let’s form the required matrices:
n_steps = 140
X = vort_all[:, 0:n_steps]
X_prime = vort_all[:, 1:n_steps+1]
Let’s write code that finds the DMD modes and eigenvalues.
import numpy as np
def dmd(X, Xprime, r):
"""Compute DMD.
Parameters
----------
X : ndarray of shape (n, m)
Data matrix.
Xprime : ndarray of shape (n, m)
Data matrix.
r : int
Rank truncation.
Returns
-------
Phi : ndarray of shape (n, r)
DMD modes.
Lambda : ndarray of shape (r, r)
Eigenvalues.
b : ndarray of shape (r,)
Initial amplitude.
"""
# Step 1 - SVD of X
U, s, VT = np.linalg.svd(X, full_matrices=False)
Ur = U[:, :r]
Sr = np.diag(s[:r])
VTr = VT[:r, :]
# Step 2 - compute Atilde
Atilde = np.linalg.solve(Sr.T, (Ur.T @ Xprime @ VTr.T).T).T
# Step 3 - Eigenvalues and eigenvectors of Atilde
Lambda, W = np.linalg.eig(Atilde)
Lambda = np.diag(Lambda)
# Step 4 - Compute first r DMD modes
Phi = Xprime @ np.linalg.solve(Sr.T, VTr).T @ W
# Step 5 - Compute initial amplitude
x1tilde = Sr @ VTr[:, 0]
b = np.linalg.solve(W @ Lambda, x1tilde)
return Phi, Lambda, b
Let’s try it out on the vorticity:
Phi, Lambda, b = dmd(X, X_prime, 50)
We want to sort the modes according to their original amplitudes.
idx = np.argsort(np.abs(b))[::-1]
b = b[idx]
Phi = Phi[:, idx]
Lambda = Lambda[idx][:, idx]
Let’s start first with the DMD modes. Note that the DMD modes are imaginary. We are going to plot their magnitude.
for j in range(4):
fig, ax = plt.subplots(figsize=(5, 1.8))
cylinder = plt.Circle((52, 100), 24, color='gray')
c = ax.contourf(np.real(Phi[:, j]).reshape((nx, ny), order='F'),
cmap='jet')
ax.add_patch(cylinder)
ax.set(title=f"DMD mode {j} (real part)")
fig, ax = plt.subplots(figsize=(5, 1.8))
cylinder = plt.Circle((52, 100), 24, color='gray')
c = ax.contourf(np.imag(Phi[:, j]).reshape((nx, ny), order='F'),
cmap='jet')
ax.add_patch(cylinder)
ax.set(title=f"DMD mode {j} (imaginary part)")
plt.show()
Now, let’s look at the eigenvalues:
fig, ax = plt.subplots()
ax.plot(Lambda.real, Lambda.imag, 'xr')
# plot the unit circle
theta = np.linspace(0, 2*np.pi, 100)
ax.plot(np.cos(theta), np.sin(theta), 'k-')
# eigenvalues with magnitude less than 1
ax.set(xlabel='Real', ylabel='Imaginary', title='Eigenvalues')
# set aspect ratio to be equal
ax.set_aspect('equal')
sns.despine(trim=True)
The eigenvalues that are on the unit circle correspond to oscillatory modes. The eigenvalues that are inside the unit circle correspond to decaying modes.
Another thing that we can do is look at the evolution of the DMD amplitudes over time. Recall, that \(\mathbf{b}\) is the DMD amplitude initially. It evolves as:
where \(\omega_j = \log(\lambda_j) / \Delta t\).
Dt = 1.0
Omega = np.log(np.diag(Lambda)) / Dt
import jax.numpy as jnp
from jax import vmap, jit
predict_b = jit(vmap(lambda t: jnp.exp(Omega * t) * b))
ts = np.linspace(0, 200, 500)
bs = predict_b(ts)
fig, ax = plt.subplots()
ax.plot(ts, bs.real, 'r', lw=0.5)
ax.set(xlabel='Time', ylabel='Real part of $b_j$', title='DMD prediction')
sns.despine(trim=True);
And for each one of these times, you can reconstruct the flow field as:
vort_pred = np.einsum('ij,tj->it', Phi, bs).real
Here is the comparison at a few snapshots. Note that the data after time step \(100\) have not been used in the DMD calculation. So, we are just doing pure predictions.
for i in [10, 50, 100, 140, 145, 150]:
fig, ax = plt.subplots(figsize=(5, 1.8))
cylinder = plt.Circle((50, 99), 27, color='gray')
c = ax.contourf(vort_all[:,i].reshape((nx, ny), order='F'), cmap='jet',
levels=np.linspace(-5, 5, 100))
ax.add_patch(cylinder)
plt.colorbar(c)
plt.title(f"Vorticity at time step {i}")
fig, ax = plt.subplots(figsize=(5, 1.8))
cylinder = plt.Circle((50, 99), 27, color='gray')
c = ax.contourf(vort_pred[:,i].reshape((nx, ny), order='F'), cmap='jet',
levels=np.linspace(-5, 5, 100))
ax.add_patch(cylinder)
plt.colorbar(c)
plt.title(f"Predicted vorticity at time step {i}")
