Changeset 7860 for anuga_work/development/2010-projects/anuga_1d/sww/run_dry_dam.py

Timestamp:

Jun 18, 2010, 5:49:57 PM (14 years ago)

Author:

steve

Message:

Continuing to numpy the for loops

File:

• anuga_work/development/2010-projects/anuga_1d/sww/run_dry_dam.py (modified) (4 diffs)

anuga_work/development/2010-projects/anuga_1d/sww/run_dry_dam.py

@@ -1 +1 @@
 import os
 from math import sqrt, pi
-from shallow_water_vel_domain import *
-from Numeric import allclose, array, zeros, ones, Float, take, sqrt
-from config import g, epsilon
+import numpy
+import time
+#from Numeric import allclose, array, zeros, ones, Float, take, sqrt
 
+
+
+from anuga_1d.sww.sww_domain import *
+from anuga_1d.config import g, epsilon
+from anuga_1d.base.generic_mesh import uniform_mesh
 
 h1 = 10.0
@@ -49 +54 @@
     return h , u*h, u
 
-print "TEST 1D-SOLUTION III -- DRY BED"
+
 
 def stage(x):
-    y = zeros(len(x),Float)
-    for i in range(len(x)):
-        if x[i]<=L/4.0:
-            y[i] = h0
-        elif x[i]<=3*L/4.0:
-            y[i] = h1
-        else:
-            y[i] = h0
+    import numpy
+
+    y = numpy.where( (x>=L/4.0) & (x<=3*L/4.0), h1 , h0)
+
+#    for i in range(len(x)):
+#        if x[i]<=L/4.0:
+#            y[i] = h0
+#        elif x[i]<=3*L/4.0:
+#            y[i] = h1
+#        else:
+#            y[i] = h0
     return y
 
 
-import time
 
-finaltime = 2.0
+print "TEST 1D-SOLUTION III -- DRY BED"
+
+
+finaltime = 4.0
 yieldstep = 0.1
 L = 2000.0     # Length of channel (m)
@@ -73 +83 @@
 N = 800
 print "Evaluating domain with %d cells" %N
-cell_len = L/N # Origin = 0.0
-points = zeros(N+1,Float)
-
-for j in range(N+1):
-    points[j] = j*cell_len
-        
-domain = Domain(points)
+domain = Domain(*uniform_mesh(N))
     
 domain.set_quantity('stage', stage)
-domain.set_boundary({'exterior': Reflective_boundary(domain)})
+
+Br = Reflective_boundary(domain)
+
+domain.set_boundary({'left': Br, 'right' : Br})
 domain.order = 2
 domain.set_timestepping_method('euler')
 domain.set_CFL(1.0)
-domain.set_limiter("vanleer")
+domain.set_limiter("minmod")
 #domain.h0=0.0001
 
@@ -94 +101 @@
     domain.write_time()
 
-print 'end'
+print 'That took %.2f seconds' %(time.time()-t0)
 
 
+N = float(N)
+HeightC    = domain.quantities['height'].centroid_values
+StageC     = domain.quantities['stage'].centroid_values
+BedC       = domain.quantities['elevation'].centroid_values
+C          = domain.centroids
+
+HeightV    = domain.quantities['height'].vertex_values
+StageV     = domain.quantities['stage'].vertex_values
+BedV       = domain.quantities['elevation'].vertex_values
+VelocityV  = domain.quantities['velocity'].vertex_values
+X          = domain.vertices
+
+
+import pylab
+
+pylab.hold(False)
+
+plot1 = pylab.subplot(211)
+
+pylab.plot(X.flat,BedV.flat,X.flat,StageV.flat)
+
+plot1.set_ylim([-1,11])
+
+pylab.xlabel('Position')
+pylab.ylabel('Stage')
+pylab.legend(('Analytical Solution', 'Numerical Solution'),
+           'upper right', shadow=True)
+
+
+plot2 = pylab.subplot(212)
+
+pylab.plot(X.flat,VelocityV.flat)
+plot2.set_ylim([-20,10])
+
+pylab.xlabel('Position')
+pylab.ylabel('Velocity')
+
+pylab.show()
+