Changeset 6453

Timestamp:

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

Author:

steve

Message:

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

Location:

anuga_work/development/anuga_1d

Files:

• backup.py (added)
• channel_domain.py (modified) (3 diffs)
• channel_domain2.py (added)
• channel_domain_Ab.py (added)
• channel_flow_1_padarn.py (added)
• channel_flow_2_padarn.py (added)
• compile.py (added)
• dam_padarn.py (modified) (7 diffs)
• domain.py (modified) (1 diff)
• dry_dam_vel.py (modified) (1 diff)
• old_compile_all.py (added)
• parabolic_cannal_padarn.py (added)
• quantity.py (modified) (2 diffs)
• radial_dam_break.py (added)
• shallow_water_domain_steve.py (added)
• stead_state_varying_width.py (added)
• straight_padarn.py (added)
• test_padarn.py (added)
• test_shallow_water_domain_suggestion1.py (added)
• var_depth_only.py (added)
• var_depth_only2.py (added)
• var_width_depth.py (added)
• var_width_only.py (added)
• var_width_only_random.py (added)

anuga_work/development/anuga_1d/channel_domain.py

@@ -6 +6 @@
 consisting of methods specific to the Shallow Water Wave Equation
 
+This particular modification of the Domain class implements the ability to
+vary the width of the 1D channel that the water flows in. As a result the
+conserved variables are different than previous implementations and so are the
+equations.
 
 U_t + E_x = S
 
 where
-
-U = [w, uh]
-E = [uh, u^2h + gh^2/2]
+------------!!!! NOTE THIS NEEDS UPDATING !!!!------------------
+U = [A, Q]
+E = [Q, Q^2/A + gh^2/2]
 S represents source terms forcing the system
 (e.g. gravity, friction, wind stress, ...)
@@ -32 +36 @@
 
 The conserved quantities are w, uh
-
+--------------------------------------------------------------------------
 For details see e.g.
 Christopher Zoppou and Stephen Roberts,
@@ -39 +43 @@
 
 
-John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
-Geoscience Australia, 2006
+John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou,
+Padarn Wilson, Geoscience Australia, 2008
 """
 

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
+

anuga_work/development/anuga_1d/domain.py

@@ -916 +916 @@
         # Update ghosts
         self.update_ghosts()
-
+        
         # Initial update of vertex and edge values
-        self.distribute_to_vertices_and_edges()
-
+        self.distribute_to_vertices_and_edges()  
+        
         # Update extrema if necessary (for reporting)
         self.update_extrema()

anuga_work/development/anuga_1d/dry_dam_vel.py

@@ -78 +78 @@
 for j in range(N+1):
     points[j] = j*cell_len
-
-boundary = {}
-boundary[0,0] = 'left'
-boundary[N-1,1] = 'right'
-
-domain = Domain(points, boundary = boundary)
-
+        
+domain = Domain(points)
+    
 domain.set_quantity('stage', stage)
-
-Br =  Reflective_boundary(domain)
-domain.set_boundary({'left':  Br, 'right': Br})
+domain.set_boundary({'exterior': Reflective_boundary(domain)})
 domain.order = 2
 domain.set_timestepping_method('euler')

anuga_work/development/anuga_1d/quantity.py

@@ -135 +135 @@
             #Use function specific method
             self.set_function_values(X, location)
+          
         elif type(X) in [types.FloatType, types.IntType, types.LongType]:
             if location == 'centroids':
@@ -163 +164 @@
         
         if location == 'centroids':
-            
+         
             P = self.domain.centroids
             self.set_values(f(P), location)
         else:
             #Vertices
+            
             P = self.domain.get_vertices()
             
             for i in range(2):
-                self.vertex_values[:,i] = f(P[:,i])        
-
+               
+                self.vertex_values[:,i] = f(P[:,i])
+               
     def set_array_values(self, values, location='vertices'):
         """Set values for quantity