source: trunk/anuga_work/development/Hinwood_2008/run_dam.py @ 7885

Last change on this file since 7885 was 5714, checked in by ole, 16 years ago

Handover from Duncan in regard the Hinwood study and the validation paper.
Added comments about locations of data
File size: 12.4 KB
Line 
1"""
2
3Script for running a breaking wave simulation of Jon Hinwoods wave tank.
4
5
6Duncan Gray, GA - 2007
7
8
9
10"""
11
12
13#----------------------------------------------------------------------------
14# Import necessary modules
15#----------------------------------------------------------------------------
16
17# Standard modules
18import time
19from time import localtime, strftime
20import sys
21from shutil import copy
22from os import path, sep
23from os.path import dirname, join  #, basename
24from Numeric import zeros, size, Float
25
26# Related major packages
27from anuga.shallow_water import Domain, Reflective_boundary, \
28                            Dirichlet_boundary,  Time_boundary, \
29                            File_boundary, \
30                            Transmissive_Momentum_Set_Stage_boundary
31from anuga.fit_interpolate.interpolate import interpolate_sww2csv
32
33from anuga.abstract_2d_finite_volumes.util import file_function
34from anuga.shallow_water.data_manager import copy_code_files, \
35     start_screen_catcher
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 anuga.utilities.interp import interp
44
45
46class Elevation_function:
47    """
48   
49    Create a callable instance that returns the bed surface of the
50    flume given the scenario metadata.
51   
52    """
53    def __init__(self, slope):
54        self.xslope_position = [slope['xleft'][0],slope['xtoe'][0],
55                  slope['xbeach'][0],slope['xright'][0]]
56        self.yslope_height = [slope['xleft'][1],slope['xtoe'][1],
57                  slope['xbeach'][1],slope['xright'][1]]
58       
59    def __call__(self, x,y):
60       
61        z = interp(self.yslope_height, self.xslope_position, x)
62        return z
63
64def main(boundary_file,
65         metadata_dic,
66         maximum_triangle_area,
67         yieldstep,
68         boundary_path=None,
69         friction=0.012,  # planed wood. http://www.lmnoeng.com/manningn.htm
70         outputdir_name=None,
71         run_type=1,
72         width=1.0,
73         use_limits=True,
74         end_tag = '_limiterD'):
75
76   
77    basename = 'zz_' + metadata_dic['scenario_id']
78   
79    if run_type == 0:
80        yieldstep = 1.0
81        finaltime = 15.
82        maximum_triangle_area=0.1
83        outputdir_name += '_test'
84         
85    elif run_type == 1:
86        finaltime = None 
87       
88    elif run_type == 2:
89        yieldstep = 0.5
90        finaltime = None
91        maximum_triangle_area=0.01
92        outputdir_name += '_test_long_time'
93       
94    elif run_type == 3:
95        yieldstep = 0.1
96        finaltime = None       
97        maximum_triangle_area=0.01
98        #outputdir_name += '_yieldstep_0.1'
99       
100    elif run_type == 4:
101        # this is not a test
102        # Output will go to a file
103        # The sww file will be interpolated
104        yieldstep = 0.01
105        finaltime = None       
106        maximum_triangle_area=0.01
107       
108    elif run_type == 5:
109        # this is not a test
110        # Output will go to a file
111        # The sww file will be interpolated
112        yieldstep = 0.01
113        finaltime = None       
114        maximum_triangle_area=0.001
115        #outputdir_name += '_good'
116       
117    elif run_type == 6:
118        # this is not a test
119        # Output will go to a file
120        # The sww file will be interpolated
121        yieldstep = 0.01
122        finaltime = None       
123        maximum_triangle_area=0.0001
124        #outputdir_name += '_good'
125       
126    elif run_type == 7:
127        # this is not a test
128        # Output will go to a file
129        # The sww file will be interpolated
130        yieldstep = 0.01
131        finaltime = None       
132        maximum_triangle_area=0.00001
133        #outputdir_name += '_good'
134       
135    elif run_type == 8:
136        # this is not a test
137        # Output will go to a file
138        # The sww file will be interpolated
139        yieldstep = 0.05
140        finaltime = None       
141        maximum_triangle_area=0.00001
142        #outputdir_name += '_good'
143    #outputdir_name += '_no_velocity'
144
145    outputdir_name = paras2outputdir_tag(mta=maximum_triangle_area,
146                        yieldstep=yieldstep,
147                        width=width,
148                        use_limits=use_limits,
149                        friction=friction,
150                        end_tag=end_tag,
151                        outputdir_name=outputdir_name
152                        )
153
154    metadata_dic = set_z_origin_to_water_depth(metadata_dic)   
155       
156    pro_instance = project.Project(['data','flumes','Hinwood_2008'],
157                                   outputdir_name=outputdir_name)
158    print "The output dir is", pro_instance.outputdir
159    copy_code_files(pro_instance.outputdir,__file__,
160                    dirname(project.__file__) \
161                    + sep + project.__name__+'.py')
162    copy (pro_instance.codedir + 'run_dam.py',
163          pro_instance.outputdir + 'run_dam.py')
164    copy (pro_instance.codedir + 'create_mesh.py',
165          pro_instance.outputdir + 'create_mesh.py')
166
167    boundary_final_time = prepare_time_boundary(metadata_dic,
168                                       pro_instance.raw_data_dir,
169                                       pro_instance.boundarydir)
170    #return pro_instance
171    if finaltime is None:
172        finaltime = boundary_final_time - 0.1 # Edge boundary problems
173    # Boundary file manipulation
174    if boundary_path is None:
175        boundary_path = pro_instance.boundarydir
176    boundary_file_path = join(boundary_path, boundary_file)
177   #  # Convert the boundary file, .csv to .tsm
178#     try:
179#         temp = open(boundary_file_path)
180#         temp.close()
181#     except IOError:
182#         prepare_time_boundary(boundary_file_path)
183   
184    mesh_filename = pro_instance.meshdir + basename + '.msh'
185
186    #--------------------------------------------------------------------------
187    # Copy scripts to output directory and capture screen
188    # output to file
189    #--------------------------------------------------------------------------
190
191    # creates copy of code in output dir
192    if run_type >= 1:
193        #start_screen_catcher(pro_instance.outputdir, rank, pypar.size())
194        start_screen_catcher(pro_instance.outputdir)
195
196    print 'USER:    ', pro_instance.user
197    #-------------------------------------------------------------------------
198    # Create the triangular mesh
199    #-------------------------------------------------------------------------
200
201    # this creates the mesh
202    #gate_position = 12.0
203    create_mesh.generate(mesh_filename, metadata_dic, width=width,
204                         maximum_triangle_area=maximum_triangle_area)
205
206    head,tail = path.split(mesh_filename)
207    copy (mesh_filename,
208          pro_instance.outputdir + tail )
209    #-------------------------------------------------------------------------
210    # Setup computational domain
211    #-------------------------------------------------------------------------
212    domain = Domain(mesh_filename, use_cache = False, verbose = True)
213   
214
215    print 'Number of triangles = ', len(domain)
216    print 'The extent is ', domain.get_extent()
217    print domain.statistics()
218
219   
220    domain.set_name(basename)
221    domain.set_datadir(pro_instance.outputdir)
222    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
223    domain.set_minimum_storable_height(0.0001)
224
225    if use_limits is True:
226        domain.set_default_order(2) # Use second order spatial scheme
227        domain.set_timestepping_method('rk2')
228        domain.use_edge_limiter = True
229        domain.tight_slope_limiters = True
230       
231        domain.beta_w      = 0.6
232        domain.beta_uh     = 0.6
233        domain.beta_vh     = 0.6
234   
235
236    #-------------------------------------------------------------------------
237    # Setup initial conditions
238    #-------------------------------------------------------------------------
239
240    domain.set_quantity('stage', 0.) #the origin is the still water level
241    domain.set_quantity('friction', friction)
242    elevation_function = Elevation_function(metadata_dic)
243    domain.set_quantity('elevation', elevation_function)
244
245   
246    print 'Available boundary tags', domain.get_boundary_tags()
247
248    # Create boundary function from timeseries provided in file
249    #function = file_function(project.boundary_file, domain, verbose=True)
250    #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
251    try:
252        function = file_function(boundary_file_path, domain,
253                                 verbose=True)
254    except IOError:
255        msg = 'Run prepare_time_boundary.py. File "%s" could not be opened.'\
256                  %(pro_instance.boundary_file)
257        raise msg
258       
259    Br = Reflective_boundary(domain)
260    Bd = Dirichlet_boundary([0.3,0,0]) 
261    Bts = Time_boundary(domain, function)
262    #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
263    domain.set_boundary( {'wall': Br, 'wave': Bts} )
264    #domain.set_boundary( {'wall': Br, 'wave': Bd} )
265
266    #-------------------------------------------------------------------------
267    # Evolve system through time
268    #-------------------------------------------------------------------------
269    t0 = time.time()
270
271    # It seems that ANUGA can't handle a starttime that is >0.
272    #domain.starttime = 1.0 #!!! what was this doing?
273    for t in domain.evolve(yieldstep, finaltime):   
274        domain.write_time()
275    print 'That took %.2f seconds' %(time.time()-t0)
276    print 'finished'
277
278    flume_y_middle = 0.0
279    points = []
280    for gauge_x in metadata_dic['gauge_x']:
281        points.append([gauge_x, flume_y_middle])
282
283
284    #-------------------------------------------------------------------------
285    # Calculate gauge info
286    #-------------------------------------------------------------------------
287
288    if run_type >= 1:
289        id = metadata_dic['scenario_id'] + ".csv"
290        interpolate_sww2csv(pro_instance.outputdir + basename +".sww",
291                            points,
292                            pro_instance.outputdir + "depth_" + id,
293                            pro_instance.outputdir + "velocity_x_" + id,
294                            pro_instance.outputdir + "velocity_y_" + id,
295                            pro_instance.outputdir + "stage_" + id,
296                            pro_instance.outputdir + "froude_" + id)
297 
298    return pro_instance
299
300
301
302def paras2outputdir_tag(mta,
303                        yieldstep,
304                        width,
305                        use_limits,
306                        friction,
307                        end_tag,
308                        outputdir_name=None
309                        ):
310    """
311    Given run parameters return a string that can be used to name files.
312    """
313
314    if outputdir_name is None:
315        outputdir_name = ''
316       
317    if use_limits is True:
318        outputdir_name += '_lmts'
319    else:
320        outputdir_name += '_nolmts'
321    outputdir_name += '_wdth_' + str(width)
322    outputdir_name += '_z_' + str(friction)
323    outputdir_name += '_ys_' + str(yieldstep)
324    outputdir_name += '_mta_' + str(mta)
325    outputdir_name += end_tag
326
327    return outputdir_name
328
329   
330def set_z_origin_to_water_depth(seabed_coords):
331    offset = seabed_coords['offshore_water_depth']
332    keys = ['xleft', 'xtoe', 'xbeach', 'xright']
333    for x in keys:
334            seabed_coords[x][1] -= offset
335    return seabed_coords
336#-------------------------------------------------------------
337if __name__ == "__main__":
338   
339    from scenarios import scenarios
340    from slope import gauges_for_slope
341    #from plot import plot
342
343    # 1 is fast and dirty
344    # 4 is 0.01
345    # 5 is 0.001
346    # 6 is 0.0001
347    # 7 is 0.00001
348    # 8 is 0.00001  yieldstep = 0.05
349   
350    run_type = 1
351    #run_type = 0  # for testing
352   
353    #for run_data in [scenarios[5]]:
354    #scenarios = scenarios[3:]
355    #scenarios = [scenarios[0]]
356   
357    #width = 1.0
358    width = 0.1
359    #width = 0.01
360
361    maximum_triangle_area=0.01
362    yieldstep = 0.01
363    friction=0.0
364   
365    for run_data in scenarios:
366        pro_instance = main( run_data['scenario_id'] + '_boundary.tsm'  ,
367                             run_data,
368                             maximum_triangle_area=maximum_triangle_area,
369                             yieldstep=yieldstep,
370                             width=width,
371                             run_type=run_type,
372                             outputdir_name=run_data['scenario_id'],
373                             use_limits=False,
374                             friction=friction,
375                             end_tag='_Z')
376        #gauges_for_slope(pro_instance.outputdir,[run_data])
Note: See TracBrowser for help on using the repository browser.