Changeset 8241 for trunk/anuga_work/development/2010-projects/anuga_1d/channel/var_width_depth.py

Timestamp:

Oct 27, 2011, 8:23:34 PM (14 years ago)

Author:

steve

Message:

Got the ctac figures working

File:

• trunk/anuga_work/development/2010-projects/anuga_1d/channel/var_width_depth.py (modified) (8 diffs)

trunk/anuga_work/development/2010-projects/anuga_1d/channel/var_width_depth.py

@@ -1 +1 @@
 import os
 import random
-from math import sqrt, pow, pi 
-from channel_domain_Ab import *
-from Numeric import allclose, array, zeros, ones, Float, take, sqrt
-from config import g, epsilon
+from math import sqrt, pow, pi
+from channel_domain import *
+from numpy import allclose, array, zeros, ones, take, sqrt
+from anuga_1d.config import g, epsilon
 
+
+random.seed(111)
 
 print "Variable Width Only Test"
@@ -11 +13 @@
 # Define functions for initial quantities
 def initial_area(x):
-    y = zeros(len(x),Float)
+    y = zeros(len(x),'f')
     for i in range(len(x)):
         y[i]=(10-randomarray[i])
@@ -20 +22 @@
 def width(x):
 
-    y = zeros(len(x),Float)
+    y = zeros(len(x),'f')
     for i in range(len(x)):
         y[i]=randomarray[i]
@@ -28 +30 @@
 def bed(x):
 
-    y = zeros(len(x),Float)
+    y = zeros(len(x),'f')
     for i in range(len(x)):
         y[i]=randomarray2[i]
@@ -38 +40 @@
 def stage(x):
 
-    y = zeros(len(x),Float)
+    y = zeros(len(x),'f')
     for i in range(len(x)):
         if x[i]<200 and x[i]>150:
@@ -64 +66 @@
 print "Evaluating domain with %d cells" %N
 cell_len = L/N # Origin = 0.0
-points = zeros(N+1,Float)
+points = zeros(N+1,'f')
 
 # Define the centroid points
@@ -74 +76 @@
 
 # Define random array for width
-randomarray=zeros(len(points),Float)
+randomarray=zeros(len(points),'f')
 for j in range(N+1):
     randomarray[j]=random.normalvariate(3,1)
-randomarray2=zeros(len(points),Float)
+randomarray2=zeros(len(points),'f')
 for j in range(N+1):
     randomarray2[j]=random.normalvariate(3,1)
@@ -98 +100 @@
 t0 = time.time()
 
-print domain.quantities['elevation'].vertex_values
+
+AreaC = domain.quantities['area'].centroid_values
+BedC = domain.quantities['elevation'].centroid_values
+WidthC = domain.quantities['width'].centroid_values
+#
+AreaC[:] = (10.0 - BedC)* WidthC
+
+
+#print domain.quantities['elevation'].vertex_values
+
+finaltime = 100.0
+yieldstep = 1.0
+
+domain.initialize_plotting(stage_lim = [-1.0, 20.0],
+                           width_lim = [0.0, 6.0],
+                           velocity_lim = [-1.0, 1.0])
+
 
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     domain.write_time()
 
-print domain.quantities['elevation'].vertex_values
+    domain.update_plotting()
 
-N = float(N)
-HeightC = domain.quantities['height'].centroid_values
-DischargeC = domain.quantities['discharge'].centroid_values
-C = domain.centroids
-print 'That took %.2f seconds' %(time.time()-t0)
-X = domain.vertices
-HeightQ = domain.quantities['height'].vertex_values
-VelocityQ = domain.quantities['velocity'].vertex_values
-x = X.flat
-z = domain.quantities['elevation'].vertex_values.flat
-b = domain.quantities['width'].vertex_values.flat
-stage=HeightQ.flat+z
+#print domain.quantities['elevation'].vertex_values
 
+domain.hold_plotting()
 
-
-from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
-
-#hold(False)
-   
-plot1=subplot(211)
-
-plot(x,stage,x,z,x,b)
- 
-plot1.set_ylim([-5,15])
-xlabel('Position')
-ylabel('Stage')
-legend(('Solution', 'Bed','Width'),
-            'upper right', shadow=True)
-
-plot2=subplot(212)
-plot(x,VelocityQ.flat)
-plot2.set_ylim([-10,10])
-xlabel('Position')
-ylabel('Velocity')
-   
-show()
-