From GeoMod
from visual import *
from raster_map import *
from visual.graph import *
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 = 100.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 = 1000 * 365.25 * 24.0 * 60.0 * 60.0
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
dz_crust = dz * ones((nrows,ncols)) #thickness of crust
'''Temperature array'''
T = AT * ones((nrows, ncols)) #intial conditons
'''boundary temp'''
T[:,0]=AT # bottom row
T[:,nrows-1]=AT # top row
T[0,:]=AT #left col
T[ncols- 1,:]=AT #right col
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*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), 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 and cr_base'''
embayment = layer_raster(Z, 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)
'''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:
tt += dt
nt += 1
M = dx * dx * dz_crust[1:nrows-1, 1:ncols-1] * density[1:nrows-1, 1:ncols-1]
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
'''update z.data'''
Z.data[1:nrows-1,1:ncols-1] = Z.data[1:nrows-1,1:ncols-1] + Exp
'''Erosion'''
Z_new = Z.data[1:nrows-1, 1:ncols-1] + S * (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
'''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
'''update'''
Z.data[1:nrows-1,1:ncols-1] = Z_new
#Z.update_line_3d()
'''output'''
topo_cross = Z.data[21,:]-(32000. * ones((ncols)))
plot_line(topo_cross)
if nt % 100 == 0:
print "exp=", dz_crust[21,:]
print "position", ball.pos, r, c
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()
