[Numpy-discussion] Some help on matlab to numpy translation
Bruce Southey
bsouthey at gmail.com
Sun Mar 14 10:55:07 EDT 2010
On Sun, Mar 14, 2010 at 7:14 AM, Nicolas Rougier
<Nicolas.Rougier at 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 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
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 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
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
More information about the NumPy-Discussion
mailing list