[7578] | 1 | """Simple water flow example using ANUGA |
---|
| 2 | """ |
---|
| 3 | |
---|
| 4 | #------------------------------------------------------------------------------ |
---|
| 5 | # Import necessary modules |
---|
| 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 * |
---|
| 12 | import numpy |
---|
| 13 | from matplotlib import * |
---|
| 14 | from pylab import * |
---|
| 15 | import csv |
---|
| 16 | |
---|
| 17 | #Parameters |
---|
| 18 | |
---|
| 19 | filename = "WORKING-RIP-LAB_Expt-Geometry_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 |
---|
| 24 | halfchannelwidth = 5 |
---|
| 25 | bank_slope = 0.1 |
---|
| 26 | simulation_length = 1 |
---|
| 27 | timestep = 1 |
---|
| 28 | |
---|
| 29 | #------------------------------------------------------------------------------ |
---|
| 30 | # Setup computational domain |
---|
| 31 | #------------------------------------------------------------------------------ |
---|
| 32 | length = 120 |
---|
| 33 | width = 170 |
---|
| 34 | seafloor_resolution = 60.0 # Resolution: Max area of triangles in the mesh for seafloor |
---|
| 35 | feature_resolution = 1.0 |
---|
| 36 | beach_resolution = 10.0 |
---|
| 37 | |
---|
| 38 | sea_boundary_polygon = [[0,0],[length,0],[length,width],[0,width]] |
---|
| 39 | feature_boundary_polygon = [[0,100],[length,100],[length,150],[0,150]] |
---|
| 40 | beach_interior_polygon = [[0,150],[length,150],[length,width],[0,width]] |
---|
| 41 | |
---|
| 42 | meshname = str(filename)+'.msh' |
---|
| 43 | |
---|
| 44 | #interior regions |
---|
| 45 | feature_regions = [[feature_boundary_polygon, feature_resolution], |
---|
| 46 | [beach_interior_polygon, beach_resolution]] |
---|
| 47 | |
---|
| 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) |
---|
| 59 | domain.set_name(filename) # Output name |
---|
| 60 | print domain.statistics() |
---|
| 61 | |
---|
| 62 | #------------------------------------------------------------------------------ |
---|
| 63 | # Setup initial conditions |
---|
| 64 | #------------------------------------------------------------------------------ |
---|
| 65 | |
---|
| 66 | def topography(x,y): |
---|
| 67 | """Complex topography defined by a function of vectors x and y.""" |
---|
| 68 | |
---|
| 69 | #general slope and buildings |
---|
| 70 | z=0.05*(y-(location_of_shore)) |
---|
| 71 | |
---|
| 72 | N = len(x) |
---|
| 73 | for i in range(N): |
---|
| 74 | if y[i] < 25: |
---|
| 75 | z[i] = (0.2*(y[i]-25)) + 0.05*(y[i]-(location_of_shore)) |
---|
| 76 | for i in range(N): |
---|
| 77 | if y[i]>150: |
---|
| 78 | z[i] = (0.1*(y[i]-150)) + 0.05*(y[i] - (location_of_shore)) |
---|
| 79 | |
---|
| 80 | return z |
---|
| 81 | |
---|
| 82 | |
---|
| 83 | def topography3(x,y): |
---|
| 84 | z=0*x |
---|
| 85 | |
---|
| 86 | 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): |
---|
| 89 | 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: |
---|
| 92 | z[i] += sandbar*(cos((y[i]-118)/steepness)) |
---|
| 93 | return z |
---|
| 94 | |
---|
| 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 |
---|
| 101 | |
---|
| 102 | #------------------------------------------------------------------------------ |
---|
| 103 | # Setup boundary conditions |
---|
| 104 | #------------------------------------------------------------------------------ |
---|
| 105 | Bi = Dirichlet_boundary([0.4, 0, 0]) # Inflow |
---|
| 106 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
| 107 | Bo = Dirichlet_boundary([-5, 0, 0]) # Outflow |
---|
| 108 | |
---|
| 109 | |
---|
| 110 | amplitude = 0.4 #amplitude of wave (wave height m) |
---|
| 111 | period = 5 #wave period (sec) |
---|
| 112 | |
---|
| 113 | def wave(t): |
---|
| 114 | |
---|
| 115 | A = amplitude # Amplitude [m] (Wave height) |
---|
| 116 | T = period # Wave period [s] |
---|
| 117 | |
---|
| 118 | if t < 30000000000: |
---|
| 119 | return [A*sin(2*pi*t/T) + 1, 0, 0] |
---|
| 120 | else: |
---|
| 121 | return [0.0, 0, 0] |
---|
| 122 | |
---|
| 123 | Bt = Time_boundary(domain, f=wave) |
---|
| 124 | |
---|
| 125 | |
---|
| 126 | domain.set_boundary({'left': Br, 'right': Br, 'top': Bo, 'bottom': Bt}) |
---|
| 127 | |
---|
| 128 | #------------------------------------------------------------------------------ |
---|
| 129 | # Evolve system through time |
---|
| 130 | #------------------------------------------------------------------------------ |
---|
| 131 | |
---|
| 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)): |
---|
| 141 | 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) |
---|
| 144 | depth = S-E |
---|
| 145 | |
---|
| 146 | uh = domain.get_quantity('xmomentum').get_values(interpolation_points=list1) |
---|
| 147 | vh = domain.get_quantity('ymomentum').get_values(interpolation_points=list1) |
---|
| 148 | u += uh/depth |
---|
| 149 | 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=",") ] |
---|
| 173 | U = u_average.tolist() |
---|
| 174 | V = v_average.tolist() |
---|
| 175 | |
---|
| 176 | print "U = ", U |
---|
| 177 | print "U has type", type(U) |
---|
| 178 | |
---|
| 179 | |
---|
| 180 | from matplotlib import * |
---|
| 181 | from pylab import * |
---|
| 182 | |
---|
| 183 | figure() |
---|
| 184 | quiver(X,Y,U,V) |
---|
| 185 | show() |
---|
| 186 | |
---|
| 187 | |
---|