Gw-2d-map
From GeoMod
Requires the following 3 files to be saved in the same directory:
- Raster Class: direct link raster_map.py
- Scene flying controls: direct link uFly.py
- User controls: direct link uControls.py
Then copy and past the following code or directly download the file: gw-2d-map.py
gw-2d-map.py
from visual import *
from raster_map import *
from uControls import *
from uFly import *
class well:
def __init__(self, pos, bot_elev, Qrate=0.0, screen_top=None,screen_bot=None, max_pump=-2500., max_rech=2500):
'''Draw well cylinder
pos vector is the position of the top of the well
'''
self.pos = pos
self.Qrate = Qrate
self.max_pump = max_pump
self.max_rech = max_rech
self.isactive = False
self.slider_scale = 10.0
self.casing = cylinder(pos=self.pos, axis=(0,0, bot_elev-self.pos.z), color = (0,1,0))
#self.button = sphere(pos=self.pos, radius = 2.*self.casing.radius, color=(0,1,0))
#self.offbutton = uSwitch
def move_well(self, pos): #move well to new position
self.pos = pos
self.casing.pos = pos
def activate(self, well_slider, val=None):
self.isactive = True
self.casing.color = vector(1,0,0)
'''set control slider values'''
well_slider.title.text = "Pumping Rate \n"+ str(self.pos)
well_slider.min_val = self.max_pump
well_slider.max_val = self.max_rech
if val:
well_slider.set_value(val)
else:
well_slider.set_value(self.Qrate)
def deactivate(self, well_slider):
self.isactive = False
self.casing.color = vector(0,1,0)
well_slider.title.text = "Pumping Rate \n"+"No well selected"
well_slider.min_val = -1
well_slider.max_val = 1
well_slider.value = 0
scene.background = color.white #vector(1,0.9,0.7)
scene.lights[0] = vector(0,0,1)
scene.lights[1] = -scene.lights[1]
#set up windows
w = 704 +4+4
h = 576 +24+4
scene.width=w
scene.height=h
scene.x = 0
scene.y = 0
'''Note that the boundary conditions are by default constant head '''
''' set dimensions of grid'''
nrows = 40
ncols = 40
dx = 10.0
''' wells '''
wells = []
wells.append(well(vector(dx*20., -dx*20., 70.0), bot_elev = 0.0, Qrate=3000.0))
wells.append(well(vector(dx*10.0, -dx*20., 70.0), 0.0, -3000.0))
Q = zeros((nrows,ncols)) #array used in head calculations
#Q[20,20] = -2500.0
'''set parameters'''
K = 0.01
porosity = 0.25
dt = 0.1
S = K * dt / (2.0*porosity)
'''Create arrays and raster_map's for saturated aquifer with the top of
the aquifer being the water table (wt), and the bottom of the aquifer
being at 0.0 (aq_base).'''
wt = raster_map( 60.0*ones((nrows,ncols)), dx=dx)
aq_base = raster_map(0.01* ones((nrows,ncols)), dx=dx)
'''Create aquifer as a layer_raster with the wt and aq_base'''
aquifer = layer_raster(wt, aq_base)
'''Draw map of aquifer'''
#color_scale = color_map(40.0, 60.0, 1.0)
color_scale = color_map(55.0, 65.0, 1.0, colmin=vector(1,0.,0), colmax=vector(0.,0.,1))
aquifer.bot.draw_boundary()
aquifer.top.draw_boundary()
'''MAIN SWITCHES'''
l_stop = False
l_pause = False
l_line_3d = False
l_draw_layer = False
l_draw_vectors = True
l_draw_surface = False
l_draw_particles = True
l_draw_trails = True
l_show_trails = True
'''draw lines for aquifer top with the color_map'''
if l_line_3d:
aquifer.top.line_3d(color_map=color_scale, radius=1.0)
#aquifer.top.draw_map(color_map=color_scale, bzmin=0.0)
#aquifer.bot.line_3d()
#aquifer.bot.surface_3d()
if l_draw_surface:
aquifer.top.surface_3d(color_map=color_scale)
'''draw layer surfaces with the color_map (may not be stable'''
#aquifer.draw_layer_surfaces(color_map_top=color_scale)
'''draw whole layer with colors determined by the elevation of the top layer'''
if l_draw_layer:
aquifer.draw_layer(color_map=color_scale)
if l_draw_vectors:
aquifer.top.draw_vectors(vscale=100.)
'''create control window - lurbano'''
cwin = display(title="Controls", width=200, height=h,
x=w, y=0)
'''create a slider to control wells'''
print "length", len(wells)
active_well = None
well_slider = slider(init_val=0.0, title="Pumping Rate",scale=1.0,
min_val = int(-1),
max_val = int(1) )
if len(wells) > 0:
active_well = wells[0]
active_well.activate(well_slider)
## self.slider = slider(pos=self.pos+vector(0,0,10.), init_val=Qrate, \
## min_val=self.max_pump, max_val=self.max_rech, scale=self.slider_scale)
## self.slider.val_label.pos = self.slider.cyl.axis * 0.5
pick = None
l_cwin = True #control window is initially active
l_scene = False #scene window is intitally inactive
if l_line_3d:
init = 1
else:
init = 0
line_3d_but = uSwitch(title="Lines", pos=vector(0,-0.9), radius=0.1, init_val=init)
if l_draw_layer:
init = 1
else:
init = 0
draw_layer_but = uSwitch(title="Blocks", pos=vector(-0.2,-0.9), radius=0.1, init_val=init)
if l_draw_surface:
init = 1
else:
init = 0
draw_surface_but = uSwitch(title="Surface", pos=vector(0.2, -0.9), radius=0.1, init_val=init)
if l_draw_vectors:
init = 1
else:
init = 0
draw_vectors_but = uSwitch(title="Vectors", pos=vector(0.2, -1.3), radius=0.1, init_val=init)
if l_draw_particles:
init = 1
else:
init = 0
draw_particles_but = uSwitch(title="Particles", pos=vector(-0.2, -1.3), radius=0.1, init_val=init)
if l_show_trails:
init = 1
else:
init = 0
draw_trails_but = uSwitch(title="Trails", pos=vector(0.0, -1.3), radius=0.1, init_val=init)
reset_particles_but = uSwitch(title="Reset\nParticles", pos=vector(-0.2, -1.7), radius=0.1, init_val=0)
reset_trails_but = uSwitch(title="Reset\nTrails", pos=vector(0.0, -1.7), radius=0.1, init_val=0)
mapview_but = uSwitch(title="Map View", pos=vector(0,-0.5), radius = 0.1)
map_view = vector(0,0,-1)
if l_pause:
init = 1
else:
init = 0
pause_but = uSwitch(title="Pause", pos=vector(-0.2,-0.5), radius = 0.1, init_val=init)
#print "display", display._selected()
'''initial particles'''
if l_draw_particles:
scene.select()
for i in range(1, nrows-1, 2):
for j in range(1, ncols-1, 2):
aquifer.top.add_particle(pos=vector(dx*i, -dx*j), l_trail=l_draw_trails)
if not l_show_trails:
aquifer.top.particles[-1].hide_trail()
'''Run model'''
# infinite loop
ct = 0
while 1:
'''Mouse events'''
if scene.mouse.events:
m1 = scene.mouse.getevent()
l_scene = True
l_cwin = False
for i in wells:
if m1.pick == i.casing and m1.release:
if active_well:
active_well.deactivate(well_slider)
active_well = i
active_well.activate(well_slider)
if m1.pick == i.casing and m1.drag:
l_stop = True
if active_well:
active_well.deactivate(well_slider)
active_well = i
active_well.activate(well_slider)
drag_pos = scene.mouse.project(normal=vector(0,0,1), point=active_well.pos) #m1.pickpos
pick = m1.pick
scene.cursor.visible = 0
if m1.drop:
if pick == active_well.casing:
(r, c) = aquifer.top.get_rc(pick.pos)
active_well.move_well(vector(aquifer.top.get_x(c), aquifer.top.get_y(r), pick.pos.z))
pick = None # end dragging
scene.cursor.visible = 1 # cursor visible
l_stop = False
if cwin.mouse.events:
m1 = cwin.mouse.getevent() # obtain drag or drop event
l_cwin = True
l_scene = False
if m1.drag and (m1.pick == well_slider.but):
drag_pos = m1.pickpos
pick = m1.pick
#print m1.pick.pos
cwin.cursor.visible = 0 # make cursor invisible
l_stop = True
elif m1.drop:
pick = None # end dragging
cwin.cursor.visible = 1 # cursor visible
l_stop = False
if (m1.pick == line_3d_but.button) and m1.release:
line_3d_but.switch()
if line_3d_but.value == 1 and not l_line_3d:
scene.select()
aquifer.top.line_3d(color_map=color_scale)
l_line_3d = True
if line_3d_but.value == 1:
aquifer.top.line_3d_show()
else:
aquifer.top.line_3d_hide()
if (m1.pick == draw_layer_but.button) and m1.release:
draw_layer_but.switch()
if draw_layer_but.value == 1 and not l_draw_layer:
scene.select()
aquifer.draw_layer(color_map=color_scale)
l_draw_layer = True
if draw_layer_but.value == 1:
aquifer.layer_show()
else:
aquifer.layer_hide()
if (m1.pick == draw_vectors_but.button) and m1.release:
draw_vectors_but.switch()
if draw_vectors_but.value == 1 and not l_draw_vectors:
scene.select()
aquifer.top.draw_vectors(color_map=color_scale)
l_draw_vectors = True
if draw_vectors_but.value == 1:
aquifer.top.vectors_show()
else:
aquifer.top.vectors_hide()
if (m1.pick == draw_particles_but.button) and m1.release:
draw_particles_but.switch()
if draw_vectors_but.value == 1 and not l_draw_particles: #initialize if necessary
scene.select()
aquifer.top.add_particles_at_nodes(2)
l_draw_particles = True
l_draw_trails = False
if draw_particles_but.value == 1:
aquifer.top.particles_show()
else:
aquifer.top.particles_hide()
if (m1.pick == draw_trails_but.button) and m1.release:
draw_trails_but.switch()
if draw_trails_but.value == 1 and not l_draw_trails:
scene.select()
aquifer.top.draw_trails() # create trails
l_draw_trails = True
if draw_trails_but.value == 1:
aquifer.top.trails_show()
else:
#fish = hens
aquifer.top.trails_hide()
if (m1.pick == reset_trails_but.button) and m1.release:
reset_trails_but.switch()
aquifer.top.reset_trails()
if l_draw_trails and draw_trails_but.value==1:
draw_trails_but.switch()
#l_draw_trails = False
l_draw_trails = False
reset_trails_but.switch()
if (m1.pick == reset_particles_but.button) and m1.release:
reset_particles_but.switch()
scene.select()
aquifer.top.delete_all_particles()
if l_draw_particles:
draw_particles_but.switch()
#l_draw_trails = False
l_draw_particles = False
reset_particles_but.switch()
if (m1.pick == pause_but.button) and m1.release:
pause_but.switch()
if pause_but.value == 1:
l_pause = True
else:
l_pause = False
if (m1.pick == mapview_but.button) and m1.release:
'''Rotate scene to top down view'''
mapview_but.switch()
print "scene.forward", scene.forward
af = vector(scene.forward.x, scene.forward.y)
print "af", af
angle = degrees(af.diff_angle(vector(0,1,0)))
if angle != 0.0:
ac = cross(af, vector(0,1,0))
print "ac", ac
rot_sign = ac.z / abs(ac.z)
print "rot_sign", rot_sign
print "angle", angle
#drate = 100.0
#nsteps = angle / drate
rotate_view(rot_sign*angle, vector(0,0,1), int(angle*2.0))
print "scene.forward", scene.forward
af = vector(0.0, scene.forward.y, scene.forward.z)
print "af", af
print "map_view", map_view
angle = degrees(af.diff_angle(map_view))
if angle != 0.0:
print "angle", angle
ac = cross(af, map_view)
print "ac", ac
rot_sign = ac.x / abs(ac.x)
print "rot_sign", rot_sign
rotate_view(rot_sign*angle, vector(1,0,0), int(angle*2.0))
#scene.forward = map_view
mapview_but.switch()
if pick and l_cwin:
new_pos = cwin.mouse.project(normal=(0,0,1))
if new_pos != drag_pos:
pick.pos += new_pos - drag_pos
drag_pos = new_pos
if pick == well_slider.but:
pick.pos = well_slider.but_move(pick.pos)
active_well.Qrate = well_slider.value
if pick and l_scene:
if pick == active_well.casing:
new_pos = scene.mouse.project(normal=vector(0,0,1), point=active_well.pos)
if new_pos != drag_pos:
pick.pos += new_pos - drag_pos
drag_pos = new_pos
print pick.pos
if pick == active_well.casing:
#pick.pos = well_slider.but_move(pick.pos)
active_well.move_well(pick.pos)
'''finite difference model'''
if not l_stop and not l_pause:
ct += 1
print "timestep = ", ct
'''Add well pumping rates'''
Q = Q * 0.0
for i in wells:
(r,c) = wt.get_rc(i.pos)
Q[r,c] = Q[r,c] + i.Qrate
print r, c, i.Qrate
# Note: wt.data = head_old
head_new = wt.data[1:nrows-1,1:ncols-1] + \
S * (power(wt.data[:nrows-2,1:ncols-1],2) + \
power(wt.data[2:,1:ncols-1],2) + \
power(wt.data[1:nrows-1,:ncols-2],2) + \
power(wt.data[1:nrows-1,2:],2) - \
4 * power(wt.data[1:nrows-1,1:ncols-1],2) + \
Q[1:nrows-1,1:ncols-1] *dt/porosity)
dh = abs(head_new - wt.data[1:nrows-1,1:ncols-1])
#print max(max(dh))
#print head_new[19,19]
#print head_new[18:21,18:21]
wt.data[1:nrows-1,1:ncols-1] = head_new
'''update image'''
if l_draw_vectors:
aquifer.top.update_vectors()
if line_3d_but .value == 1:
aquifer.top.update_line_3d()
#aquifer.top.update_draw_map()
if draw_surface_but.value == 1:
aquifer.top.update_surface_3d()
#aquifer.update_layer_surfaces()
if draw_layer_but.value == 1:
aquifer.update_layer()
'''move particles'''
if len(aquifer.top.particles) > 0:
if not l_draw_vectors: #don't want to repeat this if we don't have to.
aquifer.top.calc_surface_vectors()
for n, i in enumerate(aquifer.top.particles):
vx = aquifer.top.interpolate_onarray(i.particle.pos, aquifer.top.vx)
vy = aquifer.top.interpolate_onarray(i.particle.pos, aquifer.top.vy)
dx = vx * dt *500.
dy = vy * dt * 500.
#print "particle dx,dy", dx, dy
new_pos = i.particle.pos + vector(dx, dy, 0.0)
'''remove particle if out of bounds'''
#print n, new_pos
i.particle.color = vector(0,1,0)
if not aquifer.top.inbounds(new_pos, padding=aquifer.top.dx):
#print "out of bounds"
i.particle.color = vector(0,1,0)
aquifer.top.delete_particle(n)
#i.visible = 0
#del_particle = aquifer.top.particles.pop(n)
else:
#print "in bounds"
new_pos.z = aquifer.top.interpolate_linear(new_pos)
i.move_particle(new_pos)
i.particle.color = vector(1,0,0)
'''remove particle if it hits a pumping well'''
for j in wells:
#print "distance", j.pos, i.pos, mag(j.pos-i.pos)
#print i.pos
#print j.pos
if j.Qrate < 0.0 and mag(vector(j.pos.x, j.pos.y)-vector(i.particle.pos.x, i.particle.pos.y)) < aquifer.top.dx/2.0:
#i.particle.visible = 0
#print "magnitude=", mag(vector(j.pos.x, j.pos.y)-vector(i.particle.pos.x, i.particle.pos.y))
aquifer.top.delete_particle(n)
#del_particle = aquifer.top.particles.pop(n)
break
