In [1]:
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
plt.style.use('bmh')
In [2]:
# Implementation of Eq. (1) in the exam set
def doublegyre(x, y, t, A, e, w):
a = e * np.sin(w*t)
b = 1 - 2*e*np.sin(w*t)
f = a*x**2 + b*x
return np.array([
-np.pi*A*np.sin(np.pi*f) * np.cos(np.pi*y), # x component of velocity
np.pi*A*np.cos(np.pi*f) * np.sin(np.pi*y) * (2*a*x + b) # y component of velocity
])
# Wrapper function to pass to integrator
# X0 is a two-component vector [x, y]
def f(X, t):
# Parameters of the velocity field
A = 0.1
e = 0.25 # epsilon
w = 1 # omega
return doublegyre(X[0], X[1], t, A, e, w)
# Forward Euler integrator
# X0 is a two-component vector [x, y]
def euler(X, t, dt, f):
k1 = f(X, t)
return X + dt*k1
# 4th order Runge-Kutta integrator
# X0 is a two-component vector [x, y]
def rk4(X, t, dt, f):
k1 = f(X, t)
k2 = f(X + k1*dt/2, t + dt/2)
k3 = f(X + k2*dt/2, t + dt/2)
k4 = f(X + k3*dt, t + dt)
return X + dt*(k1 + 2*k2 + 2*k3 + k4) / 6
# Function to calculate a trajectory from an
# initial position X0 at t = 0, moving forward
# until t = tmax, using the given timestep and
# integrator
def trajectory(X0, tmax, dt, integrator, f):
t = 0
# Number of timesteps
Nx = int(tmax / dt)
# Array to hold the entire trajectory
PX = np.zeros((2, Nx+1))
# Initial position
PX[:,0] = X0
# Loop over all timesteps
for i in range(1, Nx+1):
PX[:,i] = integrator(PX[:,i-1], t, dt, f)
t += dt
# Return entire trajectory
return PX
In [3]:
# Create a grid of points
x = np.linspace(0, 2, 50)
y = np.linspace(0, 1, 25)
X, Y = np.meshgrid(x, y)
# Time and parameters of current field
t = 0.0
A = 0.1
e = 0.25
w = 1
# Calculate field
V = doublegyre(X, Y, t, A, e, w)
# Vector plot
fig = plt.figure(figsize = (12,6))
plt.quiver(X, Y, V[0], V[1], linewidths = -1.0, scale = 10, alpha = 0.6)
# Add text showing time, and set plot limits
plt.text(1.65, 0.9, '$t = %s$' % t, size = 36)
plt.xlim(0,2)
plt.ylim(0,1)
Out[3]:
In [4]:
# Using forward Euler
X0 = [1.12, 0.6]
tmax = 50
# Set figure size
fig = plt.figure(figsize = (12,6))
# Plot trajectory for dt = 0.1
dt = 0.1
T = trajectory(X0, tmax, dt, euler, f)
plt.plot(T[0,:], T[1,:], label = 'dt = %s' % dt)
# Plot trajectory for dt = 0.01
dt = 0.01
T = trajectory(X0, tmax, dt, euler, f)
plt.plot(T[0,:], T[1,:], '--', label = 'dt = %s' % dt)
# Add legend and set limits
plt.legend()
plt.xlim(0,2)
plt.ylim(0,1)
Out[4]:
In [5]:
# Using Runge-Kutta
X0 = [1.12, 0.6]
tmax = 50
# Set figure size
fig = plt.figure(figsize = (12,6))
# Plot trajectory for dt = 0.1
dt = 0.1
T = trajectory(X0, tmax, dt, rk4, f)
plt.plot(T[0,:], T[1,:], label = 'dt = %s' % dt)
# Plot trajectory for dt = 0.01
dt = 0.01
T = trajectory(X0, tmax, dt, rk4, f)
plt.plot(T[0,:], T[1,:], '--', label = 'dt = %s' % dt)
# Add legend and set limits
plt.legend()
plt.xlim(0,2)
plt.ylim(0,1)
Out[5]:
In [10]:
#manual trajectory
X0 = [1.12, 0.6]
tmax = 1
# Plot trajectory for dt = 0.1
dt = 0.1
#T = trajectory(X0, tmax, dt, rk4, f)
#def trajectory(X0, tmax, dt, integrator, f):
t = 0
# Number of timesteps
Nx = int(tmax / dt)
# Array to hold the entire trajectory
PX = np.zeros((2, Nx+1))
# Initial position
PX[:,0] = X0
In [11]:
PX
Out[11]:
In [12]:
# Loop over all timesteps
for i in range(1, Nx+1):
PX[:,i] = rk4(PX[:,i-1], t, dt, f)
t += dt
# Return entire trajectory
#return PX
In [13]:
PX
Out[13]:
In [14]:
# Set figure size
fig = plt.figure(figsize = (12,6))
plt.plot(PX[0,:], PX[1,:], '--', label = 'dt = %s' % dt)
# Add legend and set limits
plt.legend()
plt.xlim(0,2)
plt.ylim(0,1)
Out[14]:
In [48]:
# Specify resolution
Nx = 10
Ny = 5
n = Nx*Ny
# Specify limit of particles
particles_x_min = 0.9
particles_x_max = 1.1
particles_y_min = 0.4
particles_y_max = 0.6
# Create lists of all x and y positions
X = np.linspace(particles_x_min, particles_x_max, Nx)
Y = np.linspace(particles_y_min, particles_y_max, Ny)
In [49]:
# Array to hold all grid points
X0 = np.zeros((2, Ny, Nx))
# Use meshgrid to turn lists into rank 2 arrays
# of x and y positions
X0[:] = np.meshgrid(X, Y)
# Make scatter plot to show all grid points
fig = plt.figure(figsize = (12,6))
plt.scatter(X0[0,:], X0[1,:], marker = '.', c = '#348ABD')
plt.xlim(0,2)
plt.ylim(0,1)
Out[49]:
In [50]:
# Array to hold all grid points after transport
X1 = np.zeros((2, Ny, Nx))
# Transport parameters
tmax = 2.0
dt = 0.5
# Loop over grid and update all positions
# This is where parallelisation would happen, since
# each position is independent of all the others
for i in range(Nx):
for j in range(Ny):
# Keep only the last position, not the entire trajectory
X1[:,j,i] = trajectory(X0[:,j,i], tmax, dt, rk4, f)[:,-1]
# Make scatter plot to show all grid points
fig = plt.figure(figsize = (12,6))
plt.scatter(X1[0,:], X1[1,:], marker = '.', c = '#348ABD')
plt.xlim(0,2)
plt.ylim(0,1)
Out[50]:
In [51]:
### test of coordinate arrays instead:
#particle_x:
x = np.hstack(X0[0][:])
x
Out[51]:
In [52]:
# particle_y
y = np.hstack(X0[1][:])
y
Out[52]:
In [64]:
xy1 = np.zeros((2, n))
xy1
Out[64]:
In [66]:
print(X0)
In [68]:
m = [x,y]
print(m)
In [ ]:
In [30]:
# Array to hold all grid points after transport
x1 = np.zeros(n)
y1 = np.zeros(n)
xy1 = np.zeros((2, n))
# Transport parameters
tmax = 2.0
dt = 0.5
# Loop over grid and update all positions
# This is where parallelisation would happen, since
# each position is independent of all the others
for i in range(n):
# Keep only the last position, not the entire trajectory
xy1[:,j,i] = trajectory(X0[:,j,i], tmax, dt, rk4, f)[:,-1]
# Make scatter plot to show all grid points
fig = plt.figure(figsize = (12,6))
plt.scatter(X1[0,:], X1[1,:], marker = '.', c = '#348ABD')
plt.xlim(0,2)
plt.ylim(0,1)
Out[30]:
In [31]:
Out[31]:
In [33]:
# Array to hold all grid points after transport
X1 = np.zeros((2, Ny, Nx))
# Transport parameters
tmax = 4.0
dt = 0.5
# Loop over grid and update all positions
# This is where parallelisation would happen, since
# each position is independent of all the others
for i in range(Nx):
for j in range(Ny):
# Keep only the last position, not the entire trajectory
X1[:,j,i] = trajectory(X0[:,j,i], tmax, dt, rk4, f)[:,-1]
# Make scatter plot to show all grid points
fig = plt.figure(figsize = (12,6))
plt.scatter(X1[0,:], X1[1,:], marker = '.', c = '#348ABD')
plt.xlim(0,2)
plt.ylim(0,1)
Out[33]:
In [34]:
# Array to hold all grid points after transport
X1 = np.zeros((2, Ny, Nx))
# Transport parameters
tmax = 8.0
dt = 0.5
# Loop over grid and update all positions
# This is where parallelisation would happen, since
# each position is independent of all the others
for i in range(Nx):
for j in range(Ny):
# Keep only the last position, not the entire trajectory
X1[:,j,i] = trajectory(X0[:,j,i], tmax, dt, rk4, f)[:,-1]
# Make scatter plot to show all grid points
fig = plt.figure(figsize = (12,6))
plt.scatter(X1[0,:], X1[1,:], marker = '.', c = '#348ABD')
plt.xlim(0,2)
plt.ylim(0,1)
Out[34]:
In [35]:
# Array to hold all grid points after transport
X1 = np.zeros((2, Ny, Nx))
# Transport parameters
tmax = 16.0
dt = 0.5
# Loop over grid and update all positions
# This is where parallelisation would happen, since
# each position is independent of all the others
for i in range(Nx):
for j in range(Ny):
# Keep only the last position, not the entire trajectory
X1[:,j,i] = trajectory(X0[:,j,i], tmax, dt, rk4, f)[:,-1]
# Make scatter plot to show all grid points
fig = plt.figure(figsize = (12,6))
plt.scatter(X1[0,:], X1[1,:], marker = '.', c = '#348ABD')
plt.xlim(0,2)
plt.ylim(0,1)
Out[35]:
In [36]:
# Array to hold all grid points after transport
X1 = np.zeros((2, Ny, Nx))
# Transport parameters
tmax = 32.0
dt = 0.5
# Loop over grid and update all positions
# This is where parallelisation would happen, since
# each position is independent of all the others
for i in range(Nx):
for j in range(Ny):
# Keep only the last position, not the entire trajectory
X1[:,j,i] = trajectory(X0[:,j,i], tmax, dt, rk4, f)[:,-1]
# Make scatter plot to show all grid points
fig = plt.figure(figsize = (12,6))
plt.scatter(X1[0,:], X1[1,:], marker = '.', c = '#348ABD')
plt.xlim(0,2)
plt.ylim(0,1)
Out[36]:
In [ ]:
# TODO:
# calculate the deformation for each of these points, and make a 2D plot showing log(γ − 1) as a funct
#γ =
# a/b if a ≥ b
# b/a if a < b
In [9]:
X1[0,:]
Out[9]:
In [10]:
X1
Out[10]:
In [12]:
np.zeros((2, 10, 5))
Out[12]:
In [17]:
np.meshgrid(np.linspace(0, 0.9, 10), np.linspace(0, 0.8, 5))
Out[17]: