[Numpy-discussion] Some help on matlab to numpy translation

[Nicolas Rougier]

Thanks a lot for the translation.

I think the "opp" relates to the opp array declared at the top and the circshift can be made using numpy roll. I modified your translation to include them. 

The results should be something like http://www.lbmethod.org/openlb/gif/karman.avi (I think) but unfortunately, it does not seem to do anything at the moment, I need to investigate further where is the problem.

Nicolas



New version:


'''
Channel flow past a cylindrical obstacle, using a LB method
Copyright (C) 2006 Jonas Latt
Address: Rue General Dufour 24,  1211 Geneva 4, Switzerland
E-mail: Jonas.Latt at cui.unige.ch
'''
import numpy
import matplotlib
matplotlib.use('MacOSX')
import matplotlib.pyplot as plt

# General flow constants
lx = 250
ly = 51
obst_x = lx/5.+1                  # position of the cylinder; (exact
obst_y = ly/2.+1                  # y-symmetry is avoided)
obst_r = ly/10.+1                 # radius of the cylinder
uMax   = 0.02                     # maximum velocity of Poiseuille inflow
Re     = 100                     # Reynolds number
nu     = uMax * 2.*obst_r / Re    # kinematic viscosity
omega  = 1. / (3*nu+1./2.)       # relaxation parameter
maxT   = 4000                    # total number of iterations
tPlot  = 5                       # cycles

# D2Q9 Lattice constants
t  = numpy.array([4/9., 1/9.,1/9.,1/9.,1/9., 1/36.,1/36.,1/36.,1/36.])
cx = numpy.array([  0,   1,  0, -1,  0,    1,  -1,  -1,   1])
cy = numpy.array([  0,   0,  1,  0, -1,    1,   1,  -1,  -1])
opp = numpy.array([ 0,   3,  4,  1,  2,    7,   8,   5,   6])
col = numpy.arange(1, ly - 1)

y,x = numpy.meshgrid(numpy.arange(ly),numpy.arange(lx))
obst = (x-obst_x)**2 + (y-obst_y)**2 <= obst_r**2
obst[:,0] = obst[:,ly-1] = 1
bbRegion = numpy.nonzero(obst)

# Initial condition: (rho=0, u=0) ==> fIn(i) = t(i)
fIn = t[:, numpy.newaxis, numpy.newaxis].\
		repeat(lx, axis = 1).\
		repeat(ly, axis = 2)

# Main loop (time cycles)
for cycle in range(maxT):

   # Macroscopic variables
   rho = fIn.sum(axis = 0)
   ux = (cx[:,numpy.newaxis,numpy.newaxis] * fIn).sum(axis = 0) / rho
   uy = (cy[:,numpy.newaxis,numpy.newaxis] * fIn).sum(axis = 0) / rho

   # Macroscopic (Dirichlet) boundary condition
   L = ly-2.0
   y = col-1.5
   ux[:, 1, col] = 4 * uMax / (L ** 2) * (y * L - y ** 2)
   uy[:, 1, col] = 0
   rho[0, col] = 1 / (1 - ux[0, col]) * \
		   (fIn[[1, 3, 5]][:, 0][:, col].sum(axis = 0) + \
		    2 * fIn[[4, 7, 8]][:, 1][:, col].sum(axis = 0))
   rho[lx - 1, col] = rho[lx - 2, col]
   uy[lx - 1, col] = 0
   ux[lx - 1, col] = ux[:, lx - 2, col]

   fEq = numpy.zeros((9, lx, ly))
   fOut = numpy.zeros((9, lx, ly))
   for i in xrange(0, 9):
	   cu = 3 * (cx[i] * ux + cy[i] * uy)
	   fEq[i] = rho * t[i] * (1 + cu + 0.5 * cu ** 2 - \
			   1.5 * (ux ** 2 + uy ** 2))
	   fOut[i] = fIn[i] - omega * (fIn[i] - fIn[i])

   # Microscopic boundary conditions
   for i in xrange(0, 9):
      # Left boundary:
      fOut[i, 1, col] = fEq[i,0,col] + 18 * t[i] * cx[i] * cy[i] * \
          (fIn[7,0,col] - fIn[6,0,col] - fEq[7,0,col] + fEq[6,0,col])
      fOut[i,lx-1,col] = fEq[i,lx-1,col] + 18 * t[i] * cx[i] * cy[i] * \
          (fIn[5,lx-1,col] - fIn[8,lx-1,col] - fEq[5,lx-1,col] + fEq[8,lx-1,col])
      fOut[i,bbRegion] = fIn[opp[i],bbRegion]

   # Streaming step
   for i in xrange(0,9):
      fIn = numpy.roll(fIn, cx[i], 1)
      fIn = numpy.roll(fIn, cy[i], 2)
   
   if not cycle%tPlot:
      u = numpy.sqrt(ux**2+uy**2)
      #u[bbRegion] = numpy.nan
      print u.min(), u.max()
      #plt.imshow(u)
      #plt.show()




On Mar 13, 2010, at 16:59 , Friedrich Romstedt wrote:

> 2010/3/13 Nicolas Rougier <Nicolas.Rougier at loria.fr>:
>> I'm trying to translate a small matlab program for the simulation in a 2D flow in a channel past a cylinder and since I do not have matlab access, I would like to know if someone can help me, especially on array indexing. The matlab source code is available at: http://www.lbmethod.org/openlb/lb.examples.html and below is what I've done so far in my translation effort.
> 
>> In the matlab code, there is a "ux" array of shape (1,lx,ly) and I do not understand syntax: "ux(:,1,col)" with "col = [2:(ly-1)]". If someone knows, that would help me a lot...
> 
> It means that you select all in the 0-axis all indices, in the 1-axis
> the index 0 (matlab: index 1), and in the 2-axis the indices given by
> the list {col}.  {col} is in our case an ndarray of .ndim = 1.
> 
> I attach a modified version of your script which is running, computing
> *something*.  If you could provide information about matlab functions
> opp() and circshift() that would be helpful.  I marked sections I
> changed with "CHANGED", todos with TODO and lonely notes with NOTE and
> so on.
> 
> Friedrich
> <lbmethod.py>_______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion