source: anuga_validation/Hinwood_2008/run_dam.py @ 5710

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

comments

File size: 7.1 KB
Line 
1"""
2
3Script for running a breaking wave simulation of Jon Hinwood's wave tank.
4
5Duncan Gray, GA - 2007
6"""
7
8
9#----------------------------------------------------------------------------
10# Import necessary modules
11#----------------------------------------------------------------------------
12
13# Standard modules
14import time
15from os.path import join
16
17# Related major packages
18from anuga.shallow_water import Domain, Reflective_boundary, Time_boundary
19from anuga.fit_interpolate.interpolate import interpolate_sww2csv
20from anuga.abstract_2d_finite_volumes.util import file_function
21from anuga.utilities.interp import interp
22
23# Scenario specific imports
24import project
25import create_mesh
26from prepare_time_boundary import csv2tms
27
28
29def main(metadata_dic,
30         maximum_triangle_area,
31         yieldstep,
32         boundary_path=None,
33         friction=0.,  # planed wood 0.012 http://www.lmnoeng.com/manningn.htm
34         outputdir_name=None,
35         isTest=False,
36         width=1.0,
37         use_limits=True,
38         end_tag = ''):
39    """
40    Run a simulation.
41    """
42
43    id = metadata_dic['scenario_id']
44   
45    if isTest is True:
46        yieldstep = 1.0
47        finaltime = 15.
48        maximum_triangle_area=0.1
49        outputdir_name += '_test'         
50    else:
51        finaltime = None
52       
53    outputdir_name = paras2outputdir_tag(mta=maximum_triangle_area,
54                        yieldstep=yieldstep,
55                        width=width,
56                        use_limits=use_limits,
57                        friction=friction,
58                        end_tag=end_tag,
59                        outputdir_name=outputdir_name)
60    basename = outputdir_name       
61
62    if finaltime is None:
63        finaltime = metadata_dic['wave_times'][1] + 0.1 
64       
65    boundary_file_path = id + '_boundary.tsm'
66       
67    # Convert the boundary file, .csv to .tsm
68    csv2tms(boundary_file_path, metadata_dic['xleft'][1])
69   
70    mesh_filename =  join(project.output_dir,basename + '.msh')
71
72    #-------------------------------------------------------------------------
73    # Create the triangular mesh
74    #-------------------------------------------------------------------------
75
76    create_mesh.generate(mesh_filename, metadata_dic, width=width,
77                         maximum_triangle_area=maximum_triangle_area)
78   
79    #-------------------------------------------------------------------------
80    # Setup computational domain
81    #-------------------------------------------------------------------------
82    domain = Domain(mesh_filename, use_cache = False, verbose = True) 
83
84    print 'Number of triangles = ', len(domain)
85    print 'The extent is ', domain.get_extent()
86    print domain.statistics()
87   
88    domain.set_name(basename)
89    domain.set_datadir(project.output_dir)
90    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
91    domain.set_minimum_storable_height(0.0001)
92
93    if use_limits is True:
94        domain.set_default_order(2) # Use second order spatial scheme
95        domain.set_timestepping_method('rk2')
96        domain.use_edge_limiter = True
97        domain.tight_slope_limiters = True
98       
99        domain.beta_w      = 0.6
100        domain.beta_uh     = 0.6
101        domain.beta_vh     = 0.6
102   
103
104    #-------------------------------------------------------------------------
105    # Setup initial conditions
106    #-------------------------------------------------------------------------
107
108    domain.set_quantity('stage', 0.) #the z origin is the still water level
109    domain.set_quantity('friction', friction)
110    elevation_function = Elevation_function(metadata_dic)
111    domain.set_quantity('elevation', elevation_function)
112
113    function = file_function(boundary_file_path, domain, verbose=True)
114       
115    Br = Reflective_boundary(domain)
116    Bts = Time_boundary(domain, function)
117    domain.set_boundary( {'wall': Br, 'wave': Bts} )
118
119    #-------------------------------------------------------------------------
120    # Evolve system through time
121    #-------------------------------------------------------------------------
122    t0 = time.time()
123
124    for t in domain.evolve(yieldstep, finaltime):   
125        domain.write_time()
126    print 'That took %.2f seconds' %(time.time()-t0)
127    print 'finished'
128
129
130    #-------------------------------------------------------------------------
131    # Calculate gauge info
132    #-------------------------------------------------------------------------
133
134    if isTest is not True:
135        points = []
136        for gauge_x in metadata_dic['gauge_x']:
137            points.append([gauge_x, 0.0])
138
139        interpolate_sww2csv(join(project.output_dir,basename+".sww"),
140                            points,
141                            join(project.output_dir,basename + "_depth.csv"),
142                            join(project.output_dir,basename + \
143                                 "_velocit_x.csv"),
144                           join(project.output_dir, basename + \
145                                "_velocity_y.csv"),
146                            join(project.output_dir,basename + "_stage.csv"))
147 
148
149def paras2outputdir_tag(mta,
150                        yieldstep,
151                        width,
152                        use_limits,
153                        friction,
154                        end_tag,
155                        outputdir_name=None):
156    """
157    Given run parameters return a string that can be used to name files.
158    """
159    if outputdir_name is None:
160        outputdir_name = ''
161       
162    if use_limits is True:
163        outputdir_name += '_lmts'
164    else:
165        outputdir_name += '_nolmts'
166    outputdir_name += '_wdth_' + str(width)
167    outputdir_name += '_z_' + str(friction)
168    outputdir_name += '_ys_' + str(yieldstep)
169    outputdir_name += '_mta_' + str(mta)
170    outputdir_name += end_tag
171
172    return outputdir_name
173
174
175class Elevation_function:
176    """
177   
178    Create a callable instance that returns the bed surface of the
179    flume given the scenario metadata.
180   
181    """
182    def __init__(self, metadata_dic):
183        self.xslope_position = [metadata_dic['xleft'][0],
184                                metadata_dic['xtoe'][0],
185                                metadata_dic['xbeach'][0],
186                                metadata_dic['xright'][0]]
187        self.yslope_height = [metadata_dic['xleft'][1],
188                              metadata_dic['xtoe'][1],
189                              metadata_dic['xbeach'][1],
190                              metadata_dic['xright'][1]]
191       
192    def __call__(self, x,y):       
193        z = interp(self.yslope_height, self.xslope_position, x)
194        return z
195
196#-------------------------------------------------------------
197if __name__ == "__main__":
198   
199    from scenarios import scenarios
200
201    #scenarios = [scenarios[0]]
202   
203    width = 0.1
204    maximum_triangle_area=0.01
205    yieldstep = 0.01
206    #yieldstep = 0.5
207    friction=0.0
208    isTest=True
209    isTest=False
210   
211    for run_data in scenarios:
212        main(run_data,
213             maximum_triangle_area=maximum_triangle_area,
214             yieldstep=yieldstep,
215             width=width,
216             isTest=isTest,
217             outputdir_name=run_data['scenario_id'],
218             use_limits=False,
219             friction=friction,
220             end_tag='_A')
Note: See TracBrowser for help on using the repository browser.