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

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

Current Hinwood scenario

File size: 11.3 KB
Line 
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
24from os.path import dirname, join  #, basename
25from Numeric import zeros, size, Float
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, \
34      file_function
35from anuga.shallow_water.data_manager import copy_code_files
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
42from prepare_time_boundary import prepare_time_boundary
43from interp import interp
44
45
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
57
58def main(boundary_file,
59         metadata_dic,
60         boundary_path=None,
61         friction=0.012,  # planed wood. http://www.lmnoeng.com/manningn.htm
62         outputdir_name=None,
63         run_type=0,
64         width=1.0,
65         use_limits=True,
66         end_tag = '_limiterD'):
67
68   
69    basename = 'zz_' + metadata_dic['scenario_id']
70   
71    if run_type == 1:
72        yieldstep = 1.0
73        finaltime = 15.
74        maximum_triangle_area=0.1
75        outputdir_name += '_test'
76       
77    elif run_type == 2:
78        yieldstep = 0.5
79        finaltime = None
80        maximum_triangle_area=0.01
81        outputdir_name += '_test_long_time'
82       
83    elif run_type == 3:
84        yieldstep = 0.1
85        finaltime = None       
86        maximum_triangle_area=0.01
87        #outputdir_name += '_yieldstep_0.1'
88       
89    elif run_type == 4:
90        # this is not a test
91        # Output will go to a file
92        # The sww file will be interpolated
93        yieldstep = 0.01
94        finaltime = None       
95        maximum_triangle_area=0.01
96        #outputdir_name += '_no_velocity'
97       
98    elif run_type == 5:
99        # this is not a test
100        # Output will go to a file
101        # The sww file will be interpolated
102        yieldstep = 0.01
103        finaltime = None       
104        maximum_triangle_area=0.001
105        #outputdir_name += '_good'
106       
107    elif run_type == 6:
108        # this is not a test
109        # Output will go to a file
110        # The sww file will be interpolated
111        yieldstep = 0.01
112        finaltime = None       
113        maximum_triangle_area=0.0001
114        #outputdir_name += '_good'
115       
116    elif run_type == 7:
117        # this is not a test
118        # Output will go to a file
119        # The sww file will be interpolated
120        yieldstep = 0.01
121        finaltime = None       
122        maximum_triangle_area=0.00001
123        #outputdir_name += '_good'
124       
125    elif run_type == 8:
126        # this is not a test
127        # Output will go to a file
128        # The sww file will be interpolated
129        yieldstep = 0.05
130        finaltime = None       
131        maximum_triangle_area=0.00001
132        #outputdir_name += '_good'
133       
134    if use_limits is True:
135        outputdir_name += '_lmts'
136    else:
137        outputdir_name += '_nolmts'
138    outputdir_name += '_wdth_' + str(width)
139    outputdir_name += '_z_' + str(friction)
140    outputdir_name += '_ys_' + str(yieldstep)
141    outputdir_name += '_mta_' + str(maximum_triangle_area)
142    outputdir_name += end_tag
143   
144    metadata_dic = set_z_origin_to_water_depth(metadata_dic)   
145       
146    pro_instance = project.Project(['data','flumes','Hinwood_2008'],
147                                   outputdir_name=outputdir_name)
148    print "The output dir is", pro_instance.outputdir
149    copy_code_files(pro_instance.outputdir,__file__,
150                    dirname(project.__file__) \
151                    + sep + project.__name__+'.py')
152    copy (pro_instance.codedir + 'run_dam.py',
153          pro_instance.outputdir + 'run_dam.py')
154    copy (pro_instance.codedir + 'create_mesh.py',
155          pro_instance.outputdir + 'create_mesh.py')
156
157    boundary_final_time = prepare_time_boundary(metadata_dic,
158                                       pro_instance.raw_data_dir,
159                                       pro_instance.boundarydir)
160    #return pro_instance
161    if finaltime is None:
162        finaltime = boundary_final_time - 0.1 # Edge boundary problems
163    # Boundary file manipulation
164    if boundary_path is None:
165        boundary_path = pro_instance.boundarydir
166    boundary_file_path = join(boundary_path, boundary_file)
167   #  # Convert the boundary file, .csv to .tsm
168#     try:
169#         temp = open(boundary_file_path)
170#         temp.close()
171#     except IOError:
172#         prepare_time_boundary(boundary_file_path)
173   
174    mesh_filename = pro_instance.meshdir + basename + '.msh'
175
176    #--------------------------------------------------------------------------
177    # Copy scripts to output directory and capture screen
178    # output to file
179    #--------------------------------------------------------------------------
180
181    # creates copy of code in output dir
182    if run_type >= 2:
183        #start_screen_catcher(pro_instance.outputdir, rank, pypar.size())
184        start_screen_catcher(pro_instance.outputdir)
185
186    print 'USER:    ', pro_instance.user
187    #-------------------------------------------------------------------------
188    # Create the triangular mesh
189    #-------------------------------------------------------------------------
190
191    # this creates the mesh
192    #gate_position = 12.0
193    create_mesh.generate(mesh_filename, metadata_dic, width=width,
194                         maximum_triangle_area=maximum_triangle_area)
195
196    head,tail = path.split(mesh_filename)
197    copy (mesh_filename,
198          pro_instance.outputdir + tail )
199    #-------------------------------------------------------------------------
200    # Setup computational domain
201    #-------------------------------------------------------------------------
202    domain = Domain(mesh_filename, use_cache = False, verbose = True)
203   
204
205    print 'Number of triangles = ', len(domain)
206    print 'The extent is ', domain.get_extent()
207    print domain.statistics()
208
209   
210    domain.set_name(basename)
211    domain.set_datadir(pro_instance.outputdir)
212    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
213    domain.set_minimum_storable_height(0.0001)
214
215    if use_limits is True:
216        domain.set_default_order(2) # Use second order spatial scheme
217        domain.set_timestepping_method('rk2')
218        domain.use_edge_limiter = True
219        domain.tight_slope_limiters = True
220       
221        domain.beta_w      = 0.6
222        domain.beta_uh     = 0.6
223        domain.beta_vh     = 0.6
224   
225
226    #-------------------------------------------------------------------------
227    # Setup initial conditions
228    #-------------------------------------------------------------------------
229
230    domain.set_quantity('stage', 0.) #the origin is the still water level
231    domain.set_quantity('friction', friction)
232    elevation_function = Elevation_function(metadata_dic)
233    domain.set_quantity('elevation', elevation_function)
234
235   
236    print 'Available boundary tags', domain.get_boundary_tags()
237
238    # Create boundary function from timeseries provided in file
239    #function = file_function(project.boundary_file, domain, verbose=True)
240    #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
241    try:
242        function = file_function(boundary_file_path, domain,
243                                 verbose=True)
244    except IOError:
245        msg = 'Run prepare_time_boundary.py. File "%s" could not be opened.'\
246                  %(pro_instance.boundary_file)
247        raise msg
248       
249    Br = Reflective_boundary(domain)
250    Bd = Dirichlet_boundary([0.3,0,0]) 
251    Bts = Time_boundary(domain, function)
252    #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
253    domain.set_boundary( {'wall': Br, 'wave': Bts} )
254    #domain.set_boundary( {'wall': Br, 'wave': Bd} )
255
256    #-------------------------------------------------------------------------
257    # Evolve system through time
258    #-------------------------------------------------------------------------
259    t0 = time.time()
260
261    # It seems that ANUGA can't handle a starttime that is >0.
262    #domain.starttime = 1.0 #!!! what was this doing?
263    for t in domain.evolve(yieldstep, finaltime):   
264        domain.write_time()
265    print 'That took %.2f seconds' %(time.time()-t0)
266    print 'finished'
267
268    flume_y_middle = 0.0
269    points = []
270    for gauge_x in metadata_dic['gauge_x']:
271        points.append([gauge_x, flume_y_middle])
272    print "points",points
273
274
275    #-------------------------------------------------------------------------
276    # Calculate gauge info
277    #-------------------------------------------------------------------------
278
279    if run_type >= 1:
280        id = metadata_dic['scenario_id'] + ".csv"
281        interpolate_sww2csv(pro_instance.outputdir + basename +".sww",
282                            points,
283                            pro_instance.outputdir + "depth_" + id,
284                            pro_instance.outputdir + "velocity_x_" + id,
285                            pro_instance.outputdir + "velocity_y_" + id,
286                            pro_instance.outputdir + "stage_" + id,
287                            pro_instance.outputdir + "froude_" + id)
288 
289    return pro_instance
290
291def set_z_origin_to_water_depth(seabed_coords):
292    offset = seabed_coords['offshore_water_depth']
293    keys = ['xleft', 'xtoe', 'xbeach', 'xright']
294    for x in keys:
295            seabed_coords[x][1] -= offset
296    return seabed_coords
297#-------------------------------------------------------------
298if __name__ == "__main__":
299   
300    from scenarios import scenarios
301    from slope import gauges_for_slope
302    #from plot import plot
303
304    # 1 is fast and dirty
305    # 4 is 0.01
306    # 5 is 0.001
307    # 6 is 0.0001
308    # 7 is 0.00001
309    # 8 is 0.00001  yieldstep = 0.05
310   
311    #run_type = 1
312    run_type = 4
313    #for run_data in [scenarios[5]]:
314    #scenarios = scenarios[2:]
315    #scenarios = [scenarios[0]]
316    width = 1.0
317    width = 0.1
318    #width = 0.01
319    for run_data in scenarios:
320        pro_instance = main( run_data['scenario_id'] + '_boundary.tsm'  ,
321                             run_data,
322                             width=width,
323                             run_type=run_type,
324                             outputdir_name=run_data['scenario_id'],
325                             use_limits=False,
326                             friction=0.012,
327                             end_tag='_G')
328        #gauges_for_slope(pro_instance.outputdir,[run_data])
Note: See TracBrowser for help on using the repository browser.