source: anuga_work/development/dam_2006/run_dam.py @ 4993

Last change on this file since 4993 was 4044, checked in by duncan, 17 years ago

make friction a parameter

File size: 4.8 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, sep
17from os.path import dirname  #, basename
18
19# Related major packages
20from anuga.shallow_water import Domain, Reflective_boundary, \
21                            Dirichlet_boundary, Time_boundary, File_boundary
22from anuga.abstract_2d_finite_volumes.region import Set_region
23from anuga.fit_interpolate.interpolate import interpolate_sww2csv
24from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \
25     copy_code_files
26
27
28# Scenario specific imports
29import project                 # Definition of file names and polygons
30import create_mesh
31
32
33def main(friction):
34
35    #--------------------------------------------------------------------------
36    # Copy scripts to time stamped output directory and capture screen
37    # output to file
38    #--------------------------------------------------------------------------
39
40    # creates copy of code in output dir
41
42    print "The output dir is", project.outputtimedir
43    copy_code_files(project.outputtimedir,__file__,
44                    dirname(project.__file__)+sep+ project.__name__+'.py' )
45    copy (project.codedir + 'create_mesh.py',
46          project.outputtimedir + 'create_mesh.py')
47    myid = 0
48    numprocs = 1
49    #start_screen_catcher(project.outputtimedir, myid, numprocs)
50
51    print 'USER:    ', project.user
52    #-------------------------------------------------------------------------
53    # Create the triangular mesh
54    #-------------------------------------------------------------------------
55
56    gate_position = 0.85
57    create_mesh.generate(project.mesh_filename,
58                         gate_position,
59                         #is_coarse=True) # this creates the mesh
60                         is_coarse=False) # this creates the mesh
61
62    head,tail = path.split(project.mesh_filename)
63    copy (project.mesh_filename,
64          project.outputtimedir + tail )
65    #-------------------------------------------------------------------------
66    # Setup computational domain
67    #-------------------------------------------------------------------------
68    domain = Domain(project.mesh_filename, use_cache = False, verbose = True)
69   
70
71    print 'Number of triangles = ', len(domain)
72    print 'The extent is ', domain.get_extent()
73    print domain.statistics()
74
75   
76    domain.set_name(project.basename)
77    domain.set_datadir(project.outputtimedir)
78    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
79    domain.set_minimum_storable_height(0.01)
80    #domain.set_store_vertices_uniquely(True)  # for writting to sww
81
82    #-------------------------------------------------------------------------
83    # Setup initial conditions
84    #-------------------------------------------------------------------------
85
86    slope = 0.05
87
88    def elevation_tilt(x, y):
89        return x*slope
90       
91    domain.set_quantity('stage', elevation_tilt)
92    domain.set_quantity('friction', friction) 
93    domain.set_quantity('elevation',elevation_tilt)
94
95    print 'Available boundary tags', domain.get_boundary_tags()
96    domain.set_region('dam','stage',0.20,
97                                 location = 'unique vertices') 
98    #domain.set_region(Set_region('dam','stage',0.2,
99    #                             location = 'unique vertices'))
100
101    Br = Reflective_boundary(domain)
102    Bd = Dirichlet_boundary([0,0,0])  # to drain the water out.
103    domain.set_boundary( {'wall': Br, 'edge': Bd} )
104
105    #-------------------------------------------------------------------------
106    # Evolve system through time
107    #-------------------------------------------------------------------------
108    t0 = time.time()
109
110    for t in domain.evolve(yieldstep = 0.01, finaltime = 30):
111        domain.write_time()
112   
113        print 'That took %.2f seconds' %(time.time()-t0)
114        print 'finished'
115
116    points = [[gate_position - 0.65,0.2],
117              [gate_position - 0.55,0.2],
118              [gate_position - 0.45,0.2],
119              [gate_position - 0.35,0.2],
120              [gate_position - 0.25,0.2]
121              ]
122
123    #-------------------------------------------------------------------------
124    # Calculate gauge info
125    #-------------------------------------------------------------------------
126    interpolate_sww2csv(project.outputtimedir + project.basename +".sww",
127                        points,
128                        project.depth_filename,
129                        project.velocity_x_filename,
130                        project.velocity_y_filename)
131
132#-------------------------------------------------------------
133if __name__ == "__main__":
134    main(0.000)
135   
Note: See TracBrowser for help on using the repository browser.