source: anuga_validation/circular_island/sqrt_table_run_circular.py @ 7077

Last change on this file since 7077 was 5442, checked in by ole, 17 years ago

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.
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
227sww_file=output_dir+sep+model_name+'.sww'
228
229
230#-------------------------
231# Boundary Conditions
232#-------------------------
233
234# Create boundary function from timeseries provided in file
235
236boundary_filename="ts2cnew1_input_20_80sec_new"
237prepare_timeboundary(boundary_filename+'.txt')
238
239function = file_function(boundary_filename+'.tms',
240                         domain, verbose=True)
241
242# Create and assign boundary objects
243Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
244
245Br = Reflective_boundary(domain)
246Bd = Dirichlet_boundary([water_height,0,0])
247#domain.set_boundary({'wave': Bts, 'wall': Br, 'buffer':Bd})
248domain.set_boundary({'wave': Bts, 'wall': Bd, 'buffer':Bd})
249#domain.set_time(20.0)
250#-------------------------
251# Evolve through time
252#-------------------------
253import time
254t0 = time.time()
255
256for t in domain.evolve(yieldstep = 1, finaltime = 28):
257    domain.write_time()
258#    domain.write_time(track_speeds=False)
259
260for t in domain.evolve(yieldstep = 0.1, finaltime = 36
261#for t in domain.evolve(yieldstep = 0.5, finaltime = 36
262                       ,skip_initial_step = True): 
263    domain.write_time()
264   
265for t in domain.evolve(yieldstep = 1, finaltime = 45
266                       ,skip_initial_step = True): 
267    domain.write_time()
268
269print 'That took %.2f seconds' %(time.time()-t0)
270
271
272
273#define the segments to find the run-up height to compare to wave tank study
274angles = (0.0, 22.5, 45.0, 67.5, 75.0, 80.0, 85.0, 87.5, 90.0, 92.5,
275          95.0, 100.0, 105.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0,
276          247.5, 270.0, 292.5, 315.0, 337.5)
277r1 = 1.5
278r2 = 2.6
279#degrees to check
280d=1
281
282poly_segment = []
283
284for i,angle in enumerate(angles):
285    angle = ((angle-180)/180.0)*pi
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, location = get_maximum_inundation_data(filename=sww_file,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, location[0],location[1]) 
298
299#sww_file=anuga_dir+'outputs'+sep+'20080220_010628_res_0.05_nbartzis'+sep+'wavetank_model.sww'
300#sww_file=anuga_dir+'outputs'+sep+'20080219_233611_res_0.05_nbartzis'+sep+'wavetank_model.sww'
301sww2csv_gauges(sww_file=sww_file,
302                   gauge_file='gauges.csv',
303                   quantities = ['stage','depth', 'elevation', 'xmomentum', 'ymomentum'],
304                   verbose=True,
305                   use_cache = True)
306
307#csv2timeseries_graphs(directories_dic={anuga_dir+sep+'boundaries'+sep:['wavetank',0, 0],
308#                                       anuga_dir+'outputs'+sep+'20080219_233611_res_0.05_nbartzis'+sep:['Fixed Wave',0,0]},
309#                            output_dir=anuga_dir+'outputs'+sep+'20080219_233611_res_0.05_nbartzis'+sep,
310#                            base_name='gauge_',
311#                            plot_numbers='',
312#                            quantities=['stage'],
313#                            extra_plot_name='',
314#                            assess_all_csv_files=True,                           
315#                            create_latex=False,
316#                            verbose=True)
317
318
319
320
321
Note: See TracBrowser for help on using the repository browser.