 Timestamp:
 Jun 18, 2010, 5:49:57 PM (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/2010projects/anuga_1d/sww/run_dry_dam.py
r7857 r7860 1 1 import os 2 2 from math import sqrt, pi 3 from shallow_water_vel_domain import * 4 from Numeric import allclose, array, zeros, ones, Float, take, sqrt 5 from config import g, epsilon 3 import numpy 4 import time 5 #from Numeric import allclose, array, zeros, ones, Float, take, sqrt 6 6 7 8 9 from anuga_1d.sww.sww_domain import * 10 from anuga_1d.config import g, epsilon 11 from anuga_1d.base.generic_mesh import uniform_mesh 7 12 8 13 h1 = 10.0 … … 49 54 return h , u*h, u 50 55 51 print "TEST 1DSOLUTION III  DRY BED" 56 52 57 53 58 def stage(x): 54 y = zeros(len(x),Float) 55 for i in range(len(x)): 56 if x[i]<=L/4.0: 57 y[i] = h0 58 elif x[i]<=3*L/4.0: 59 y[i] = h1 60 else: 61 y[i] = h0 59 import numpy 60 61 y = numpy.where( (x>=L/4.0) & (x<=3*L/4.0), h1 , h0) 62 63 # for i in range(len(x)): 64 # if x[i]<=L/4.0: 65 # y[i] = h0 66 # elif x[i]<=3*L/4.0: 67 # y[i] = h1 68 # else: 69 # y[i] = h0 62 70 return y 63 71 64 72 65 import time66 73 67 finaltime = 2.0 74 print "TEST 1DSOLUTION III  DRY BED" 75 76 77 finaltime = 4.0 68 78 yieldstep = 0.1 69 79 L = 2000.0 # Length of channel (m) … … 73 83 N = 800 74 84 print "Evaluating domain with %d cells" %N 75 cell_len = L/N # Origin = 0.0 76 points = zeros(N+1,Float) 77 78 for j in range(N+1): 79 points[j] = j*cell_len 80 81 domain = Domain(points) 85 domain = Domain(*uniform_mesh(N)) 82 86 83 87 domain.set_quantity('stage', stage) 84 domain.set_boundary({'exterior': Reflective_boundary(domain)}) 88 89 Br = Reflective_boundary(domain) 90 91 domain.set_boundary({'left': Br, 'right' : Br}) 85 92 domain.order = 2 86 93 domain.set_timestepping_method('euler') 87 94 domain.set_CFL(1.0) 88 domain.set_limiter(" vanleer")95 domain.set_limiter("minmod") 89 96 #domain.h0=0.0001 90 97 … … 94 101 domain.write_time() 95 102 96 print ' end'103 print 'That took %.2f seconds' %(time.time()t0) 97 104 98 105 106 N = float(N) 107 HeightC = domain.quantities['height'].centroid_values 108 StageC = domain.quantities['stage'].centroid_values 109 BedC = domain.quantities['elevation'].centroid_values 110 C = domain.centroids 111 112 HeightV = domain.quantities['height'].vertex_values 113 StageV = domain.quantities['stage'].vertex_values 114 BedV = domain.quantities['elevation'].vertex_values 115 VelocityV = domain.quantities['velocity'].vertex_values 116 X = domain.vertices 117 118 119 import pylab 120 121 pylab.hold(False) 122 123 plot1 = pylab.subplot(211) 124 125 pylab.plot(X.flat,BedV.flat,X.flat,StageV.flat) 126 127 plot1.set_ylim([1,11]) 128 129 pylab.xlabel('Position') 130 pylab.ylabel('Stage') 131 pylab.legend(('Analytical Solution', 'Numerical Solution'), 132 'upper right', shadow=True) 133 134 135 plot2 = pylab.subplot(212) 136 137 pylab.plot(X.flat,VelocityV.flat) 138 plot2.set_ylim([20,10]) 139 140 pylab.xlabel('Position') 141 pylab.ylabel('Velocity') 142 143 pylab.show() 144
Note: See TracChangeset
for help on using the changeset viewer.