3rd movement in model: Bronson

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

'''molding a ball'''
ball = (sphere_pos.x(20,0), color, yellow, radius, 4.0)

'''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
g = -9.8
vel = 0.0
xvel = .005


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


'''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.)

nsteps = 1000

for nt in range(1, nsteps+1):

       vel = g * dt
       sphere_pos.x = vel * dt

       sphere_pos.z < .125=(sphere_pos.z, .125)
       xvel = xvel * .8

       sphere_pos.x > 3.2=(vel, -vel * .83)

       sphere_pos.x < 3.1=(vel, -vel * .853)

    
###    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, :])
    
                   
'''update'''
T[1:nrows-1,1:ncols-1] = temp_new
print nt, T[21, :]
sphere_pos.x = xvel
Tmap.update_line_3d()

plot_line(T[21,:])
    
'''update image'''
Tmap.update_line_3d()