Irrotational random flow¶
Last updated: March 3, 2022
Setup
In [1]:
Copied!
import numpy as np
import skimage.filters
import scipy.ndimage
import scipy.sparse
# Plotting
import matplotlib as mpl
import matplotlib.pyplot as plt
import proplot as pplt
# Lima
import lima
import numpy as np
import skimage.filters
import scipy.ndimage
import scipy.sparse
# Plotting
import matplotlib as mpl
import matplotlib.pyplot as plt
import proplot as pplt
# Lima
import lima
In [2]:
Copied!
lima.plot.init_mpl_style()
lima.plot.init_mpl_style()
In [3]:
Copied!
debug = True
debug = True
In [4]:
Copied!
subplotgrid = [
[1, 1, 2, 2, 3, 3],
[0, 4, 4, 5, 5, 0],
]
subplotgrid = [
[1, 1, 2, 2, 3, 3],
[0, 4, 4, 5, 5, 0],
]
Functions
In [5]:
Copied!
def vortZ(u, v):
dudy, dudx = np.gradient(u)
dvdy, dvdx = np.gradient(v)
return dvdx - dudy
def divUV(u, v):
dudy, dudx = np.gradient(u)
dvdy, dvdx = np.gradient(v)
return dudx + dvdy
def random_dataset(vscale=1.0, height=256, width=256, filter_size=5, seed=None):
if seed:
np.random.seed(seed)
# Coordinates
x, y = np.meshgrid(np.arange(height, dtype='float32'),
np.arange(width, dtype='float32'))
#u, v = np.random.rand(2, height, width).astype('float32') - 0.5
u, v = np.random.uniform(-1, 1, size=(2, height, width)).astype('float32')
u = skimage.filters.gaussian(u, filter_size, mode='reflect')
v = skimage.filters.gaussian(v, filter_size, mode='reflect')
u = u/np.abs(u).max() * vscale
v = v/np.abs(v).max() * vscale
# Parameters
#print(np.max(u), np.max(v))
return x, y, u, v
def quiver(ax, x, y, u, v, skip=1, *args, **kwargs):
return ax.quiver(
x[::skip, ::skip], y[::skip, ::skip], u[::skip, ::skip], v[::skip, ::skip],
*args, **kwargs
)
norm = lambda u,v: np.sqrt(u**2 + v**2)
def vortZ(u, v):
dudy, dudx = np.gradient(u)
dvdy, dvdx = np.gradient(v)
return dvdx - dudy
def divUV(u, v):
dudy, dudx = np.gradient(u)
dvdy, dvdx = np.gradient(v)
return dudx + dvdy
def random_dataset(vscale=1.0, height=256, width=256, filter_size=5, seed=None):
if seed:
np.random.seed(seed)
# Coordinates
x, y = np.meshgrid(np.arange(height, dtype='float32'),
np.arange(width, dtype='float32'))
#u, v = np.random.rand(2, height, width).astype('float32') - 0.5
u, v = np.random.uniform(-1, 1, size=(2, height, width)).astype('float32')
u = skimage.filters.gaussian(u, filter_size, mode='reflect')
v = skimage.filters.gaussian(v, filter_size, mode='reflect')
u = u/np.abs(u).max() * vscale
v = v/np.abs(v).max() * vscale
# Parameters
#print(np.max(u), np.max(v))
return x, y, u, v
def quiver(ax, x, y, u, v, skip=1, *args, **kwargs):
return ax.quiver(
x[::skip, ::skip], y[::skip, ::skip], u[::skip, ::skip], v[::skip, ::skip],
*args, **kwargs
)
norm = lambda u,v: np.sqrt(u**2 + v**2)
In [6]:
Copied!
def construct_LHS(nx, ny, h=1):
n = nx * ny
B = np.zeros((5, n))
for i in range(ny):
for j in range(nx):
k = j*nx + i
if i == 0:
if j == 0:
B[2, k] = 2*h/(h + h) + 2*h/(h + h)
B[3, k+1] = -2*h/(h + h)
B[4, k+nx] = -2*h/(h + h)
elif j == nx-1:
B[0, k-nx] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h)
B[3, k+1] = -2*h/(h + h)
else:
B[0, k-nx] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h) + 2*h/(h + h)
B[3, k+1] = -2*h/(h + h)
B[4, k+nx] = -2*h/(h + h)
elif i == nx-1:
if j == 0:
B[1, k-1] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h)
B[4, k+nx] = -2*h/(h + h)
elif j == nx-1:
B[0, k-nx] = -2*h/(h + h)
B[1, k-1] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h)
else:
B[0, k-nx] = -2*h/(h + h)
B[1, k-1] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h) + 2*h/(h + h)
B[4, k+nx] = -2*h/(h + h)
elif j == 0:
if ( i > 0 and i < nx-1 ):
B[1, k-1] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h) + 2*h/(h + h)
B[3, k+1] = -2*h/(h + h)
B[4, k+nx] = -2*h/(h + h)
elif j == nx-1:
if ( i > 0 and i < nx-1 ):
B[0, k-nx] = -2*h/(h + h)
B[1, k-1] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h) + 2*h/(h + h)
B[3, k+1] = -2*h/(h + h)
else:
B[0, k-nx] = -2*h/(h + h)
B[1, k-1] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h) + 2*h/(h + h) + 2*h/(h + h)
B[3, k+1] = -2*h/(h + h)
B[4, k+nx] = -2*h/(h + h)
# Diagonal indices
diags = np.array([-nx, -1, 0, 1, nx])
# Construct sparse diagonal matrix
L = scipy.sparse.spdiags(B, diags, n, n).tocsr()
return L
def construct_RHS(u, v):
f = -divUV(u, v)
return f.ravel()
def solve_poisson_equation(u, v, L=None, debug=False, maxiter=100):
if L is None:
L = construct_LHS(u.shape[0], u.shape[1])
f = construct_RHS(u, v)
# Solve using gmres
phi, info = scipy.sparse.linalg.gmres(L, f, maxiter=maxiter)
if debug:
print(f"Residual: {np.sum(np.abs(L.dot(phi) - f))}")
phi = phi.reshape(u.shape[0], u.shape[1])
return phi
def calculate_irrotational_flow(u, v, L=None, debug=False, maxiter=100):
# Solve poisson equation: $\nabla^2 \phi = -\nabla \cdot u$
phi = solve_poisson_equation(u, v, L, debug, maxiter)
# Calculate irrotation flow
v_phi, u_phi = np.gradient(phi)
return u_phi, v_phi
def construct_LHS(nx, ny, h=1):
n = nx * ny
B = np.zeros((5, n))
for i in range(ny):
for j in range(nx):
k = j*nx + i
if i == 0:
if j == 0:
B[2, k] = 2*h/(h + h) + 2*h/(h + h)
B[3, k+1] = -2*h/(h + h)
B[4, k+nx] = -2*h/(h + h)
elif j == nx-1:
B[0, k-nx] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h)
B[3, k+1] = -2*h/(h + h)
else:
B[0, k-nx] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h) + 2*h/(h + h)
B[3, k+1] = -2*h/(h + h)
B[4, k+nx] = -2*h/(h + h)
elif i == nx-1:
if j == 0:
B[1, k-1] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h)
B[4, k+nx] = -2*h/(h + h)
elif j == nx-1:
B[0, k-nx] = -2*h/(h + h)
B[1, k-1] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h)
else:
B[0, k-nx] = -2*h/(h + h)
B[1, k-1] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h) + 2*h/(h + h)
B[4, k+nx] = -2*h/(h + h)
elif j == 0:
if ( i > 0 and i < nx-1 ):
B[1, k-1] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h) + 2*h/(h + h)
B[3, k+1] = -2*h/(h + h)
B[4, k+nx] = -2*h/(h + h)
elif j == nx-1:
if ( i > 0 and i < nx-1 ):
B[0, k-nx] = -2*h/(h + h)
B[1, k-1] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h) + 2*h/(h + h)
B[3, k+1] = -2*h/(h + h)
else:
B[0, k-nx] = -2*h/(h + h)
B[1, k-1] = -2*h/(h + h)
B[2, k] = 2*h/(h + h) + 2*h/(h + h) + 2*h/(h + h) + 2*h/(h + h)
B[3, k+1] = -2*h/(h + h)
B[4, k+nx] = -2*h/(h + h)
# Diagonal indices
diags = np.array([-nx, -1, 0, 1, nx])
# Construct sparse diagonal matrix
L = scipy.sparse.spdiags(B, diags, n, n).tocsr()
return L
def construct_RHS(u, v):
f = -divUV(u, v)
return f.ravel()
def solve_poisson_equation(u, v, L=None, debug=False, maxiter=100):
if L is None:
L = construct_LHS(u.shape[0], u.shape[1])
f = construct_RHS(u, v)
# Solve using gmres
phi, info = scipy.sparse.linalg.gmres(L, f, maxiter=maxiter)
if debug:
print(f"Residual: {np.sum(np.abs(L.dot(phi) - f))}")
phi = phi.reshape(u.shape[0], u.shape[1])
return phi
def calculate_irrotational_flow(u, v, L=None, debug=False, maxiter=100):
# Solve poisson equation: $\nabla^2 \phi = -\nabla \cdot u$
phi = solve_poisson_equation(u, v, L, debug, maxiter)
# Calculate irrotation flow
v_phi, u_phi = np.gradient(phi)
return u_phi, v_phi
Pure-random flow¶
In [7]:
Copied!
filter_size = 5
vscale = 1
height = width = 64
filter_size = 5
vscale = 1
height = width = 64
In [8]:
Copied!
# Random dist.
np.random.seed(234)
x, y, _,_ = random_dataset(vscale=1, height=height, width=width, filter_size=filter_size)
u, v = np.random.uniform(-1, 1, size=(2, height, width)).astype('float32')
# Filtered
uf = skimage.filters.gaussian(u, filter_size, mode="reflect")#, preserve_range=True)
vf = skimage.filters.gaussian(v, filter_size, mode="reflect")#, preserve_range=True)
uf = uf/np.abs(uf).max() * vscale
vf = vf/np.abs(vf).max() * vscale
# Scaled
us = uf
vs = vf
# Random dist.
np.random.seed(234)
x, y, _,_ = random_dataset(vscale=1, height=height, width=width, filter_size=filter_size)
u, v = np.random.uniform(-1, 1, size=(2, height, width)).astype('float32')
# Filtered
uf = skimage.filters.gaussian(u, filter_size, mode="reflect")#, preserve_range=True)
vf = skimage.filters.gaussian(v, filter_size, mode="reflect")#, preserve_range=True)
uf = uf/np.abs(uf).max() * vscale
vf = vf/np.abs(vf).max() * vscale
# Scaled
us = uf
vs = vf
Rotational and irrotational component¶
In [9]:
Copied!
fig = pplt.figure(refwidth='1.8', span=False)
axes = fig.subplots(subplotgrid)
s = 5
# u
ax = axes[0]
im = ax.imshow(us, interpolation='none', cmap='turbo')
cb = ax.colorbar(im, width='0.75em')
ax.format(title=r"U-Velocity: $u$")
# v
ax = axes[1]
im = ax.imshow(vs, interpolation='none', cmap='turbo')
cb = ax.colorbar(im, width='0.75em')
ax.format(title=r"V-Velocity: $v$")
# |u|
ax = axes[2]
im = ax.imshow((us**2 + vs**2)**0.5, interpolation='none', cmap='turbo')
cb = ax.colorbar(im, width='0.75em')
quiver(ax, x, y, us, vs, skip=2, color='k', scale=s)
ax.format(title=r"Velocity mag.: $|u|$")
# vorticity |u|
ax = axes[3]
im = ax.imshow(np.abs(vortZ(us, vs)), interpolation='none', cmap='turbo')
cb = ax.colorbar(im, width='0.75em')
ax.format(title=fr"Vorticity mag.: $|\nabla \times u|$ (sum={np.abs(vortZ(us, vs)).sum():2g})")
ax.streamplot(x, y, us, vs, c='k', lw=0.5)
ax.axis([0, height, 0, width])
# divergence |u|
ax = axes[4]
im = ax.imshow(divUV(us, vs), interpolation='none', cmap='turbo')
cb = ax.colorbar(im, width='0.75em')
ax.format(title=r"Divergence: $\nabla \cdot u$")
# Quiver
# # Format
axes.format(
abc='(a)', abcloc='ul', abcbbox=True,
xlabel='$x$ (px)', ylabel='$y$ (px)',
suptitle='Pure random vector field',
)
fig = pplt.figure(refwidth='1.8', span=False)
axes = fig.subplots(subplotgrid)
s = 5
# u
ax = axes[0]
im = ax.imshow(us, interpolation='none', cmap='turbo')
cb = ax.colorbar(im, width='0.75em')
ax.format(title=r"U-Velocity: $u$")
# v
ax = axes[1]
im = ax.imshow(vs, interpolation='none', cmap='turbo')
cb = ax.colorbar(im, width='0.75em')
ax.format(title=r"V-Velocity: $v$")
# |u|
ax = axes[2]
im = ax.imshow((us**2 + vs**2)**0.5, interpolation='none', cmap='turbo')
cb = ax.colorbar(im, width='0.75em')
quiver(ax, x, y, us, vs, skip=2, color='k', scale=s)
ax.format(title=r"Velocity mag.: $|u|$")
# vorticity |u|
ax = axes[3]
im = ax.imshow(np.abs(vortZ(us, vs)), interpolation='none', cmap='turbo')
cb = ax.colorbar(im, width='0.75em')
ax.format(title=fr"Vorticity mag.: $|\nabla \times u|$ (sum={np.abs(vortZ(us, vs)).sum():2g})")
ax.streamplot(x, y, us, vs, c='k', lw=0.5)
ax.axis([0, height, 0, width])
# divergence |u|
ax = axes[4]
im = ax.imshow(divUV(us, vs), interpolation='none', cmap='turbo')
cb = ax.colorbar(im, width='0.75em')
ax.format(title=r"Divergence: $\nabla \cdot u$")
# Quiver
# # Format
axes.format(
abc='(a)', abcloc='ul', abcbbox=True,
xlabel='$x$ (px)', ylabel='$y$ (px)',
suptitle='Pure random vector field',
)