1d groundwater model

From GeoMod

1-dimensional groundwater model with 11 nodes where:

• initial condition: at t = 0; h(for the entire domain) = 5.0
• boundary conditions: at t > 0: h(node=1) = 30; h(node=11)
• convergence criteria: Model is at steady-state when dh/dt < dhdt_limit (dhdt_limit is set to 0.001 in the code below)

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
porosity = 0.3
dt = 3

S = K * dt / porosity

nsteps = 10 #number of timesteps
dhdt_limit = 0.001

'''initial conditions'''
h = ones((11,)) * 5.0
hold = ones((11,)) * 0.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] = 30.0
    h[10] = 5.0
    #print h

    '''solve for heads'''
    #print h[1:10]
    h[1:10] = h[1:10] + (-2.0*h[1:10]*h[1:10] + h[0:9]*h[0:9] + h[2:11]*h[2:11]) * S
    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