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

[Nicolas Rougier]

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.Latt at 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)]);