Changeset 6167
 Timestamp:
 Jan 14, 2009, 4:10:50 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/anuga_1d/dam_padarn.py
r6038 r6167 1 1 import os 2 2 from math import sqrt, pi 3 from shallow_water_vel_domainimport *3 from channel_domain_Ab import * 4 4 from Numeric import allclose, array, zeros, ones, Float, take, sqrt 5 5 from config import g, epsilon 6 6 7 7 8 h1 = 10.08 h1 = 5.0 9 9 h0 = 0.0 10 10 … … 65 65 # g.plot(plot1,plot2) 66 66 67 67 h2=5.0 68 k=1 68 69 69 70 print "TEST 1DSOLUTION III  DRY BED" … … 73 74 for i in range(len(x)): 74 75 if x[i]<=L/4.0: 75 y[i] = h0 76 y[i] = h0*width([x[i]]) 76 77 elif x[i]<=3*L/4.0: 77 y[i] = h 178 y[i] = h2*width([x[i]]) 78 79 else: 79 y[i] = h0 80 y[i] = h0*width([x[i]]) 80 81 return y 81 82 83 def width(x): 84 return k 82 85 86 83 87 import time 84 88 … … 86 90 yieldstep = finaltime 87 91 L = 2000.0 # Length of channel (m) 88 number_of_cells = [810]#,200,500,1000,2000,5000,10000,20000] 89 h_error = zeros(len(number_of_cells),Float) 90 uh_error = zeros(len(number_of_cells),Float) 92 number_of_cells = [200] 91 93 k = 0 92 for i in range(len(number_of_cells)): 93 N = int(number_of_cells[i]) 94 print "Evaluating domain with %d cells" %N 95 cell_len = L/N # Origin = 0.0 96 points = zeros(N+1,Float) 97 for j in range(N+1): 98 points[j] = j*cell_len 94 widths = [1] 95 heights= [] 96 for i in range(len(widths)): 97 k=widths[i] 98 for i in range(len(number_of_cells)): 99 N = int(number_of_cells[i]) 100 print "Evaluating domain with %d cells" %N 101 cell_len = L/N # Origin = 0.0 102 points = zeros(N+1,Float) 103 for j in range(N+1): 104 points[j] = j*cell_len 99 105 100 domain = Domain(points)106 domain = Domain(points) 101 107 102 domain.set_quantity('stage', stage) 103 domain.set_boundary({'exterior': Reflective_boundary(domain)}) 104 domain.order = 2 105 domain.set_timestepping_method('rk2') 106 domain.set_CFL(1.0) 107 domain.set_limiter("vanleer") 108 #domain.h0=0.0001 108 domain.set_quantity('area', stage) 109 domain.set_quantity('width',width) 110 print "width in cell 1",domain.quantities['width'].vertex_values[1] 111 domain.set_boundary({'exterior': Reflective_boundary(domain)}) 112 domain.order = 2 113 domain.set_timestepping_method('rk2') 114 domain.set_CFL(1.0) 115 domain.set_limiter("vanleer") 116 #domain.h0=0.0001 117 118 t0 = time.time() 109 119 110 t0 = time.time() 120 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): 121 domain.write_time() 111 122 112 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): 113 domain.write_time() 114 115 N = float(N) 116 StageC = domain.quantities['stage'].centroid_values 117 XmomC = domain.quantities['xmomentum'].centroid_values 118 C = domain.centroids 119 h, uh, u = analytical_sol(C,domain.time) 120 h_error[k] = 1.0/(N)*sum(abs(hStageC)) 121 uh_error[k] = 1.0/(N)*sum(abs(uhXmomC)) 122 print "h_error %.10f" %(h_error[k]) 123 print "uh_error %.10f"% (uh_error[k]) 124 k = k+1 125 print 'That took %.2f seconds' %(time.time()t0) 126 X = domain.vertices 127 StageQ = domain.quantities['stage'].vertex_values 128 XmomQ = domain.quantities['xmomentum'].vertex_values 129 velQ = domain.quantities['velocity'].vertex_values 123 N = float(N) 124 HeightC = domain.quantities['height'].centroid_values 125 DischargeC = domain.quantities['discharge'].centroid_values 126 C = domain.centroids 127 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)) 130 #print "h_error %.10f" %(h_error[k]) 131 #print "uh_error %.10f"% (uh_error[k]) 132 k = k+1 133 print 'That took %.2f seconds' %(time.time()t0) 134 X = domain.vertices 135 heights.append(domain.quantities['height'].vertex_values) 136 VelocityQ = domain.quantities['velocity'].vertex_values 137 #stage = domain.quantities['stage'].vertex_values 138 h, uh, u = analytical_sol(X.flat,domain.time) 139 x = X.flat 130 140 131 h, uh, u = analytical_sol(X.flat,domain.time) 132 x = X.flat 141 #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot 142 from pylab import * 143 import pylab as p 144 import matplotlib.axes3d as p3 145 print 'test 2' 146 hold(False) 147 print 'test 3' 148 plot1 = subplot(211) 149 print 'test 4' 150 plot(x,h,x,heights[0].flat) 151 print 'test 5' 152 plot1.set_ylim([1,11]) 153 xlabel('Position') 154 ylabel('Stage') 155 legend(('Analytical Solution', 'Numerical Solution'), 156 'upper right', shadow=True) 157 plot2 = subplot(212) 158 plot(x,u,x,VelocityQ.flat) 159 plot2.set_ylim([35,35]) 133 160 134 from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot 135 print 'test 2' 136 hold(False) 137 print 'test 3' 138 plot1 = subplot(211) 139 print 'test 4' 140 plot(x,h,x,StageQ.flat) 141 print 'test 5' 142 plot1.set_ylim([1,11]) 143 xlabel('Position') 144 ylabel('Stage') 145 legend(('Analytical Solution', 'Numerical Solution'), 146 'upper right', shadow=True) 147 plot2 = subplot(212) 148 plot(x,u,x,velQ.flat) 149 plot2.set_ylim([35,35]) 150 151 xlabel('Position') 152 ylabel('Velocity') 161 xlabel('Position') 162 ylabel('Velocity') 153 163 154 155 156 157 158 164 file = "dry_bed_" 165 file += str(number_of_cells[i]) 166 file += ".eps" 167 #savefig(file) 168 show() 159 169 160 170 161 print "Error in height", h_error162 print "Error in xmom", uh_error171 #print "Error in height", h_error 172 #print "Error in xmom", uh_error
Note: See TracChangeset
for help on using the changeset viewer.