Changeset 6167

Timestamp:

Jan 14, 2009, 4:10:50 PM (16 years ago)

Author:

wilson

Message:

update 14/01/09

File:

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

anuga_work/development/anuga_1d/dam_padarn.py

@@ -1 +1 @@
 import os
 from math import sqrt, pi
-from shallow_water_vel_domain import *
+from channel_domain_Ab import *
 from Numeric import allclose, array, zeros, ones, Float, take, sqrt
 from config import g, epsilon
 
 
-h1 = 10.0
+h1 = 5.0
 h0 = 0.0
 
@@ -65 +65 @@
 #    g.plot(plot1,plot2)
 
-
+h2=5.0
+k=1
 
 print "TEST 1D-SOLUTION III -- DRY BED"
@@ -73 +74 @@
     for i in range(len(x)):
         if x[i]<=L/4.0:
-            y[i] = h0
+            y[i] = h0*width([x[i]])
         elif x[i]<=3*L/4.0:
-            y[i] = h1
+            y[i] = h2*width([x[i]])
         else:
-            y[i] = h0
+            y[i] = h0*width([x[i]])
     return y
 
+def width(x):
+    return k
 
+ 
 import time
 
@@ -86 +90 @@
 yieldstep = finaltime
 L = 2000.0     # Length of channel (m)
-number_of_cells = [810]#,200,500,1000,2000,5000,10000,20000]
-h_error = zeros(len(number_of_cells),Float)
-uh_error = zeros(len(number_of_cells),Float)
+number_of_cells = [200]
 k = 0
-for i in range(len(number_of_cells)):
-    N = int(number_of_cells[i])
-    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
+widths = [1]
+heights= []
+for i in range(len(widths)):
+    k=widths[i]
+    for i in range(len(number_of_cells)):
+        N = int(number_of_cells[i])
+        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(points)
     
-    domain.set_quantity('stage', stage)
-    domain.set_boundary({'exterior': Reflective_boundary(domain)})
-    domain.order = 2
-    domain.set_timestepping_method('rk2')
-    domain.set_CFL(1.0)
-    domain.set_limiter("vanleer")
-    #domain.h0=0.0001
+        domain.set_quantity('area', stage)
+        domain.set_quantity('width',width)
+        print "width in cell 1",domain.quantities['width'].vertex_values[1]
+        domain.set_boundary({'exterior': Reflective_boundary(domain)})
+        domain.order = 2
+        domain.set_timestepping_method('rk2')
+        domain.set_CFL(1.0)
+        domain.set_limiter("vanleer")
+        #domain.h0=0.0001
+ 
+        t0 = time.time()
 
-    t0 = time.time()
+        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+            domain.write_time()
 
-    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
-        domain.write_time()
-
-    N = float(N)
-    StageC = domain.quantities['stage'].centroid_values
-    XmomC = domain.quantities['xmomentum'].centroid_values
-    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))
-    print "h_error %.10f" %(h_error[k])
-    print "uh_error %.10f"% (uh_error[k])
-    k = k+1
-    print 'That took %.2f seconds' %(time.time()-t0)
-    X = domain.vertices
-    StageQ = domain.quantities['stage'].vertex_values
-    XmomQ = domain.quantities['xmomentum'].vertex_values
-    velQ = domain.quantities['velocity'].vertex_values
+        N = float(N)
+        HeightC = domain.quantities['height'].centroid_values 
+        DischargeC = domain.quantities['discharge'].centroid_values
+        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))
+        #print "h_error %.10f" %(h_error[k])
+        #print "uh_error %.10f"% (uh_error[k])
+        k = k+1
+        print 'That took %.2f seconds' %(time.time()-t0)
+        X = domain.vertices
+        heights.append(domain.quantities['height'].vertex_values)
+        VelocityQ = domain.quantities['velocity'].vertex_values
+        #stage = domain.quantities['stage'].vertex_values
+        h, uh, u = analytical_sol(X.flat,domain.time)
+        x = X.flat
     
-    h, uh, u = analytical_sol(X.flat,domain.time)
-    x = X.flat
+        #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+from pylab import *
+import pylab as p
+import matplotlib.axes3d as p3
+print 'test 2'
+hold(False)
+print 'test 3'
+plot1 = subplot(211)
+print 'test 4'
+plot(x,h,x,heights[0].flat)
+print 'test 5'
+plot1.set_ylim([-1,11])
+xlabel('Position')
+ylabel('Stage')
+legend(('Analytical Solution', 'Numerical Solution'),
+           'upper right', shadow=True)
+plot2 = subplot(212)
+plot(x,u,x,VelocityQ.flat)
+plot2.set_ylim([-35,35])
     
-    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
-    print 'test 2'
-    hold(False)
-    print 'test 3'
-    plot1 = subplot(211)
-    print 'test 4'
-    plot(x,h,x,StageQ.flat)
-    print 'test 5'
-    plot1.set_ylim([-1,11])
-    xlabel('Position')
-    ylabel('Stage')
-    legend(('Analytical Solution', 'Numerical Solution'),
-           'upper right', shadow=True)
-    plot2 = subplot(212)
-    plot(x,u,x,velQ.flat)
-    plot2.set_ylim([-35,35])
-    
-    xlabel('Position')
-    ylabel('Velocity')
+xlabel('Position')
+ylabel('Velocity')
    
-    file = "dry_bed_"
-    file += str(number_of_cells[i])
-    file += ".eps"
-    #savefig(file)
-    show()
+file = "dry_bed_"
+file += str(number_of_cells[i])
+file += ".eps"
+#savefig(file)
+show()
     
     
-print "Error in height", h_error
-print "Error in xmom", uh_error
+#print "Error in height", h_error
+#print "Error in xmom", uh_error