source: development/dam_2006/run_dam.py @ 3447

Last change on this file since 3447 was 3447, checked in by duncan, 18 years ago

Making set_region easier to use.

File size: 4.9 KB
Line 
1"""Script for running a tsunami inundation scenario for Onslow, WA, Australia.
2
3Source data such as elevation and boundary data is assumed to be available in
4directories specified by project.py
5The output sww file is stored in project.outputtimedir
6
7The scenario is defined by a triangular mesh created from project.polygon,
8the elevation data and a simulated submarine landslide.
9
10Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
11"""
12
13
14#----------------------------------------------------------------------------
15# Import necessary modules
16#----------------------------------------------------------------------------
17
18# Standard modules
19import time
20import sys
21from shutil import copy
22from os import mkdir, access, F_OK, path
23
24# Related major packages
25from pyvolution.shallow_water import Domain, Reflective_boundary, \
26                            Dirichlet_boundary, Time_boundary, File_boundary
27from pyvolution.util import Screen_Catcher
28from pyvolution.region import Set_region
29from fit_interpolate.interpolate import interpolate_sww2csv
30import create_mesh
31
32# Application specific imports
33import project                 # Definition of file names and polygons
34
35def main():
36    #-------------------------------------------------------------------------
37    # Setup archiving of simulations
38    #-------------------------------------------------------------------------
39
40    copy (project.codedirname, project.outputtimedir + 'project.py')
41    copy (project.codedir + 'run_dam.py', project.outputtimedir + 'run_dam.py')
42    copy (project.codedir + 'create_mesh.py',
43          project.outputtimedir + 'create_mesh.py')
44    print'output dir', project.outputtimedir
45
46    #FIXME this isn't working
47    #normal screen output is stored in
48    screen_output_name = project.outputtimedir + "screen_output.txt"
49    screen_error_name = project.outputtimedir + "screen_error.txt"
50
51    #-------------------------------------------------------------------------
52    # Create the triangular mesh
53    #-------------------------------------------------------------------------
54
55    gate_position = 0.85
56    create_mesh.generate(project.mesh_filename,
57                         gate_position,
58                         is_course=True) # this creates the mesh
59                         #is_course=False) # this creates the mesh
60
61    head,tail = path.split(project.mesh_filename)
62    copy (project.mesh_filename,
63          project.outputtimedir + tail )
64    #-------------------------------------------------------------------------
65    # Setup computational domain
66    #-------------------------------------------------------------------------
67    domain = Domain(project.mesh_filename, use_cache = False, verbose = True)
68   
69
70    print 'Number of triangles = ', len(domain)
71    print 'The extent is ', domain.get_extent()
72    print domain.statistics()
73
74    domain.set_name(project.basename)
75    domain.set_datadir(project.outputtimedir)
76    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
77
78    #-------------------------------------------------------------------------
79    # Setup initial conditions
80    #-------------------------------------------------------------------------
81
82    tide = 0.0
83    slope = 0.05
84
85    def elevation_tilt(x, y):
86        return x*slope
87       
88    domain.set_quantity('stage', elevation_tilt)
89    domain.set_quantity('friction', 0.03) 
90    domain.set_quantity('elevation',elevation_tilt)
91
92    print 'Available boundary tags', domain.get_boundary_tags()
93    domain.set_region('dam','stage',0.2,
94                                 location = 'unique vertices') 
95    domain.set_region(Set_region('dam','stage',0.2,
96                                 location = 'unique vertices'))
97    domain.set_store_vertices_uniquely(True)  # for writting to sww
98
99    Br = Reflective_boundary(domain)
100    Bd = Dirichlet_boundary([0,0,0])
101    domain.set_boundary( {'wall': Br, 'edge': Bd} )
102
103    #-------------------------------------------------------------------------
104    # Evolve system through time
105    #-------------------------------------------------------------------------
106    t0 = time.time()
107
108    for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
109        domain.write_time()
110   
111        print 'That took %.2f seconds' %(time.time()-t0)
112        print 'finished'
113
114    points = [[gate_position - 0.65,0.2],
115              [gate_position - 0.55,0.2],
116              [gate_position - 0.45,0.2],
117              [gate_position - 0.35,0.2],
118              [gate_position - 0.25,0.2],
119              [0.55,0.2]
120              ]
121
122    #-------------------------------------------------------------------------
123    # Calc gauge info
124    #-------------------------------------------------------------------------
125    interpolate_sww2csv(project.outputtimedir + project.basename +".sww",
126                        points,
127                        project.depth_filename,
128                        project.velocity_filename)
129
130#-------------------------------------------------------------
131if __name__ == "__main__":
132    main()
133   
Note: See TracBrowser for help on using the repository browser.