Show code cell source Hide 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 Hide 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 Hide 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:

\[ \omega = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}. \]

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:

\[ b_j(t) = e^{\omega_j t} b_j, \]

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:

\[ \mathbf{x}(t) = \sum_{j=1}^{r} b_j e^{\omega_j t} \boldsymbol{\phi}_j. \]

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}")