1 | """Validation of the AnuGA implementation of the shallow water wave equation. |
---|
2 | |
---|
3 | This script sets up 2D circular island..... |
---|
4 | |
---|
5 | """ |
---|
6 | |
---|
7 | # Module imports |
---|
8 | from anuga.shallow_water import Domain |
---|
9 | from anuga.shallow_water import Reflective_boundary |
---|
10 | from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary |
---|
11 | from anuga.shallow_water import Dirichlet_boundary, Time_boundary |
---|
12 | from anuga.shallow_water.data_manager import get_maximum_inundation_data, start_screen_catcher |
---|
13 | from anuga.abstract_2d_finite_volumes.util import file_function, sww2csv_gauges |
---|
14 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
15 | from anuga.utilities.polygon import read_polygon#, plot_polygons |
---|
16 | #from cmath import cos,sin,cosh,pi |
---|
17 | from math import cos,pi,sin,tan,sqrt |
---|
18 | from Numeric import array, zeros, Float, allclose,resize |
---|
19 | from anuga.shallow_water.data_manager import csv2dict |
---|
20 | |
---|
21 | #------------------------- |
---|
22 | # Create Domain from mesh |
---|
23 | #------------------------- |
---|
24 | |
---|
25 | #start_screen_catcher() |
---|
26 | |
---|
27 | def get_xy(x=0,y=0,r=1,angle=1): |
---|
28 | #return x and y co-ordinates |
---|
29 | x1 = x + cos(angle)*r |
---|
30 | y1 = y + sin(-angle)*r |
---|
31 | # print "get",x,y,r,x1,y1,"angle",angle, cos(angle),sin(angle) |
---|
32 | return x1,y1, |
---|
33 | |
---|
34 | def prepare_timeboundary(filename): |
---|
35 | """Converting benchmark 2 time series to NetCDF tms file. |
---|
36 | This is a 'throw-away' code taylor made for files like |
---|
37 | 'Benchmark_2_input.txt' from the LWRU2 benchmark |
---|
38 | """ |
---|
39 | |
---|
40 | textversion = filename[:-4] + '.txt' |
---|
41 | |
---|
42 | print 'Preparing time boundary from %s' %textversion |
---|
43 | from Numeric import array |
---|
44 | |
---|
45 | fid = open(textversion) |
---|
46 | |
---|
47 | #Skip first line |
---|
48 | line = fid.readline() |
---|
49 | |
---|
50 | #Read remaining lines |
---|
51 | lines = fid.readlines() |
---|
52 | fid.close() |
---|
53 | |
---|
54 | |
---|
55 | N = len(lines) |
---|
56 | T = zeros(N, Float) #Time |
---|
57 | Q = zeros(N, Float) #Values |
---|
58 | |
---|
59 | for i, line in enumerate(lines): |
---|
60 | fields = line.split() |
---|
61 | # print 'fields',i,line |
---|
62 | |
---|
63 | T[i] = float(fields[0]) |
---|
64 | Q[i] = float(fields[1]) |
---|
65 | |
---|
66 | |
---|
67 | #Create tms file |
---|
68 | from Scientific.IO.NetCDF import NetCDFFile |
---|
69 | |
---|
70 | print 'Writing to', filename |
---|
71 | fid = NetCDFFile(filename[:-4] + '.tms', 'w') |
---|
72 | |
---|
73 | fid.institution = 'Geoscience Australia' |
---|
74 | fid.description = 'Input wave for Benchmark 2' |
---|
75 | fid.starttime = 0.0 |
---|
76 | fid.createDimension('number_of_timesteps', len(T)) |
---|
77 | fid.createVariable('time', Float, ('number_of_timesteps',)) |
---|
78 | fid.variables['time'][:] = T |
---|
79 | |
---|
80 | fid.createVariable('stage', Float, ('number_of_timesteps',)) |
---|
81 | fid.variables['stage'][:] = Q[:] |
---|
82 | |
---|
83 | fid.createVariable('xmomentum', Float, ('number_of_timesteps',)) |
---|
84 | fid.variables['xmomentum'][:] = 0.0 |
---|
85 | |
---|
86 | fid.createVariable('ymomentum', Float, ('number_of_timesteps',)) |
---|
87 | fid.variables['ymomentum'][:] = 0.0 |
---|
88 | |
---|
89 | fid.close() |
---|
90 | |
---|
91 | |
---|
92 | |
---|
93 | angles1 = (0.0, 22.5, 45.0, 67.5, 90.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0, |
---|
94 | 247.5, 270.0, 292.5, 315.0, 337.5) |
---|
95 | |
---|
96 | #x = 12.96 |
---|
97 | #y = 13.80 |
---|
98 | center_x = 12.96 |
---|
99 | center_y = 13.80 |
---|
100 | |
---|
101 | #wavelenght = |
---|
102 | |
---|
103 | r1 = 1.5 |
---|
104 | r2 = 2.5 |
---|
105 | d=0.75 |
---|
106 | res=.05 |
---|
107 | #res=1. |
---|
108 | L=4.11 |
---|
109 | |
---|
110 | |
---|
111 | #create polygons (circles) to have higher resolution |
---|
112 | poly = [] |
---|
113 | poly1 = [] |
---|
114 | internal_poly=[] |
---|
115 | for i,angle in enumerate(angles1): |
---|
116 | #convert Degs to Rads |
---|
117 | angle = ((angle-180)/180.0)*pi |
---|
118 | |
---|
119 | p = get_xy(center_x,center_y,3.6,angle) |
---|
120 | poly1.append([p[0],p[1]]) |
---|
121 | |
---|
122 | poly.append(poly1) |
---|
123 | internal_poly.append([poly1,.1*res]) |
---|
124 | |
---|
125 | poly2 = [] |
---|
126 | for i,angle in enumerate(angles1): |
---|
127 | #convert Degs to Rads |
---|
128 | angle = ((angle-180)/180.0)*pi |
---|
129 | p = get_xy(center_x,center_y,2.4,angle) |
---|
130 | poly2.append([p[0],p[1]]) |
---|
131 | |
---|
132 | poly.append(poly2) |
---|
133 | internal_poly.append([poly2,.01*res]) |
---|
134 | |
---|
135 | poly3 = [] |
---|
136 | for i,angle in enumerate(angles1): |
---|
137 | #convert Degs to Rads |
---|
138 | angle = ((angle-180)/180.0)*pi |
---|
139 | p = get_xy(center_x,center_y,1.6,angle) |
---|
140 | poly3.append([p[0],p[1]]) |
---|
141 | |
---|
142 | poly.append(poly3) |
---|
143 | internal_poly.append([poly3,1*res]) |
---|
144 | |
---|
145 | poly4 = [] |
---|
146 | for i,angle in enumerate(angles1): |
---|
147 | #convert Degs to Rads |
---|
148 | angle = ((angle-180)/180.0)*pi |
---|
149 | p = get_xy(center_x,center_y,1.1,angle) |
---|
150 | poly4.append([p[0],p[1]]) |
---|
151 | |
---|
152 | poly.append(poly4) |
---|
153 | internal_poly.append([poly2,.01*res]) |
---|
154 | |
---|
155 | #plot_polygons(poly,'poly.png') |
---|
156 | |
---|
157 | #poly_all= [[0,7],[0,25],[30,25],[30,7]] |
---|
158 | poly_all= [[0,L],[0,25],[30,25],[30,L]] |
---|
159 | create_mesh_from_regions(poly_all, |
---|
160 | boundary_tags={'wall': [0,2],'wave': [3], 'buffer':[1]}, |
---|
161 | maximum_triangle_area=1, |
---|
162 | interior_regions=internal_poly, |
---|
163 | filename='temp.msh', |
---|
164 | use_cache=False, |
---|
165 | verbose=True) |
---|
166 | |
---|
167 | domain = Domain('temp.msh', use_cache=True, verbose=True) |
---|
168 | print domain.statistics() |
---|
169 | |
---|
170 | #------------------------- |
---|
171 | # Initial Conditions |
---|
172 | #------------------------- |
---|
173 | domain.set_quantity('friction', 0.0) |
---|
174 | water_height = 0.322 |
---|
175 | domain.set_quantity('stage', water_height) |
---|
176 | |
---|
177 | |
---|
178 | attribute_dic, title_index_dic = csv2dict('sqrt_table.csv') |
---|
179 | #print attribute_dic |
---|
180 | print attribute_dic['number'][10] |
---|
181 | |
---|
182 | def test(x,y): |
---|
183 | # center_x = 12.96 |
---|
184 | # center_y = 13.80 |
---|
185 | # center_y = 6.80 |
---|
186 | |
---|
187 | du = (center_x-x)**2+(center_y-L-y)**2 |
---|
188 | # print 'x',min(x),max(x) |
---|
189 | # print 'y',min(y),max(y) |
---|
190 | z1=[] |
---|
191 | for d in du: |
---|
192 | # this finds the int that relates to the correct sqrt in the table |
---|
193 | nn=int(round(d,2)*100) |
---|
194 | zz=-float(attribute_dic['sqrt'][nn]) |
---|
195 | # print nn,zz |
---|
196 | |
---|
197 | # z = -square_root((center_x-x)**2+(center_y-y)**2) |
---|
198 | z = zz*.25+0.8975 |
---|
199 | |
---|
200 | if z <0: |
---|
201 | z=0 |
---|
202 | if z > .625: |
---|
203 | z=0.625 |
---|
204 | # print 'hello',z |
---|
205 | z1.append(z) |
---|
206 | return z1 |
---|
207 | |
---|
208 | domain.set_quantity('elevation',test, verbose=True) |
---|
209 | |
---|
210 | |
---|
211 | #------------------------- |
---|
212 | # Set simulation parameters |
---|
213 | #------------------------- |
---|
214 | model_name='wavetank_model' |
---|
215 | domain.set_name(model_name) # Name of output sww file |
---|
216 | domain.set_default_order(1) # Apply second order scheme |
---|
217 | #domain.set_all_limiters(0.9) # Max second order scheme (old lim) |
---|
218 | domain.set_minimum_storable_height(0.001) # Don't store w < 0.001m |
---|
219 | domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK) |
---|
220 | domain.tight_slope_limiters = 1 |
---|
221 | domain.beta_h = 0.0 |
---|
222 | |
---|
223 | |
---|
224 | #------------------------- |
---|
225 | # Boundary Conditions |
---|
226 | #------------------------- |
---|
227 | |
---|
228 | # Create boundary function from timeseries provided in file |
---|
229 | #function = file_function(project.boundary_filename, |
---|
230 | # domain, verbose=True) |
---|
231 | |
---|
232 | # Create and assign boundary objects |
---|
233 | boundary_filename="ts2cnew1_input" |
---|
234 | prepare_timeboundary(boundary_filename+'.txt') |
---|
235 | |
---|
236 | function = file_function(boundary_filename+'.tms', |
---|
237 | domain, verbose=True) |
---|
238 | |
---|
239 | # Create and assign boundary objects |
---|
240 | Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) |
---|
241 | |
---|
242 | Bw = Time_boundary(domain=domain, |
---|
243 | f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0]) |
---|
244 | Br = Reflective_boundary(domain) |
---|
245 | Bd = Dirichlet_boundary([water_height,0,0]) |
---|
246 | domain.set_boundary({'wave': Bts, 'wall': Br, 'buffer':Bd}) |
---|
247 | domain.starttime = 10 |
---|
248 | #------------------------- |
---|
249 | # Evolve through time |
---|
250 | #------------------------- |
---|
251 | import time |
---|
252 | t0 = time.time() |
---|
253 | |
---|
254 | for t in domain.evolve(yieldstep = 1, finaltime = 30): |
---|
255 | domain.write_time() |
---|
256 | # domain.write_time(track_speeds=False) |
---|
257 | |
---|
258 | for t in domain.evolve(yieldstep = 0.1, finaltime = 36 |
---|
259 | #for t in domain.evolve(yieldstep = 0.5, finaltime = 36 |
---|
260 | ,skip_initial_step = True): |
---|
261 | domain.write_time() |
---|
262 | |
---|
263 | for t in domain.evolve(yieldstep = 1, finaltime = 45 |
---|
264 | ,skip_initial_step = True): |
---|
265 | domain.write_time() |
---|
266 | |
---|
267 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
268 | |
---|
269 | |
---|
270 | |
---|
271 | #define the segments to find the run-up height to compare to wave tank study |
---|
272 | angles = (0.0, 22.5, 45.0, 67.5, 75.0, 80.0, 85.0, 87.5, 90.0, 92.5, |
---|
273 | 95.0, 100.0, 105.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0, |
---|
274 | 247.5, 270.0, 292.5, 315.0, 337.5) |
---|
275 | #x = 12.96 |
---|
276 | #y = 13.80 |
---|
277 | r1 = 1.1 |
---|
278 | r2 = 3.6 |
---|
279 | d=2 |
---|
280 | |
---|
281 | poly_segment = [] |
---|
282 | |
---|
283 | for i,angle in enumerate(angles): |
---|
284 | angle = ((angle-180)/180.0)*pi |
---|
285 | # angle = angle1*3.14157 |
---|
286 | p1 = get_xy(center_x,center_y,r1,angle-0.01745*d) |
---|
287 | p2 = get_xy(center_x,center_y,r2,angle-0.01745*d) |
---|
288 | p3 = get_xy(center_x,center_y,r2,angle+0.01745*d) |
---|
289 | p4 = get_xy(center_x,center_y,r1,angle+0.01745*d) |
---|
290 | # print p1,p2,p3,p4,angle |
---|
291 | poly_segment =[[p1[0],p1[1]],[p2[0],p2[1]],[p3[0],p3[1]],[p4[0],p4[1]]] |
---|
292 | # print i, poly_segment, angle |
---|
293 | #print i,poly[i] |
---|
294 | |
---|
295 | run_up, x_y = get_maximum_inundation_data(filename=model_name+'.sww',polygon=poly_segment, verbose=False) |
---|
296 | |
---|
297 | print 'maxd %.3f, %.6f, %.6f, %.5s, %.5s '%(((angle/pi)*180)+180, run_up*100, run_up-water_height*100, x_y[0],x_y[1]) |
---|
298 | |
---|
299 | |
---|
300 | |
---|
301 | |
---|
302 | |
---|
303 | |
---|
304 | |
---|