Changeset 5827 for anuga_work/development/anuga_1d/parabolic_cannal.py
 Timestamp:
 Oct 9, 2008, 12:54:08 PM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/anuga_1d/parabolic_cannal.py
r5535 r5827 1 1 import os 2 2 from math import sqrt, pi, sin, cos 3 from shallow_water_ 1dimport *3 from shallow_water_domain import * 4 4 from Numeric import allclose, array, zeros, ones, Float, take, sqrt 5 5 from config import g, epsilon 6 6 """ 7 7 def newLinePlot(title='Simple Plot'): 8 8 import Gnuplot … … 20 20 plot3 = Gnuplot.PlotItems.Data(x3.flat,y3.flat,with="linespoints 1") 21 21 g.plot(plot1,plot2,plot3) 22 22 """ 23 23 24 24 def analytic_cannal(C,t): … … 54 54 55 55 L_x = 2500.0 # Length of channel (m) 56 N = 8# Number of compuational cells56 N = 100 # Number of compuational cells 57 57 cell_len = 4*L_x/N # Origin = 0.0 58 58 … … 63 63 domain = Domain(points) 64 64 65 domain.order = 2 #make this unnecessary 66 domain.default_order = 2 67 domain.default_time_order = 2 65 domain.order = 2 66 domain.set_timestepping_method('rk2') 68 67 domain.cfl = 1.0 69 68 domain.beta = 1.0 … … 72 71 #domain.limiter = "vanalbada" 73 72 #domain.limiter = "superbee" 74 domain.limiter = "steve_minmod" 73 #domain.limiter = "steve_minmod" 74 domain.limiter = "pyvolution" 75 75 76 76 def stage(x): … … 113 113 domain.set_quantity('elevation',elevation) 114 114 #domain.set_quantity('height',height) 115 116 115 domain.set_boundary({'exterior': Reflective_boundary(domain)}) 117 116 118 117 import time 119 118 t0 = time.time() 120 finaltime = 1122.0*0.75121 yieldstep = finaltime 122 yieldstep = 10.0123 plot1 = newLinePlot("Stage")124 plot2 = newLinePlot("Xmomentum")119 #finaltime = 1122.0*0.75 120 yieldstep = finaltime = 10.0 121 #yieldstep = 10.0 122 #plot1 = newLinePlot("Stage") 123 #plot2 = newLinePlot("Xmomentum") 125 124 C = domain.centroids 126 125 X = domain.vertices 127 126 StageQ = domain.quantities['stage'].vertex_values 128 127 XmomQ = domain.quantities['xmomentum'].vertex_values 128 #Bed = domain.quantities['elevation'].vertex_values 129 129 import time 130 130 t0 = time.time() … … 133 133 domain.write_time() 134 134 u,h,z,z_b = analytic_cannal(X.flat,t) 135 linePlot(plot1,X,z,X,z_b,X,StageQ)136 linePlot(plot2,X,u*h,X,u*h,X,XmomQ)135 #linePlot(plot1,X,z,X,z_b,X,StageQ) 136 #linePlot(plot2,X,u*h,X,u*h,X,XmomQ) 137 137 HeightQ = domain.quantities['stage'].centroid_valuesdomain.quantities['elevation'].centroid_values 138 138 u,hc,z,z_b = analytic_cannal(C,t) … … 142 142 error = 1.0/(N)*sum(abs(hcHeightQ)) 143 143 print 'Error measured at centroids', error 144 #raw_input() 145 146 #from Gnuplot import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot 147 from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot 148 #import Gnuplot 149 X = domain.vertices 150 u,h,z,z_b = analytic_cannal(X.flat,domain.time) 151 StageQ = domain.quantities['stage'].vertex_values 152 XmomQ = domain.quantities['xmomentum'].vertex_values 153 hold(False) 154 plot1 = subplot(211) 155 plot(X,h,X,StageQ) 156 #plot1.set_ylim([4,11]) 157 #title('Free Surface Elevation of a Dry DamBreak') 158 xlabel('Position') 159 ylabel('Stage') 160 legend(('Analytical Solution', 'Numerical Solution'), 161 'upper right', shadow=True) 162 plot2 = subplot(212) 163 plot(X,u*h,X,XmomQ) 164 #plot2.set_ylim([1,25]) 165 #title('Xmomentum Profile of a Dry DamBreak') 166 xlabel('Position') 167 ylabel('Xmomentum') 168 show() 144 # #raw_input() 145 # #from Gnuplot import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot 146 from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot 147 # #import Gnuplot 148 X = domain.vertices 149 u,h,z,z_b = analytic_cannal(X.flat,domain.time) 150 StageQ = domain.quantities['stage'].vertex_values 151 XmomQ = domain.quantities['xmomentum'].vertex_values 152 hold(False) 153 plot1 = subplot(211) 154 plot(X,h,X,StageQ,X,z_b) 155 #plot1.set_ylim([4,11]) 156 #title('?????????') 157 xlabel('Position') 158 ylabel('Stage') 159 legend(('Analytical Solution', 'Numerical Solution', 'Bed'), 160 'upper right', shadow=True) 161 plot2 = subplot(212) 162 plot(X,u*h,X,XmomQ) 163 #plot2.set_ylim([1,25]) 164 #title('?????????????????????????') 165 xlabel('Position') 166 ylabel('Xmomentum') 167 show()
Note: See TracChangeset
for help on using the changeset viewer.