diffusionOmega.py
#
# Iterative iterations for steady diffusion in cooling blades
# Comparison between vectorized Jacobi and loops SOR implementations !
#
# === Gauss-Seidel : omega = 1.87 : elapsed time is 13.600645 seconds.
# T(0,Ly) = 366.099943 [K] after 47046 steps
# |Qin-Qout|/Qin = 9.999992e-04
# === Gauss-Seidel : omega = 1.50 : elapsed time is 41.657255 seconds.
# T(0,Ly) = 366.099921 [K] after 144120 steps
# |Qin-Qout|/Qin = 9.999960e-04
# === Jacobi : elapsed time is 33.492100 seconds.
# T(0,Ly) = 366.099909 [K] after 735885 steps
# |Qin-Qout|/Qin = 9.999907e-04
#
# In terms of CPU for this example, the winner is SOR with optimal omega !
#
# 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
matplotlib.rcParams['toolbar'] = 'None'
# =========================================================================
def diffusionSolveSteady(omega,tol,n,method):
#
# -1- Geometrical, material an numerical aata
#
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
dx = Lx/n; nx = n + 1; ny = 6*n + 1
beta = 0.25; ntmax = 10000000
#
# -2- SOR iterations : cannot be easily vectorized.... but still facter than Jacobi
# if you use the optimal omega = 1.87
#
tic = timer()
T = Tf * ones((nx+2,ny+2))
E = zeros(ntmax); Echeck = 2*tol; it = 0
while (abs(Echeck) >= tol and abs(Echeck) <= 2.0 and it < ntmax) :
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
Qin = qin * Lx
Qout = (sum(T[2:-2,0] - T[2:-2,2])*k/2 +
sum(T[[1,-2],0] - T[[1,-2],2])*k/4 +
sum(T[-1,2:-2] - T[-3,2:-2])*k/2 +
sum(T[-1,[1,-2]] - T[-3,[1,-2]])*k/4 )
Echeck = E[it] = abs(Qin + Qout) / Qin
if (method == "Gauss-Seidel") :
for i in range(1,nx+1):
for j in range(1,ny+1):
T[i,j] += omega * beta * (T[i+1,j] + T[i-1,j]
+ T[i,j+1] + T[i,j-1] - 4*T[i,j])
if (method == "Jacobi") :
T[1:-1,1:-1] += omega * beta*(T[2:,1:-1] + T[:-2,1:-1]
+ T[1:-1,2:] + T[1:-1,:-2] - 4*T[1:-1,1:-1])
it = it+1
toc = timer() - tic
print(' === %s : omega = %4.2f : elapsed time is %f seconds.' % (method,omega,toc))
print(' T(0,Ly) = %10.6f [K] after %d steps ' % (T[1,-2],it))
print(' |Qin-Qout|/Qin = %14.6e ' % E[it-1])
return E[0:it]
# ============================= mainProgram ===============================
plt.figure("Convergence of flux balance")
E = diffusionSolveSteady(1.87,0.001,5,"Gauss-Seidel")
plt.plot(range(len(E)),E,'-b')
E = diffusionSolveSteady(1.50,0.001,5,"Gauss-Seidel")
plt.plot(range(len(E)),E,'-g')
E = diffusionSolveSteady(1.0,0.001,5,"Jacobi")
plt.plot(range(len(E)),E,'-r')
plt.show()
# =========================================================================