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 = 1.0* 365.25 * 24.0 * 60.0 * 60.0
BT = 335.0 + 273.0 #Base Temperature of crust
q = K * (dT/dz)
C = 3.2 # W/mK specific heat capacity of granite
S = K / (dx * C)
AT = (TT + BT) / 2.0 #average temp cell: boundary
density = 2600 * ones((nrows,ncols)) # density of granite
dz_crust = dz * ones((nrows,ncols)) #thickness of crust
'''Temperature array'''
T = AT * ones((nrows, ncols)) #intial conditons
'''Initial conditions: create elevation array'''
Z = raster_map(dz*ones((nrows,ncols)), dx=dx) #elevation array
color_scale = color_map(32000.0, 33000.0, 100.0, colmin=vector(1,0,0), colmax=vector(0,1,0))
Z.data[20,20]=dz+999.0 #high elevation near center
'''Boundary conditions'''
Z.data[:,0]=dz #South
Z.data[:,nrows-1]=dz #North
Z.data[0,:]=dz #West
Z.data[ncols-1,:]=dz #East
print Z.data
Z.line_3d(color_map=color_scale,scale=100.0)
'''Draw map of elevation'''
nsteps = 1000
for nt in range(1, nsteps+1):
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])
'''update'''
Z.data[1:nrows-1,1:ncols-1] = Z_new
print nt, Z_new
'''set boundary conditions'''
Q = 0.0 * ones((nrows-2,ncols-2)) #hotspot
Q[20,20] = q * A
'''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, 5.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.)
###color_scale = color_map(50.0, 62.0, 1.0, colmin=vector(1,0,0), colmax=vector(0,1,0))
###'''draw lines for aquifer top with the color_map'''
###Tmap.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'''
###aquifer.draw_layer(color_map=color_scale)
nsteps = 1000
for nt in range(1, nsteps+1):
## for i in range(ncols):
## for j in range(nrows):
## temp_new = T[i,j] + S *(T[i+1,j] ...
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*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))
A=dx*dz
print "S=",S,A,(S*A/(density[21, :]*dz_crust[21, :]))
## temp_new = T[1:nrows-1, 1:ncols-1] + (Q*dt/(C*M))+ \
## (S*A/(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])- \
## (T[1:nrows-1, 1:ncols-1] - T[0:nrows-2, 1:ncols-1])+ \
## (T[1:nrows-1, 2:ncols]-T[1:nrows-1, 1:ncols-1])- \
## (T[1:nrows-1, 1:ncols-1]-T[1:nrows-1, 0:ncols-2]))
'''update'''
T[1:nrows-1,1:ncols-1] = temp_new
print nt, T[21, :]
Tmap.update_line_3d()
plot_line(T[21,:])
#temp_new = T[20,20]+S*(T[20+1,20]-4*T[20,20]+T[20-1,20]+T[20,20+1]+T[20,20-1])+(Q*dt)#/c
'''Thermal expansion of Rocks'''
#temp_new * 9 * .000001 #need to multiply each cell; T[:,:]???
#print temp_new
'''update image'''
Tmap.update_line_3d()
