1d thermal expansion model example

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 = 3.0 #W/mK
dx = 3200
dt = 1000
S = K * dt / (dx*dx)

nsteps = 10 #number of timesteps
dhdt_limit = 0.00001

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

'''Boundary conditions'''
Qin = 1.0e-12

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
    h[10] = h[10] + 6.952e9 * Qin + (K/dx)*(h[9]-h[10])
    #print h

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

    '''Thermal expansion of Rocks'''
    #h[1:10] * 9 * .000001
     
    #print h
    if t % 100 == 0:
        print "time =", t*dt / (60*60*24*365.25)
        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