source: anuga_work/production/karratha_2006/run_karratha.py @ 3627

Last change on this file since 3627 was 3627, checked in by ole, 18 years ago

New files for Karratha study

File size: 6.0 KB
Line 
1"""Script for running a tsunami inundation scenario for Broome, 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# Import necessary modules
14#------------------------------------------------------------------------------
15
16# Standard modules
17from os import sep
18from os.path import dirname, basename
19import time
20
21# Related major packages
22from anuga.shallow_water import Domain
23from anuga.shallow_water import Dirichlet_boundary
24from anuga.shallow_water import File_boundary
25from anuga.shallow_water import Reflective_boundary
26
27from anuga.pmesh.mesh_interface import create_mesh_from_regions
28
29from shutil import copy
30from os import mkdir, access, F_OK
31from anuga.geospatial_data.geospatial_data import *
32import sys
33from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
34
35# Application specific imports
36import project                 # Definition of file names and polygons
37
38from anuga_parallel.parallel_api import distribute, myid
39
40#------------------------------------------------------------------------------
41# Copy scripts to time stamped output directory and capture screen
42# output to file
43#------------------------------------------------------------------------------
44
45# filenames
46meshname = project.meshname+'.msh'
47source_dir = project.boundarydir
48
49
50
51if myid == 0:
52    # creates copy of code in output dir if dir doesn't exist
53    if access(project.outputtimedir,F_OK) == 0:
54        mkdir (project.outputtimedir)
55    copy (dirname(project.__file__) +sep+ project.__name__+'.py',
56          project.outputtimedir + project.__name__+'.py')
57    copy (__file__, project.outputtimedir + basename(__file__))
58    print 'project.outputtimedir',project.outputtimedir
59
60    # normal screen output is stored in
61    screen_output_name = project.outputtimedir + 'screen_output.txt'
62    screen_error_name = project.outputtimedir + 'screen_error.txt'
63
64    # used to catch screen output to file
65    #sys.stdout = Screen_Catcher(screen_output_name)
66    #sys.stderr = Screen_Catcher(screen_error_name)
67    print 'USER:    ', project.user
68
69
70    #--------------------------------------------------------------------------
71    # Create the triangular mesh based on overall clipping polygon with a tagged
72    # boundary and interior regions defined in project.py along with
73    # resolutions (maximal area of per triangle) for each polygon
74    #--------------------------------------------------------------------------
75
76
77    resolution = 4000
78    interior_regions = [[project.neil1_polygon, resolution],
79                        [project.neil2_polygon, 64000]]
80
81    print 'number of interior regions', len(interior_regions)
82
83
84    print 'start create mesh from regions'
85    from caching import cache
86
87    meshname = project.meshname + '_%d.msh' %myid
88    _ = cache(create_mesh_from_regions,
89              project.bounding_polygon,
90              {'boundary_tags': {'back': [7, 8], 'side': [0, 6],
91                                 'ocean': [1, 2, 3, 4, 5]},
92               'maximum_triangle_area': 100000,
93               'filename': meshname},
94              #'interior_regions': interior_regions},
95              verbose = True,
96              evaluate = False)
97
98    #------------------------------------------------------------------------- 
99    # Setup computational domain
100    #-------------------------------------------------------------------------
101
102    domain = Domain(meshname, use_cache = True, verbose = True)
103    print domain.statistics()
104
105
106    domain.set_name(project.basename)
107    domain.set_datadir(project.outputtimedir)
108
109
110    #-------------------------------------------------------------------------
111    # Setup initial conditions
112    #-------------------------------------------------------------------------
113
114    tide = 0.
115   
116    domain.set_quantity('stage', tide)
117    domain.set_quantity('friction', 0.0) 
118    domain.set_quantity('elevation', 
119                        filename = project.datadir + project.basename + '.pts',
120                        use_cache = False,
121                        verbose = True,
122                        alpha = 0.1
123                        )
124   
125    #-------------------------------------------------------------------------
126    # Setup boundary conditions
127    #-------------------------------------------------------------------------
128
129    print 'Available boundary tags', domain.get_boundary_tags()
130
131    Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
132                       domain, verbose = True)
133    Br = Reflective_boundary(domain)
134    Bd = Dirichlet_boundary([tide,0,0])
135    domain.set_boundary({'back': Br,
136                         'side': Bd,
137                         #'ocean': None}) # Bind this one later
138                         'ocean': Bf}) # Bind this one later   
139else:
140    domain = None
141
142
143#--------------------------------------------------------------------------
144# Create the parallel domain
145#--------------------------------------------------------------------------
146domain = distribute(domain, verbose=True)
147
148
149# Set those boundaries that can't be communicated automatically
150Bf = File_boundary(source_dir + project.boundary_basename + '.sww',
151                   domain, verbose = True)                   
152boundary_map = domain.boundary_map
153boundary_map['ocean'] = Bf
154
155print boundary_map
156domain.set_boundary(boundary_map)
157
158
159
160#----------------------------------------------------------------------------
161# Evolve system through time
162#----------------------------------------------------------------------------
163import time
164t0 = time.time()
165
166for t in domain.evolve(yieldstep = 60, finaltime = 28600):
167    domain.write_time()
168    domain.write_boundary_statistics(tags = 'ocean')     
169   
170print 'That took %.2f seconds' %(time.time()-t0)
171
Note: See TracBrowser for help on using the repository browser.