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 = 2200.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
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)
color_scale = color_map(300.0, 400.0, 1.0, colmin=vector(1,0,0), colmax=vector(0,1,0))
Tmap.line_3d(scale=1.0, color_map = color_scale)
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)
##
##'''Create embayment as a layer_raster with the cr 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))
##color_scale = color_map(50.0, 62.0, 1.0, colmin=vector(1,0,0), colmax=vector(0,1,0))
##'''draw lines for embayment top with the color_map'''
##crust.top.line_3d(color_map=color_scale)
##'''draw layer surfaces with the color_map (may not be stable'''
###aquifer.draw_layer_surfaces(color_scale)
##'''draw whole layer with colors determined by the elevation of the top layer'''
embayment.draw_layer(color_map=color_scale)
'''molding a ball'''
ball = sphere(pos=(0,-20*dx,-dz-10000), color=(0,0,1),radius=50000.0)
'''make pos Q = ball.pos'''
ball.pos = (dx,-20*dx,-dz-10000)
(r,c)=cr_top.get_rc(ball.pos)
'''moving ball'''
vel = 0.0
t = 0
xvel = .03 / (365.25 * 24.0 * 60.0 * 60.0) ##meters per year
##if ball.pos.x > - 1:
##
## vel += xvel * dt
## ball.pos.x += vel * dt
'''set boundary conditions'''
Q = 0.0 * ones((nrows,ncols)) #hotspot
Q[r,c] = K * ((NT-AT)/dz) * A
nsteps = 100000
for nt in range(1, nsteps+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
#print Exp.shape
#print dz_crust.shape
dz_crust[1:nrows-1,1:ncols-1] = dz_crust[1:nrows-1,1:ncols-1] \
+ Exp * dz_crust[1:nrows-1,1:ncols-1]
'''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'''
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, :]
Tmap.update_line_3d()
##plot_line(dz_crust[21,:])
embayment.update_layer()
'''update image'''
Tmap.update_line_3d()
