Changeset 5066 for anuga_work/development/MISG_2008
- Timestamp:
- Feb 20, 2008, 4:01:00 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/MISG_2008/sensitivity_study.py
r4937 r5066 3 3 Water driven up a linear slope and time varying boundary, 4 4 similar to a beach environment 5 6 Test case for MISG 2008, University of Wollongong 5 7 """ 6 8 … … 16 18 from anuga.shallow_water import Time_boundary 17 19 from anuga.shallow_water import Transmissive_boundary 20 from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary 21 22 from anuga.pmesh.mesh_interface import create_mesh_from_regions 18 23 19 24 from Numeric import zeros, Float … … 25 30 # Setup computational domain 26 31 #------------------------------------------------------------------------------ 27 28 points, vertices, boundary = rectangular_cross(90, 30, 29 len1=3.0, 30 len2=1.0) # Basic mesh 31 32 domain = Domain(points, vertices, boundary) # Create domain 32 y_extent=2000.0 33 bounding_polygon = [[0,0],[2000,0],[2000,2000],[0,2000]] 34 beach_polygon = [[5,900],[225,900],[225,1100],[5,1100]] 35 interior_regions = [[beach_polygon, 50]] 36 create_mesh_from_regions(bounding_polygon, 37 boundary_tags={'bottom':[0], 38 'right':[1], 39 'top':[2], 40 'left':[3]}, 41 maximum_triangle_area=2000, 42 interior_regions= interior_regions, 43 filename='misg.msh', 44 use_cache=False, 45 verbose=True) 46 domain = Domain('misg.msh',use_cache=False, verbose=True) 47 48 33 49 domain.set_name('sensitivity') # Output to file runup.sww 34 50 domain.set_datadir('.') # Use current directory for output 35 51 domain.tight_slope_limiters = 1 36 52 domain.beta_h = 0 53 domain.set_minimum_storable_height(0.01) 54 55 print 'number of elements: ', len(domain) 37 56 38 57 #------------------------------------------------------------------------------ … … 44 63 N = len(x) 45 64 z = zeros(N, Float) 46 47 65 48 66 for i in range(N): 49 # linear bed slope 50 z[i] = -x[i]/15 51 52 # IID noise 53 z[i] += normal(0.0, 0.01) 67 if x[i] <= (depth-c)/m: 68 z[i] = m*x[i]+c 69 else: 70 # linear bed slope 71 #z[i] = slope*x[i] + depth-(slope)*(depth-c)/m 72 z[i] = slope*x[i] + depth-(slope)*(depth-c)/m + 1.0*sin((x[i]-(depth-c)/m)/T) + 1.0*sin((y[i]-(depth-c)/m)/T) 73 74 # IID noise 75 #z[i] += normal(-100.0, 5.0) 54 76 55 77 return z 56 57 domain.set_quantity('elevation', topography) # Use function for elevation58 domain.set_quantity('friction', 0.01) # Constant friction59 domain.set_quantity('stage', -.4) # Constant negative initial stage60 61 78 62 79 #------------------------------------------------------------------------------ … … 69 86 Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values 70 87 Bw = Time_boundary(domain=domain, # Time dependent boundary 71 f=lambda t: [(.1*sin(t*2*pi)-0.1) , 0.0, 0.0]) 72 73 # Associate boundary tags with boundary objects 74 domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br}) 75 76 77 import time 78 t0 = time.time() 79 #------------------------------------------------------------------------------ 80 # Evolve system through time 81 #------------------------------------------------------------------------------ 82 83 for t in domain.evolve(yieldstep = 0.1, finaltime = 50.0): 84 domain.write_time() 85 86 print 'That took %.2f seconds' %(time.time()-t0) 87 88 f=lambda t: [0.0, 0.0, 0.0]) 89 90 from Numeric import zeros, Float, arange 91 start=-1.0/40.0 92 finish=-1.8/40.0 93 n_vec=arange(start,finish,-2.0/400.0) 94 N=len(n_vec) 95 #N=1 # use for 2.5% and 5% change tests 96 max0=zeros(N, Float) 97 max1=zeros(N, Float) 98 max2=zeros(N, Float) 99 max3=zeros(N, Float) 100 mmax0=zeros(N, Float) 101 mmax1=zeros(N, Float) 102 mmax2=zeros(N, Float) 103 mmax3=zeros(N, Float) 104 output_file ='misg_influence_data.csv' 105 domain.store=True 106 fid_out = open(output_file, 'w') 107 s = 'variable, h0, h1, h2, h3, v0, v1, v2, v3\n' 108 fid_out.write(s) 109 110 for i in range(N): 111 manning = 0.01*1.0 112 A=1.0*1.0 113 w=1.0/100.0 114 m=-1.0/20.0 115 depth=-10.0 116 c=10.0 117 T=100.0*1.0 118 #slope=-1.0/40.0*1.0 119 slope = n_vec[i] 120 121 var = slope 122 domain.set_quantity('elevation', topography) # Use function for elevation 123 domain.set_quantity('friction', manning) 124 domain.set_quantity('stage', 0.0) # Constant negative initial stage 125 domain.set_quantity('xmomentum', 0.0) # Use function for elevation 126 domain.set_quantity('ymomentum', 0.0) # Use function for elevation 127 swwname = 'sensitivity_2d_wavey_bottom'#+str(i) 128 domain.set_name(swwname) 129 Btd = Transmissive_Momentum_Set_Stage_boundary(domain=domain, 130 function=lambda t: [(0<t<pi/w)*A*sin(w*t), 0.0, 0.0]) 131 domain.set_boundary({'left': Br, 'right': Btd, 'top': Br, 'bottom': Br}) 132 133 # locations to record wave height and velocity 134 tri_id0 = domain.get_triangle_containing_point([-c/m,y_extent*0.5]) 135 tri_id1 = domain.get_triangle_containing_point([-c/m-10,y_extent*0.5]) 136 tri_id2 = domain.get_triangle_containing_point([-c/m-20,y_extent*0.5]) 137 tri_id3 = domain.get_triangle_containing_point([-c/m+50,y_extent*0.5]) 138 139 stage_centroid = domain.quantities['stage'].centroid_values 140 xmom_centroid = domain.quantities['xmomentum'].centroid_values 141 elevation_centroid = domain.quantities['elevation'].centroid_values 142 143 elev0 = elevation_centroid[tri_id0] 144 elev1 = elevation_centroid[tri_id1] 145 elev2 = elevation_centroid[tri_id2] 146 elev3 = elevation_centroid[tri_id3] 147 148 import time 149 t0 = time.time() 150 #------------------------------------------------------------------------------ 151 # Evolve system through time 152 #------------------------------------------------------------------------------ 153 domain.time=0.0 154 finaltime=4.0 #pi/w+150.0 155 max=0 156 for t in domain.evolve(yieldstep = 1.0, finaltime = finaltime): 157 158 if stage_centroid[tri_id0]-elev0>max0[i]: 159 max0[i]=stage_centroid[tri_id0]-elev0 160 if stage_centroid[tri_id1]-elev1>max1[i]: 161 max1[i]=stage_centroid[tri_id1]-elev1 162 if stage_centroid[tri_id2]-elev2>max2[i]: 163 max2[i]=stage_centroid[tri_id2]-elev2 164 if stage_centroid[tri_id3]-elev3>max3[i]: 165 max3[i]=stage_centroid[tri_id3]-elev3 166 167 if abs(xmom_centroid[tri_id0])>mmax0[i]: 168 mmax0[i]=abs(xmom_centroid[tri_id0]) 169 if abs(xmom_centroid[tri_id1])>mmax1[i]: 170 mmax1[i]=abs(xmom_centroid[tri_id1]) 171 if abs(xmom_centroid[tri_id2])>mmax2[i]: 172 mmax2[i]=abs(xmom_centroid[tri_id2]) 173 if abs(xmom_centroid[tri_id3])>mmax3[i]: 174 mmax3[i]=abs(xmom_centroid[tri_id3]) 175 176 print 'finished cycle: ', i 177 print 'That took %.2f seconds' %(time.time()-t0) 178 179 180 s = '%.6f, %.6f, %.6f, %.6f, %.6f, %.6f, %.6f, %.6f, %.6f\n' %(var, max0[i], max1[i], max2[i], max3[i], mmax0[i], mmax1[i], mmax2[i], mmax3[i]) 181 fid_out.write(s) 182 183 from pylab import plot, ion, savefig, xlabel, ylabel, title, close, legend,hist,subplot 184 ion() 185 plot(n_vec,max2,'bo-',n_vec,max1,'go-',n_vec,max0,'ro-',n_vec,max3,'co-') 186 title('MISG testing: varying slope') 187 xlabel('slope') 188 ylabel('depth (m)') 189 #legend(['loc2','loc1','loc0','loc3']) 190 name='bed_slope' 191 savefig(name) 192 ''' 193 #title('MISG testing: varying slope PDF') 194 subplot(2,2,1) 195 [h0,bins,patches]=hist(max0,10) 196 subplot(2,2,2) 197 [h1,bins,pa3ches]=hist(max1,10) 198 subplot(2,2,3) 199 [h2,bins,patches]=hist(max2,10) 200 subplot(2,2,4) 201 [h3,bins,patches]=hist(max3,10) 202 #plot(bins,h0,'bo-',bins,h1,'ro-',bins,h2,'go-',bins,h3,'co-') 203 xlabel('Runup') 204 ylabel('Frequency') 205 name='bed_slope_PDF' 206 savefig(name) 207 ''' 208 209 close('all')
Note: See TracChangeset
for help on using the changeset viewer.