source: anuga_validation/circular_island/run_circular.py @ 5123

Last change on this file since 5123 was 5123, checked in by nick, 16 years ago

updates to circular island model. Now uses Numeric sqrt instead of a sqrt table.

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