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

Bruce Southey bsouthey@gmail....
Sun Mar 14 09:55:07 CDT 2010

```On Sun, Mar 14, 2010 at 7:14 AM, Nicolas Rougier
<Nicolas.Rougier@loria.fr> wrote:
>
> 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@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

You probably should split these into the creation part
(cy[:,numpy.newaxis,numpy.newaxis]) and multiplication etc part ( *
fIn).sum(axis = 0) / rho). It will save the repeated memory
allocation.

>
>   # Macroscopic (Dirichlet) boundary condition
>   L = ly-2.0
>   y = col-1.5
These two variables can be moved out in the loop

>   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])

This line is probably wrong,  (fln[i]-fEq[i])?

In anycase, I do not see the need for the loop (which why it caught my
attention). If you are just indexing the same 'row' of an array using
a loop then you probably don't need the loop. So you probably should
be able to drop the indexing (assuming arrays) something like this,
which is untested:
cu-3*(cx*ux +cy*uy)
fEq= rho * t * (1 + cu + 0.5 * cu ** 2 -1.5 * (ux ** 2 + uy ** 2))
fOut = fIn - omega * (fIn - fEq)

If these are correct then you don't need the creation of fEQ and fOut
as numpy.zeros.

>
>   # 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]

You should be able to drop the indexing and loop here as well.

>
>   # Streaming step
>   for i in xrange(0,9):
>      fIn = numpy.roll(fIn, cx[i], 1)
>      fIn = numpy.roll(fIn, cy[i], 2)

Likewise, it would not surprise me if this loop is not needed.

Bruce

>
>   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@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@scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
```