Introduces the **PyNetwork** class used to set up and run simulations
An instance of this **PyNetwork** class is a conceptual representation of a water distribution network. To create an instance of this class, you will use two files to specify parameters:
You need to specify network connectivity externally from a .inp file (I suggest generating with EPANET. If the formatting is different from that of EPANET, the file may not be read properly. All other parameters--pipe length, diameter, manning coefficient, elevation, etc can be specified from python notebook and used to write new files to run simulations.
In general, before you begin:
python cleanupinpfiles.py mynetwork.inp
in a terminal in the folder where you have created myfile.
**Double check** that lengths, diameters, elevations etc in mynetwork2.0.inp have correct units (should be meters for lengths, and a manning coefficient for roughness)
Now follow steps below to change parameters and set up your own simulation.
In [3]:
from IPython.display import Image
Image("/Users/anna/Desktop/export_inp.png", width=720, height=450)
Out[3]:
The cell below loads modules you will need for creating networks, solving, plotting, etc.
In [4]:
import sys
sys.path.append("..")
from allthethings import PyNetwork, PyPipe_ps
from allthethings import PyBC_opt_dh
import numpy as np
import matplotlib.pyplot as plt
from writeit import *
%pylab inline
Now we're going to create a network based on the files we just generated using cleanupinpfile.py
In [5]:
fi = "../indata/mynetwork2.0.inp"
fc = "../indata/mynetwork2.0.config"
mtype =1 #this specifies Preissman slot model
n0 = PyNetwork(fi,fc,mtype) # an instance of the network class
Now let's look at what we've created. A network object consists of "edges" (pipes), "nodes" (junctions), runtime info, and various methods.
In [6]:
#"print" shows the memory address and the number of nodes and edge
print n0
In [7]:
#showLayout() shows the connectivity--how the nodes and edges are connected
n0.showLayout()
#you can also plot the network layout
(xs,ys,conns,ls) = getBasicConnectivity(fi)
Np = shape(conns)[0]
plotNetworkLayout(xs,ys,conns,ls,Np)
The edges have the following attributes:
In [8]:
#print network parameters--these cannot be changed once the network is instantiated
print "pipe L D Mr #grid cells"
for k in range(n0.Nedges):
print "%d %.2f %.2f %.3f %d"%(k, n0.Ls[k], n0.Ds[k],n0.Mrs[k], n0.Ns[k])
In [9]:
#print pipe initial conditions--these can be modified using setIC
#n0.showCurrentData()
Now look at junction information
In [10]:
#look at specification for boundary value types for junction1s--these cannot be changed once network is instantiated
n0.showExternalBoundaries()
Now let's check out the runtime parameters--these cannot be changed once network is instantiated
In [11]:
print "Simulation run time is %.2f s, number of time steps is %d, and pressure wavespeed is %.2f m/s"%(n0.T, n0.M, n0.a[0])
Now we're going to modify initial conditions and run a simulation
In [12]:
#PyPipe_ps?
x = np.linspace(0,n0.Ls[0],n0.Ns[0])# make a vector with as many entries as grid cells in pipe 0
Af =(n0.Ds[0]**2)*pi/4 #full cross sectional area in pipe 0
A0 = []#we're going to store a list of these initial conditions to compare with final results later
Q0 = []
p0 = PyPipe_ps(n0.Ns[0], n0.Ds[0], n0.Ls[0], n0.M, n0.a[0])
print (n0.Ds[0]**2)*pi/4.
Af =(n0.Ds[k]**2)*pi/4
print p0.AofH(3,False)
#assume pipes are 3/4 full with moving water
# reflective boundary is like we've closed valves at all exit points at time t=0
A0.append(.75*Af*np.ones(n0.Ns[0]))
Q0.append(0.007*np.ones(n0.Ns[0]))
plot(x,A0[0],label='pipe 0')
n0.setIC(0,A0[0],Q0[0]) #set the initial conditions in pipe 0
#now make other pipes 'empty' (very small cross sectional area, 0 discharge)
for k in range(1,3):
Af =(n0.Ds[k]**2)*pi/4
x = np.linspace(0,n0.Ls[k],n0.Ns[k])
A0.append(.75*Af*np.ones(n0.Ns[k]))
Q0.append(0.007*np.ones(n0.Ns[k]))
plot(x,A0[k],label='pipe %d'%k)
n0.setIC(k,A0[k],Q0[k])
ylabel('A(x,t=0) in pipe k')
legend()
#and we'll plot these initial cross-sectional areas
Out[12]:
In [13]:
#now run a simulation, keeping track of solve time and the initial and final volume of water in the network
#the simulation is run up to time T using the method runForwardProblem(dt)
import time
V0 = n0.getTotalVolume()
dt = n0.T/float(n0.M)
dx = n0.Ls/[float(nn) for nn in n0.Ns]
t0 = time.clock()
n0.runForwardProblem(dt)
tf = time.clock()
Vf = n0.getTotalVolume()
print "Solve time is %.5f s"%(tf-t0)
print "Simulated time is %.5f s"%n0.T
print "change in volume is %e m^3"%(Vf-V0)
**IMPORTANT:** once you instantiate a network, you can only run the method *runForwardProblem(dt)* **ONCE** because only enough memory is allocated to store a single run of M time steps. If you want to re-run a simulation, you need to either
For now we're just going to look at the data from this run using several attributes of the PyNetwork class.
In [14]:
#The PyNetwork function q(k) returns a list [A,Q] of the current data in pipe k
#Below we use this to retrieve and plot the final cross sectional area
# along length of each pipe and compare with our initial conditions
fig, ax1= plt.subplots(figsize=(15,10), nrows = n0.Nedges)
fig, ax2= plt.subplots(figsize=(15,10), nrows = n0.Nedges)
for k in range(n0.Nedges):
N =n0.Ns[k]
x = linspace(0,n0.Ls[k], N)
q = n0.q(k)#this is the current data (A,Q)
A = q[0:N]#first N entries are values of A in each cell)
Q = q[N:]#next N entries are values of Q in each cell
ax1[k].plot(x,A0[k],'g',label='pipe %d initial'%k)#we stored this above
ax1[k].plot(x,A,'k',label='pipe %d final'%k)
ax1[k].legend()
ax1[k].set_ylabel('A(x)')
ax1[k].set_xlabel('x')
ax2[k].plot(x,Q0[k],'g',label='pipe %d initial'%k)#we stored this above
ax2[k].plot(x,Q,'k',label='pipe %d final'%k)
ax2[k].set_ylabel('Q(x)')
ax2[k].set_xlabel('x')
ax2[k].legend()
The vector of the entire history of (A,Q) at each time step is stored and can be accessed using the PyNetwork function q_hist(i)
However, this requires some slightly complicated indexing (see water_hammer.ipynb for an example)
Below we use two shortcut functions to look at pressure head $H = \bar{p}/(\rho g)$ as a function of space and time in each pipe.
In [15]:
#first look at a time series of pressure data near the end of each pipe
t = linspace(0,n0.T, n0.M+1)
for k in range(n0.Nedges):
K = n0.Ns[k]-2 #index of cell we want to look at
Ht = n0.pressureTimeSeries(k,K)#this function returns H in cell K of pipe k, at each time step
plot(t,Ht, label = 'pipe %d'%k)
plot(t,n0.Ds[k]*np.ones(len(t)),'k:')
legend(loc = 'upper left')
xlabel('t (s)')
ylabel('H (m)')
print 'dashed line denotes pipe crown'
In [16]:
#now look at pressure head H as a function of x in pipe 2 at time step m, for m = 0,100, ...M
# import a nice colormap to do this
from matplotlib import cm
import matplotlib.colors as colors
cNorm = colors.Normalize(vmin=0, vmax=n0.M)
sMap = cm.ScalarMappable(norm=cNorm, cmap=cm.get_cmap('Blues') )
figure(figsize(10,6))
for m in range(0,n0.M,100):
Hx = n0.pressureSpaceSeries(2,m)#this returns H as a function of x in pipe 2 at time step m
x = np.linspace(0,n0.Ls[2],n0.Ns[2])
plot(x,Hx, lw = 2, color = sMap.to_rgba(m), label = 't=%.2f'%(dt*m))
xlabel('x (m)')
ylabel('H (m) in pipe 2')
xlim(0,n0.Ls[2])
legend(loc = 'upper left')
Out[16]:
Ok, we've now seen how to instantiate a network class, look at some attributes, run a simulation, and extract the data. However, it is possible (likely!) that you will want to change runtime parameters like the simulation time, the boundary condition type, the boundary values, and more. Below we will use writeit.py to create new .inp and .config files with desired parameters. Then we can instantiate another network and run it!
In [17]:
Np = n0.Nedges #number of pipes
T = 20 #run time in seconds
a = 150 #pressure wave speed in m/s (this sets the Preissman slot width)
#pipe stuff
Ls = [100,50,25] #try new pipe lengths
Ns = [int(l) for l in Ls] # set dx = 1 meter by assigning 1 grid cell per meter of length
Ds = [0.1,0.05,0.05] #try new pipe diameters
Mrs = n0.Mrs #keep Manning coeff the same
h0s = [5,5,5] #constant initial value of water height (if you want nonconstant, specify later with setIC)
q0s = [0.007, 0.007,0.007] #constant initial discharge Q (if you want nonconstant, specify later with setIC)
#runtime stuff
dx = [Ls[i]/Ns[i] for i in range(Np)]
M = int(T*a/(max(dx)*.8))#set time steps based on dx to assure CFL condition (may need to adjust up if slopes are steep)
#junction stuff
elevs = [0,0,0,0] #make it flat this time
jt = n0.nodeTypes #junction type (can also write this out as a list--in this case, [1,3,1,1]
# note that below, r, bt, and bval only matter for the nodes with one incoming pipe--nodes 0,1 and 2
r = [0,0,1,0] #reflection: we're doing specified bounary at nodes 0 and 3, reflection at node 2
bt = [1,0,0,2] #boundary type (only matters if r=0). Here we're specifying Q at node 0 and orifice outflow at node 2
bv = [0.007,0,0,n0.Ds[2]/2] #node 0: Q = 0.007; node 3: orifice opening height =D/2
In [18]:
#specify .inp file with network connectivity
oldinp = "../indata/mynetwork2.0.inp"
#new prefix for .inp and .config files that will be used for actual runtime parameters
fn = "../indata/myfile_new"
#and write the files!
(fi2, fc2) = rewritePipes(fn,oldinp, Ns, Ls, Mrs, Ds, jt, bt, bv, r, h0s, q0s, T, M, a, elevs)
In [19]:
print "name of new files is %s and %s "%(fi2, fc2)
In [20]:
# now make a new network
n1 = PyNetwork(fi2,fc2,mtype)
In [21]:
n1.showExternalBoundaries()
T = n1.T
M = n1.M
dt = T/float(M)
In [22]:
V0 = n1.getTotalVolume()
t0 = time.clock()
n1.runForwardProblem(dt)
tf = time.clock()
Vf = n1.getTotalVolume()
print "Solve time is %.5f s"%(tf-t0)
print "Simulated time is %.5f s"%n0.T
print "change in volume is %e m^3"%(Vf-V0)
print "note we don't expect volume to be constant since we have orifice and outflow boundaries"
In [28]:
t = linspace(0,n1.T, n1.M+1)
for k in range(0,n1.Nedges):
K = n1.Ns[k]/2 #index of cell we want to look at (midway in each pipe)
Ht = n1.pressureTimeSeries(k,K)#this function returns H in cell K of pipe k, at each time step
plot(t,Ht, label = 'pipe %d'%k)
plot(t,n1.Ds[k]*np.ones(len(t)),'k:')
legend(loc = 'lower right')
xlabel('t (s)')
ylabel('H (m)')
print 'dashed line denotes pipe crown'
Parting thoughts: what to try if your simulation isn't working (a partial list)