Iterating within a function using for loop

[Simplicity.]

Dear Folks,

When I iterate within the very bottom function under the comment "#Update conserved quantities" in this module using for loop (see below), I have noticed that the block of statements below it is not iterated,

""" """
from numpy import *
from math import *
from scitools.std import *

#Global constants
N = 100; M = 200; dx = 1.0/M ; cfl = 0.75 ; a_3 =   5.045; a_2 =1.68; a_1 = 1 ; dt = cfl*dx/sqrt(a_3); Gamma = 1.4 ; Eps = 10e-8
K  = linspace(0.0,1.0, M+1) ;  rho =zeros(M); u = zeros(M); p = zeros(M)

#Initialize conserved quantities
for i in range(M):			
    if K[i]>= 0 and K[i] <= 0.5:  u[i] = 0.0
    elif K[i] > 0.5 and K[i] <= 1 : u[i] = 0.0

def velocity(u):
    return u

for i in range(M):			
    if K[i]>= 0 and K[i] <= 0.5:  p[i] = 1.0
    elif K[i] > 0.5 and K[i] <= 1 : p[i] = 0.1	

def pressure(p):
    return p

#Conserved variables	
for i in range(M):		                       
    if K[i]>=0 and K[i] <= 0.5:  rho[i] = 1.0
    elif K[i] > 0.5 and K[i] <= 1: rho[i] = 0.125

def G(rho):
    return rho

m = G(rho)*velocity(u)
def momentum(m):
    return m

E = pressure(p)/(Gamma-1) + 0.5*(G(rho))*(momentum(m)/G(rho))**2
def Energy(E):
    return E

#Initialize v_rho, v_m, v_E
v_rho = momentum(m)                                   
def rho_v(v_rho):
    return v_rho

v_m = (momentum(m))**2/(G(rho)) + (Gamma -1)*(Energy(E)- (momentum(m))**2/(2*G(rho)))
def m_v(v_m):
    return v_m

v_E = (Energy(E) + (Gamma -1.0)*(Energy(E)- (momentum(m))**2/(2*G(rho))))*(momentum(m)/G(rho))
def E_v(v_E):
    return v_E

#Solve for Density
def Evolve_rho_s(rho):                                
    rho_s = zeros((M,), float32)
    for k in range (1, M-1):	
        rho_s[k] =	Eps/(Eps -dt)*(v_rho[k] - (dt/Eps)*m[k])
        rho_s[0] = rho_s[1]
        rho_s[M-1] = rho_s[M-2]
    return rho_s

def Evolve_rho1(rho):
    rho1 = zeros((M,), float32)
    rhos = Evolve_rho_s(rho)
    for k in range (1, M-1):
        rho1[k] =  rho[k]-0.5*(dt/dx)*(rhos[k+1] - rhos[k-1]) + 0.5*(dt/dx)*sqrt(a_1)*(rho[k+1] - 2*rho[k] + rho[k-1])
        rho1[0] = rho1[1]
        rho1[M-1] = rho1[M-2]
    return rho1

def Evolve_rho_v(rho):
    rho_v = zeros((M,), float32)
    rhos = Evolve_rho_s(rho)
    for k in range (1, M-1):
        rho_v[k] = rhos[k]- 0.5*(dt/dx)*a_1*(rho[k+1] - rho[k-1]) + 0.5*(dt/dx)*(a_1/sqrt(a_1))*(rhos[k+1] - 2*rhos[k] + rhos[k-1])
        rho_v[0] = rho_v[1]
        rho_v[M-1] = rho_v[M-2]
    return rho_v

#Solve for momentum, m
def Evolve_ms(m):                             
    ms = zeros((M,), float32)
    for k in range (1, M-1):
        ms[k] =  Eps/(Eps -dt)*(v_m[k] - (dt/Eps)*((m[k])**2/(rho[k]) + (Gamma -1)*(E[k]- (m[k])**2/(2*rho[k]))))
        ms[0] = ms[1]
        ms[M-1] = ms[M-2]
    return ms

def Evolve_m1(m):
    m1 = zeros((M,), float32)
    ms = Evolve_ms(m)
    for k in range (1, M-1):
        m1[k] =  m[k]-0.5*(dt/dx)*(ms[k+1] - ms[k-1]) + 0.5*(dt/dx)*sqrt(a_2)*(m[k+1] - 2*m[k] + m[k-1])	
        m1[0] = m1[1]
        m1[M-1] = m1[M-2]
    return m1

def Evolve_mv(m):
    mv = zeros((M,), float32)
    ms = Evolve_ms(m)
    for k in range (1, M-1):
        mv[k] = ms[k]- 0.5*(dt/dx)*a_2*(m[k+1] - m[k-1]) + 0.5*(dt/dx)*(a_2/sqrt(a_2))*(ms[k+1] - 2*ms[k] + ms[k-1])
        mv[0] = mv[1]
        mv[M-1] = mv[M-2]
    return mv

#Solve for Energy, E
def Evolve_Es(E):                                              
    E_s = zeros((M,), float32)
    for k in range (1, M-1):	
        E_s[k] =  Eps/(Eps -dt)*(v_E[k] - (dt/Eps)*(E[k] + (Gamma -1.0)*(E[k]- (m[k])**2/(2*rho[k])))*(m[k]/rho[k]))
        E_s[0] = E_s[1]
        E_s[M-1] = E_s[M-2]	
    return E_s

def Evolve_E1(E):
    E1 = zeros((M,), float32)
    Es = Evolve_Es(E)
    for k in range (1, M-1):	
        E1[k] =  E[k]-0.5*(dt/dx)*(Es[k+1] - Es[k-1]) + 0.5*(dt/dx)*sqrt(a_3)*(E[k+1] - 2*E[k] + E[k-1])
        E1[0] = E1[1]
        E1[M-1] = E1[M-2]	
    return E1

def Evolve_Ev(E):
    Ev = zeros((M,), float32)
    Es = Evolve_Es(E)
    for k in range (1, M-1):
        Ev[k] = Es[k] - 0.5*(dt/dx)*a_3*(E[k+1] - E[k-1]) + 0.5*(dt/dx)*(a_3/sqrt(a_3))*(Es[k+1] - 2*Es[k] + Es[k-1])
        Ev[0] = Ev[1]
        Ev[M-1] = Ev[M-2]
    return Ev

#Update conserved quantities 
def F(rho, m, E, v_rho, v_m, v_E):
    for n in xrange(N):                                  
        rho = G(Evolve_rho1(rho))
        m  = momentum(Evolve_m1(m))
        E = Energy(Evolve_E1(E))
        v_rho = rho_v(Evolve_rho_v(rho))
        v_m =   m_v(Evolve_mv(m))	
        v_E  =  E_v(Evolve_Ev(E))
    return (rho, m, E, v_rho, v_m, v_E)

Calling the function as,
density, momentum, Energy, dflux, mflux, Eflux = F(rho, m, E, v_rho, v_m, v_E), and
plot for example, density, gives a wrong profile.
But if I only use for loop as below it works fine.

#Update conserved quantities
for n in xrange(N):
    rho = G(Evolve_rho1(rho))
    m = momentum(Evolve_m1(m))
    E = Energy(Evolve_E1(E))
    v_rho = rho_v(Evolve_rho_v(rho))
    v_m = m_v(Evolve_mv(m))
    v_E = E_v(Evolve_Ev(E))

Plotting for example rho,
plot(linspace(0.0,1.0,M),rho,'r0.2'),
gives a correct profile.

Any suggestion is welcome and appreciated.

[Gribouillis]

The body of Evolve_ms(), which is called by Evolve_m1() and Evolve_mv(), uses the global E and not the local E wich is updated in each iteration in F(). A solution is to pass E as a parameter and define Evolve_ms(m, E), Evolve_m1(m, E) and Evolve_mv(m, E).

There may be other errors like this, but I couldn't see them. It would be a good idea to use a class with methods instead of a bunch of global functions and variables.

[JoshuaBurleson]

Gribouillis is right, using classes and creating instances would dramatically cut down on bugs like whether the variable your calling is within global scope or not as it would eliminate the need for several global variables. Just in case you're very new and not aware, to utilize a global variable within a function you need to declare that it is of global scope prior to utilizing it. i.e.

for n in xrange(N):
    global E
    rho = G(Evolve_rho1(rho))
    m = momentum(Evolve_m1(m))
    E = Energy(Evolve_E1(E))
    v_rho = rho_v(Evolve_rho_v(rho))
    v_m = m_v(Evolve_mv(m))
    v_E = E_v(Evolve_Ev(E))

[Simplicity.]

What Gribouillis said is true, passing rho, v_rho, m, v_m, etc., as parameters and define Evolve_rho_s(rho,v_rho,m), Evolve_rho1(rho,v_rho,m),Evolve_rho_v(rho,v_rho,m), Evolve_ms(rho, m,v_m,E) and so, worked well. Using class with methods and declaring globals (with caution since declaring globals can mess up the code) could also be other pretty alternatives. Thank you Gribouillis and pyguy62 for your constructive ideas.