1d thermal expansion model: Smith

From GeoMod

from visual 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)
    

'''1d groundwater model'''

'''set up parameters'''
K = 0.001
dx = 10
dt = 10
S = K * dt / (dx*dx)

nsteps = 10 #number of timesteps
dhdt_limit = 0.001

'''initial conditions'''
h = ones((11,)) * 15.0
hold = ones((11,)) * 15.0
#print hold[0:11]

l_loop = True
t = 0
#for t in range(1,nsteps+1):
while (l_loop):
    t += 1
    print "timestep =", t
    '''record old heads'''
    for i, n in enumerate(h):
        hold[i] = h[i]
   
    '''set boundary conditions'''
    h[0] = 15.0
    h[10] = 335.0
    #print h

    '''solve for heads'''
    #print h[1:10]
    h[1:10] = h[1:10] + (h[1:10] -2 * h[0:9] + h[2:11]) * S 

    '''Thermal expansion of Rocks'''
    h[1:10] * 9 * .000001
     
    print h
    plot_line(h)

    '''post processing'''
    '''determine if model is at steady state'''
    l_loop = False
    for i, n in enumerate(h):
        dhdt = hold[i] - h[i]
        #print hold[i], h[i], dhdt
        if abs(dhdt) > dhdt_limit:
            l_loop = True