source: anuga_work/development/Hinwood_2008/run_dam.py @ 5410

Last change on this file since 5410 was 5410, checked in by duncan, 16 years ago

Current Hinwood scenario

File size: 17.7 KB
RevLine 
[5076]1"""
2
3Script for running a breaking wave simulation of Jon Hinwoods wave tank.
4Note: this is based on the frinction_ua_flume_2006 structure.
5
6
7Duncan Gray, GA - 2007
8
9
10
11"""
12
13
14#----------------------------------------------------------------------------
15# Import necessary modules
16#----------------------------------------------------------------------------
17
18# Standard modules
19import time
20from time import localtime, strftime
21import sys
22from shutil import copy
23from os import path, sep
[5370]24from os.path import dirname, join  #, basename
25from Numeric import zeros, size, Float
[5076]26
27# Related major packages
28from anuga.shallow_water import Domain, Reflective_boundary, \
29                            Dirichlet_boundary,  Time_boundary, \
30                            File_boundary, \
31                            Transmissive_Momentum_Set_Stage_boundary
32from anuga.fit_interpolate.interpolate import interpolate_sww2csv
33from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \
[5392]34      file_function
35from anuga.shallow_water.data_manager import copy_code_files
[5076]36from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
37     import File_boundary_time
38
39# Scenario specific imports
40import project                 # Definition of file names and polygons
41import create_mesh
[5350]42from prepare_time_boundary import prepare_time_boundary
[5370]43from interp import interp
[5076]44
[5350]45
[5370]46class Elevation_function:
47    def __init__(self, slope):
48        self.xslope_position = [slope['xleft'][0],slope['xtoe'][0],
49                  slope['xbeach'][0],slope['xright'][0]]
50        self.yslope_height = [slope['xleft'][1],slope['xtoe'][1],
51                  slope['xbeach'][1],slope['xright'][1]]
52       
53    def __call__(self, x,y):
54       
55        z = interp(self.yslope_height, self.xslope_position, x)
56        return z
[5350]57
[5370]58def main(boundary_file,
[5395]59         metadata_dic,
[5370]60         boundary_path=None,
61         friction=0.01,
62         outputdir_name=None,
[5392]63         run_type=0):
[5092]64
[5350]65   
[5395]66    basename = 'zz_' + metadata_dic['scenario_id']
[5392]67    if run_type == 1:
[5076]68        outputdir_name += '_test'
69        yieldstep = 0.1
70        finaltime = 15.
71        maximum_triangle_area=0.01
[5395]72       
[5392]73    elif run_type == 2:
[5395]74        outputdir_name += '_test_long_time'
75        yieldstep = 0.5
76        finaltime = None
[5392]77        maximum_triangle_area=0.01
[5076]78       
[5392]79    elif run_type == 3:
80        outputdir_name += '_test_good_time_mesh'
81        yieldstep = 0.1
[5395]82        finaltime = None       
[5392]83        maximum_triangle_area=0.001
84    elif run_type == 4:
[5410]85        outputdir_name += '_good_tri_area_0.01_A'
[5392]86        # this is not a test
87        # Output will go to a file
88        # The sww file will be interpolated
[5395]89        yieldstep = 0.01
90        finaltime = None       
[5410]91        maximum_triangle_area=0.01
[5405]92    elif run_type == 5:
[5410]93        outputdir_name += '_good_tri_area_0.001_A'
[5405]94        # this is not a test
95        # Output will go to a file
96        # The sww file will be interpolated
97        yieldstep = 0.01
98        finaltime = None       
[5410]99        maximum_triangle_area=0.001
[5405]100    elif run_type == 6:
[5410]101        outputdir_name += '_good_tri_area_0.0001_A'
[5405]102        # this is not a test
103        # Output will go to a file
104        # The sww file will be interpolated
105        yieldstep = 0.01
106        finaltime = None       
[5410]107        maximum_triangle_area=0.0001
[5395]108     
109    metadata_dic = set_z_origin_to_water_depth(metadata_dic)   
[5076]110       
111    pro_instance = project.Project(['data','flumes','Hinwood_2008'],
[5350]112                                   outputdir_name=outputdir_name)
[5076]113    print "The output dir is", pro_instance.outputdir
114    copy_code_files(pro_instance.outputdir,__file__,
115                    dirname(project.__file__) \
116                    + sep + project.__name__+'.py')
117    copy (pro_instance.codedir + 'run_dam.py',
118          pro_instance.outputdir + 'run_dam.py')
119    copy (pro_instance.codedir + 'create_mesh.py',
120          pro_instance.outputdir + 'create_mesh.py')
[5350]121
[5395]122    boundary_final_time = prepare_time_boundary(metadata_dic,
123                                       pro_instance.raw_data_dir,
124                                       pro_instance.boundarydir)
125    if finaltime is None:
126        finaltime = boundary_final_time
[5370]127    # Boundary file manipulation
128    if boundary_path is None:
129        boundary_path = pro_instance.boundarydir
130    boundary_file_path = join(boundary_path, boundary_file)
[5392]131   #  # Convert the boundary file, .csv to .tsm
132#     try:
133#         temp = open(boundary_file_path)
134#         temp.close()
135#     except IOError:
136#         prepare_time_boundary(boundary_file_path)
[5350]137   
[5076]138    mesh_filename = pro_instance.meshdir + basename + '.msh'
139
140    #--------------------------------------------------------------------------
141    # Copy scripts to output directory and capture screen
142    # output to file
143    #--------------------------------------------------------------------------
144
145    # creates copy of code in output dir
[5395]146    if run_type >= 2:
147        #start_screen_catcher(pro_instance.outputdir, rank, pypar.size())
148        start_screen_catcher(pro_instance.outputdir)
[5076]149
150    print 'USER:    ', pro_instance.user
151    #-------------------------------------------------------------------------
152    # Create the triangular mesh
153    #-------------------------------------------------------------------------
154
155    # this creates the mesh
156    #gate_position = 12.0
[5395]157    create_mesh.generate(mesh_filename, metadata_dic,
[5076]158                         maximum_triangle_area=maximum_triangle_area)
159
160    head,tail = path.split(mesh_filename)
161    copy (mesh_filename,
162          pro_instance.outputdir + tail )
163    #-------------------------------------------------------------------------
164    # Setup computational domain
165    #-------------------------------------------------------------------------
166    domain = Domain(mesh_filename, use_cache = False, verbose = True)
167   
168
169    print 'Number of triangles = ', len(domain)
170    print 'The extent is ', domain.get_extent()
171    print domain.statistics()
172
173   
174    domain.set_name(basename)
175    domain.set_datadir(pro_instance.outputdir)
176    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
177    domain.set_minimum_storable_height(0.001)
178    #domain.set_store_vertices_uniquely(True)  # for writting to sww
179
180    #-------------------------------------------------------------------------
181    # Setup initial conditions
182    #-------------------------------------------------------------------------
183
[5395]184    domain.set_quantity('stage', 0.) #the origin is the still water level
[5370]185    domain.set_quantity('friction', friction)
[5395]186    elevation_function = Elevation_function(metadata_dic)
[5076]187    domain.set_quantity('elevation', elevation_function)
188
189   
190    print 'Available boundary tags', domain.get_boundary_tags()
191
192    # Create boundary function from timeseries provided in file
193    #function = file_function(project.boundary_file, domain, verbose=True)
194    #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
[5350]195    try:
[5370]196        function = file_function(boundary_file_path, domain,
[5350]197                                 verbose=True)
198    except IOError:
199        msg = 'Run prepare_time_boundary.py. File "%s" could not be opened.'\
200                  %(pro_instance.boundary_file)
201        raise msg
202       
[5076]203    Br = Reflective_boundary(domain)
[5370]204    Bd = Dirichlet_boundary([0.3,0,0]) 
[5076]205    Bts = Time_boundary(domain, function)
206    domain.set_boundary( {'wall': Br, 'wave': Bts} )
[5395]207    #domain.set_boundary( {'wall': Br, 'wave': Bd} )
[5076]208
209    #-------------------------------------------------------------------------
210    # Evolve system through time
211    #-------------------------------------------------------------------------
212    t0 = time.time()
213
[5370]214    # It seems that ANUGA can't handle a starttime that is >0.
215    domain.starttime = 1.0
[5078]216    for t in domain.evolve(yieldstep, finaltime):   
[5076]217        domain.write_time()
[5078]218    print 'That took %.2f seconds' %(time.time()-t0)
219    print 'finished'
[5076]220
[5405]221    flume_y_middle = 0.5
[5395]222    points = []
223    for gauge_x in metadata_dic['gauge_x']:
224        points.append([gauge_x, flume_y_middle])
225    print "points",points
[5076]226
227
228    #-------------------------------------------------------------------------
229    # Calculate gauge info
230    #-------------------------------------------------------------------------
231
[5395]232    if run_type >= 2:
233        id = metadata_dic['scenario_id']
[5076]234        interpolate_sww2csv(pro_instance.outputdir + basename +".sww",
235                            points,
[5395]236                            pro_instance.outputdir + "depth_" + id + ".csv",
237                            pro_instance.outputdir + "velocity_x_" + id + ".csv",
238                            pro_instance.outputdir + "velocity_y_" + id + ".csv")
[5076]239 
240    return pro_instance
241
[5390]242def set_z_origin_to_water_depth(seabed_coords):
243    offset = seabed_coords['offshore_water_depth']
[5392]244    keys = ['xleft', 'xtoe', 'xbeach', 'xright']
245    for x in keys:
[5390]246            seabed_coords[x][1] -= offset
247    return seabed_coords
[5076]248#-------------------------------------------------------------
249if __name__ == "__main__":
[5370]250    #slopes = [[-4.5,0.0],[0.0,0.0],[1.285,0.090],[16.1,.960]]
[5395]251    # Note, gauge A has been removed, since it is used as the
252    # boundary.
[5405]253    run_type = 6
[5395]254
255    # T1R5
256    run_data = {'xleft':[-3.106,0.0],  # Av' of ADV and Gauge A
[5390]257                     'xtoe':[0.0,0.0],
258                     'xbeach':[1.285,0.090],
259                     'xright':[16.1,.960],
[5392]260                     'offshore_water_depth':.4,
[5395]261                     'scenario_id':'T1R5',
262                     'gauge_names':['B','1','2','3','4','5','6','7','8',
263                                    '9','10','11','12','13','14'],
264                     'gauge_x':[-0.68, 1.572, 2.572, 3.572, 4.572, 5.572,
265                                6.572, 7.572, 8.572, 9.572, 10.572, 11.572,
266                                12.572, 13.572, 14.572],
267                     'gauge_bed_elevation':[-0.400000, -0.293158,
268                     -0.234473, -0.175788, -0.117104, -0.058419, 0.000266,
269                     0.058950, 0.117635, 0.176320, 0.235004, 0.293689,
270                     0.352374, 0.411058, 0.469743]
271                     }
272    main( run_data['scenario_id'] + '_boundary.tsm' , run_data,
273          run_type = run_type,
274         outputdir_name=run_data['scenario_id'])
275
276    # T1R3
[5410]277    run_data = {'xleft':[-3.106,0.0],  # Av' of ADV and Gauge A
278                     'xtoe':[0.0,0.0],
279                     'xbeach':[1.285,0.090],
280                     'xright':[16.1,.960],
281                     'offshore_water_depth':.4,
282                     'scenario_id':'T1R3',
283                     'gauge_names':['B','1','2','3','4','5','6','7','8',
284                                    '9','10','11','12','13','14'],
285                     'gauge_x':[-0.68, 1.572, 2.572, 3.572, 4.572, 5.572,
286                                6.572, 7.572, 8.572, 9.572, 10.572, 11.572,
287                                12.572, 13.572, 14.572],
288                     'gauge_bed_elevation':[-0.400000, -0.293158,
289                     -0.234473, -0.175788, -0.117104, -0.058419, 0.000266,
290                     0.058950, 0.117635, 0.176320, 0.235004, 0.293689,
291                     0.352374, 0.411058, 0.469743]
292                     }
[5395]293    main( run_data['scenario_id'] + '_boundary.tsm'  , run_data,
294          run_type = run_type,
295         outputdir_name=run_data['scenario_id'])
296
297    # #T2R7
298    # xleft is different
299    run_data = {'xleft':[-4.586,0.0],  # Av' of ADV and Gauge A
300                     'xtoe':[0.0,0.0],
301                     'xbeach':[1.285,0.090],
302                     'xright':[16.1,.960],
303                     'offshore_water_depth':.4,
304                     'scenario_id':'T2R7',
305                     'gauge_names':['B','1','2','3','4','5','6','7','8',
306                                    '9','10','11','12','13','14'],
307                     'gauge_x':[-0.68, 1.572, 2.572, 3.572, 4.572, 5.572,
308                                6.572, 7.572, 8.572, 9.572, 10.572, 11.572,
309                                12.572, 13.572, 14.572],
310                     'gauge_bed_elevation':[-0.400000, -0.293158,
311                     -0.234473, -0.175788, -0.117104, -0.058419, 0.000266,
312                     0.058950, 0.117635, 0.176320, 0.235004, 0.293689,
313                     0.352374, 0.411058, 0.469743]
314                     }
315    main( run_data['scenario_id'] + '_boundary.tsm'  , run_data,
316          run_type = run_type,
317          outputdir_name=run_data['scenario_id'])
318
319    # #T2R8
320    # xleft is different
321    run_data = {'xleft':[-4.586,0.0],  # Av' of ADV and Gauge A
322                     'xtoe':[0.0,0.0],
323                     'xbeach':[1.285,0.090],
324                     'xright':[16.1,.960],
325                     'offshore_water_depth':.4,
326                     'scenario_id':'T2R8',
327                     'gauge_names':['B','1','2','3','4','5','6','7','8',
328                                    '9','10','11','12','13','14'],
329                     'gauge_x':[-0.68, 1.572, 2.572, 3.572, 4.572, 5.572,
330                                6.572, 7.572, 8.572, 9.572, 10.572, 11.572,
331                                12.572, 13.572, 14.572],
332                     'gauge_bed_elevation':[-0.400000, -0.293158,
333                     -0.234473, -0.175788, -0.117104, -0.058419, 0.000266,
334                     0.058950, 0.117635, 0.176320, 0.235004, 0.293689,
335                     0.352374, 0.411058, 0.469743]
336                     }
337    main( run_data['scenario_id'] + '_boundary.tsm'  , run_data,
338          run_type = run_type,
339          outputdir_name=run_data['scenario_id'])
340
341    # #T3R29
342    # xleft is different
343    run_data = {'xleft':[-3.875,0.0],  # Av' of ADV and Gauge A
344                'xtoe':[0.0,0.0],
345                'xbeach':[1.285,0.090],
346                'xright':[16.1,.440],
347                'offshore_water_depth':.336,
348                'scenario_id':'T3R29',
349                'gauge_names':['1','2','3','4','5','6','7','8',
350                               '9','10','11','12','13','14','B'],
351                'gauge_x':[1.572, 2.572, 3.572, 4.572, 5.572,
352                           6.572, 7.572, 8.572, 9.572, 10.572, 11.572,
353                           12.572, 13.572, 14.572, -0.325],
354                'gauge_bed_elevation':[-0.237263, -0.213789, -0.190315,
355                                       -0.166841, -0.143368, -0.119894,
356                                       -0.096420, -0.072946, -0.049472,
357                                       -0.025998, -0.002524, 0.020949,
358                                       0.044423, 0.067897, -0.336000]
359                     }
360    main( run_data['scenario_id'] + '_boundary.tsm'  , run_data,
361          run_type = run_type,
362          outputdir_name=run_data['scenario_id'])
363
364    # #T3R28
365    run_data = {'xleft':[-3.875,0.0],  # Av' of ADV and Gauge A
366                'xtoe':[0.0,0.0],
367                'xbeach':[1.285,0.090],
368                'xright':[16.1,.440],
369                'offshore_water_depth':.336,
370                'scenario_id':'T3R28',
371                'gauge_names':['1','2','3','4','5','6','7','8',
372                               '9','10','11','12','13','14','B'],
373                'gauge_x':[1.572, 2.572, 3.572, 4.572, 5.572,
374                           6.572, 7.572, 8.572, 9.572, 10.572, 11.572,
375                           12.572, 13.572, 14.572, -0.325],
376                'gauge_bed_elevation':[-0.237263, -0.213789, -0.190315,
377                                       -0.166841, -0.143368, -0.119894,
378                                       -0.096420, -0.072946, -0.049472,
379                                       -0.025998, -0.002524, 0.020949,
380                                       0.044423, 0.067897, -0.336000]
381                     }
382    main( run_data['scenario_id'] + '_boundary.tsm'  , run_data,
383          run_type = run_type,
384          outputdir_name=run_data['scenario_id'])
385
386    # #T4R31
387    # xleft is different
388    run_data = {'xleft':[-2.43,0.0],  # Av' of ADV and Gauge A
389                'xtoe':[0.0,0.0],
390                'xbeach':[1.285,0.090],
391                'xright':[16.1,.440],
392                'offshore_water_depth':.336,
393                'scenario_id':'T4R31',
394                'gauge_names':['1','2','3','4','5','6','7','8',
395                               '9','10','11','12','13','14','B'],
396                'gauge_x':[1.572, 2.572, 3.572, 4.572, 5.572,
397                           6.572, 7.572, 8.572, 9.572, 10.572, 11.572,
398                           12.572, 13.572, 14.572, -0.325],
399                'gauge_bed_elevation':[-0.237263, -0.213789, -0.190315,
400                                       -0.166841, -0.143368, -0.119894,
401                                       -0.096420, -0.072946, -0.049472,
402                                       -0.025998, -0.002524, 0.020949,
403                                       0.044423, 0.067897, -0.336000]
404                     }
405    main( run_data['scenario_id'] + '_boundary.tsm'  , run_data,
406          run_type = run_type,
407          outputdir_name=run_data['scenario_id'])
408   
409    # #T4R32
410    run_data = {'xleft':[-2.43,0.0],  # Av' of ADV and Gauge A
411                'xtoe':[0.0,0.0],
412                'xbeach':[1.285,0.090],
413                'xright':[16.1,.440],
414                'offshore_water_depth':.336,
415                'scenario_id':'T4R32',
416                'gauge_names':['1','2','3','4','5','6','7','8',
417                               '9','10','11','12','13','14','B'],
418                'gauge_x':[1.572, 2.572, 3.572, 4.572, 5.572,
419                           6.572, 7.572, 8.572, 9.572, 10.572, 11.572,
420                           12.572, 13.572, 14.572, -0.325],
421                'gauge_bed_elevation':[-0.237263, -0.213789, -0.190315,
422                                       -0.166841, -0.143368, -0.119894,
423                                       -0.096420, -0.072946, -0.049472,
424                                       -0.025998, -0.002524, 0.020949,
425                                       0.044423, 0.067897, -0.336000]
426                     }
427    main( run_data['scenario_id'] + '_boundary.tsm'  , run_data,
428          run_type = run_type,
429          outputdir_name=run_data['scenario_id'])
Note: See TracBrowser for help on using the repository browser.