In [ ]:

    

%pylab inline


    

In [ ]:

    

import numpy as np
import time
import matplotlib.pyplot as plt
from IPython.display import clear_output
from mpl_toolkits.mplot3d.axes3d import Axes3D

# Grid spacing
dx = 0.5
dy = 0.5

# Heat diffusion constant
kappa = 0.3

# Maximum time step size (constrained by CFL)
#dt = (dx*dx)/(4*kappa) # ok when dx==dy
dt = (dx*dx*dy*dy)/(dx*dx+dy*dy) / (2*kappa) # if dx!=dy

# Number of grid cells
n = 6

# Boundary conditions
ub = np.zeros( (n, n) )
ub[0]      = np.linspace(0.5, 1.5, n) # Top row
ub[n-1]    = np.linspace(0.5, 1.5, n) # Bottom row
ub[:, 0]   = np.linspace(0.5, 0.5, n) # Left column
ub[:, n-1] = np.linspace(1.5, 1.5, n) # Right column

# Grid (fikse denne, bruke dx og dy)
x = np.outer( np.linspace(1, 1, n), np.linspace(0.0, 1.0, n) )
y = np.outer( np.linspace(1.0, 0.0, n), np.linspace(1, 1, n) )


    

In [ ]:

    

def Fx(u_l, u_r):
    return kappa*(u_l - u_r)
def Fy(u_b, u_t):
    return kappa*(u_b - u_t)


    

In [ ]:

    

# Initial temperatures
u = np.zeros( (n, n) )
u[0]      = ub[0]
u[n-1]    = ub[n-1]
u[:, 0]   = ub[:, 0] 
u[:, n-1] = ub[:, n-1]

# Fluxes
fx = zeros( (n, n+1) ) # We have n+1 interfaces for n cells 
fy = zeros( (n+1, n) ) # We have n+1 interfaces for n cells 

# print or plot?
pl = True

if (pl):
    fig = plt.figure(figsize=(10, 6))
    ax = fig.add_subplot(1, 1, 1, projection='3d')
    ax.view_init(50, -110)

for k in range(0, 10):
    # Perform one timestep

    t0 = time.time()
    
    # Loop over all faces (with normal in x-direction) and calculate the flux
    for i in range(n):
        for j in range(1, n): # This will leave the flux in [0] = [n] = 0
            fx[i, j] = Fx(u[i, j-1], u[i, j]) # fx[i, j] corresponds to the flux between cell (i, j-1) and (i, j) (F_{i-1/2})
    # Similarly for the other direction
    for j in range(n):
        for i in range(1, n):
            fy[i, j] = Fy(u[i-1, j], u[i, j]) # fy[i, j] corresponds to the flux between cell (i-1, j) and (i, j)
    
    # For each cell, loop over all faces for that cell, and sum the flux changes
    for i in range(0, n):
        for j in range(0, n):
            u[i, j] = u[i, j] + (fx[i, j] - fx[i, j+1]) * dt / (dx*dx)
            u[i, j] = u[i, j] + (fy[i, j] - fy[i+1, j]) * dt / (dy*dy)
            
    # Apply boundary conditions
    u[0]      = ub[0]
    u[n-1]    = ub[n-1]
    u[:, 0]   = ub[:, 0]
    u[:, n-1] = ub[:, n-1]
    
    comp_time = time.time()-t0
    t0 = time.time()
    
    if (pl):
        ax.cla()
        clear_output()
        ax.set_title( "Time = " + str(dt*k) )
        ax.grid(False) # When 'True', the grid "shines through". How to avoid this?
        #ax.grid(color='r', linewidth=10) # funker ikke
        #ax.set_alpha(1.0) # funker ikke
        #ax.set_axisbelow(False) # funker ikke
        #[line.set_zorder(3) for line in ax.lines] # funker heller ikke
        ax.plot_surface( x, y, u, rstride=1, cstride=1, antialiased=False, linewidth=1, shade=True );
        ax.set_xlabel("x")
        ax.set_ylabel("y");
        ax.set_zlabel("temp");
        ax.set_zlim(0.5, 1.5)
        display(fig)
        
    print "computation: ", comp_time, "plotting: ", time.time() - t0;

plt.close()


    

In [ ]: