Changeset 8240 for trunk/anuga_work/development/2010-projects/anuga_1d/channel

Timestamp:

Oct 26, 2011, 6:32:44 PM (14 years ago)

Author:

steve

Message:

Padarn's code

Location:

trunk/anuga_work/development/2010-projects/anuga_1d/channel

Files:

• channel_domain.py (modified) (2 diffs)
• var_depth_only.py (modified) (3 diffs)

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

@@ -144 +144 @@
 
     def update_derived_quantites(self):
+
         pass
+
+
+
+    def initialize_plotting(self,
+                            stage_lim = [-1.0, 40.0],
+                            height_lim = [-1.0, 10.0],
+                            velocity_lim = [-5.0, 5.0]):
+
+        import pylab
+        pylab.ion()
+
+
+        x = self.get_vertices().flatten()
+        z = self.quantities['elevation'].vertex_values.flatten()
+        w = self.quantities['stage'].vertex_values.flatten()
+        h = self.quantities['height'].vertex_values.flatten()
+        v = self.quantities['velocity'].vertex_values.flatten()
+
+        print x.shape
+        print z.shape
+
+        self.plot1 = pylab.subplot(311)
+
+        self.zplot, = pylab.plot(x, z)
+        self.wplot, = pylab.plot(x, w)
+
+        self.plot1.set_ylim(stage_lim)
+        pylab.xlabel('Position')
+        pylab.ylabel('Stage')
+
+        self.plot2 = pylab.subplot(312)
+
+        self.hplot, = pylab.plot(x, h)
+
+        self.plot2.set_ylim(height_lim)
+        pylab.xlabel('Position')
+        pylab.ylabel('Height')
+
+        self.plot3 = pylab.subplot(313)
+
+        self.vplot, = pylab.plot(x, v)
+
+        self.plot3.set_ylim(velocity_lim)
+
+        pylab.xlabel('Position')
+        pylab.ylabel('Velocity')
+
+
+    def update_plotting(self):
+
+        import pylab
+
+        #x = self.get_vertices().flatten()
+        z = self.quantities['elevation'].vertex_values.flatten()
+        w = self.quantities['stage'].vertex_values.flatten()
+        h = self.quantities['height'].vertex_values.flatten()
+        v = self.quantities['velocity'].vertex_values.flatten()
+
+
+        self.zplot.set_ydata(z)
+        self.wplot.set_ydata(w)
+        self.hplot.set_ydata(h)
+        self.vplot.set_ydata(v)
+
+        pylab.draw()
+
+
+    def hold_plotting(self):
+
+        self.update_plotting()
+        import pylab
+        
+        pylab.ioff()
+
+        pylab.show()
+
+
+
+    def finalize_plotting(self):
+
+        pass
+
+
 
 #=============== End of Channel Domain ===============================
@@ -198 +282 @@
     # (i.e. bed and width) have been accurately computed in the previous
     # timestep. 
-    h_C[:] = numpy.where(a_C > 0.0, a_C/(b_C + h0/b_C), 0.0)
-    u_C[:] = numpy.where(a_C > 0.0, d_C/(a_C + h0/a_C), 0.0)
+    h_C[:] = numpy.where(a_C > 0.0, a_C/b_C, 0.0)
+    u_C[:] = numpy.where(a_C > 0.0, d_C/a_C, 0.0)
 
     w_C[:] = h_C + z_C

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

@@ -54 +54 @@
 import time
 
-# Set final time and yield time for simulation
-finaltime = 1
-yieldstep = 0.1
+
 
 # Length of channel (m)
@@ -103 +101 @@
 
 
-#AreaC = domain.quantities['area'].centroid_values
-#BedC = domain.quantities['elevation'].centroid_values
-#WidthC = domain.quantities['width'].centroid_values
-#
-#AreaC[:] = (6.0 - BedC)* WidthC
+AreaC = domain.quantities['area'].centroid_values
+BedC = domain.quantities['elevation'].centroid_values
+WidthC = domain.quantities['width'].centroid_values
+#
+AreaC[:] = (6.0 - BedC)* WidthC
 
 #domain.set_quantity('area', initial_area)
@@ -126 +124 @@
 print domain.quantities['width'].vertex_values
 
-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).flatten()
-HeightQ = (domain.quantities['height'].vertex_values).flatten()
-VelocityQ = (domain.quantities['velocity'].vertex_values).flatten()
-Z = (domain.quantities['elevation'].vertex_values).flatten()
-stage=HeightQ + Z
-
-
-from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,savefig
-import pylab as p
-#hold(False)
-
-plot1=subplot(211)
-
-plot(X,stage,X,Z)
-
-plot1.set_ylim([-1,11])
-xlabel('Position')
-ylabel('Stage')
-legend(('Solution', 'Bed Height'),
-                'upper right', shadow=True)
-
-plot2=subplot(212)
-
-
-plot(X,VelocityQ)
-plot2.set_ylim([-5,5])
-xlabel('Position')
-ylabel('Velocity')
-
-savefig(str(i)+'.png')
-i+=1
-plot1.clear()
-plot2.clear()
-
+#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).flatten()
+#HeightQ = (domain.quantities['height'].vertex_values).flatten()
+#VelocityQ = (domain.quantities['velocity'].vertex_values).flatten()
+#Z = (domain.quantities['elevation'].vertex_values).flatten()
+#stage=HeightQ + Z
+#
+#
+#from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,savefig
+#import pylab as p
+##hold(False)
+#
+#plot1=subplot(211)
+#
+#plot(X,stage,X,Z)
+#
+#plot1.set_ylim([-1,11])
+#xlabel('Position')
+#ylabel('Stage')
+#legend(('Solution', 'Bed Height'),
+#                'upper right', shadow=True)
+#
+#plot2=subplot(212)
+#
+#
+#plot(X,VelocityQ)
+#plot2.set_ylim([-5,5])
+#xlabel('Position')
+#ylabel('Velocity')
+#
+#savefig(str(i)+'.png')
+#i+=1
+#plot1.clear()
+#plot2.clear()
+
+domain.distribute_to_vertices_and_edges()
+
+
+
+# Set final time and yield time for simulation
+finaltime = 100.0
+yieldstep = 1.0
+
+
+
+domain.initialize_plotting(stage_lim = [-2.0, 14.0],
+                           height_lim = [-2.0, 14.0],
+                           velocity_lim = [-10.0, 10.0])
 
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     domain.write_time()
 
-    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).flatten()
-    HeightQ = (domain.quantities['height'].vertex_values).flatten()
-    VelocityQ = (domain.quantities['velocity'].vertex_values).flatten()
-    Z = (domain.quantities['elevation'].vertex_values).flatten()
-    stage=HeightQ + Z
-
-
-    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,savefig
-    import pylab as p
-#hold(False)
+    domain.update_plotting()
+
+#    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).flatten()
+#    HeightQ = (domain.quantities['height'].vertex_values).flatten()
+#    VelocityQ = (domain.quantities['velocity'].vertex_values).flatten()
+#    Z = (domain.quantities['elevation'].vertex_values).flatten()
+#    stage=HeightQ + Z
+#
+#
+#    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,savefig
+#    import pylab as p
+##hold(False)
+#
+#    plot1=subplot(211)
+#
+#    plot(X,stage,X,Z)
+#
+#    plot1.set_ylim([-1,11])
+#    xlabel('Position')
+#    ylabel('Stage')
+#    legend(('Solution', 'Bed Height'),
+#                'upper right', shadow=True)
+#
+#    plot2=subplot(212)
+#
+#
+#    plot(X,VelocityQ)
+#    plot2.set_ylim([-5,5])
+#    xlabel('Position')
+#    ylabel('Velocity')
+#
+#    savefig(str(i)+'.png')
+#    i+=1
+#    plot1.clear()
+#    plot2.clear()
+
+print domain.quantities['elevation'].vertex_values
+
+
+domain.hold_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).flatten()
+#HeightQ = (domain.quantities['height'].vertex_values).flatten()
+#VelocityQ = (domain.quantities['velocity'].vertex_values).flatten()
+#Z = (domain.quantities['elevation'].vertex_values).flatten()
+#stage=HeightQ + Z
+#
+#
+#from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+#
+##hold(False)
+#
+#plot1=subplot(211)
+#
+#plot(X,stage,X,Z)
+#
+#plot1.set_ylim([-1,11])
+#xlabel('Position')
+#ylabel('Stage')
+#legend(('Solution', 'Bed Height'),
+#           'upper right', shadow=True)
+#
+#plot2=subplot(212)
+#plot(X,VelocityQ)
+#plot2.set_ylim([-5,5])
+#xlabel('Position')
+#ylabel('Velocity')
    
-    plot1=subplot(211)
-
-    plot(X,stage,X,Z)
- 
-    plot1.set_ylim([-1,11])
-    xlabel('Position')
-    ylabel('Stage')
-    legend(('Solution', 'Bed Height'),
-                'upper right', shadow=True)
-
-    plot2=subplot(212)
-
-    
-    plot(X,VelocityQ)
-    plot2.set_ylim([-5,5])
-    xlabel('Position') 
-    ylabel('Velocity')
-
-    savefig(str(i)+'.png')
-    i+=1
-    plot1.clear()
-    plot2.clear()
-
-print domain.quantities['elevation'].vertex_values
-
-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).flatten()
-HeightQ = (domain.quantities['height'].vertex_values).flatten()
-VelocityQ = (domain.quantities['velocity'].vertex_values).flatten()
-Z = (domain.quantities['elevation'].vertex_values).flatten()
-stage=HeightQ + Z
-
-
-from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
-
-#hold(False)
-   
-plot1=subplot(211)
-
-plot(X,stage,X,Z)
- 
-plot1.set_ylim([-1,11])
-xlabel('Position')
-ylabel('Stage')
-legend(('Solution', 'Bed Height'),
-           'upper right', shadow=True)
-
-plot2=subplot(212)
-plot(X,VelocityQ)
-plot2.set_ylim([-5,5])
-xlabel('Position')
-ylabel('Velocity')
-   
-show()
-
+#show()
+