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

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

update broome to incorporate new copy files and screen catcher

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