Spring2008:Smith
From GeoMod
Contents |
[edit]
01/28/2008
from visual import *
from random import uniform, seed
scene.forward = (0, -0.5,1)
scene.up = (0,0,1)
scene.forward = (0, 1,0)
scene.lights = [vector(0, -10, 0)]
xaxis = curve(pos=[ (10, 0, 0), (0, 0,0)], color=(1,0,0))
yaxis = curve(pos=[ (0, 10, 0), (0, 0,0)], color=(0,1,0))
zaxis = curve(pos=[ (0, 0, 10), (0, 0,0)], color=(0,0,1))
seed (2)
nballs = 5
vx = zeros((nballs,))
vz = zeros((nballs,))
print vx
print vz
balls = []
for i in range(nballs):
rad = uniform(.5, 2)
vx[i] = uniform(.005, .01)
balls.append (sphere(radius = rad, pos = (1,0,10)))
g = -9.8
dt = 0.01
t = 0
while 1:
rate #(100)
t += dt
vz = vz + g * dt
for n,i in enumerate(balls):
i.pos.z += vz[n] * dt
if i.pos.z < i.radius:
i.pos.z = i.radius
vz[n] = -vz[n] * 0.75
vx[n] = vx[n] * 0.95
i.pos.x = i.pos.x + vx[n]
print dt, vz[3:5]
[edit]
02/01/2008
from visual import *
from random import uniform,seed
def contact(ball1,ball2):
hit=false
#calculate distance
a=ball1.pos.z-ball2.pos.z
b=ball1.pos.x-ball2.pos.x
d=ball1.pos.y-ball1.pos.y
c=((a*a)+(b*b)+(d*d))**0.5
#sum radii
r=ball1.radius+ball2.radius
#determine if there is a hit
if c <= r:
hit=true
return hit
scene.forward = (0, -0.5,1)
scene.up = (0,0,1)
scene.forward = (0, 1,0)
scene.lights = [vector(0, -10, 0)]
xaxis = curve(pos=[ (20, 0, 0), (0, 0,0)], color=(1,0,0))
yaxis = curve(pos=[ (0, 10, 0), (0, 0,0)], color=(0,1,0))
zaxis = curve(pos=[ (0, 0, 10), (0, 0,0)], color=(0,0,1))
balls=[]
seed(4)
nballs=5
vx=zeros((nballs,))
vz=zeros((nballs,))
for i in range(nballs):
rad=uniform(0.5,2)
vx[i]=uniform(1,2)
balls.append(sphere(radius=rad,pos=(1,0,10)))
'''gravity constant'''
g = -9.8
'''change in time'''
dt = 0.01
'''initial z velocity'''
velz = 0.0
'''initial time'''
t = 0.0
print vx,vz
while 1:
rate(60)
'''time'''
t = t + dt
vz=vz+g*dt
for n,i in enumerate(balls):
i.pos.z = i.pos.z + vz[n] * dt #+ 0.5 * g * dt * dt
i.pos.x = i.pos.x + vx[n] * dt
if i.pos.z < i.radius:
vz[n] = -vz[n] * 0.75
vx[n] = vx[n] * 0.75
i.pos.z = i.radius
for n in range(nballs-1):
for m in range (n+1,nballs):
ball1=balls[n]
ball2=balls[m]
print 'contact',n,m
if contact(ball1,ball2):
print n,m
[edit]
2/1/09 end of class
from visual import *
from random import uniform,seed
def contact(ball1,ball2):
hit=false
#calculate distance
a=ball1.pos.z-ball2.pos.z
b=ball1.pos.x-ball2.pos.x
d=ball1.pos.y-ball1.pos.y
c=((a*a)+(b*b)+(d*d))**0.5
#sum radii
r=ball1.radius+ball2.radius
#determine if there is a hit
if c <= r:
hit=true
return hit
scene.forward = (0, -0.5,1)
scene.up = (0,0,1)
scene.forward = (0, 1,0)
scene.lights = [vector(0, -10, 0)]
xaxis = curve(pos=[ (20, 0, 0), (0, 0,0)], color=(1,0,0))
yaxis = curve(pos=[ (0, 10, 0), (0, 0,0)], color=(0,1,0))
zaxis = curve(pos=[ (0, 0, 10), (0, 0,0)], color=(0,0,1))
balls=[]
seed(4)
nballs=5
vx=zeros((nballs,))
vz=zeros((nballs,))
'''gravity constant'''
g = -9.8
'''change in time'''
dt = 0.01
'''initial z velocity'''
velz = 0.0
'''initial time'''
t = 0.0
b = -1
print vx,vz
while 1:
if t == 0.00 or int(t * 100) or int(t * 200) or int(t * 100) == 100:
b = b + 1
nballs = b
rad=uniform(0.5,2)
print vx, t
vx[b]=uniform(1,2)
balls.append(sphere(radius=rad,pos=(1,0,10)))
rate(60)
'''time'''
t = t + dt
vz=vz+g*dt
vx=vx
for n,i in enumerate(balls):
i.pos.z = i.pos.z + vz[n] * dt #+ 0.5 * g * dt * dt
i.pos.x = i.pos.x + vx[n] * dt
if i.pos.z < i.radius:
vz[n] = -vz[n] * 0.75
vx[n] = vx[n]
i.pos.z = i.radius
for n in range(nballs-1):
for m in range (n+1,nballs):
ball1=balls[n]
ball2=balls[m]
[edit]
02/03/08
from visual import *
from random import uniform,seed
def contact(ball1,ball2):
hit=false
#calculate distance
a=ball1.pos.z-ball2.pos.z
b=ball1.pos.x-ball2.pos.x
d=ball1.pos.y-ball1.pos.y
c=((a*a)+(b*b)+(d*d))**0.5
#sum radii
r=ball1.radius+ball2.radius
#determine if there is a hit
if c <= r:
hit=true
return hit
scene.forward = (0, -0.5,1)
scene.up = (0,0,1)
scene.forward = (0, 1,0)
scene.lights = [vector(0, -10, 0)]
zaxis = curve(pos=[ (0, 0, 15), (0, 0,0)], color=(0,0,1))
xaxis = curve(pos=[ (15, 0, 0), (0, 0, 0)], color=(1,0,0))
zaxis2 = curve(pos=[ (15, 0, 0), (15, 0,15)], color=(1,0,0))
xaxis2 = curve(pos=[ (15, 0, 15), (0, 0,15)], color=(0,1,0))
balls=[]
seed(4)
nballs=5
vx=zeros((nballs,))
vz=zeros((nballs,))
'''gravity constant'''
g = -9.8
'''change in time'''
dt = 0.01
'''initial z velocity'''
velz = 0.0
'''initial time'''
t = 0.0
b = -1
print vx,vz
while 1:
if t == 0.00 or int(t * 100) == 100:
b = b + 1
nballs = b
rad=uniform(0.5,2)
print vx, t
vx[b]=uniform(1,3)
balls.append(sphere(radius=rad,pos=(1,0,10)))
if t == 0.00 or int(t * 100) == 100:
b = b + 1
nballs = b
rad=uniform(0.5,2)
print vx, t
vx[b]=uniform(-1,-3)
balls.append(sphere(radius=rad,pos=(15,0,10)))
rate(60)
'''time'''
t = t + dt
vz=vz+g*dt
vx=vx
a=1
e=1
c=1
r=1
d=1
for n,i in enumerate(balls):
a=i.pos.z-i.pos.z
e=i.pos.x-i.pos.x
d=i.pos.y-i.pos.y
c=((a*a)+(e*e)+(d*d))**0.5
r=i.radius+i.radius
i.pos.z = i.pos.z + vz[n] * dt #+ 0.5 * g * dt * dt
i.pos.x = i.pos.x + vx[n] * dt
if i.pos.z < i.radius:
vz[n] = -vz[n] * 0.75
vx[n] = vx[n]
i.pos.z = i.radius
if c <= r:
vx==-vx
for n in range(nballs-1):
for m in range (n+1,nballs):
ball1=balls[n]
ball2=balls[m]
[edit]
02/17/08
[edit]
02/17/08
Today we went over the Mississippi Embayment with Dr. Cox and broke down how we were going to model it. We also worked on the equations getting closer to having a final answer.
[edit]
02/20/02
Today we went over both implicit and explicit equations for groundwater flow and modeled these equations in excell.
[edit]
02/22/08
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)
'''id groundwatter model'''
'''set up parameters'''
K = 0.001
porosity = 0.3
dt = 3
S = K * dt / porosity
nsteps = 350 #number of timesteps
dhdt_limit = 0.001
'''intial conditions'''
h = ones((11,)) * 5.0
hold = zeros((11,))
l_loop = True
t=0
for t in range(1, nsteps+1):
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
'''solve for heads'''
print h[1:10]
h[1:10] = h[1:10] + (-2*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]
if abs(dhdt) > dhdt_limit:
l_loop = True
[edit]
02/25/08
Today we went over the implicit and explicit for heat transfer. The implicit uses the new timestep and the explicit uses the old timestep.
