1040.sagews

from numpy import matrix, multiply, power, abs, transpose, sign, sum
from numpy.linalg import inv
L_pipe=matrix([3,2,3,2,3,2,3]) #[km]
D_pipe=matrix([30,20,20,20,20,20,20]) #[cm]
CHW_pipe=matrix([100,100,100,100,100,100,100]) # Friction coefficient of Hazen Williams
_=multiply(power(D_pipe,4.87),power(CHW_pipe,1.852))
r_pipe =1.526e7*L_pipe/_ 
Q0_pipe =matrix([110,30,110,110,80,60,140]) #[m^3/hr]
Qd = matrix([-220,0,0,0,20,200]).transpose()  #[m^3/hr]
H1_assumed = 100.0 #[m]
Q_pipe=Q0_pipe
H=0.0
OLD_H=1.0
N_iter = 100
# start the loop
while sum(abs(OLD_H-H))> 0.01:
    M = -1.0/multiply(r_pipe,power(abs(Q_pipe),0.85))  
    # construct the LHS
    LHS=matrix([[-(M[0,0]+M[0,3]),M[0,0],0,M[0,3],0,0],
    [M[0,0],-(M[0,0]+M[0,1]+M[0,4]),M[0,1],0,M[0,4],0],
    [0,M[0,1],-(M[0,1]+M[0,2]+M[0,6]),M[0,2],0,M[0,6]],
    [M[0,3],0,M[0,2],-(M[0,2]+M[0,3]),0,0],
    [0,M[0,4],0,0,-(M[0,4]+M[0,5]),M[0,5]],
    [0,0,M[0,6],0,M[0,5],-(M[0,5]+M[0,6])]])
    LHS[0,0]=10**10
    RHS = -Qd
    RHS[0,0]=H1_assumed*10**10
    OLD_H=H
    H=inv(LHS)*RHS
    hloss = matrix([H[0,0]-H[1,0],  H[1,0]-H[2,0],  H[3,0]-H[2,0] , H[0,0]-H[3,0] , H[1,0]-H[4,0],  H[4,0]-H[5,0] , H[2,0]-H[5,0]])
    _ = power(abs(hloss)/r_pipe,1.0/1.852)
    Q_pipe=multiply(sign(abs(hloss)),_)
    
    
print "Heads:\n", H
print ""
print "Discharges: \n", Q_pipe.transpose()

Heads: [[ 100.00000032] [ 92.89818533] [ 87.71815683] [ 95.08726293] [ 70.52116406] [ 60.46294669]] Discharges: [[ 161.24691038] [ 58.28155829] [ 56.63794734] [ 56.63794734] [ 103.17533775] [ 83.39513087] [ 114.76823954]]