Changeset 5827 for anuga_work/development/anuga_1d/steady_flow.py
- Timestamp:
- Oct 9, 2008, 12:54:08 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/anuga_1d/steady_flow.py
r5535 r5827 1 1 import os 2 import Gnuplot2 #import Gnuplot 3 3 from math import sqrt, pi 4 from shallow_water_ 1dimport *4 from shallow_water_domain import * 5 5 from Numeric import allclose, array, zeros, ones, Float, take, sqrt 6 6 from config import g, epsilon … … 21 21 h_1 = 0.5 22 22 u_1 = 0.6 23 g = 9.81 23 24 24 25 25 #h_c = (q**2/g)**(1/3) … … 103 103 104 104 L = 50.0 105 N = 50 #800105 N = 400 106 106 cell_len = L/N 107 107 points = zeros(N+1,Float) … … 134 134 wc, uc, hc = analytical_sol(C) 135 135 136 from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ylim,xlim, rc137 hold(False)138 rc('text', usetex=True)136 #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ylim,xlim, rc 137 #hold(False) 138 #rc('text', usetex=True) 139 139 h_error = zeros(1,Float) 140 140 uh_error = zeros(1,Float) 141 141 k = 0 142 yieldstep = 25.0143 finaltime = 25.0142 yieldstep = 1.0 143 finaltime = 1.0 144 144 domain.limiter = "vanleer" 145 #domain.limiter = "steve_minmod" 146 domain.default_order = 2 147 domain.default_time_order = 2#check order is not working 145 domain.order = 2 146 domain.set_timestepping_method('rk2') 147 domain.cfl = 1.0 148 print 'THE DOMAIN LIMITER is', domain.limiter 149 import time 150 t0=time.time() 148 151 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): 149 152 domain.write_time() 150 print 'Test 1' 151 plot1 = subplot(211) 152 print 'Test 2' 153 plot(X,w,X,StageQ,X,ElevationQ) 154 print 'Test 3' 155 #xlabel('Position') 156 #ylabel('Stage (m)') 157 #legend(('Analytical Solution', 'Numerical Solution', 'Channel Bed'), 158 # 'upper right', shadow=True) 159 plot1.set_ylim([-0.1,1.0]) 160 print 'Test 4' 161 plot2 = subplot(212) 162 print 'Test 5' 163 plot(X,u*h,X,XmomQ)#/(StageQ-ElevationQ)) 164 print 'Test 6' 165 #xlabel('x (m)') 166 #ylabel(r'X-momentum ($m^2/s$)') 167 h_error[k] = 1.0/(N)*sum(abs(wc-StageC)) 168 uh_error[k] = 1.0/(N)*sum(abs(uc*hc-XmomC)) 169 print "h_error %.10f" %(h_error[k]) 170 print "uh_error %.10f"% (uh_error[k]) 171 print "h_max %.10f"%max(wc-StageC) 172 print "uh max %.10f"%max(uc*hc-XmomC) 173 filename = "steady_flow_" 174 filename += domain.limiter 175 filename += str(400) 176 filename += ".ps" 177 #savefig(filename) 178 show() 179 153 print "integral", domain.quantities['stage'].get_integral() 154 if t>0.0: 155 print 'Test 1' 156 from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion 157 #ion() 158 hold(False) 159 #rc('text', usetex=True) 160 clf() 161 162 plot1 = subplot(211) 163 print 'Test 2' 164 plot(X,w,X,StageQ,X,ElevationQ) 165 print 'Test 3' 166 xlabel('Position') 167 ylabel('Stage (m)') 168 legend(('Analytical Solution', 'Numerical Solution', 'Channel Bed'), 169 'upper right', shadow=True) 170 plot1.set_ylim([-0.1,1.0]) 171 print 'Test 4' 172 plot2 = subplot(212) 173 print 'Test 5' 174 plot(X,u*h,X,XmomQ)#/(StageQ-ElevationQ)) 175 print 'Test 6' 176 xlabel('x (m)') 177 ylabel(r'X-momentum (m^2/s)') 178 plot2.set_ylim([0.297,0.301]) 179 h_error[k] = 1.0/(N)*sum(abs(wc-StageC)) 180 uh_error[k] = 1.0/(N)*sum(abs(uc*hc-XmomC)) 181 print "h_error %.10f" %(h_error[k]) 182 print "uh_error %.10f"% (uh_error[k]) 183 print "h_max %.10f"%max(wc-StageC) 184 print "uh max %.10f"%max(uc*hc-XmomC) 185 #filename = "steady_flow_" 186 #filename += domain.limiter 187 #filename += str(400) 188 #filename += ".ps" 189 #savefig(filename) 190 show() 191 print 'That took %.2f seconds'%(time.time()-t0) 192 raw_input("Press return key!")
Note: See TracChangeset
for help on using the changeset viewer.