Changeset 7579
 Timestamp:
 Dec 7, 2009, 6:18:23 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/classroom/WORKINGRipCurrent.py
r7578 r7579 5 5 # Import necessary modules 6 6 # 7 from anuga. abstract_2d_finite_volumes.mesh_factory import *8 from anuga.shallow_water import *9 from anuga.shallow_water.shallow_water_domain import *10 from math import *11 from numpy import *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 12 12 import numpy 13 from matplotlib import *14 from pylab import *15 13 import csv 16 17 #Parameters 18 14 import time 15 16 17 # 18 # Parameters 19 # 19 20 filename = "WORKINGRIPLAB_ExptGeometry_Triangular_Mesh" 20 location_of_shore = 140 #the position along the y axis of the shorefront 21 sandbar = 1.2 #height of sandbar 22 sealevel = 0 #height of coast above sea level 23 steepness = 8000 #period of sandbar  larger number gives smoother slope  longer period 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 24 26 halfchannelwidth = 5 25 27 bank_slope = 0.1 … … 27 29 timestep = 1 28 30 31 29 32 # 30 33 # Setup computational domain … … 32 35 length = 120 33 36 width = 170 34 seafloor_resolution = 60.0 # Resolution: Max area of triangles in the mesh for seafloor37 seafloor_resolution = 60.0 # Resolution: Max area of triangles in the mesh 35 38 feature_resolution = 1.0 36 39 beach_resolution = 10.0 … … 42 45 meshname = str(filename)+'.msh' 43 46 44 # interior regions47 # Interior regions 45 48 feature_regions = [[feature_boundary_polygon, feature_resolution], 46 49 [beach_interior_polygon, beach_resolution]] 47 50 48 create_mesh_from_regions(sea_boundary_polygon,49 boundary_tags={'bottom': [0],50 'right' : [1],51 'top' : [2],52 'left': [3]},53 maximum_triangle_area = seafloor_resolution,54 filename = meshname,55 interior_regions = feature_regions,56 use_cache = False,57 verbose = False)58 domain = Domain(meshname, use_cache=True, verbose=True) 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 59 62 domain.set_name(filename) # Output name 60 63 print domain.statistics() 61 64 65 62 66 # 63 67 # Setup initial conditions 64 68 # 65 66 69 def topography(x,y): 67 70 """Complex topography defined by a function of vectors x and y.""" 68 71 69 # general slope and buildings72 # General slope and buildings 70 73 z=0.05*(y(location_of_shore)) 71 74 … … 76 79 for i in range(N): 77 80 if y[i]>150: 78 z[i] = (0.1*(y[i]150)) + 0.05*(y[i] (location_of_shore))81 z[i] = (0.1*(y[i]150)) + 0.05*(y[i](location_of_shore)) 79 82 80 83 return z … … 85 88 86 89 N = len(x) 87 for i in range(N): 88 if 1*(bank_slope)*x[i] + 112 < y[i] < 1*(bank_slope)*x[i] + 124 and 0<x[i]<((length/2)halfchannelwidth): 90 91 # It would be great with a comment about what this does 92 for i in range(N): 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: 89 98 z[i] += sandbar*((cos((y[i]118)/steepness))) 90 for i in range(N): 91 if 1*(bank_slope)*(x[i](length/2)) + (1*(bank_slope)*(length/2)+112) < y[i] < 1*(bank_slope)*(x[i](length/2)) + (1*(bank_slope)*(length/2)+124) and ((length/2)+halfchannelwidth)<x[i]<183: 99 100 # It would be great with a comment about what this does 101 for i in range(N): 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: 92 107 z[i] += sandbar*(cos((y[i]118)/steepness)) 108 93 109 return z 94 110 95 domain.set_quantity('elevation', topography) # elevation is a function 96 domain.set_quantity('friction', 0.01) # Constant friction 97 98 domain.add_quantity('elevation', topography3) 99 100 domain.set_quantity('stage', 0) # Dry initial condition 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 115 101 116 102 117 # … … 107 122 Bo = Dirichlet_boundary([5, 0, 0]) # Outflow 108 123 109 110 amplitude = 0.4 #amplitude of wave (wave height m)111 period = 5 #wave period (sec)112 113 124 def wave(t): 114 115 A = amplitude # Amplitude [m] (Wave height) 116 T = period # Wave period [s] 125 """Define wave driving the system 126 """ 127 128 A = 0.4 # Amplitude of wave [m] (wave height) 129 T = 5 # Wave period [s] 117 130 118 131 if t < 30000000000: … … 126 139 domain.set_boundary({'left': Br, 'right': Br, 'top': Bo, 'bottom': Bt}) 127 140 141 128 142 # 129 143 # Evolve system through time 130 144 # 131 145 132 list1 = [x for x in csv.reader(open('New_gauges.csv','r'),dialect='excel',delimiter=",") ] 133 134 print list1 135 136 u = numpy.zeros(len(list1), 'd') 137 v = numpy.zeros(len(list1), 'd') 138 139 140 for t in domain.evolve(yieldstep = (timestep), finaltime = (simulation_length)): 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 154 155 # Allocate space for velocity values 156 u = numpy.zeros(number_of_gauges) 157 v = numpy.zeros(number_of_gauges) 158 159 t0 = time.time() 160 for t in domain.evolve(yieldstep = timestep, finaltime = simulation_length): 141 161 print domain.timestepping_statistics() 142 S = domain.get_quantity('stage').get_values(interpolation_points= list1)143 E = domain.get_quantity('elevation').get_values(interpolation_points= list1)162 S = domain.get_quantity('stage').get_values(interpolation_points=gauges) 163 E = domain.get_quantity('elevation').get_values(interpolation_points=gauges) 144 164 depth = SE 145 165 146 uh = domain.get_quantity('xmomentum').get_values(interpolation_points= list1)147 vh = domain.get_quantity('ymomentum').get_values(interpolation_points= list1)166 uh = domain.get_quantity('xmomentum').get_values(interpolation_points=gauges) 167 vh = domain.get_quantity('ymomentum').get_values(interpolation_points=gauges) 148 168 u += uh/depth 149 169 v += vh/depth 150 151 n_time_intervals = (simulation_length)/(timestep) 152 u_average = u / (n_time_intervals) 153 v_average = v / (n_time_intervals) 154 155 print "there were", n_time_intervals, "time steps" 156 157 print "sum y velocity", v 158 print "average y velocity", v_average 159 print "sum x velocity", u 160 print "average x velocity", u_average 161 162 x_output = file("x_velocity.txt", 'w') 163 y_output = file("y_velocity.txt", 'w') 164 165 print >> x_output, " " 166 print >> y_output, " " 167 168 print >> x_output, u_average 169 print >> y_output, v_average 170 171 X = [x[0] for x in csv.reader(open('New_gauges.csv','r'),dialect='excel',delimiter=",") ] 172 Y = [x[1] for x in csv.reader(open('New_gauges.csv','r'),dialect='excel',delimiter=",") ] 170 print 'Evolution took %.2f seconds' % (time.time()t0) 171 172 173 # 174 # Post processing 175 # 176 177 # Use this 178 179 # from anuga.fit_interpolate.interpolate import Interpolation_function 180 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 186 187 # # Find all intersections and associated triangles. 188 # segments = mesh.get_intersecting_segments(polyline, verbose=verbose) 189 190 # # Get midpoints 191 # interpolation_points = segment_midpoints(segments) 192 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 173 237 U = u_average.tolist() 174 238 V = v_average.tolist() 175 239 176 print "U = ", U 177 print "U has type", type(U) 178 179 180 from matplotlib import * 181 from pylab import * 240 #print "U = ", U 241 #print "U has type", type(U) 242 182 243 183 244 figure()
Note: See TracChangeset
for help on using the changeset viewer.