source: production/broome_2006/run_broome.py @ 3514

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

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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 simulated submarine landslide.
9
10Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
11"""
12#-------------------------------------------------------------------------------# Import necessary modules
13#-------------------------------------------------------------------------------
14
15# Standard modules
16from os import sep
17from os.path import dirname, basename
18import time
19
20# Related major packages
21from anuga.pyvolution.shallow_water import Domain, Reflective_boundary, \
22                            Dirichlet_boundary, Time_boundary, File_boundary
23from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
24from anuga.pyvolution.combine_pts import combine_rectangular_points_files
25from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
26from shutil import copy
27from os import mkdir, access, F_OK
28from anuga.geospatial_data.geospatial_data import *
29import sys
30from anuga.pyvolution.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 (dirname(project.__file__) +sep+ project.__name__+'.py', project.outputtimedir + project.__name__+'.py')
44copy (__file__, project.outputtimedir + basename(__file__))
45print 'project.outputtimedir',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)
54print 'USER:    ', project.user
55
56#-------------------------------------------------------------------------------
57# Preparation of topographic data
58#
59# Convert ASC 2 DEM 2 PTS using source data and store result in source data
60# Do for coarse and fine data
61# Fine pts file to be clipped to area of interest
62#-------------------------------------------------------------------------------
63
64# filenames
65meshname = project.meshname+'.msh'
66source_dir = project.boundarydir
67
68# creates DEM from asc data
69convert_dem_from_ascii2netcdf(project.onshore_dem_name, use_cache=True, verbose=True)
70
71#creates pts file from DEM
72dem2pts(project.onshore_dem_name,
73        easting_min=project.eastingmin,
74        easting_max=project.eastingmax,
75        northing_min=project.northingmin,
76        northing_max= project.northingmax,
77        use_cache=True,
78        verbose=True)
79
80print 'create G1'
81G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')
82print 'create G2'
83G2 = Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')
84print 'create G3'
85G3 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
86print 'create G4'
87G4 = Geospatial_data(file_name = project.coast_dem_name + '.xya')
88print 'add G1+G2+G3+G4'
89G = G1 + G2 + G3 + G4
90print 'export G'
91G.export_points_file(project.combined_dem_name + '.pts')
92
93#-------------------------------------------------------------------------------                                 
94# Create the triangular mesh based on overall clipping polygon with a tagged
95# boundary and interior regions defined in project.py along with
96# resolutions (maximal area of per triangle) for each polygon
97#-------------------------------------------------------------------------------
98
99from pmesh.mesh_interface import create_mesh_from_regions
100
101region_res = 500000
102coast_res = 500
103broome_res = 5000
104interior_regions = [[project.poly_broome, broome_res],
105                    [project.poly_region, region_res]]
106
107print 'number of interior regions', len(interior_regions)
108
109from anuga.utilities.polygon import plot_polygons
110if sys.platform == 'win32':
111    #figname = project.outputtimedir + 'pt_hedland_polys'
112    figname = 'broome_polys_test'
113    plot_polygons([project.polyAll,project.poly_broome,project.poly_region],
114              figname,
115              verbose = True)   
116
117print 'start create mesh from regions'
118from caching import cache
119_ = cache(create_mesh_from_regions,
120          project.polyAll,
121          {'boundary_tags': {'topright': [0], 'topleft': [1],
122                             'left': [2], 'bottom0': [3],
123                             'bottom1': [4], 'bottom2': [5],
124                             'bottom3': [6], 'right': [7]},
125           'maximum_triangle_area': 250000,
126           'filename': meshname,           
127           'interior_regions': interior_regions},
128          verbose = True, evaluate=True)
129
130#-------------------------------------------------------------------------------                                 
131# Setup computational domain
132#-------------------------------------------------------------------------------                                 
133domain = Domain(meshname, use_cache = False, verbose = True)
134
135print domain.statistics()
136print 'Number of triangles = ', len(domain)
137print 'The extent is ', domain.get_extent()
138print domain.statistics()
139
140domain.set_name(project.basename)
141domain.set_datadir(project.outputtimedir)
142domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
143
144#-------------------------------------------------------------------------------                                 
145# Setup initial conditions
146#-------------------------------------------------------------------------------
147
148tide = 0.
149#high
150#tide =
151#low
152#tide =
153
154domain.set_quantity('stage', tide)
155domain.set_quantity('friction', 0.0) 
156print 'hi and file',project.combined_dem_name + '.pts'
157
158domain.set_quantity('elevation', 
159                    filename = project.combined_dem_name + '.pts',
160                    use_cache = True,
161                    verbose = True,
162                    alpha = 0.1
163                    )
164
165#-------------------------------------------------------------------------------                                 
166# Setup boundary conditions
167#-------------------------------------------------------------------------------
168print 'start ferret2sww'
169# skipped as results in file SU-AU_clipped is correct for all WA
170
171from anuga.pyvolution.data_manager import ferret2sww
172
173south = project.south
174north = project.north
175west = project.west
176east = project.east
177
178#note only need to do when an SWW file for the MOST boundary doesn't exist
179cache(ferret2sww,
180      (source_dir + project.boundary_basename,
181       source_dir + project.boundary_basename+'_'+project.basename), 
182      {'verbose': True,
183       'minlat': south,
184       'maxlat': north,
185       'minlon': west,
186       'maxlon': east,
187#       'origin': project.mesh_origin,
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
197print 'Available boundary tags', domain.get_boundary_tags()
198
199Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
200                    domain, verbose = True)
201Br = Reflective_boundary(domain)
202Bd = Dirichlet_boundary([tide,0,0])
203domain.set_boundary( {'topright': Bf,'topleft': Bf, 'left':  Bd, 'bottom0': Bd,
204                      'bottom1': Bd, 'bottom2': Bd, 'bottom3': Bd, 
205                        'right': Bd})
206
207#-------------------------------------------------------------------------------                                 
208# Evolve system through time
209#-------------------------------------------------------------------------------
210import time
211t0 = time.time()
212
213for t in domain.evolve(yieldstep = 240, finaltime = 10800): 
214    domain.write_time()
215    domain.write_boundary_statistics(tags = 'topright')     
216
217for t in domain.evolve(yieldstep = 120, finaltime = 16200
218                       ,skip_initial_step = True): 
219    domain.write_time()
220    domain.write_boundary_statistics(tags = 'topright')     
221
222for t in domain.evolve(yieldstep = 60, finaltime = 21600
223                       ,skip_initial_step = True): 
224    domain.write_time()
225    domain.write_boundary_statistics(tags = 'topright')     
226   
227for t in domain.evolve(yieldstep = 120, finaltime = 27000
228                       ,skip_initial_step = True): 
229    domain.write_time()
230    domain.write_boundary_statistics(tags = 'topright')     
231
232for t in domain.evolve(yieldstep = 240, finaltime = 36000
233                       ,skip_initial_step = True): 
234    domain.write_time()
235    domain.write_boundary_statistics(tags = 'topright')   
236 
237print 'That took %.2f seconds' %(time.time()-t0)
238
239print 'finished'
Note: See TracBrowser for help on using the repository browser.