1-D Finite difference code for Ground Water

From GeoMod

from Numeric import *
from visual import *
from random import uniform

n = 10
dx = 5.0
Acx = 1.0
S = 0.01
dt = 100
K = 0.0001
n_steps = 500

h_init = 20.0

h = h_init * ones((n,),Float)
h_old = h

print h


#boundary conditions
h[0] = 30.0
h[9] = 20.0

cells = []
for i in range(n):
    cells.append(box(pos=(i*dx,h[i]/2.0,0), length=dx, height = h[i], color=(uniform(0,1),uniform(0,1), 1)))

for t in range(n_steps):
    rate(1)
    for i in range(1, n-1):
        #print i
        h[i] = ( (S*h_old[i]/dt) - (K*Acx*(h_old[i]-h_old[i-1])/dx) - (K*Acx*(h_old[i]-h_old[i+1])/dx) ) / (S/dt)
        cells[i].pos.y = h[i]/2.0
        cells[i].height = h[i]
    h_old = h
    print h
#print h