Changeset 6453 for anuga_work/development/anuga_1d/dam_padarn.py

Timestamp:

Mar 4, 2009, 3:30:29 PM (14 years ago)

Author:

steve

Message:

Added Padarn's (2008/2009 summer student) files

File:

• anuga_work/development/anuga_1d/dam_padarn.py (modified) (7 diffs)

anuga_work/development/anuga_1d/dam_padarn.py

@@ -10 +10 @@
 
 def analytical_sol(C,t):
-    
+    if t==0:
+        t=0.0001
     #t  = 0.0     # time (s)
     # gravity (m/s^2)
@@ -90 +91 @@
 yieldstep = finaltime
 L = 2000.0     # Length of channel (m)
-number_of_cells = [200]
+number_of_cells = [810]
 k = 0
-widths = [1]
+widths = [1,2,5]
 heights= []
+velocities = []
+
 for i in range(len(widths)):
     k=widths[i]
@@ -118 +121 @@
         t0 = time.time()
 
-        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):            
             domain.write_time()
 
@@ -126 +129 @@
         C = domain.centroids
         h, uh, u = analytical_sol(C,domain.time)
-        #h_error[k] = 1.0/(N)*sum(abs(h-StageC))
-        #uh_error[k] = 1.0/(N)*sum(abs(uh-XmomC))
+        h_error = 1.0/(N)*sum(abs(h-HeightC))
+        u_error = 1.0/(N)*sum(abs(uh-DischargeC))
         #print "h_error %.10f" %(h_error[k])
         #print "uh_error %.10f"% (uh_error[k])
@@ -134 +137 @@
         X = domain.vertices
         heights.append(domain.quantities['height'].vertex_values)
-        VelocityQ = domain.quantities['velocity'].vertex_values
+        velocities.append( domain.quantities['velocity'].vertex_values)
         #stage = domain.quantities['stage'].vertex_values
         h, uh, u = analytical_sol(X.flat,domain.time)
         x = X.flat
-    
+            
+        print "Error in height", h_error
+        print "Error in xmom", u_error
         #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
 from pylab import *
@@ -144 +149 @@
 import matplotlib.axes3d as p3
 print 'test 2'
-hold(False)
+#hold(False)
 print 'test 3'
 plot1 = subplot(211)
 print 'test 4'
-plot(x,h,x,heights[0].flat)
+plot(x,heights[0].flat,x,h)
 print 'test 5'
 plot1.set_ylim([-1,11])
 xlabel('Position')
 ylabel('Stage')
-legend(('Analytical Solution', 'Numerical Solution'),
-           'upper right', shadow=True)
+legend(('Numerical Solution','Analytic Soltuion'), 'upper right', shadow=True)
 plot2 = subplot(212)
-plot(x,u,x,VelocityQ.flat)
+
+
+plot(x,velocities[0].flat,x,u)
 plot2.set_ylim([-35,35])
-    
 xlabel('Position')
 ylabel('Velocity')
+
+print heights[0].flat-heights[1].flat
    
 file = "dry_bed_"
@@ -168 +175 @@
 show()
     
-    
-#print "Error in height", h_error
-#print "Error in xmom", uh_error
+