[Numpy-discussion] Some help on matlab to numpy translation
Hello, 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... Nicolas ''' 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.l...@cui.unige.ch ''' import numpy # 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 # matlab: t = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36]; # matlab: cx = [ 0, 1, 0, -1, 0,1, -1, -1, 1]; # matlab: cy = [ 0, 0, 1, 0, -1,1, 1, -1, -1]; # matlab: opp = [ 1, 4, 5, 2, 3,8, 9, 6, 7]; # matlab: col = [2:(ly-1)]; 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([ 1, 4, 5, 2, 3,8, 9, 6, 7]) col = numpy.arange(2,(ly-1)) # matlab: [y,x] = meshgrid(1:ly,1:lx); # matlab: obst = (x-obst_x).^2 + (y-obst_y).^2 = obst_r.^2; # matlab: obst(:,[1,ly]) = 1; # matlab: bbRegion = find(obst); 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) # matlab: fIn = reshape( t' * ones(1,lx*ly), 9, lx, ly); fIn = numpy.ones((lx,ly,9)) fIn [:,:] = t # MAIN LOOP (TIME CYCLES) # matlab: for cycle = 1:maxT for cycle in range(maxT): # MACROSCOPIC VARIABLES # matlab: rho = sum(fIn); # matlab: ux = reshape ( ... # matlab: (cx * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho; # matlab: uy = reshape ( ... # matlab: (cy * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho; rho = fIn.sum(-1) ux = (cx*fIn).sum(-1)/rho uy = (cx*fIn).sum(-1)/rho # MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS # Inlet: Poiseuille profile # matlab: L = ly-2; y = col-1.5; # matlab: ux(:,1,col) = 4 * uMax / (L*L) * (y.*L-y.*y); # matlab: uy(:,1,col) = 0; # matlab: rho(:,1,col) = 1 ./ (1-ux(:,1,col)) .* ( ... # matlab: sum(fIn([1,3,5],1,col)) + ... # matlab: 2*sum(fIn([4,7,8],1,col)) ); L = ly-2.0 y = col-1.5 # Is that right ? ux[0,1:-1] = 4 * uMax / (L*L) * (y.*L-y.*y) # Is that right ? uy[0,1:-1] = 0 rho[0,1:-1] = 1./(1-ux[1,1:-1])*(sum(fIn[ ([1,3,5],1,col)) + 2*sum(fIn([4,7,8],1,col))) # Outlet: Zero gradient on rho/ux # matlab: rho(:,lx,col) = rho(:,lx-1,col); # matlab: uy(:,lx,col) = 0; # matlab: ux(:,lx,col) = ux(:,lx-1,col); # % Outlet: Zero gradient on rho/ux # rho(:,lx,col) = rho(:,lx-1,col); # uy(:,lx,col) = 0; # ux(:,lx,col) = ux(:,lx-1,col); # % COLLISION STEP # for i=1:9 #cu = 3*(cx(i)*ux+cy(i)*uy); #fEq(i,:,:) = rho .* t(i) .* ( 1 + cu + 1/2*(cu.*cu) - 3/2*(ux.^2+uy.^2) ); #fOut(i,:,:) = fIn(i,:,:) - omega .* (fIn(i,:,:)-fEq(i,:,:)); # end # % MICROSCOPIC BOUNDARY CONDITIONS # for i=1:9 # % Left boundary # fOut(i,1,col) = fEq(i,1,col) + 18*t(i)*cx(i)*cy(i)* ( fIn(8,1,col) - fIn(7,1,col)-fEq(8,1,col)+fEq(7,1,col) ); # % Right boundary # fOut(i,lx,col) = fEq(i,lx,col) + 18*t(i)*cx(i)*cy(i)* ( fIn(6,lx,col) - fIn(9,lx,col)-fEq(6,lx,col)+fEq(9,lx,col) ); # % Bounce back region # fOut(i,bbRegion) = fIn(opp(i),bbRegion); # % STREAMING STEP # for i=1:9 #fIn(i,:,:) = circshift(fOut(i,:,:), [0,cx(i),cy(i)]); ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Some help on matlab to numpy translation
Le samedi 13 mars 2010 à 10:20 +0100, Nicolas Rougier a écrit : Hello, 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... As ux 's shape is (1,lx,ly), ux(:,1,col) is equal to ux(1,1,col) which is a vector with the elements [ux(1,1,2), ... ux(1,1,ly-1)]. Using : juste after the reshape seems a lit bit silly... -- Fabrice Silva LMA UPR CNRS 7051 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Some help on matlab to numpy translation
Thanks. I agree that the use of ':' is a bit weird. Nicolas On Mar 13, 2010, at 11:45 , Fabrice Silva wrote: Le samedi 13 mars 2010 à 10:20 +0100, Nicolas Rougier a écrit : Hello, 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... As ux 's shape is (1,lx,ly), ux(:,1,col) is equal to ux(1,1,col) which is a vector with the elements [ux(1,1,2), ... ux(1,1,ly-1)]. Using : juste after the reshape seems a lit bit silly... -- Fabrice Silva LMA UPR CNRS 7051 ___ 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
Re: [Numpy-discussion] Some help on matlab to numpy translation
On Sat, Mar 13, 2010 at 4:45 AM, Fabrice Silva si...@lma.cnrs-mrs.fr wrote: Le samedi 13 mars 2010 à 10:20 +0100, Nicolas Rougier a écrit : Hello, 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... As ux 's shape is (1,lx,ly), ux(:,1,col) is equal to ux(1,1,col) which is a vector with the elements [ux(1,1,2), ... ux(1,1,ly-1)]. Using : juste after the reshape seems a lit bit silly... Except that python uses 0-based indexing and does not include the last number in a slice, while Matlab uses 1-based indexing and includes the last number, so really: ux(:,1,col) becomes: ux(0, 0, col) # or ux(:, 0, col) And if col is col = [2:(ly-1)] This needs to be: col = np.arange([1, ly - 1) Ryan -- Ryan May Graduate Research Assistant School of Meteorology University of Oklahoma ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Some help on matlab to numpy translation
2010/3/13 Nicolas Rougier nicolas.roug...@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
Description: Binary data
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Some help on matlab to numpy translation
As ux 's shape is (1,lx,ly), ux(:,1,col) is equal to ux(1,1,col) which is a vector with the elements [ux(1,1,2), ... ux(1,1,ly-1)]. Using : juste after the reshape seems a lit bit silly... Except that python uses 0-based indexing and does not include the last number in a slice, while Matlab uses 1-based indexing and includes the last number, so really: ux(:,1,col) becomes: ux(0, 0, col) # or ux(:, 0, col) And if col is col = [2:(ly-1)] This needs to be: col = np.arange([1, ly - 1) You are right about the 0 or 1 based indexing argument, but I was speaking matlab language as visible in the symbols used for indexing ( () and not [] )... :) -- Fabrice Silva LMA UPR CNRS 7051 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] problem with applying patch
I have installed latest numpy from svn on ubuntu karmic(9.10) and compiled the code in my sand box. I downloaded two patches (test_umath_complex.py.patch, build_clib.py.patch) from the review patch directory on http://projects.scipy.org/numpy/report/12 page for testing. The following patch command gives me patch -p0 -i test_umath_complex.py.patch patch unexpectedly ends in middle of line patch: Only garbage was found in the patch input. The patch files are placed in the top level numpy sandbox directory and are in xml format rather than a diff file. I used save link as on the link to download the patch on the Attachments section of the link page. Checked the test_umath_complex.py in the numpy/core/test directory and found significant differences with the listing on the patch link (both old, and the new versions). Don't think this is limited to just these patches any patch that I tried gives the same error message. Thank you regards alex ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] problem with applying patch
Sat, 13 Mar 2010 12:06:24 -0700, z99719 z99719 wrote: [clip] The patch files are placed in the top level numpy sandbox directory and are in xml format rather than a diff file. I used save link as on the link to download the patch on the Attachments section of the link page. That gives you the HTML page. Click on the link, and choose Download in other formats: Original Format from the bottom of the page. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] integer division rule?
Francesc Altet once provided an example that for integer division, numpy uses the C99 rule: round towards 0. This is different than Python's rule for integer division: round towards negative infinity. But I cannot reproduce his example. (NumPy 1.3.) Did this behavior change at some point? Thanks, Alan Isaac ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] long integers in genfromtxt
I was trying to find out what the helpful message
TypeError: expected a readable buffer object means
and it seems genfromtxt has problems identifying long integers (at
least on Windows 32)
np.array(416068,int)
Traceback (most recent call last):
File pyshell#4, line 1, in module
np.array(416068,int)
OverflowError: long int too large to convert to int
s = '''Date,Open,High,Low,Close,Volume,Adj Close
... 2010-02-12,1075.95,1077.81,1062.97,1075.51,416068,1075.51
... 2010-02-11,1067.10,1080.04,1060.59,1078.47,440087,1078.47'''
sh = StringIO(s)
data = np.genfromtxt(sh, delimiter=,, dtype=None, names=True)
Traceback (most recent call last):
File C:\Programs\Python25\Lib\site-packages\numpy\lib\io.py, line
1367, in genfromtxt
output = np.array(data, dtype=ddtype)
TypeError: expected a readable buffer object
sh = StringIO(s)
data = ts.tsfromtxt(sh, delimiter=,, dtype=None, datecols=(0,),
names=True)
Traceback (most recent call last):
File \Programs\Python25\Lib\site-packages\scikits\timeseries\extras.py,
line 504, in tsfromtxt
File \Programs\Python25\Lib\site-packages\scikits\timeseries\_preview.py,
line 1312, in genfromtxt
TypeError: expected a readable buffer object
dt=
[('','S10'),('',float),('',float),('',float),('',float),('',int),('',float)]
sh = StringIO(s)
data = np.genfromtxt(sh, delimiter=,, dtype=dt, names=True)
Traceback (most recent call last):
File C:\Programs\Python25\Lib\site-packages\numpy\lib\io.py, line
1388, in genfromtxt
rows = np.array(data, dtype=[('', _) for _ in dtype_flat])
TypeError: expected a readable buffer object
sh = StringIO(s)
data = ts.tsfromtxt(sh, delimiter=,, dtype=dt, datecols=(0,), names=True)
Traceback (most recent call last):
File \Programs\Python25\Lib\site-packages\scikits\timeseries\extras.py,
line 451, in tsfromtxt
File \Programs\Python25\Lib\site-packages\scikits\timeseries\_preview.py,
line 798, in easy_dtype
File \Programs\Python25\Lib\site-packages\scikits\timeseries\_preview.py,
line 364, in __call__
File \Programs\Python25\Lib\site-packages\scikits\timeseries\_preview.py,
line 329, in validate
TypeError: object of type 'bool' has no len()
using float works, but I needed the debugger to figure out what the
problem with my data is (or whether I just make mistakes)
dt=
[('','S10'),('',float),('',float),('',float),('',float),('',float),('',float)]
sh = StringIO(s)
data = np.genfromtxt(sh, delimiter=,, dtype=dt, names=True)
data
array([ ('2010-02-12', 1075.95, 1077.80999, 1062.97, 1075.51,
416068.0, 1075.51),
('2010-02-11', 1067.0, 1080.04, 1060.58999,
1078.47, 440087.0, 1078.47)],
dtype=[('Date', '|S10'), ('Open', 'f8'), ('High', 'f8'),
('Low', 'f8'), ('Close', 'f8'), ('Volume', 'f8'), ('Adj_Close',
'f8')])
ts.version.version
'0.91.3'
Josef
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] integer division rule?
la, 2010-03-13 kello 15:32 -0500, Alan G Isaac kirjoitti: Francesc Altet once provided an example that for integer division, numpy uses the C99 rule: round towards 0. This is different than Python's rule for integer division: round towards negative infinity. But I cannot reproduce his example. (NumPy 1.3.) Did this behavior change at some point? It was changed in r5888. What the rationale was is not clear from the commit message. Pauli ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] integer division rule?
On Sat, Mar 13, 2010 at 3:04 PM, Pauli Virtanen p...@iki.fi wrote: la, 2010-03-13 kello 15:32 -0500, Alan G Isaac kirjoitti: Francesc Altet once provided an example that for integer division, numpy uses the C99 rule: round towards 0. This is different than Python's rule for integer division: round towards negative infinity. But I cannot reproduce his example. (NumPy 1.3.) Did this behavior change at some point? It was changed in r5888. What the rationale was is not clear from the commit message. The change was before that, the logic of the loop after r5888 is the same as before. I suspect the change was made, whenever that was, in order to conform to python. Chuck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
