<jupyter_start><jupyter_text>Assignment 2.1 : 2D Heat Equation (Euler Explicit)<jupyter_code>import numpy as np
import scipy.linalg as la
from scipy import stats
from scipy.interpolate import interp1d
import scipy
import time
import scipy.io
from scipy.sparse import diags
from matplotlib import cm
from matplotlib.ticker import LinearLocator
import copy
import cv2
from matplotlib.patches import Circle


import matplotlib.pyplot as plt
import matplotlib
font = {'family' : 'serif',
        'size'   : 13}

matplotlib.rc('font', **font)
matplotlib.rc('xtick', labelsize=15) 
matplotlib.rc('ytick', labelsize=15)<jupyter_output><empty_output><jupyter_text>$$\frac{\partial T}{\partial t} = \alpha \left(\frac{\partial^2T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2}\right)$$where $\alpha$ is the thermal diffusivity. This can be solved numerically using the Forward in Time Cenetred in Space (FTCS) method. The stability requirement is $\alpha \frac{\Delta t}{\Delta x^2} + \alpha \frac{\Delta t}{\Delta y^2} \leq \frac{1}{2}$. $$T_{i, j}^{n+1} = T_{i, j}^{n} + \alpha \frac{\Delta t}{\Delta x^2} \left(T_{i-1, j}^{n} - 2T_{i, j}^{n} + T_{i+1, j}^{n} \right) + \alpha \frac{\Delta t}{\Delta y^2} \left(T_{i, j-1}^{n} - 2T_{i, j}^{n} + T_{i, j+1}^{n} \right)$$<jupyter_code>mat = scipy.io.loadmat('matlab_files/gotritons.mat')
# The temperature initial condition.    
T = mat["T"]

# The grid coordinates.
# The domain is assumed to be periodic. 
x = mat["xx"]
y = mat["yy"]

# Delta x and y.
dx = x[1, 0] - x[0, 0]
dy = y[0, 1] - y[0, 0]

# thermal diffusivity 
alpha = 2

# final time 
tf = 0.001
fig, ax = plt.subplots(nrows=1, sharex=True,  figsize=(8, 4))
pos = ax.pcolor(x, y, T,  shading='auto')
fig.colorbar(pos, ax=ax)

_ = ax.set_xlabel("x")
_ = ax.set_ylabel("y")
_ = ax.set_title(r"$T(x, y, t=0)$")

plt.tight_layout()
plt.savefig("diffusion/initial_condition_tritons.png", dpi=600)
t = np.linspace(0,tf,int(1e3))
dt = t[1] - t[0]
print("CFL = ", alpha*dt/(dx**2) + alpha*dt/(dy**2))<jupyter_output>CFL =  0.15896146146146078<jupyter_text>$$T_{i, j}^{n+1} = T_{i, j}^{n} + \alpha \frac{\Delta t}{\Delta x^2} \left(T_{i-1, j}^{n} - 2T_{i, j}^{n} + T_{i+1, j}^{n} \right) + \alpha \frac{\Delta t}{\Delta y^2} \left(T_{i, j-1}^{n} - 2T_{i, j}^{n} + T_{i, j+1}^{n} \right)$$<jupyter_code>def advance_diffusion_eq(Tn, alpha, dx, dy, dt):
    Tnext = np.zeros(Tn.shape)
    for ii in range(Tn.shape[0]):
        for jj in range(Tn.shape[1]):  
            # periodic boundary conditions
            if ii == Tn.shape[0]-1 and jj == Tn.shape[1] -1:
                Tnext[ii, jj] = Tn[ii, jj] + alpha*(dt/dx**2)*(Tn[ii-1, jj] -2*Tn[ii, jj] + Tn[0, jj]) + alpha*(dt/dy**2)*(Tn[ii, jj-1] - 2*Tn[ii, jj] + Tn[ii, 0])           
            elif ii == Tn.shape[0]-1:
                Tnext[ii, jj] = Tn[ii, jj] + alpha*(dt/dx**2)*(Tn[ii-1, jj] -2*Tn[ii, jj] + Tn[0, jj]) + alpha*(dt/dy**2)*(Tn[ii, jj-1] - 2*Tn[ii, jj] + Tn[ii, jj+1])           
            elif jj == Tn.shape[1] -1:
                Tnext[ii, jj] = Tn[ii, jj] + alpha*(dt/dx**2)*(Tn[ii-1, jj] -2*Tn[ii, jj] + Tn[ii+1, jj]) + alpha*(dt/dy**2)*(Tn[ii, jj-1] - 2*Tn[ii, jj] + Tn[ii, 0])
            else:
                Tnext[ii, jj] = Tn[ii, jj] + alpha*(dt/dx**2)*(Tn[ii-1, jj] -2*Tn[ii, jj] + Tn[ii+1, jj]) + alpha*(dt/dy**2)*(Tn[ii, jj-1] - 2*Tn[ii, jj] + Tn[ii, jj+1])
    return Tnext
T_sol = np.zeros((T.shape[0], T.shape[1], len(t)))
T_sol[:, :, 0] = T

for ii in range(len(t) - 1):
    T_sol[:, :, ii+1] = advance_diffusion_eq(T_sol[:, :, ii], alpha, dx, dy, dt)
fig, ax = plt.subplots(nrows=1, sharex=True,  figsize=(8, 4))
pos = ax.pcolor(x, y, T_sol[:, :, -1],  shading='auto')
fig.colorbar(pos, ax=ax)

_ = ax.set_xlabel("x")
_ = ax.set_ylabel("y")
_ = ax.set_title(r"$T(x, y, t=$" + str(tf) + ")")

plt.tight_layout()
plt.savefig("diffusion/final_time_diffusion.png", dpi=600)<jupyter_output><empty_output><jupyter_text>Assignment 2.2: Euler Explicit Advection $$\frac{\partial T}{\partial t} = - \bar{c} \cdot \nabla T = - c_{x} \frac{\partial T}{\partial x} - c_{y} \frac{\partial T}{\partial y}$$where $\bar{c}$ is the bulk velocity vector.<jupyter_code>mat = scipy.io.loadmat('matlab_files/gotritons.mat')
# The temperature initial condition.    
T = mat["T"]

# The grid coordinates.
# The domain is assumed to be periodic. 
x = mat["xx"]
y = mat["yy"]

# Delta x and y.
dx = x[1, 0] - x[0, 0]
dy = y[0, 1] - y[0, 0]<jupyter_output><empty_output><jupyter_text>First-order backward differences (up to t=2)<jupyter_code>tf = 2
t = np.linspace(0,tf,int(1e3))
dt = t[1] - t[0]

cx = 1
cy = 1

print("CFL = ", cx*dt/dx + cy*dt/dy)
def advance_diffusion_eq_bkwd(Tn, cx, cy, dx, dy, dt):
    Tnext = np.zeros(Tn.shape)
    for ii in range(Tn.shape[0]):
        for jj in range(Tn.shape[1]):  
            Tnext[ii, jj] = Tn[ii, jj] -cx*(dt/dx)*(Tn[ii, jj] - Tn[ii-1, jj]) -cy*(dt/dy)*(Tn[ii, jj] - Tn[ii, jj-1])
    return Tnext
T_sol = np.zeros((T.shape[0], T.shape[1], len(t)))
T_sol[:, :, 0] = T

for ii in range(len(t) - 1):
    T_sol[:, :, ii+1] = advance_diffusion_eq_bkwd(T_sol[:, :, ii], 1, 1, dx, dy, dt)
fig, ax = plt.subplots(nrows=1, sharex=True,  figsize=(8, 4))
pos = ax.pcolor(x, y, T_sol[:, :, -1], vmin=0, vmax=1, shading='auto')
fig.colorbar(pos, ax=ax)

_ = ax.set_xlabel("x")
_ = ax.set_ylabel("y")
_ = ax.set_title(r"$T(x, y, t=$" + str(tf) + ")")

plt.tight_layout()
plt.savefig("diffusion/final_time_advection_backwards.png", dpi=600)<jupyter_output><empty_output><jupyter_text>Second-order central differences (only up to t=0.25)<jupyter_code>tf = 0.25

t = np.linspace(0,tf,int(5e2))
dt = t[1] - t[0]

cx = 1
cy = 1

print("CFL = ", cx*dt/dx + cy*dt/dy)
def advance_diffusion_eq_central(Tn, cx, cy, dx, dy, dt):
    Tnext = np.zeros(Tn.shape)
    for ii in range(Tn.shape[0]):
        for jj in range(Tn.shape[1]):  
            # periodic boundary conditions
            if ii == Tn.shape[0]-1 and jj == Tn.shape[1]-1:
                Tnext[ii, jj] = Tn[ii, jj] -cx*(dt/dx)*(Tn[0, jj] - Tn[ii-1, jj]) -cy*(dt/dy)*(Tn[ii, 0] - Tn[ii, jj-1])         
            elif ii == Tn.shape[0]-1:
                Tnext[ii, jj] = Tn[ii, jj] -cx*(dt/dx)*(Tn[0, jj] - Tn[ii-1, jj]) -cy*(dt/dy)*(Tn[ii, jj+1] - Tn[ii, jj-1])           
            elif jj == Tn.shape[1]-1:
                Tnext[ii, jj] = Tn[ii, jj] -cx*(dt/dx)*(Tn[ii+1, jj] - Tn[ii-1, jj]) -cy*(dt/dy)*(Tn[ii, 0] - Tn[ii, jj-1])
            else:
                Tnext[ii, jj] = Tn[ii, jj] -cx*(dt/dx)*(Tn[ii+1, jj] - Tn[ii-1, jj]) -cy*(dt/dy)*(Tn[ii, jj+1] - Tn[ii, jj-1])
    return Tnext
T_sol = np.zeros((T.shape[0], T.shape[1], len(t)))
T_sol[:, :, 0] = T

for ii in range(len(t) - 1):
    T_sol[:, :, ii+1] = advance_diffusion_eq_central(T_sol[:, :, ii], 1, 1, dx, dy, dt)
fig, ax = plt.subplots(nrows=1, sharex=True,  figsize=(8, 4))
pos = ax.pcolor(x, y, T_sol[:, :, -1], vmin=0, vmax=1, shading='auto')
fig.colorbar(pos, ax=ax)

_ = ax.set_xlabel("x")
_ = ax.set_ylabel("y")
_ = ax.set_title(r"$T(x, y, t=$" + str(tf) + ")")

plt.tight_layout()
plt.savefig("diffusion/final_time_advection_central.png", dpi=600)<jupyter_output><empty_output><jupyter_text>Assignment 2.3: MacCormack Advection<jupyter_code>tf = 2
t = np.linspace(0,tf,int(1e3))
dt = t[1] - t[0]

cx = 1
cy = 1

print("CFL = ", cx*dt/dx + cy*dt/dy)
print(dx)
print(dy)
def advance_diffusion_eq_maccormack(Tn, cx, cy, dx, dy, dt):
    Tnext = np.zeros(Tn.shape)
    
    for ii in range(Tn.shape[0]):
        for jj in range(Tn.shape[1]): 
            if ii == Tn.shape[0] - 1 and jj == Tn.shape[1] -1:
                Ttemp1 = Tn[ii, jj] - cx*(dt/dx)*(Tn[0, jj] - Tn[ii, jj])- cy*(dt/dy)*(Tn[ii, 0] - Tn[ii, jj])
                Ttemp2 = Tn[ii-1, jj] - cx*(dt/dx)*(Tn[ii, jj] - Tn[ii-1, jj]) - cy*(dt/dy)*(Tn[ii-1, 0] - Tn[ii-1, jj])
                Ttemp3 = Tn[ii, jj-1] - cx*(dt/dx)*(Tn[0, jj-1] - Tn[ii, jj-1]) - cy*(dt/dy)*(Tn[ii, jj] - Tn[ii, jj-1])
            
            elif ii == Tn.shape[0] - 1:
                Ttemp1 = Tn[ii, jj] - cx*(dt/dx)*(Tn[0, jj] - Tn[ii, jj])- cy*(dt/dy)*(Tn[ii, jj+1] - Tn[ii, jj])
                Ttemp2 = Tn[ii-1, jj] - cx*(dt/dx)*(Tn[ii, jj] - Tn[ii-1, jj]) - cy*(dt/dy)*(Tn[ii-1, jj+1] - Tn[ii-1, jj])
                Ttemp3 = Tn[ii, jj-1] - cx*(dt/dx)*(Tn[0, jj-1] - Tn[ii, jj-1]) - cy*(dt/dy)*(Tn[ii, jj] - Tn[ii, jj-1])
            
            elif jj == Tn.shape[1] -1:
                Ttemp1 = Tn[ii, jj] - cx*(dt/dx)*(Tn[ii+1, jj] - Tn[ii, jj])- cy*(dt/dy)*(Tn[ii, 0] - Tn[ii, jj])
                Ttemp2 = Tn[ii-1, jj] - cx*(dt/dx)*(Tn[ii, jj] - Tn[ii-1, jj]) - cy*(dt/dy)*(Tn[ii-1, 0] - Tn[ii-1, jj])
                Ttemp3 = Tn[ii, jj-1] - cx*(dt/dx)*(Tn[ii+1, jj-1] - Tn[ii, jj-1]) - cy*(dt/dy)*(Tn[ii, jj] - Tn[ii, jj-1])

            else:
                Ttemp1 = Tn[ii, jj] - cx*(dt/dx)*(Tn[ii+1, jj] - Tn[ii, jj])- cy*(dt/dy)*(Tn[ii, jj+1] - Tn[ii, jj])
                Ttemp2 = Tn[ii-1, jj] - cx*(dt/dx)*(Tn[ii, jj] - Tn[ii-1, jj]) - cy*(dt/dy)*(Tn[ii-1, jj+1] - Tn[ii-1, jj])
                Ttemp3 = Tn[ii, jj-1] - cx*(dt/dx)*(Tn[ii+1, jj-1] - Tn[ii, jj-1]) - cy*(dt/dy)*(Tn[ii, jj] - Tn[ii, jj-1])
            
            Tnext[ii, jj] = 0.5*(Tn[ii, jj] + Ttemp1) -cx*dt/(2*dx)*(Ttemp1-Ttemp2) -cy*dt/(2*dy)*(Ttemp1-Ttemp3)
    return Tnext
T_sol = np.zeros((T.shape[0], T.shape[1], len(t)))
T_sol[:, :, 0] = T

for ii in range(len(t) - 1):
    T_sol[:, :, ii+1] = advance_diffusion_eq_maccormack(T_sol[:, :, ii], 1, 1, dx, dy, dt)
fig, ax = plt.subplots(nrows=1, sharex=True,  figsize=(8, 4))
pos = ax.pcolor(x, y, T_sol[:, :, -1], vmin=0, vmax=1, shading='auto')
fig.colorbar(pos, ax=ax)

_ = ax.set_xlabel("x")
_ = ax.set_ylabel("y")
_ = ax.set_title(r"$T(x, y, t=$" + str(tf) + ")")

plt.tight_layout()
plt.savefig("diffusion/final_time_advection_maccormack.png", dpi=600)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.23483226266952176, 'fraction_numerical': 0.03415927398796778, 'mean_word_length': 3.1175482787573467, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '653374', 'license': 'MIT', 'n_tokens_mistral': 4861, 'n_tokens_neox': 4673, 'n_words': 915, 'path': 'Diffusion_Advection_Equation.ipynb', 'repo_name': 'opaliss/Computational_Fluid_Dynamics', 'size': '9807'}