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

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

updates (i) export to xya for Broome (ii) gauge investigation for Hobart (compare MOST and ANUGA) (iii) script to compare onslow outputs with change in parameters

File size: 9.1 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')
106G.export_points_file(project.combined_dem_name + '.xya')
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.