source: anuga_work/production/broome_2006/run_broome.py @ 3951

Last change on this file since 3951 was 3951, checked in by sexton, 17 years ago

update broome script to incorporate additional data delivered today

File size: 9.3 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 tsunami wave generated by MOST.
9
10Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
11"""
12
13#-------------------------------------------------------------------------------
14# Import necessary modules
15#-------------------------------------------------------------------------------
16
17# Standard modules
18import os
19import time
20from shutil import copy
21from os import mkdir, access, F_OK
22import sys
23
24# Related major packages
25from anuga.shallow_water import Domain, Reflective_boundary, \
26                            Dirichlet_boundary, Time_boundary, File_boundary
27from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
28from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files
29from anuga.geospatial_data.geospatial_data import *
30from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
31
32# Application specific imports
33import project                 # Definition of file names and polygons
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 if dir doesn't exist
41if access(project.outputtimedir,F_OK) == 0 :
42    mkdir (project.outputtimedir)
43copy (project.codedirname, project.outputtimedir + project.codename)
44copy (project.codedir + 'run_broome.py', project.outputtimedir + 'run_broome.py')
45print'output dir', project.outputtimedir
46
47#normal screen output is stored in
48screen_output_name = project.outputtimedir + "screen_output.txt"
49screen_error_name = project.outputtimedir + "screen_error.txt"
50
51#used to catch screen output to file
52sys.stdout = Screen_Catcher(screen_output_name)
53sys.stderr = Screen_Catcher(screen_error_name)
54
55print 'USER:    ', project.user
56
57#-------------------------------------------------------------------------------
58# Preparation of topographic data
59#
60# Convert ASC 2 DEM 2 PTS using source data and store result in source data
61#-------------------------------------------------------------------------------
62
63# filenames
64onshore_dem_name = project.onshore_dem_name
65offshore_interp_dem_name = project.offshore_interp_dem_name
66coast_points = project.coast_dem_name
67meshname = project.meshname+'.msh'
68
69# creates DEM from asc data
70convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
71
72#creates pts file for onshore DEM
73dem2pts(onshore_dem_name, use_cache=True, verbose=True)
74
75# creates DEM from asc data
76convert_dem_from_ascii2netcdf(offshore_interp_dem_name, use_cache=True, verbose=True)
77
78#creates pts file for offshore interpolated DEM
79dem2pts(offshore_interp_dem_name, use_cache=True, verbose=True)
80
81print 'create offshore'
82G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')+\
83     Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')+\
84     Geospatial_data(file_name = project.offshore_dem_name3 + '.xya')+\
85     Geospatial_data(file_name = project.offshore_dem_name4 + '.xya')+\
86     Geospatial_data(file_name = project.offshore_dem_name5 + '.xya')+\
87     Geospatial_data(file_name = project.offshore_dem_name6 + '.xya')+\
88     Geospatial_data(file_name = project.offshore_dem_name7 + '.xya')+\
89     Geospatial_data(file_name = project.offshore_dem_name8 + '.xya')+\
90     Geospatial_data(file_name = project.offshore_dem_name9 + '.xya')+\
91     Geospatial_data(file_name = project.offshore_dem_name10 + '.xya')+\
92     Geospatial_data(file_name = project.offshore_dem_name11 + '.xya')+\
93     Geospatial_data(file_name = project.offshore_dem_name12 + '.xya')+\
94     Geospatial_data(file_name = project.offshore_dem_name13 + '.xya')+\
95     Geospatial_data(file_name = project.offshore_dem_name14 + '.xya')+\
96     Geospatial_data(file_name = project.offshore_dem_name15 + '.xya')+\
97     Geospatial_data(file_name = project.offshore_dem_name16 + '.xya')+\
98     Geospatial_data(file_name = project.offshore_dem_name17 + '.xya')+\
99     Geospatial_data(file_name = project.offshore_dem_name18 + '.xya')+\
100     Geospatial_data(file_name = project.offshore_dem_name19 + '.xya')+\
101     Geospatial_data(file_name = project.offshore_dem_name20 + '.xya')+\
102     Geospatial_data(file_name = project.offshore_dem_name21 + '.xya')+\
103     Geospatial_data(file_name = project.offshore_dem_name22 + '.xya')+\
104     Geospatial_data(file_name = project.offshore_interp_dem_name + '.pts')
105print 'create onshore'
106G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
107print 'create coast'
108G3 = Geospatial_data(file_name = project.coast_dem_name + '.xya')
109print 'add'
110G = G1 + G2 + G3
111print 'export points'
112G.export_points_file(project.combined_dem_name + '.pts')
113
114#----------------------------------------------------------------------------
115# Create the triangular mesh based on overall clipping polygon with a tagged
116# boundary and interior regions defined in project.py along with
117# resolutions (maximal area of per triangle) for each polygon
118#-------------------------------------------------------------------------------
119
120from anuga.pmesh.mesh_interface import create_mesh_from_regions
121remainder_res = 750000
122local_res = 25000
123broome_res = 5000
124coast_res = 500
125interior_regions = [[project.poly_broome1, local_res],
126                    [project.poly_broome2, broome_res],
127                    [project.poly_broome3, coast_res]]
128
129from project import number_mesh_triangles
130print 'estimate of number of triangles', \
131      number_mesh_triangles(interior_regions, project.polyAll, remainder_res)
132
133from caching import cache
134_ = cache(create_mesh_from_regions,
135          project.polyAll,
136           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
137                              'e3': [3], 'e4':[4], 'e5': [5],
138                              'e6': [6]},
139           'maximum_triangle_area': remainder_res,
140           'filename': meshname,
141           'interior_regions': interior_regions},
142          verbose = True, evaluate=False)
143print 'created mesh'
144
145#-------------------------------------------------------------------------------                                 
146# Setup computational domain
147#-------------------------------------------------------------------------------                                 
148domain = Domain(meshname, use_cache = True, verbose = True)
149
150print 'Number of triangles = ', len(domain)
151print 'The extent is ', domain.get_extent()
152print domain.statistics()
153
154domain.set_name(project.basename)
155domain.set_datadir(project.outputtimedir)
156domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
157domain.set_minimum_storable_height(0.01)
158
159#-------------------------------------------------------------------------------                                 
160# Setup initial conditions
161#-------------------------------------------------------------------------------
162
163tide = 0.0
164domain.set_quantity('stage', tide)
165domain.set_quantity('friction', 0.0) 
166domain.set_quantity('elevation', 
167                    filename = project.combined_dem_name + '.pts',
168                    use_cache = True,
169                    verbose = True,
170                    alpha = 0.1
171                    )
172
173#-------------------------------------------------------------------------------                                 
174# Setup boundary conditions
175#-------------------------------------------------------------------------------
176'''
177print 'start urs2sww'
178print '', project.boundary_basename
179from anuga.shallow_water.data_manager import urs2sww
180
181south = project.south
182north = project.north
183west  = project.west
184east  = project.east
185
186#note only need to do when an SWW file for the MOST boundary doesn't exist
187cache(urs2sww,
188      (source_dir + project.boundary_basename,
189       source_dir + project.boundary_basename),
190      {'verbose': True,
191       'minlat': south,
192       'maxlat': north,
193       'minlon': west,
194       'maxlon': east,
195       #'origin': domain.geo_reference.get_origin(),
196       'mean_stage': tide,
197       'zscale': 1,                 #Enhance tsunami
198       'fail_on_NaN': False,
199       'inverted_bathymetry': True},
200      #evaluate = True,
201       verbose = True,
202       dependencies = source_dir + project.boundary_basename + '.sww')
203
204'''
205print 'Available boundary tags', domain.get_boundary_tags()
206
207#Bf = File_boundary(source_dir + project.boundary_basename + '.sww',
208#                    domain, verbose = True)
209Br = Reflective_boundary(domain)
210Bd = Dirichlet_boundary([tide,0,0])
211
212# 7 min square wave starting at 1 min, 6m high
213Bw = Time_boundary(domain = domain,
214                   f=lambda t: [(60<t<480)*10, 0, 0])
215
216domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd,
217                      'e5': Bd,  'e6': Bd} )
218
219
220#-------------------------------------------------------------------------------                                 
221# Evolve system through time
222#-------------------------------------------------------------------------------
223import time
224t0 = time.time()
225
226for t in domain.evolve(yieldstep = 240, finaltime = 480): 
227    domain.write_time()
228    domain.write_boundary_statistics(tags = 'e14')     
229   
230print 'That took %.2f seconds' %(time.time()-t0)
231
232print 'finished'
Note: See TracBrowser for help on using the repository browser.