Changeset 6453 for anuga_work/development/anuga_1d/dam_padarn.py
 Timestamp:
 Mar 4, 2009, 3:30:29 PM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/anuga_1d/dam_padarn.py
r6167 r6453 10 10 11 11 def analytical_sol(C,t): 12 12 if t==0: 13 t=0.0001 13 14 #t = 0.0 # time (s) 14 15 # gravity (m/s^2) … … 90 91 yieldstep = finaltime 91 92 L = 2000.0 # Length of channel (m) 92 number_of_cells = [ 200]93 number_of_cells = [810] 93 94 k = 0 94 widths = [1 ]95 widths = [1,2,5] 95 96 heights= [] 97 velocities = [] 98 96 99 for i in range(len(widths)): 97 100 k=widths[i] … … 118 121 t0 = time.time() 119 122 120 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): 123 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): 121 124 domain.write_time() 122 125 … … 126 129 C = domain.centroids 127 130 h, uh, u = analytical_sol(C,domain.time) 128 #h_error[k] = 1.0/(N)*sum(abs(hStageC))129 #uh_error[k] = 1.0/(N)*sum(abs(uhXmomC))131 h_error = 1.0/(N)*sum(abs(hHeightC)) 132 u_error = 1.0/(N)*sum(abs(uhDischargeC)) 130 133 #print "h_error %.10f" %(h_error[k]) 131 134 #print "uh_error %.10f"% (uh_error[k]) … … 134 137 X = domain.vertices 135 138 heights.append(domain.quantities['height'].vertex_values) 136 VelocityQ = domain.quantities['velocity'].vertex_values139 velocities.append( domain.quantities['velocity'].vertex_values) 137 140 #stage = domain.quantities['stage'].vertex_values 138 141 h, uh, u = analytical_sol(X.flat,domain.time) 139 142 x = X.flat 140 143 144 print "Error in height", h_error 145 print "Error in xmom", u_error 141 146 #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot 142 147 from pylab import * … … 144 149 import matplotlib.axes3d as p3 145 150 print 'test 2' 146 hold(False)151 #hold(False) 147 152 print 'test 3' 148 153 plot1 = subplot(211) 149 154 print 'test 4' 150 plot(x,h ,x,heights[0].flat)155 plot(x,heights[0].flat,x,h) 151 156 print 'test 5' 152 157 plot1.set_ylim([1,11]) 153 158 xlabel('Position') 154 159 ylabel('Stage') 155 legend(('Analytical Solution', 'Numerical Solution'), 156 'upper right', shadow=True) 160 legend(('Numerical Solution','Analytic Soltuion'), 'upper right', shadow=True) 157 161 plot2 = subplot(212) 158 plot(x,u,x,VelocityQ.flat) 162 163 164 plot(x,velocities[0].flat,x,u) 159 165 plot2.set_ylim([35,35]) 160 161 166 xlabel('Position') 162 167 ylabel('Velocity') 168 169 print heights[0].flatheights[1].flat 163 170 164 171 file = "dry_bed_" … … 168 175 show() 169 176 170 171 #print "Error in height", h_error 172 #print "Error in xmom", uh_error 177
Note: See TracChangeset
for help on using the changeset viewer.