Final Code
From GeoMod
My main program goes:
from visual import *
from raster_import5a import *
from Classes import *
w = 704+4+4
h = 576+24+4
scene.x = 0
scene.y = 0
scene.width = w
scene.height = h
#import topography
topog = raster_map("priver100_cut.txt")
#topog.draw_boundary()
topog.draw_map(0, 17, 0, 31, 3.0, 1)
#initialize water levels
topog.init_waterlevels(61.5)
#timestep
dt = 0.00000005
for i in range(1500):
if (i % 10 == 0):
print "TIMESTEP = ", i
#add water to river (boundary condition)
topog.wateradd(0.0008*dt, 2, 17)
#propogate water downstream (diffusion)
topog.diffuse2(dt)
topog.wateradd(-0.8*dt, 14, 15)
if (i % 10 == 0):
topog.show_change()
My classes including the diffusion equations reads as:
from Numeric import *
from string import *
from visual import *
class raster_map:
def __init__(self, rasterfile): #datarr, xmin, ymax, dx, nrows, ncols):
infile = open(rasterfile,"r")
ct = 0
line = infile.readline()
while line:
ct += 1
l = split(strip(line))
if ct == 1:
ncols = int(l[1])
if ct == 2:
nrows = int(l[1])
if ct == 3:
xmin = float(l[1])
if ct == 4:
ymin = float(l[1])
if ct == 5:
dx = float(l[1])
if ct == 7:
nodata = l[1]
if ct >= 7:
if ct == 7:
topo = zeros((nrows*ncols,), Float)
t_ct = -1
for i in l:
t_ct += 1
if i == nodata: #set no-data values to zero
topo[t_ct] = 0.0
else:
topo[t_ct] = float(i)
line = infile.readline()
infile.close
ymax = ymin + (dx * (nrows-1)) # top left corner
'''datarr is a 2d array '''
self.xmin = xmin
self.ymax = ymax
self.dx = dx
self.nrows = nrows
self.ncols = ncols
self.data = resize(topo, (nrows, ncols))
def get_val(self, pos):
r, c = self.get_rc(pos)
return self.data[r,c]
def get_rc(self, pos):
r = int(((self.ymax - pos.y + self.dx/2.0) ) /self.dx )
c = int(((pos.x - self.xmin + self.dx/2.0) ) /self.dx)
return (r, c)
def draw_map(self, imin=0, imax=-99999, jmin=0, jmax=-99999, scale = 1.0, center = 0):
''' draw a specific region of the raster map'''
self.scale = scale
default_val = -99999
if imax == default_val:
imax = self.nrows
if jmax == default_val:
jmax = self.ncols
tdata = self.data[imin:imax,jmin:jmax]
#draw in topography
tmax = max(max(tdata))
tmin = min(min(tdata))
for i in range(imin,imax):
for j in range(jmin,jmax):
ipos = self.xmin + j* self.dx
jpos = self.ymax - i* self.dx
kpos = self.data[i,j]
colval = (kpos - tmin) / (tmax - tmin)
box(pos=(ipos,jpos, 0.5*kpos*scale),
length = self.dx, height = self.dx, width = self.data[i,j]*scale,
color = (colval,colval,colval))
#center scene on area
if center == 1:
xc = self.xmin + (self.dx * ((jmin+jmax)/2))
yc = self.ymax - (self.dx * ((imin+imax)/2))
scene.center = vector(xc, yc, 0)
def draw_boundary(self): # Drawing boundary outside the map
bound = curve(color=color.blue)
bound.append(pos=(self.xmin, self.ymax,self.data[0,0]))
bound.append(pos=(self.xmin+self.dx*self.ncols, self.ymax,self.data[0,0]))
bound.append(pos=(self.xmin+self.dx*self.ncols,
self.ymax-self.dx*self.nrows,self.data[0,0]))
bound.append(pos=(self.xmin,
self.ymax-self.dx*self.nrows,self.data[0,0]))
bound.append(pos=(self.xmin, self.ymax,self.data[0,0]))
def init_waterlevels(self, initial_elevation):
'''intitial_elevations is a single value: assuming initial water level is flat'''
self.waterlevels = float(initial_elevation) * ones([self.nrows, self.ncols])
self.waterdepths = self.waterlevels - self.data
for i in range (self.nrows):
for j in range (self.ncols):
if self.waterdepths[i,j] < 0.0:
self.waterdepths[i,j] = 0.0
self.waterlevels[i,j] = self.data[i,j]
# draw water where the surface is covered with water
self.water = []
col = 0
for i in range(self.nrows):
for j in range(self.ncols):
if self.waterlevels[i,j] > self.data[i,j]:
ipos = self.xmin + j* self.dx
jpos = self.ymax - i* self.dx
kpos = self.waterlevels[i,j]
col=self.data[i,j] / self.waterlevels[i,j]
self.water.append(box(pos=(ipos,jpos, 0.5*kpos*self.scale),
length = self.dx, height = self.dx, width = self.waterlevels[i,j]*self.scale,
color = (0,col,1)))
else:
ipos = self.xmin + j* self.dx
jpos = self.ymax - i* self.dx
kpos = 0.0
self.water.append(box(pos=(ipos,jpos, 0.5*kpos*self.scale),
length = self.dx, height = self.dx, width = 0,
color = (0,col,1)))
def wateradd (self, quant, i, j): # Adding water to one cell
self.waterlevels[i,j] += quant / (self.dx*self.dx)
def show_change (self): #draw change in level of water
c = -1
for i in range(self.nrows):
for j in range(self.ncols):
c += 1
if self.waterlevels[i,j] > self.data[i,j]:
kpos = (self.waterlevels[i,j]+self.data[i,j])/2.0
self.water[c].pos.z = kpos * self.scale
self.water[c].width = (self.waterlevels[i,j]-self.data[i,j])*self.scale
else:
kpos = 0.0
self.water[c].pos.z = 0.0
self.water[c].width = 0.0
def diffuse2(self, dt, mann=0.035, threshold=75.0):
self.mann = mann
self.dt = dt
new_wl = self.waterlevels
for i in range(1, self.nrows-1):
for j in range(1, self.ncols-1):
V_ym = 0.0
V_yp = 0.0
'''The equations would be followed only if the waterlevel of present cell is lower than the cell in next row
and if the adjacent cell has water. Similar conditions are followed everytime as we see the forward face, left and right sides'''
#For trailing face
if self.waterlevels[i,j]<self.waterlevels[i-1,j]:
if self.waterlevels[i-1,j] > self.data[i-1,j]:
slope = (self.waterlevels[i-1,j] - self.waterlevels[i,j]) / self.dx #calculating water surface elevation
self.ws_elev = self.waterdepths + self.data
hmax = self.hmax_define(float(self.ws_elev[i,j]), float(self.ws_elev[i-1,j]))
zmax = max(self.data[i,j], self.data[i-1,j])
dh = self.equations(zmax, hmax, slope) #Call function equations
if self.waterdepths[i,j] < dh:
dh = self.waterdepths[i,j]
new_wl[i,j] = self.waterlevels[i,j] + dh
kpos = new_wl[i,j]
else:
new_wl[i,j] = self.waterlevels[i,j]
else:
slope = (self.waterlevels[i,j] - self.waterlevels[i-1,j]) / self.dx #calculating water surface elevation
self.ws_elev = self.waterdepths + self.data
hmax = self.hmax_define(float(self.ws_elev[i,j]), float(self.ws_elev[i-1,j]))
zmax = max(self.data[i,j], self.data[i-1,j])
dh = self.equations(zmax, hmax, slope) #Call function equations
if self.waterdepths[i-1,j] < dh:
dh = self.waterdepths[i-1,j]
new_wl[i-1,j]=self.waterlevels[i-1,j] + dh
kpos = new_wl[i-1,j]
#For Forward Face
if self.waterlevels[i,j]<self.waterlevels[i+1,j]:
if self.waterlevels[i+1,j] > self.data[i+1,j]:
slope = (self.waterlevels[i+1,j] - self.waterlevels[i,j]) / self.dx #calculating water surface elevation
self.ws_elev = self.waterdepths + self.data
hmax = self.hmax_define(float(self.ws_elev[i,j]), float(self.ws_elev[i+1,j]))
zmax = max(self.data[i,j], self.data[i+1,j])
dh = self.equations(zmax, hmax, slope) #Call function equations
if self.waterdepths[i,j] < dh:
dh = self.waterdepths[i,j]
new_wl[i,j] = self.waterlevels[i,j] + dh
kpos = new_wl[i,j]
else:
new_wl[i,j]=self.waterlevels[i,j]
else:
slope = (self.waterlevels[i,j] - self.waterlevels[i+1,j]) / self.dx #calculating water surface elevation
self.ws_elev = self.waterdepths + self.data
hmax = self.hmax_define(float(self.ws_elev[i,j]), float(self.ws_elev[i+1,j]))
zmax = max(self.data[i,j], self.data[i+1,j])
dh = self.equations(zmax, hmax, slope) #Call function equations
if self.waterdepths[i+1,j] < dh:
dh = self.waterdepths[i+1,j]
new_wl[i+1,j]=self.waterlevels[i+1,j] + dh
kpos = new_wl[i+1,j]
#For Right Sides
if self.waterlevels[i,j]>self.waterlevels[i,j+1]:
if self.waterlevels[i,j+1] > self.data[i,j+1]:
slope = (self.waterlevels[i,j] - self.waterlevels[i,j+1]) / self.dx
self.ws_elev = self.waterdepths + self.data
hmax = self.hmax_define(float(self.ws_elev[i,j]), float(self.ws_elev[i,j+1]))
zmax = max(self.data[i,j], self.data[i,j+1])
dh = self.equations(zmax, hmax, slope) #Call function equations
if self.waterdepths[i,j+1] < dh:
dh = self.waterdepths[i,j+1]
new_wl[i,j+1]=self.waterlevels[i,j+1] + dh
kpos = new_wl[i,j+1]
else:
new_wl[i,j+1]=self.waterlevels[i,j+1]
else:
slope = (self.waterlevels[i,j] - self.waterlevels[i,j+1]) / self.dx #calculating water surface elevation
self.ws_elev = self.waterdepths + self.data
hmax = self.hmax_define(float(self.ws_elev[i,j]), float(self.ws_elev[i,j+1]))
zmax = max(self.data[i,j], self.data[i,j+1])
dh = self.equations(zmax, hmax, slope) #Call function equations
if self.waterdepths[i,j+1] < dh:
dh = self.waterdepths[i,j+1]
new_wl[i,j]=self.waterlevels[i,j] + dh
kpos = new_wl[i,j]
#For Left Side
if self.waterlevels[i,j] > self.waterlevels[i,j-1]:
if self.waterlevels[i,j-1] > self.data[i,j-1]:
slope = (self.waterlevels[i,j] - self.data[i,j-1]) / self.dx
self.ws_elev = self.waterdepths + self.data
hmax = self.hmax_define(float(self.ws_elev[i,j]), float(self.ws_elev[i,j-1]))
zmax = max(self.data[i,j], self.data[i,j-1])
dh = self.equations(zmax, hmax, slope) #Call function equations
if self.waterdepths[i,j-1] < dh:
dh = self.waterdepths[i,j-1]
new_wl[i,j-1]=self.waterlevels[i,j-1] + dh
kpos = new_wl[i,j-1]
else:
new_wl[i,j-1]=self.waterlevels[i,j-1]
else:
slope = (self.waterlevels[i,j] - self.waterlevels[i,j-1]) / self.dx #calculating water surface elevation
self.ws_elev = self.waterdepths + self.data
hmax = self.hmax_define(float(self.ws_elev[i,j]), float(self.ws_elev[i,j-1]))
zmax = max(self.data[i,j], self.data[i,j-1])
dh = self.equations(zmax, hmax, slope) #Call function equations
if self.waterdepths[i,j-1] < dh:
dh = self.waterdepths[i,j-1]
new_wl[i,j]=self.waterlevels[i,j] + dh
kpos = new_wl[i,j]
self.waterlevels = new_wl
for i in range(self.nrows):
for j in range (self.ncols):
if (self.waterlevels[i,j] < self.data[i,j]):
self.waterlevels[i,j] = self.data[i,j]
## self.waterlevels[i,j] = new_wl[i,j]
## self.waterlevels[i-1,j] = new_wl[i-1,j]
## self.waterlevels[i+1,j] = new_wl[i+1,j]
## self.waterlevels[i,j-1] = new_wl[i,j-1]
## self.waterlevels[i,j+1] = new_wl[i,j+1]
self.waterdepths = new_wl - self.data
if self.waterlevels[i,j] > threshold:
self.waterlevels[i,j] = self.waterlevels[i,j] - (self.waterlevels[i,j]-threshold)
#self.show_change()
#print "water depths"
#print self.waterdepths
def equations (self, zmax, hmax, slope): # Defining all equations
diff = (hmax - zmax) #zmax is the maximum bed surface elevation between two cells
cs_area = diff * self.dx #Calculating Cr0ss-sectional area
w_perimeter = self.dx + (2*diff) #Calculating Wetted Perimeter
h_radius = (cs_area * self.dx ) / w_perimeter #Calculating Hydraulic Radius
#print "h_radius = ", h_radius
flowrate_u = (cs_area * pow(h_radius, 2.0/3.0) * pow(abs(slope), 0.5)) / self.mann
V_ym = flowrate_u * self.dt #Calculating volume of water coming into cell
dh = V_ym / self.dx *self.dx #Calculating the change in water level
return dh
def hmax_define (self, a, b): #Defining Hmz value
hmax = 0.0 #hmax is the maximum water surface elevation
if a > b:
hmax = a
else:
hmax = b
return hmax
I have used the following DEM of Pearl River: Image:Priver100 cut.txt
The law of conservation of mass/matter ( The Lomonosov-Lavoisier law ) states that the mass of a system of substances is constant, regardless of the processes acting inside the system. An equivalent statement is that matter changes form, but cannot be created or destroyed. This implies that for any chemical process in a closed system, the mass of the reactants must equal the mass of the products. However my model does not deal with this important aspect. For future work, conservation of mass should be included in the model.
