source: anuga_validation/circular_island/run_circular.py @ 7746

Last change on this file since 7746 was 7746, checked in by hudson, 13 years ago

Upgraded validation to ANUGA 1.2

File size: 9.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.utilities.file_utils import copy_code_files
15from anuga.shallow_water.data_manager import get_maximum_inundation_data
16from anuga.abstract_2d_finite_volumes.util import file_function, \
17                            sww2csv_gauges, csv2timeseries_graphs
18from anuga.pmesh.mesh_interface import create_mesh_from_regions
19from anuga.utilities.polygon import read_polygon
20from math import cos,pi,sin,tan
21import numpy as num
22from anuga.shallow_water.data_manager import load_csv_as_dict as csv2dict
23from time import localtime, strftime, gmtime
24from anuga.utilities.system_tools import get_user_name, get_host_name
25from os import sep, environ, getenv
26from anuga.config import netcdf_float
27from Scientific.IO.NetCDF import NetCDFFile
28
29#-------------------------
30# Create Domain from mesh
31#-------------------------
32
33user = get_user_name()
34home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent dir 
35time = strftime('%Y%m%d_%H%M%S',localtime())
36
37res=0.05
38
39dir_comment = time+'_res_'+str(res)+'_'+str(user)
40
41anuga_dir = home+'anuga_validation'+sep+'circular_island_tsunami_benchmark'+sep+'anuga'+sep
42output_dir = anuga_dir+'outputs'+sep+dir_comment+sep
43#output_dir = anuga_dir+'outputs'+sep+'20080222_062917_res_0.0025_nbartzis'+sep
44copy_code_files(output_dir,__file__)
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#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.5,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,3.2,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#-------------------------
171#domain.set_quantity('friction', 0.0)
172domain.set_quantity('friction', 0.01)#for concrete
173#water_height = 0.322
174water_height = 0.0
175
176domain.set_quantity('stage', water_height)
177
178
179#attribute_dic, title_index_dic = csv2dict('sqrt_table.csv')
180
181def circular_island_elevation(x,y):
182    water_depth = 0.32
183    list_xy = num.sqrt((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    for distance in list_xy:
189
190#    The cone has a 1/4 slope and a 7.2 diameter this make the height 0.9m
191#    to get the stage at 0 elevation + cone height and - water depth
192        z = -distance*.25+0.9-water_depth
193
194        if z <0-water_depth:
195            z=0-water_depth
196        if z > .625-water_depth:
197            z=0.625-water_depth
198        list_z.append(z)
199    return list_z
200
201domain.set_quantity('elevation', circular_island_elevation, verbose=True)
202
203
204##-------------------------
205## Set simulation parameters
206##-------------------------
207model_name='wavetank_model'
208domain.set_name(model_name)  # Name of output sww file
209domain.set_datadir(output_dir)
210domain.set_default_order(2)               # Apply second order scheme
211#domain.set_all_limiters(0.9)              # Max second order scheme (old lim)
212domain.set_minimum_storable_height(0.001) # Don't store w < 0.001m
213#domain.set_maximum_allowed_speed(0.1)     # Allow a little runoff (0.1 is OK)
214domain.tight_slope_limiters = 1
215
216sww_file=output_dir+sep+model_name+'.sww'
217
218
219#-------------------------
220# Boundary Conditions
221#-------------------------
222
223# Create boundary function from timeseries provided in file
224
225boundary_filename="ts2cnew1_input_20_80sec_new"
226prepare_timeboundary(boundary_filename+'.txt')
227
228function = file_function(boundary_filename+'.tms',
229                         domain, verbose=True)
230
231# Create and assign boundary objects
232Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
233
234Br = Reflective_boundary(domain)
235Bd = Dirichlet_boundary([water_height,0,0])
236#domain.set_boundary({'wave': Bts, 'wall': Br, 'buffer':Bd})
237domain.set_boundary({'wave': Bts, 'wall': Bd, 'buffer':Bd})
238#domain.set_time(20.0)
239#-------------------------
240# Evolve through time
241#-------------------------
242import time
243t0 = time.time()
244
245for t in domain.evolve(yieldstep = 1, finaltime = 28):
246    domain.write_time()
247#    domain.write_time(track_speeds=False)
248
249for t in domain.evolve(yieldstep = 0.1, finaltime = 36
250#for t in domain.evolve(yieldstep = 0.5, finaltime = 36
251                       ,skip_initial_step = True): 
252    domain.write_time()
253   
254for t in domain.evolve(yieldstep = 1, finaltime = 45
255                       ,skip_initial_step = True): 
256    domain.write_time()
257
258print 'That took %.2f seconds' %(time.time()-t0)
259
260
261
262#define the segments to find the run-up height to compare to wave tank study
263angles = (0.0, 22.5, 45.0, 67.5, 75.0, 80.0, 85.0, 87.5, 90.0, 92.5,
264          95.0, 100.0, 105.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0,
265          247.5, 270.0, 292.5, 315.0, 337.5)
266r1 = 1.5
267r2 = 2.6
268#degrees to check
269d=1
270
271poly_segment = []
272
273for i,angle in enumerate(angles):
274    angle = ((angle-180)/180.0)*pi
275    p1 = get_xy(center_x,center_y,r1,angle-0.01745*d)
276    p2 = get_xy(center_x,center_y,r2,angle-0.01745*d)
277    p3 = get_xy(center_x,center_y,r2,angle+0.01745*d)
278    p4 = get_xy(center_x,center_y,r1,angle+0.01745*d)
279#    print p1,p2,p3,p4,angle
280    poly_segment =[[p1[0],p1[1]],[p2[0],p2[1]],[p3[0],p3[1]],[p4[0],p4[1]]]
281   
282    run_up, location = get_maximum_inundation_data(filename=sww_file,polygon=poly_segment, verbose=False)
283   
284    print 'maxd %.3f, %.6f, %.6f, %.5s, %.5s '%(((angle/pi)*180)+180, run_up*100, (run_up-water_height)*100, location[0],location[1]) 
285
286#sww_file=anuga_dir+'outputs'+sep+'20080220_010628_res_0.05_nbartzis'+sep+'wavetank_model.sww'
287#sww_file=anuga_dir+'outputs'+sep+'20080219_233611_res_0.05_nbartzis'+sep+'wavetank_model.sww'
288sww2csv_gauges(sww_file=sww_file,
289                   gauge_file=anuga_dir+'gauges'+sep+'gauges.csv',
290#                   gauge_file='gauges.csv',
291                   quantities = ['stage','depth', 'elevation', 'xmomentum', 'ymomentum'],
292                   verbose=True,
293                   use_cache = True)
294
295
296
297
298
Note: See TracBrowser for help on using the repository browser.