diffusionCpu.py

#
# Unsteady diffusion in cooling blades 
# Comparison between vectorized and loops-implementations !
#
# === Vectorized version : elapsed time is 2.662588 seconds.
#   T(0,Ly) = 335.505622 [K] after 80744 steps 
# === Loops version : elapsed time is 25.453349 seconds.
#   T(0,Ly) = 335.505622 [K] after 80744 steps 
#
# on my laptop 25/11/2018 : birthday of my Mother : 91 years :-)
#
# Vincent Legat - 2018
# Ecole Polytechnique de Louvain
#

from numpy import *
import matplotlib 
import matplotlib.pyplot as plt 
from timeit import default_timer as timer

# =========================================================================

def diffusionSolve(betaRequested,n):

# 
# -1- Geometrical and Material Data
# 
  
  rho = 2707; c = 896; k = 204
  alpha = k / (rho * c)
  Ly = 0.015; Lx = Ly/6; Lt = 60
  qin = 30000; Tf = 300; hx = 60; hy = 100
  
# 
# -2- Discretization Parameters
# 
  
  dx = Lx/n; nx = n + 1; ny = 6*n + 1
  dt = (dx * dx) * betaRequested / alpha
  nt = int(ceil(Lt/dt)); dt = Lt/nt
  beta = alpha * dt / (dx * dx)
  
# 
# -3- Smart vectorized computation 
# 
  tic = timer()
  T = Tf * ones((nx+2,ny+2))  
  for it in range(nt):
    T[1:-1,0]  = T[1:-1,2]  - 2*dx*hy*(T[1:-1,1]-Tf) /k
    T[1:-1,-1] = T[1:-1,-3] + 2*dx*qin /k
    T[0,1:-1]  = T[2,1:-1]
    T[-1,1:-1] = T[-3,1:-1] - 2*dx*hx*(T[-2,1:-1]-Tf) /k
    T[1:-1,1:-1] += beta*(T[2:,1:-1] + T[:-2,1:-1]   
                        + T[1:-1,2:] + T[1:-1,:-2] - 4*T[1:-1,1:-1])
  toc = timer() - tic
  print(' === Vectorized version : elapsed time is %f seconds.' % toc)
  print('   T(0,Ly) = %10.6f [K] after %d steps ' % (T[1,-2],nt))
  
# 
# -3- Computation with loops 
# 
  tic = timer()
  T = Tf * ones((nx+2,ny+2)) 
  Tnew = zeros((nx+2,ny+2))
  for it in range(nt):
    for i in range(1,nx+1):
      T[i,0]  = T[i,2]  - 2*dx*hy*(T[i,1]-Tf) /k
      T[i,-1] = T[i,-3] + 2*dx*qin /k
    for j in range(1,ny+1):
      T[0,j]  = T[2,j]
      T[-1,j] = T[-3,j] - 2*dx*hx*(T[-2,j]-Tf) /k
    for i in range(1,nx+1):
      for j in range(1,ny+1):
        Tnew[i,j] = T[i,j] + beta*(T[i+1,j] + T[i-1,j]   
                           + T[i,j+1] + T[i,j-1] - 4*T[i,j])
    T[:] = Tnew[:]
  toc = timer() - tic
  print(' === Loops version : elapsed time is %f seconds.' % toc)
  print('   T(0,Ly) = %10.6f [K] after %d steps ' % (T[1,-2],nt))

  
# ============================= mainProgram ===============================


diffusionSolve(0.25,5)


# =========================================================================

• Download
• Back to EPL1104