source: anuga_validation/circular_island/sqrt_table_run_circular.py @ 7276

Last change on this file since 7276 was 7276, checked in by ole, 15 years ago

Merged numpy branch back into the trunk.

In ~/sandpit/anuga/anuga_core/source
svn merge -r 6246:HEAD ../../branches/numpy .

In ~/sandpit/anuga/anuga_validation
svn merge -r 6417:HEAD ../branches/numpy_anuga_validation .

In ~/sandpit/anuga/misc
svn merge -r 6809:HEAD ../branches/numpy_misc .

For all merges, I used numpy version where conflicts existed

The suites test_all.py (in source/anuga) and validate_all.py passed using Python2.5 with numpy on my Ubuntu Linux box.

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
19import numpy as num
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
24from anuga.config import netcdf_float
25from Scientific.IO.NetCDF import NetCDFFile
26
27
28#-------------------------
29# Create Domain from mesh
30#-------------------------
31
32user = get_user_name()
33home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent dir 
34time = strftime('%Y%m%d_%H%M%S',localtime())
35
36res=0.05
37
38dir_comment = time+'_res_'+str(res)+'_'+str(user)
39
40anuga_dir = home+'anuga_validation'+sep+'circular_island_tsunami_benchmark'+sep+'anuga'+sep
41output_dir = anuga_dir+'outputs'+sep+dir_comment+sep
42copy_code_files(output_dir,__file__)
43
44start_screen_catcher(output_dir)
45
46def get_xy(x=0,y=0,r=1,angle=1):
47    #return x and y co-ordinates
48    x1 = x + cos(angle)*r
49    y1 = y + sin(-angle)*r
50#    print "get",x,y,r,x1,y1,"angle",angle, cos(angle),sin(angle)
51    return x1,y1,
52
53def prepare_timeboundary(filename):
54    """Converting benchmark 2 time series to NetCDF tms file.
55    This is a 'throw-away' code taylor made for files like
56    'Benchmark_2_input.txt' from the LWRU2 benchmark
57    """
58
59    textversion = filename[:-4] + '.txt'
60
61    print 'Preparing time boundary from %s' %textversion
62
63    fid = open(textversion)
64
65    #Skip first line
66    line = fid.readline()
67
68    #Read remaining lines
69    lines = fid.readlines()
70    fid.close()
71
72
73    N = len(lines)
74    T = num.zeros(N, num.float)  #Time
75    Q = num.zeros(N, num.float)  #Values
76
77    for i, line in enumerate(lines):
78        fields = line.split()
79#        print 'fields',i,line
80
81        T[i] = float(fields[0])
82        Q[i] = float(fields[1])
83
84
85    #Create tms file
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', netcdf_float, ('number_of_timesteps',))
94    fid.variables['time'][:] = T
95
96    fid.createVariable('stage', netcdf_float, ('number_of_timesteps',))
97    fid.variables['stage'][:] = Q[:]
98
99    fid.createVariable('xmomentum', netcdf_float, ('number_of_timesteps',))
100    fid.variables['xmomentum'][:] = 0.0
101
102    fid.createVariable('ymomentum', netcdf_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#-------------------------
172domain.set_quantity('friction', 0.0)
173#water_height = 0.322
174water_height = 0.0
175domain.set_quantity('stage', water_height)
176
177
178attribute_dic, title_index_dic = csv2dict('sqrt_table.csv')
179#print attribute_dic
180#print attribute_dic['number'][10]
181
182def circular_island_elevation(x,y):
183
184    list_xy_square = (center_x-x)**2+(center_y-L-y)**2
185    print 'x',min(x),max(x)
186    print 'y',min(y),max(y)
187    list_z=[]
188    for xy in list_xy_square:
189        # this finds the int that relates to the correct sqrt in the table
190        nn=int(round(xy,2)*100)
191        #the square root is to find the distance from the center to the current xy
192        distance = float(attribute_dic['sqrt'][nn])
193#        zz=-float(attribute_dic['sqrt'][nn])
194#        print nn,zz
195       
196        # slope of cone is 1/4 and has a radii of 3.6 therefore the cone is 0.9 high.
197        # So to have the elevation positive 0.90 has been added. Furthermore the
198        #wavetank is submerged .32m so 0.9 - 0.32=.0578
199#        z = zz*.25+0.578
200
201#        z = -distance*.25+0.9
202#        z = -distance*.25+0.578
203        z = -distance*.25+0.9-0.322
204#        z = zz*.25+0.8975
205
206        if z <0-0.322:
207            z=0-0.322
208        if z > .625-0.322:
209            z=0.625-0.322
210#        print 'hello',z
211        list_z.append(z)
212    return list_z
213
214domain.set_quantity('elevation', circular_island_elevation, verbose=True)
215
216
217#-------------------------
218# Set simulation parameters
219#-------------------------
220model_name='wavetank_model'
221domain.set_name(model_name)  # Name of output sww file
222domain.set_datadir(output_dir)
223domain.set_default_order(1)               # Apply second order scheme
224#domain.set_all_limiters(0.9)              # Max second order scheme (old lim)
225domain.set_minimum_storable_height(0.001) # Don't store w < 0.001m
226domain.set_maximum_allowed_speed(0.1)     # Allow a little runoff (0.1 is OK)
227domain.tight_slope_limiters = 1
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.