[7578] | 1 | """Simple water flow example using ANUGA |
---|
| 2 | """ |
---|
| 3 | |
---|
| 4 | #------------------------------------------------------------------------------ |
---|
| 5 | # Import necessary modules |
---|
| 6 | #------------------------------------------------------------------------------ |
---|
[7579] | 7 | from anuga.interface import create_domain_from_regions |
---|
| 8 | from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary |
---|
| 9 | from anuga.shallow_water.shallow_water_domain import Reflective_boundary |
---|
| 10 | from anuga.shallow_water.shallow_water_domain import Time_boundary |
---|
| 11 | from pylab import figure, quiver, show, cos, sin, pi |
---|
[7578] | 12 | import numpy |
---|
| 13 | import csv |
---|
[7579] | 14 | import time |
---|
[7578] | 15 | |
---|
| 16 | |
---|
[7579] | 17 | #------------------------------------------------------------------------------ |
---|
| 18 | # Parameters |
---|
| 19 | #------------------------------------------------------------------------------ |
---|
[7578] | 20 | filename = "WORKING-RIP-LAB_Expt-Geometry_Triangular_Mesh" |
---|
[7579] | 21 | location_of_shore = 140 # The position along the y axis of the shorefront |
---|
| 22 | sandbar = 1.2 # Height of sandbar |
---|
| 23 | sealevel = 0 # Height of coast above sea level |
---|
| 24 | steepness = 8000 # Period of sandbar - |
---|
| 25 | # larger number gives smoother slope - longer period |
---|
[7578] | 26 | halfchannelwidth = 5 |
---|
| 27 | bank_slope = 0.1 |
---|
| 28 | simulation_length = 1 |
---|
| 29 | timestep = 1 |
---|
| 30 | |
---|
[7579] | 31 | |
---|
[7578] | 32 | #------------------------------------------------------------------------------ |
---|
| 33 | # Setup computational domain |
---|
| 34 | #------------------------------------------------------------------------------ |
---|
| 35 | length = 120 |
---|
| 36 | width = 170 |
---|
[7579] | 37 | seafloor_resolution = 60.0 # Resolution: Max area of triangles in the mesh |
---|
[7578] | 38 | feature_resolution = 1.0 |
---|
| 39 | beach_resolution = 10.0 |
---|
| 40 | |
---|
| 41 | sea_boundary_polygon = [[0,0],[length,0],[length,width],[0,width]] |
---|
| 42 | feature_boundary_polygon = [[0,100],[length,100],[length,150],[0,150]] |
---|
| 43 | beach_interior_polygon = [[0,150],[length,150],[length,width],[0,width]] |
---|
| 44 | |
---|
| 45 | meshname = str(filename)+'.msh' |
---|
| 46 | |
---|
[7579] | 47 | # Interior regions |
---|
[7578] | 48 | feature_regions = [[feature_boundary_polygon, feature_resolution], |
---|
| 49 | [beach_interior_polygon, beach_resolution]] |
---|
| 50 | |
---|
[7579] | 51 | domain = create_domain_from_regions(sea_boundary_polygon, |
---|
| 52 | boundary_tags={'bottom': [0], |
---|
| 53 | 'right' : [1], |
---|
| 54 | 'top' : [2], |
---|
| 55 | 'left': [3]}, |
---|
| 56 | maximum_triangle_area = seafloor_resolution, |
---|
| 57 | mesh_filename = meshname, |
---|
| 58 | interior_regions = feature_regions, |
---|
| 59 | use_cache = True, |
---|
| 60 | verbose = True) |
---|
| 61 | |
---|
[7578] | 62 | domain.set_name(filename) # Output name |
---|
| 63 | print domain.statistics() |
---|
| 64 | |
---|
[7579] | 65 | |
---|
[7578] | 66 | #------------------------------------------------------------------------------ |
---|
| 67 | # Setup initial conditions |
---|
| 68 | #------------------------------------------------------------------------------ |
---|
| 69 | def topography(x,y): |
---|
| 70 | """Complex topography defined by a function of vectors x and y.""" |
---|
| 71 | |
---|
[7579] | 72 | # General slope and buildings |
---|
[7578] | 73 | z=0.05*(y-(location_of_shore)) |
---|
| 74 | |
---|
| 75 | N = len(x) |
---|
| 76 | for i in range(N): |
---|
| 77 | if y[i] < 25: |
---|
| 78 | z[i] = (0.2*(y[i]-25)) + 0.05*(y[i]-(location_of_shore)) |
---|
| 79 | for i in range(N): |
---|
| 80 | if y[i]>150: |
---|
[7579] | 81 | z[i] = (0.1*(y[i]-150)) + 0.05*(y[i]-(location_of_shore)) |
---|
[7578] | 82 | |
---|
| 83 | return z |
---|
| 84 | |
---|
| 85 | |
---|
| 86 | def topography3(x,y): |
---|
| 87 | z=0*x |
---|
| 88 | |
---|
| 89 | N = len(x) |
---|
[7579] | 90 | |
---|
| 91 | # It would be great with a comment about what this does |
---|
[7578] | 92 | for i in range(N): |
---|
[7579] | 93 | ymin = -1*(bank_slope)*x[i] + 112 |
---|
| 94 | ymax = -1*(bank_slope)*x[i] + 124 |
---|
| 95 | xmin = 0 |
---|
| 96 | xmax = (length/2)-halfchannelwidth |
---|
| 97 | if ymin < y[i] < ymax and xmin < x[i]< xmax: |
---|
[7578] | 98 | z[i] += sandbar*((cos((y[i]-118)/steepness))) |
---|
[7579] | 99 | |
---|
| 100 | # It would be great with a comment about what this does |
---|
[7578] | 101 | for i in range(N): |
---|
[7579] | 102 | ymin = -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+112) |
---|
| 103 | ymax = -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+124) |
---|
| 104 | xmin = (length/2)+halfchannelwidth |
---|
| 105 | xmax = 183 |
---|
| 106 | if ymin < y[i] < ymax and xmin < x[i] < xmax: |
---|
[7578] | 107 | z[i] += sandbar*(cos((y[i]-118)/steepness)) |
---|
[7579] | 108 | |
---|
[7578] | 109 | return z |
---|
| 110 | |
---|
[7579] | 111 | domain.set_quantity('elevation', topography) # Apply base elevation function |
---|
| 112 | domain.add_quantity('elevation', topography3) # Add elevation modification |
---|
| 113 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
| 114 | domain.set_quantity('stage', 0) # Constant initial condition at mean sea level |
---|
[7578] | 115 | |
---|
| 116 | |
---|
| 117 | #------------------------------------------------------------------------------ |
---|
| 118 | # Setup boundary conditions |
---|
| 119 | #------------------------------------------------------------------------------ |
---|
| 120 | Bi = Dirichlet_boundary([0.4, 0, 0]) # Inflow |
---|
| 121 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
| 122 | Bo = Dirichlet_boundary([-5, 0, 0]) # Outflow |
---|
| 123 | |
---|
| 124 | def wave(t): |
---|
[7579] | 125 | """Define wave driving the system |
---|
| 126 | """ |
---|
| 127 | |
---|
| 128 | A = 0.4 # Amplitude of wave [m] (wave height) |
---|
| 129 | T = 5 # Wave period [s] |
---|
[7578] | 130 | |
---|
| 131 | if t < 30000000000: |
---|
| 132 | return [A*sin(2*pi*t/T) + 1, 0, 0] |
---|
| 133 | else: |
---|
| 134 | return [0.0, 0, 0] |
---|
| 135 | |
---|
| 136 | Bt = Time_boundary(domain, f=wave) |
---|
| 137 | |
---|
| 138 | |
---|
| 139 | domain.set_boundary({'left': Br, 'right': Br, 'top': Bo, 'bottom': Bt}) |
---|
| 140 | |
---|
[7579] | 141 | |
---|
[7578] | 142 | #------------------------------------------------------------------------------ |
---|
| 143 | # Evolve system through time |
---|
| 144 | #------------------------------------------------------------------------------ |
---|
| 145 | |
---|
[7579] | 146 | # Read in gauge locations for interpolation and convert them to floats |
---|
| 147 | gauge_file = open('New_gauges.csv', 'r') |
---|
| 148 | G = [(float(x[0]), float(x[1])) for x in csv.reader(gauge_file, |
---|
| 149 | dialect='excel', |
---|
| 150 | delimiter=',')] |
---|
| 151 | gauges = numpy.array(G) |
---|
| 152 | number_of_gauges = len(gauges) |
---|
| 153 | print gauges |
---|
[7578] | 154 | |
---|
[7579] | 155 | # Allocate space for velocity values |
---|
| 156 | u = numpy.zeros(number_of_gauges) |
---|
| 157 | v = numpy.zeros(number_of_gauges) |
---|
[7578] | 158 | |
---|
[7579] | 159 | t0 = time.time() |
---|
| 160 | for t in domain.evolve(yieldstep = timestep, finaltime = simulation_length): |
---|
[7578] | 161 | print domain.timestepping_statistics() |
---|
[7579] | 162 | S = domain.get_quantity('stage').get_values(interpolation_points=gauges) |
---|
| 163 | E = domain.get_quantity('elevation').get_values(interpolation_points=gauges) |
---|
[7578] | 164 | depth = S-E |
---|
| 165 | |
---|
[7579] | 166 | uh = domain.get_quantity('xmomentum').get_values(interpolation_points=gauges) |
---|
| 167 | vh = domain.get_quantity('ymomentum').get_values(interpolation_points=gauges) |
---|
[7578] | 168 | u += uh/depth |
---|
| 169 | v += vh/depth |
---|
[7579] | 170 | print 'Evolution took %.2f seconds' % (time.time()-t0) |
---|
[7578] | 171 | |
---|
[7579] | 172 | |
---|
| 173 | #------------------------------------------------------------------------------ |
---|
| 174 | # Post processing |
---|
| 175 | #------------------------------------------------------------------------------ |
---|
[7578] | 176 | |
---|
[7579] | 177 | # Use this |
---|
[7578] | 178 | |
---|
[7579] | 179 | # from anuga.fit_interpolate.interpolate import Interpolation_function |
---|
[7578] | 180 | |
---|
[7579] | 181 | # # Get mesh and quantities from sww file |
---|
| 182 | # X = get_mesh_and_quantities_from_file(filename, |
---|
| 183 | # quantities=quantity_names, |
---|
| 184 | # verbose=verbose) |
---|
| 185 | # mesh, quantities, time = X |
---|
[7578] | 186 | |
---|
[7579] | 187 | # # Find all intersections and associated triangles. |
---|
| 188 | # segments = mesh.get_intersecting_segments(polyline, verbose=verbose) |
---|
[7578] | 189 | |
---|
[7579] | 190 | # # Get midpoints |
---|
| 191 | # interpolation_points = segment_midpoints(segments) |
---|
[7578] | 192 | |
---|
[7579] | 193 | # # Interpolate |
---|
| 194 | # if verbose: |
---|
| 195 | # log.critical('Interpolating - total number of interpolation points = %d' |
---|
| 196 | # % len(interpolation_points)) |
---|
| 197 | |
---|
| 198 | # I = Interpolation_function(time, |
---|
| 199 | # quantities, |
---|
| 200 | # quantity_names=quantity_names, |
---|
| 201 | # vertex_coordinates=mesh.nodes, |
---|
| 202 | # triangles=mesh.triangles, |
---|
| 203 | # interpolation_points=interpolation_points, |
---|
| 204 | # verbose=verbose) |
---|
| 205 | |
---|
| 206 | # return segments, I |
---|
| 207 | |
---|
| 208 | |
---|
| 209 | |
---|
| 210 | |
---|
| 211 | |
---|
| 212 | |
---|
| 213 | n_time_intervals = simulation_length/timestep |
---|
| 214 | u_average = u/n_time_intervals |
---|
| 215 | v_average = v/n_time_intervals |
---|
| 216 | |
---|
| 217 | #print "there were", n_time_intervals, "time steps" |
---|
| 218 | |
---|
| 219 | #print "sum y velocity", v |
---|
| 220 | #print "average y velocity", v_average |
---|
| 221 | #print "sum x velocity", u |
---|
| 222 | #print "average x velocity", u_average |
---|
| 223 | |
---|
| 224 | x_output = file('x_velocity.txt', 'w') |
---|
| 225 | y_output = file('y_velocity.txt', 'w') |
---|
| 226 | |
---|
| 227 | #print >> x_output, " " |
---|
| 228 | #print >> y_output, " " |
---|
| 229 | |
---|
| 230 | #print >> x_output, u_average |
---|
| 231 | #print >> y_output, v_average |
---|
| 232 | |
---|
| 233 | |
---|
| 234 | X = gauges[:,0] |
---|
| 235 | Y = gauges[:,1] |
---|
| 236 | |
---|
[7578] | 237 | U = u_average.tolist() |
---|
| 238 | V = v_average.tolist() |
---|
| 239 | |
---|
[7579] | 240 | #print "U = ", U |
---|
| 241 | #print "U has type", type(U) |
---|
[7578] | 242 | |
---|
| 243 | |
---|
| 244 | figure() |
---|
| 245 | quiver(X,Y,U,V) |
---|
| 246 | show() |
---|
| 247 | |
---|
| 248 | |
---|