source: anuga_validation/Hinwood_2008/run_dam.py @ 5693

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

Adding a header to the boundary files

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