source: anuga_validation/UQ_runup_2006/run_dam.py @ 4748

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

changes to matts script

File size: 5.5 KB
Line 
1"""Script for running a dam break simulation of UQ's dam break tank.
2
3
4Ole Nielsen and Duncan Gray, GA - 2006
5"""
6
7
8#----------------------------------------------------------------------------
9# Import necessary modules
10#----------------------------------------------------------------------------
11
12# Standard modules
13import time
14import sys
15from shutil import copy
16from os import path
17
18# Related major packages
19from anuga.shallow_water import Domain, Reflective_boundary, \
20                            Dirichlet_boundary, Time_boundary, File_boundary
21from anuga.abstract_2d_finite_volumes.region import Set_region
22from anuga.fit_interpolate.interpolate import interpolate_sww2csv
23from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \
24     copy_code_files
25import create_mesh
26
27# Application specific imports
28
29
30def main():
31     
32    slopes = [0]
33    frictions = [0.01] #can enter multiple frictions eg [0.00, 0.01]
34    inital_depths = [0.2] #can enter multiple initial depths eg [0.1, 0.2]
35    gate_positions = [0.75]
36    flume_lengths = [3] #can enter multiple flume lengths eg [3.0, 5.0]
37    # NOTE flume_lengths is not used.  It's hardcoded into create_mesh
38   
39    for slope in slopes:
40        for friction in frictions:
41            for inital_depth in inital_depths:
42                for flume_length in flume_lengths:
43                    for gate_position in gate_positions:
44                        scenario(slope, friction, inital_depth, gate_position, flume_length)
45           
46
47def scenario(slope, friction, inital_depth, gate_position, flume_length):
48   
49    import project                 # Definition of file names and polygons
50
51    #-------------------------------------------------------------------------
52    # Setup archiving of simulations
53    #-------------------------------------------------------------------------
54
55    id = 'd'+ str(inital_depth) + '_s'+str(slope)+'_g'+ str(gate_position)+ '_n'+ str(friction)
56    copy (project.codedirname, project.outputtimedir + 'project.py')
57    run_name = 'run_dam.py'
58    run_name_out = 'run_dam'+id+'.py'
59
60    copy (project.codedir + run_name, project.outputtimedir + run_name_out)
61    copy (project.codedir + 'create_mesh.py',
62          project.outputtimedir + 'create_mesh.py')
63    print'output dir', project.outputtimedir
64
65    #-------------------------------------------------------------------------
66    # Create the triangular mesh
67    #-------------------------------------------------------------------------
68
69    create_mesh.generate(project.mesh_filename,
70                         gate_position, flume_length,
71                         #is_course=True) # this creates the mesh
72                         is_course=False) # this creates the mesh
73
74    head,tail = path.split(project.mesh_filename)
75    copy (project.mesh_filename,
76          project.outputtimedir + tail )
77    #-------------------------------------------------------------------------
78    # Setup computational domain
79    #-------------------------------------------------------------------------
80    domain = Domain(project.mesh_filename, use_cache = False, verbose = True)
81   
82
83    print 'Number of triangles = ', len(domain)
84    print 'The extent is ', domain.get_extent()
85    print domain.statistics()
86
87   
88    domain.set_name(project.basename + id)
89    domain.set_datadir(project.outputtimedir)
90    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
91    domain.set_minimum_storable_height(0.001)
92    domain.set_store_vertices_uniquely(True)  # for writting to sww
93
94    #-------------------------------------------------------------------------
95    # Setup initial conditions
96    #-------------------------------------------------------------------------
97
98
99    def elevation_tilt(x, y):
100        return x*slope
101       
102    domain.set_quantity('stage', elevation_tilt)
103    domain.set_quantity('friction', friction) 
104    domain.set_quantity('elevation',elevation_tilt)
105
106    print 'Available boundary tags', domain.get_boundary_tags()
107    domain.set_region('dam','stage',inital_depth,
108                                 location = 'unique vertices') 
109
110    Br = Reflective_boundary(domain)
111    Bd = Dirichlet_boundary([0,0,0])  # to drain the water out.
112    domain.set_boundary( {'wall': Br, 'edge': Bd} )
113
114    #-------------------------------------------------------------------------
115    # Evolve system through time
116    #-------------------------------------------------------------------------
117    t0 = time.time()
118
119    for t in domain.evolve(yieldstep = 0.1, finaltime = 10): #enter timestep and final time
120        domain.write_time()
121   
122        print 'That took %.2f seconds' %(time.time()-t0)
123        print 'finished'
124       
125    points = [[gate_position - 0.65,0.2],
126              [gate_position - 0.55,0.2],
127              [gate_position - 0.45,0.2],
128              [gate_position - 0.35,0.2],
129              [0.45,0.2],
130              [gate_position - 0.25,0.2]
131              ]
132    #-------------------------------------------------------------------------
133    # Calculate gauge info
134    #-------------------------------------------------------------------------
135    print "name", project.outputtimedir + project.basename \
136                   + id+".sww"
137    interpolate_sww2csv(project.outputtimedir + project.basename \
138                   + id+".sww",
139                        points,
140                        project.depth_filename + id + '.csv',
141                        project.velocity_x_filename + id + '.csv',
142                        project.velocity_y_filename + id + '.csv')
143
144
145#-------------------------------------------------------------
146if __name__ == "__main__":
147    main()
148   
Note: See TracBrowser for help on using the repository browser.