source: anuga_validation/circular_island/sqrt_table_run_circular.py @ 5390

Last change on this file since 5390 was 5131, checked in by nick, 17 years ago

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