MEP with controls
From GeoMod
In the following code there are 3 sections marked with "CONTROLS" comments. To use:
- copy and paste these section into your code, and
- add the file uControls.py (from the User controls page) to your directory.
from visual import *
from raster_map import *
from visual.graph import *
'''CONTROLS 1/3'''
from uControls import *
'''END CONTROLS 1/3'''
def plot_line(data_array):
glist = []
ct = -1
for i in data_array:
ct+=1
glist.append((ct,i))
c = gcurve(pos=glist)
'''Note that the boundary conditions are by default constant head '''
''' set dimensions of grid'''
nrows = 41
ncols = 41
'''set parameters'''
K = 3.0 # thermal conductivty (W/mK)
dx = 100000.0 #meters
TT = 15.0 + 273.0#assumption: may need to take average of cell
NT = 8400.0 + 273.0# Hot spot temp. Kelvin
dz =32000.0 #meters
A = (dx * dx)
dT = NT-TT
dt = 500.0* 365.25 * 24.0 * 60.0 * 60.0
BT = 335.0 + 273.0 #Base Temperature of crust
#q = K * (dT/dz)
C = 790 # W/mK specific heat capacity of granite
S = K*dt / (dx * dx * C)
AT = (TT + BT) / 2.0 #average temp cell: boundary
'''Erosion'''
dt_E = dt
K_E = 50000.0/ (365.25 * 24.0 * 60.0 * 60.0)
S_E = K_E * dt_E / (dx*dx)
'''ball'''
vel = 0.0
t = 0
xvel = .3 / (365.25 * 24.0 * 60.0 * 60.0) ##meters per year
density = 2600.0 * ones((nrows,ncols)) # density of granite
density_mantle = 3300 * ones((nrows,ncols)) #density of mantle
dz_crust = dz * ones((nrows,ncols)) #thickness of crust
'''Temperature array'''
T = AT * ones((nrows, ncols)) #intial conditons
#print T
Tmap_base = raster_map(zeros((nrows,ncols)), dx=dx)
Tmap = raster_map(T, dx=dx)
tcolor_scale = color_map(300.0, 1000.0, 1.0, colmin=vector(1,0,0), colmax=vector(0,1,0))
Tmap.line_3d(scale=1.0, color_map = tcolor_scale, center=1)
'''Initial conditions: create elevation array'''
Z = raster_map(dz_crust, dx=dx)
#Z = raster_map(dz*ones((nrows,ncols)), dx=dx) #elevation array
scene.center = vector(dx*20.0,dx*20.0,4500.)
'''Create arrays and raster_map's for embayment with the top of
the model being the crust (cr_top), and the bottom of the model
being at 0.0 (cr_base).'''
cr_top = raster_map(dz_crust, dx=dx)
cr_base = raster_map(0.01 * ones((nrows,ncols)), dx=dx)
'''Set Z = cr_top'''
Z = cr_top
'''Create embayment as a layer_raster with the cr_top and cr_base'''
embayment = layer_raster(cr_top, cr_base)
'''Draw map of embayment'''
color_scale = color_map(31990.0, 32050.0, 5.0, colmin=vector(1,0,0), colmax=vector(0,1,0))
embayment.draw_layer(color_map=color_scale)
'''molding a ball'''
ball = sphere(pos=(dx,-20*dx,dz), color=(0,0,1),radius=dx*.25)
'''make pos Q = ball.pos'''
##ball.pos = (dx,-20*dx,-dz-10000)
(r,c)=cr_top.get_rc(ball.pos)
''' CONTROLS 2/3'''
'''create control window - lurbano'''
cwin = display(title="Controls")
stop_button = uSwitch(title="Pause", pos=vector(0.2, -0.9), radius=0.1, init_val=0)
Tslider = slider(init_val=NT, title="Hot Spot \n Temperature \n(K)",scale=1.0,
min_val = 0.,
max_val = 10000. )
TKslider = slider(init_val=K, title="Thermal \n conductivity",scale=1.0,
min_val = 0.1,
max_val = 3. ,
pos = vector(-1,0))
Vslider = slider(init_val=xvel*(365.25 * 24.0 * 60.0 * 60.0) , title="Plume \nVelocity",scale=1.0,
min_val = -.5,
max_val = .5 ,
pos = vector(1,0))
pick = None
'''END CONTROLS 2/3'''
'''set boundary conditions'''
Q = 0.0 * ones((nrows,ncols)) #hotspot
Q[r,c] = K * ((NT-AT)/dz) * A
nt = 0
runtime = 20000000.0 * 365.25 * 24.0 * 60.0 * 60.0
tt = 0.0
##for nt in range(1, nsteps+1):
while tt < runtime:
''' CONTROLS 3/3'''
if cwin.mouse.events:
m1 = cwin.mouse.getevent() # obtain drag or drop event
if (m1.pick == stop_button.button) and m1.release:
stop_button.switch()
if m1.drag and (m1.pick == Tslider.but):
drag_pos = m1.pickpos
pick = m1.pick
cwin.cursor.visible = 0 # make cursor invisible
elif m1.drag and (m1.pick == TKslider.but):
drag_pos = m1.pickpos
pick = m1.pick
cwin.cursor.visible = 0 # make cursor invisible
elif m1.drag and (m1.pick == Vslider.but):
drag_pos = m1.pickpos
pick = m1.pick
cwin.cursor.visible = 0 # make cursor invisible
elif m1.drop:
pick = None # end dragging
cwin.cursor.visible = 1 # cursor visible
l_stop = False
if pick:
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 == Tslider.but:
pick.pos = Tslider.but_move(pick.pos)
NT = Tslider.value
if pick == TKslider.but:
pick.pos = TKslider.but_move(pick.pos)
K = TKslider.value
S = K*dt / (dx * dx * C)
if pick == Vslider.but:
pick.pos = Vslider.but_move(pick.pos)
#print "V", Vslider.value
xvel = Vslider.value/ (365.25 * 24.0 * 60.0 * 60.0)
if stop_button.value == 1: #pause simulation by not doing the rest of model
continue
'''END CONTROLS 3/3'''
tt += dt
nt += 1
M = dx * dx * dz_crust[1:nrows-1, 1:ncols-1] * density[1:nrows-1, 1:ncols-1]
'''boundary temp'''
T[:,0]=T[:,1] # Left column
T[:,ncols-1]=T[:,ncols-2] # Right column
T[0,:]=AT #top row
T[nrows- 1,:]=AT#bottom col
'''Temperature Diffusion Equation'''
temp_new = T[1:nrows-1, 1:ncols-1] + (Q[1:nrows-1, 1:ncols-1]*dt/(C*M))+ \
(S/(density[1:nrows-1, 1:ncols-1]*dz_crust[1:nrows-1, 1:ncols-1]))* \
((T[2:nrows, 1:ncols-1]- T[1:nrows-1, 1:ncols-1])*(dx*(dz_crust[1:nrows-1, 1:ncols-1]+dz_crust[2:nrows, 1:ncols-1])/2.0)- \
(T[1:nrows-1, 1:ncols-1] - T[0:nrows-2, 1:ncols-1])*(dx*(dz_crust[1:nrows-1, 1:ncols-1]+dz_crust[0:nrows-2, 1:ncols-1])/2.0)+ \
(T[1:nrows-1, 2:ncols]-T[1:nrows-1, 1:ncols-1])*(dx*(dz_crust[1:nrows-1, 1:ncols-1]+dz_crust[1:nrows-1, 2:ncols])/2.0)- \
(T[1:nrows-1, 1:ncols-1]-T[1:nrows-1, 0:ncols-2])*(dx*(dz_crust[1:nrows-1, 1:ncols-1]+dz_crust[1:nrows-1, 0:ncols-2])/2.0))
'''Thermal Expansion'''
Exp = (temp_new - T[1:nrows-1,1:ncols-1]) * 9. * 0.0000001 * dz_crust[1:nrows-1,1:ncols-1]
#print Exp.shape
#print dz_crust.shape
dz_crust[1:nrows-1,1:ncols-1] = dz_crust[1:nrows-1,1:ncols-1] + Exp
'''Erosion'''
Z_new = Z.data[1:nrows-1, 1:ncols-1] + S_E * (Z.data[2:nrows, 1:ncols-1] - (4 * Z.data[1:nrows-1, 1:ncols-1]) +
Z.data[0:nrows-2, 1:ncols-1] + Z.data[1:nrows-1, 2:ncols] +
Z.data[1:nrows-1, 0:ncols-2])
ze = Z.data[1:nrows-1,1:ncols-1] - Z_new
'''isostatic rebound'''
dz_i = Z_new * (1.0 - density[1:nrows-1, 1:ncols-1]/density_mantle[1:nrows-1, 1:ncols-1]) - \
Z.data[1:nrows-1,1:ncols-1] * (1.0 - density[1:nrows-1, 1:ncols-1]/density_mantle[1:nrows-1, 1:ncols-1])
'''update z.data'''
Z.data[1:nrows-1,1:ncols-1] = Z_new + dz_i
Z.data[1:nrows-1,1:ncols-1] = Z.data[1:nrows-1,1:ncols-1] + Exp
'''update'''
#Z.data[1:nrows-1,1:ncols-1] = Z_new
#Z.update_line_3d()
'''update'''
T[1:nrows-1,1:ncols-1] = temp_new
'''move plume'''
ball.pos = ball.pos + vector(xvel*dt,0.0,0.0)
(r, c) = cr_top.get_rc(ball.pos)
'''update plume heat flux'''
Q = 0.0 * Q
Q[r,c] = K * ((NT-T[r,c])/dz_crust[r,c]) * dx * dx *5000.0
'''output'''
topo_cross = cr_top.data[21,:]-(32000. * ones((ncols)))
#plot_line(topo_cross)
if nt % 100 == 0:
print #Exp, ze
print #"exp=", dz_crust[21,:]
print "position", ball.pos, r, c
print cr_top.data[21,:]
print #"Q", NT, T[r,c], Q[r,c]*dt/(C*M[r,c])
print #nt, T[21, :], ze[20, :]
Tmap.update_line_3d()
##plot_line(dz_crust[21,:])
embayment.update_layer()
'''update image'''
Tmap.update_line_3d()
