Changeset 5827 for anuga_work/development/anuga_1d/parabolic_cannal.py

Timestamp:

Oct 9, 2008, 12:54:08 PM (15 years ago)

Author:

sudi

Message:

Adding new versions of shallow_water_domain

File:

• anuga_work/development/anuga_1d/parabolic_cannal.py (modified) (8 diffs)

anuga_work/development/anuga_1d/parabolic_cannal.py

@@ -1 +1 @@
 import os
 from math import sqrt, pi, sin, cos
-from shallow_water_1d import *
+from shallow_water_domain import *
 from Numeric import allclose, array, zeros, ones, Float, take, sqrt
 from config import g, epsilon
-
+"""
 def newLinePlot(title='Simple Plot'):
     import Gnuplot
@@ -20 +20 @@
     plot3 = Gnuplot.PlotItems.Data(x3.flat,y3.flat,with="linespoints 1")
     g.plot(plot1,plot2,plot3)
-
+"""
 
 def analytic_cannal(C,t):
@@ -54 +54 @@
 
 L_x = 2500.0     # Length of channel (m)
-N = 8         # Number of compuational cells
+N = 100         # Number of compuational cells
 cell_len = 4*L_x/N   # Origin = 0.0
 
@@ -63 +63 @@
 domain = Domain(points)
 
-domain.order = 2 #make this unnecessary
-domain.default_order = 2
-domain.default_time_order = 2
+domain.order = 2
+domain.set_timestepping_method('rk2')
 domain.cfl = 1.0
 domain.beta = 1.0
@@ -72 +71 @@
 #domain.limiter = "vanalbada"
 #domain.limiter = "superbee"
-domain.limiter = "steve_minmod"
+#domain.limiter = "steve_minmod"
+domain.limiter = "pyvolution"
 
 def stage(x):
@@ -113 +113 @@
 domain.set_quantity('elevation',elevation)
 #domain.set_quantity('height',height)
-
 domain.set_boundary({'exterior': Reflective_boundary(domain)})
  
 import time
 t0 = time.time()
-finaltime = 1122.0*0.75
-yieldstep = finaltime
-yieldstep = 10.0
-plot1 = newLinePlot("Stage")
-plot2 = newLinePlot("Xmomentum")
+#finaltime = 1122.0*0.75
+yieldstep = finaltime = 10.0
+#yieldstep = 10.0
+#plot1 = newLinePlot("Stage")
+#plot2 = newLinePlot("Xmomentum")
 C = domain.centroids
 X = domain.vertices
 StageQ = domain.quantities['stage'].vertex_values
 XmomQ = domain.quantities['xmomentum'].vertex_values
+#Bed = domain.quantities['elevation'].vertex_values
 import time
 t0 = time.time()
@@ -133 +133 @@
     domain.write_time()
     u,h,z,z_b = analytic_cannal(X.flat,t)
-    linePlot(plot1,X,z,X,z_b,X,StageQ)
-    linePlot(plot2,X,u*h,X,u*h,X,XmomQ)
+    #linePlot(plot1,X,z,X,z_b,X,StageQ)
+    #linePlot(plot2,X,u*h,X,u*h,X,XmomQ)
     HeightQ = domain.quantities['stage'].centroid_values-domain.quantities['elevation'].centroid_values
     u,hc,z,z_b = analytic_cannal(C,t)
@@ -142 +142 @@
     error = 1.0/(N)*sum(abs(hc-HeightQ))
     print 'Error measured at centroids', error
-    #raw_input()         
-
-#from Gnuplot import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
-from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
-#import Gnuplot  
-X = domain.vertices
-u,h,z,z_b = analytic_cannal(X.flat,domain.time)
-StageQ = domain.quantities['stage'].vertex_values
-XmomQ = domain.quantities['xmomentum'].vertex_values
-hold(False)
-plot1 = subplot(211)
-plot(X,h,X,StageQ)
-#plot1.set_ylim([4,11])
-#title('Free Surface Elevation of a Dry Dam-Break')
-xlabel('Position')
-ylabel('Stage')
-legend(('Analytical Solution', 'Numerical Solution'),
-       'upper right', shadow=True)
-plot2 = subplot(212)
-plot(X,u*h,X,XmomQ)
-#plot2.set_ylim([-1,25])
-#title('Xmomentum Profile of a Dry Dam-Break')
-xlabel('Position')
-ylabel('Xmomentum')
-show()
+    # #raw_input()         
+    # #from Gnuplot import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+    # #import Gnuplot  
+    X = domain.vertices
+    u,h,z,z_b = analytic_cannal(X.flat,domain.time)
+    StageQ = domain.quantities['stage'].vertex_values
+    XmomQ = domain.quantities['xmomentum'].vertex_values
+    hold(False)
+    plot1 = subplot(211)
+    plot(X,h,X,StageQ,X,z_b)
+    #plot1.set_ylim([4,11])
+    #title('?????????')
+    xlabel('Position')
+    ylabel('Stage')
+    legend(('Analytical Solution', 'Numerical Solution', 'Bed'),
+           'upper right', shadow=True)
+    plot2 = subplot(212)
+    plot(X,u*h,X,XmomQ)
+    #plot2.set_ylim([-1,25])
+    #title('?????????????????????????')
+    xlabel('Position')
+    ylabel('Xmomentum')
+    show()