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.