Changeset 3693
 Timestamp:
 Oct 4, 2006, 2:36:27 PM (18 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/examples/runup_convergence.py
r3689 r3693 4 4 similar to a beach environment. 5 5 6 The study area is discretised as a regular triangular grid 100m x 100m7 6 """ 8 7 … … 23 22 from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary 24 23 from abstract_2d_finite_volumes.util import file_function 25 from pylab import plot, xlabel, ylabel, title, ion, close, savefig, figure, axis, legend, grid, hold 24 from pylab import plot, xlabel, ylabel, title, ion, close, savefig,\ 25 figure, axis, legend, grid, hold 26 26 27 27 … … 30 30 # Model constants 31 31 32 slope = 0.02 # 1:50 Slope, reaches h=20m 1000m from western bndry, and h=0 (coast) at 300m 32 slope = 0.02 # 1:50 Slope, reaches h=20m 1000m from western bndry, 33 # and h=0 (coast) at 300m 33 34 highest_point = 6 # Highest elevation (m) 34 35 sea_level = 0 # Mean sea level 35 36 min_elevation = 20 # Lowest elevation (elevation of offshore flat part) 36 37 offshore_depth = sea_levelmin_elevation # offshore water depth 37 amplitude = 0.5 # Solitary wave height H 38 39 amplitude = 0.5 # Solitary wave height H 38 40 normalized_amplitude = amplitude/offshore_depth 39 41 simulation_name = 'runup_convergence' … … 57 59 length = eastwest 58 60 width = northsouth 59 points, vertices, boundary = rectangular_cross(length/dx, width/dy, len1=length, len2=width, 61 points, vertices, boundary = rectangular_cross(length/dx, width/dy, 62 len1=length, len2=width, 60 63 origin = (west, south)) 61 64 … … 63 66 64 67 65 66 68 # Unstructured mesh 67 69 polygon = [[east,north],[west,north],[west,south],[east,south]] 68 interior_polygon = [[400,north10],[west+10,north10],[west+10,south+10],[400,south+10]] 70 interior_polygon = [[400,north10],[west+10,north10], 71 [west+10,south+10],[400,south+10]] 69 72 meshname = simulation_name + '.msh' 70 73 create_mesh_from_regions(polygon, 71 boundary_tags={'top': [0], 'left': [1], 'bottom': [2], 'right': [3]}, 72 maximum_triangle_area=dx*dy/4, # Triangle area commensurate with structured mesh 74 boundary_tags={'top': [0], 'left': [1], 75 'bottom': [2], 'right': [3]}, 76 maximum_triangle_area=dx*dy/4, 73 77 filename=meshname, 74 78 interior_regions=[[interior_polygon,dx*dy/32]]) 79 75 80 domain = Domain(meshname, use_cache=True, verbose = True) 76 81 domain.set_minimum_storable_height(0.01) 77 78 82 79 83 domain.set_name(simulation_name) … … 85 89 86 90 #def topography(x,y): 87 # return slope*x+highest_point # Return linear bed slope bathymetry as vector91 # return slope*x+highest_point # Return linear bed slope (vector) 88 92 89 93 def topography(x,y): … … 121 125 122 126 def waveform(t): 123 return sea_level + amplitude/cosh(((t50)/offshore_depth)*(0.75*g*amplitude)**0.5)**2 124 125 Bw = Time_boundary(domain=domain, # Time dependent boundary 126 f=lambda t: [waveform(t), 0.0, 0.0]) 127 return sea_level +\ 128 amplitude/cosh(((t50)/offshore_depth)*(0.75*g*amplitude)**0.5)**2 127 129 128 130 # Time dependent boundary for stage, where momentum is set automatically 129 131 Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform) 130 131 132 132 133 # Associate boundary tags with boundary objects … … 134 135 135 136 137 # Find initial runup location and height (coastline) 136 138 w0 = domain.get_maximum_inundation_elevation() 137 139 x0, y0 = domain.get_maximum_inundation_location() 138 140 print 'Coastline elevation = %.2f at (%.2f, %.2f)' %(w0, x0, y0) 141 142 # Sanity check 139 143 w_i = domain.get_quantity('stage').get_values(interpolation_points=[[x0,y0]]) 140 144 print 'Interpolated elevation at (%.2f, %.2f) is %.2f' %(x0, y0, w_i) 145 141 146 142 147 # … … 145 150 146 151 w_max = w0 147 stagestep = []148 152 for t in domain.evolve(yieldstep = 1, finaltime = 300): 149 153 domain.write_time() … … 151 155 w = domain.get_maximum_inundation_elevation() 152 156 x, y = domain.get_maximum_inundation_location() 153 print ' Coastline elevation = %.2f at (%.2f, %.2f)' %(w, x, y) 154 #w_i = domain.get_quantity('stage').get_values(interpolation_points=[[x,y]]) 155 #print ' Interpolated elevation at (%.2f, %.2f) is %.2f' %(x, y, w_i) 156 157 print ' Coastline elevation = %.2f at (x,y)=(%.2f, %.2f)' %(w, x, y) 157 158 158 159 if w > w_max: … … 161 162 y_max = y 162 163 163 # Let's find the maximum runup here working directly with the quantities,164 # and stop when it has been detected.165 166 # 1 Find coastline as x where z==0167 # 2 Workout travel time to coastline168 # 3 Find min x where h>0 over all t.169 # 4 Perhaps do this across a number of ys170 164 171 165 print 'Max coastline elevation = %.2f at (%.2f, %.2f)' %(w_max, x_max, y_max) 172 173 print 'Run up distance  %.2f' %sqrt( (xx0)**2 + (yy0)**2 ) 166 print 'Run up distance = %.2f' %sqrt( (x_maxx0)**2 + (y_maxy0)**2 ) 174 167 175 168 176 169 177 170 # 178 # Interrogate solution179 # 171 # Interrogate further 172 # 180 173 181 '''182 # Define line of gauges through center of domain183 def gauge_line(west,east,north,south):184 from Numeric import arange185 x_vector = arange(west,600, 10) # Gauges every 1 meter from west to 600m from western bdry186 y = (north+south)/2.187 188 gauges = []189 for x in x_vector:190 gauges.append([x,y])191 192 return gauges, x_vector193 194 gauges, x_vector = gauge_line(west,east,north,south)195 196 # Obtain interpolated timeseries at gauges197 f = file_function(domain.get_name()+'.sww',198 quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'],199 interpolation_points = gauges,200 verbose = True,201 use_cache = True)202 203 204 # Find runup distance from western boundary through a linear search205 max_stage = []206 min_stage = []207 runup_point = west208 coastline = east209 for k, g in enumerate(gauges):210 z = f(0, point_id=k)[1] # Elevation211 212 min_w = sys.maxint213 max_w = min_w214 for i, t in enumerate(f.get_time()):215 w = f(t, point_id = k)[0]216 if w > max_w: max_w = w217 if w < min_w: min_w = w218 219 if max_wz <= 0.01: # Find first gauge where max depth > eps (runup)220 runup_point = g[0]221 222 if min_wz <= 0.01: # Find first gauge where min depth > eps (coastline)223 coastline = g[0]224 225 max_stage.append(max_w)226 min_stage.append(min_w)227 228 229 # Print230 print 'wave height [m]: ', amplitude231 runup_height = topography([runup_point], [(north+south)/2.])[0]232 print 'run up height [m]: ', runup_height233 234 runup_distance = runup_pointcoastline235 print 'run up distance from coastline [m]: ', runup_distance236 237 print 'Coastline (meters form west): ', coastline238 239 '''240 174 # Generate time series of "gauge" situated at right hand boundary 241 175 from anuga.abstract_2d_finite_volumes.util import sww2timeseries … … 259 193 verbose = True) 260 194 261 # Stop here262 import sys; sys.exit()263 264 # Take snapshots and plot265 ion()266 figure(1)267 plot(x_vector, topography(x_vector,(north+south)/2.), 'r')268 xlabel('x')269 ylabel('Elevation')270 #legend(('Max stage', 'Min stage', 'Elevation'), shadow=True, loc='upper right')271 title('Stage snapshots (t=0, 10, ...) for gauge line')272 grid()273 hold(True)274 275 for i, t in enumerate(f.get_time()):276 if i % 10 == 0:277 # Take only some timesteps to avoid clutter278 stages = []279 for k, g in enumerate(gauges):280 w = f(t, point_id = k)[0]281 stages.append(w)282 283 plot(x_vector, stages, 'b')284 285 savefig('snapshots')286 195 287 196 288 197 289 # Store290 filename = 'maxrunup'+str(amplitude)+'.csv'291 fid = open(filename,'w')292 s = 'Waveheight,Runup distance,Runup height\n'293 fid.write(s)294 295 s = '%.2f,%.2f,%.2f\n' %(amplitude, runup_distance, runup_height)296 fid.write(s)297 298 fid.close()299 300 # Plot max runup etc301 ion()302 figure(1)303 plot(x_vector, max_stage, 'g+',304 x_vector, min_stage, 'b+',305 x_vector, topography(x_vector,(north+south)/2.), 'r')306 xlabel('x')307 ylabel('stage')308 legend(('Max stage', 'Min stage', 'Elevation'), shadow=True, loc='upper right')309 title('Maximum stage for gauge line')310 grid()311 #axis([33000, 47000, 1000, 3000])312 savefig('max_stage')313 314 close('all')315 316 317
Note: See TracChangeset
for help on using the changeset viewer.