Python Thermal lattice Boltzmann model ignoring source term -- python code help please

[Gwen]

LBM model for phase change- relevant equations found here. Also here.

#Thermal LBM
#solves 1D 1 phase phase-change
#D2Q5 Lattice

nx=100                                   # the number of nodes in x direction lattice direction
ny=5                                    # the number of nodes in y direction lattice direction
alpha=.5/3                     # heat diffusion coefficient                                   # the dimension of the problem
mstep=1000                               # the number of time step
tau=3.*alpha+0.5

Tleft=0.0                                 # left wall temperature
Tright=1.0                                # right wall temperature
k=5 # k=0,1,2,3,4,5,6,8,9

x=numpy.linspace(0,1,nx+1) #start at zero, end at 1, fill with nx+1 even spaced intervals
y=numpy.linspace(0,1,ny+1)
t=np.zeros(mstep)
s=np.zeros(mstep)
w=numpy.ones(k)                              # witghting factor
T=numpy.ones((ny+1,nx+1) )         # Temperature matrix
f= numpy.ones((k, ny+1,nx+1))                # distribution function

Hl=1
Hs=0.5
H=numpy.ones((ny+1,nx+1) )                   # Enthalpy matrix
Fl=numpy.ones((ny+1,nx+1) )                  # Liquid fraction matrix (Fl=1 for liquid, Fl=0 for solid)

##================ Initial boundary condition
w[0]=1./3. #0.0
w[1:5]=1./6. #1./4.

##================== Initial value

T[0:ny+1,0:nx+1]=1.0   #temperature in the whole region (including bottom wall)
T[0:ny+1,0]=0        #temperature on the left wall
T[0:ny+1,nx]=1.0       #temperature one node in from the right wall
T[ny,1:nx]=1.0         #temp one node in from the top wall (and one node in from left and right sides)

for i in range(nx+1):
    for j in range(ny+1):
        for l in range (k): #k=0,1,2,3,4    
            f[l,j,i]=w[l]*T[j,i]
  
##   Main loop  : comprised two parts :collision and streaming
##=====================
for n in range(mstep) :
    t[n]=n  #track the time
    time=t[n]

    epsilon=1e-8
    error=1
    Fl_old=Fl
    while error>epsilon:
        Fl_old_iter=Fl
        T_old_iter=T

# collision process
# ==========================
        for i in range(nx+1):
            for j in range(ny+1):
                    for l in range (k):
                        feq=w[l]*T[j,i]  
                        f[l,j,i]=(1.-1/tau)*f[l,j,i]+(1/tau)*feq-w[l]*(Fl[j,i]-Fl_old[j,i])
 
 #streaming process
# ==========================
        for i in range(nx):
            for j in range(ny,0,-1):  #backwards from top to bottom
                f[2,j,i]=f[2,j-1,i]

        for i in range(nx,0,-1):   #backwards from right to left
            for j in range(ny,0,-1):  #backwards from top to bottom
                f[1,j,i]=f[1,j,i-1]

        for i in range(nx,0,-1):   #backwards from right to left
            for j in range(0,ny):     #forward from bottom to second-to-top lattice node
                f[4,j,i]=f[4,j+1,i]

        for i in range(0,nx):      #forward from left to second-to-right lattice node
            for j in range(0,ny):     #forward from bottom to second-to-top lattice node
                f[3,j,i]=f[3,j,i+1]

# Boundary conditions
#  =============================
        for j in range(0,ny+1) :               #left Boundary. Dirichlet boundary condition: constant temperature.
            f[1,j,0]=( Tleft*(w[1]+w[3]) )-f[3,j,0]

        for j in range(0,ny+1):                #right Boundary. adiabatic
            f[3,j,nx]=f[1,j,nx]

        for i in range(0,nx+1):                # bottom and top Boundary
            f[4,ny,i]=f[2,ny,i]                  #adiabatic
 #================================ #calculate temperature
        for i in range(nx+1):
            for j in range(ny+1):
                sum=0.0
                for l in range (k):
                    sum=sum+f[l,j,i]
                T[j,i]=sum
        T[0:ny+1,0]=Tleft           #Dirichlet BC      
        T[0:ny+1,nx]=T[0:ny+1,nx]   #adiabatic BC        
        T[ny,1:nx]=T[ny-1,1:nx]     #adiabatic BC        
        T[0,1:nx]=T[1,1:nx]         #adiabatic BC       
#==============================   #calculate nodal enthalpy and liquid fraction
        for i in range(nx+1):
            for j in range(ny+1):
                H[j,i]=0.5*T[j,i]+0.5*Fl[j,i]
        for i in range(nx+1):
            for j in range(ny+1):
                if H[j,i]<=Hs:
                    Fl[j,i]=0
                elif H[j,i]>Hs and H[j,i] < Hl:
                    Fl[j,i]=(H[j,i]-Hs)/(Hl-Hs)
                else:
                    Fl[j,i]=1
#==============================   #convergence? If no, go back
        for i in range(nx+1):
            for j in range(ny+1):
                error_Fl=abs(np.max(np.max((Fl[j,i]-Fl_old_iter[j,i])/Fl[j,i])))
                error_T=abs(np.max(np.max((T[j,i]-T_old_iter[j,i])/T[j,i])))
                error=np.max([error_Fl, error_T])      
#find position of phase change boundary (where Fl<=0.5)  
    Fl_col=Fl[3,:]<=0.5
    max = Fl_col[0]
    index = 0
    for i in range(1,len(Fl_col)):
        if Fl_col[I] >= 0.5:
            max = Fl_col
            index = i[/I][/I]
 
        s[n]=index/nx   
#==============================[/I][/I]

 

[Arjan82]

I get a division by zero error at line 119