source: production/pt_hedland_2006/run_pt_hedland.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.1 KB
Line 
1"""Script for running a tsunami inundation scenario for Onslow, 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, \
24     dem2pts
25from anuga.pyvolution.combine_pts import combine_rectangular_points_files
26from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
27from shutil import copy
28from os import mkdir, access, F_OK
29from anuga.geospatial_data.geospatial_data import *
30import sys
31from anuga.pyvolution.util import Screen_Catcher
32
33# Application specific imports
34import project                 # Definition of file names and polygons
35
36#-------------------------------------------------------------------------------
37# Copy scripts to time stamped output directory and capture screen
38# output to file
39#-------------------------------------------------------------------------------
40
41# creates copy of code in output dir if dir doesn't exist
42if access(project.outputtimedir,F_OK) == 0 :
43    mkdir (project.outputtimedir)
44copy (dirname(project.__file__) +sep+ project.__name__+'.py', project.outputtimedir + project.__name__+'.py')
45copy (__file__, project.outputtimedir + basename(__file__))
46print 'project.outputtimedir',project.outputtimedir
47
48# normal screen output is stored in
49screen_output_name = project.outputtimedir + "screen_output.txt"
50screen_error_name = project.outputtimedir + "screen_error.txt"
51
52# used to catch screen output to file
53sys.stdout = Screen_Catcher(screen_output_name)
54sys.stderr = Screen_Catcher(screen_error_name)
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# Do for coarse and fine data
62# Fine pts file to be clipped to area of interest
63#-------------------------------------------------------------------------------
64
65# filenames
66meshname = project.meshname+'.msh'
67source_dir = project.boundarydir
68
69# fine data (clipping the points file to smaller area)
70# creates DEM from asc data
71convert_dem_from_ascii2netcdf(project.onshore_dem_name, use_cache=True, verbose=True)
72
73#creates pts file from DEM
74dem2pts(project.onshore_dem_name,
75        easting_min=project.eastingmin,
76        easting_max=project.eastingmax,
77        northing_min=project.northingmin,
78        northing_max= project.northingmax,
79        use_cache=True,
80        verbose=True)
81
82print 'create G1'
83G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')
84print 'create G2'
85G2 = Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')
86print 'create G3'
87G3 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
88print 'create G4'
89G4 = Geospatial_data(file_name = project.coast_dem_name + '.xya')
90print 'add G1+G2+G3+G4'
91G = G1 + G2 + G3 + G4
92print 'export G'
93G.export_points_file(project.combined_dem_name + '.pts')
94
95#-------------------------------------------------------------------------------                                 
96# Create the triangular mesh based on overall clipping polygon with a tagged
97# boundary and interior regions defined in project.py along with
98# resolutions (maximal area of per triangle) for each polygon
99#-------------------------------------------------------------------------------
100
101from pmesh.mesh_interface import create_mesh_from_regions
102
103region_res = 500000
104coast_res = 500
105pt_hedland_res = 5000
106interior_regions = [[project.poly_pt_hedland, pt_hedland_res],
107                    [project.poly_region, region_res]]
108
109print 'number of interior regions', len(interior_regions)
110
111from anuga.utilities.polygon import plot_polygons
112if sys.platform == 'win32':
113    #figname = project.outputtimedir + 'pt_hedland_polys'
114    figname = 'pt_hedland_polys_test'
115    plot_polygons([project.polyAll,project.poly_pt_hedland,project.poly_region],
116              figname,
117              verbose = True)   
118
119print 'start create mesh from regions'
120from caching import cache
121_ = cache(create_mesh_from_regions,
122          project.polyAll,
123          {'boundary_tags': {'topright': [0], 'topleft': [1],
124                             'left': [2], 'bottom0': [3],
125                             'bottom1': [4], 'bottom2': [5],
126                             'bottom3': [6], 'right': [7]},
127           'maximum_triangle_area': 250000,
128           'filename': meshname,           
129           'interior_regions': interior_regions},
130          verbose = True, evaluate=True)
131
132#-------------------------------------------------------------------------------                                 
133# Setup computational domain
134#-------------------------------------------------------------------------------                                 
135domain = Domain(meshname, use_cache = False, verbose = True)
136
137print domain.statistics()
138print 'Number of triangles = ', len(domain)
139print 'The extent is ', domain.get_extent()
140print domain.statistics()
141
142domain.set_name(project.basename)
143domain.set_datadir(project.outputtimedir)
144domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
145
146#-------------------------------------------------------------------------------                                 
147# Setup initial conditions
148#-------------------------------------------------------------------------------
149
150tide = 0.
151#high
152#tide = 3.6
153#low
154#tide = -3.9
155
156domain.set_quantity('stage', tide)
157domain.set_quantity('friction', 0.0) 
158print 'hi and file',project.combined_dem_name + '.pts'
159
160domain.set_quantity('elevation', 
161                    filename = project.combined_dem_name + '.pts',
162                    use_cache = True,
163                    verbose = True,
164                    alpha = 0.1
165                    )
166
167#-------------------------------------------------------------------------------                                 
168# Setup boundary conditions (all reflective)
169#-------------------------------------------------------------------------------
170print 'start ferret2sww'
171# skipped as results in file SU-AU_clipped is correct for all WA
172
173from anuga.pyvolution.data_manager import ferret2sww
174
175south = project.south
176north = project.north
177west = project.west
178east = project.east
179
180#note only need to do when an SWW file for the MOST boundary doesn't exist
181cache(ferret2sww,
182      (source_dir + project.boundary_basename,
183       source_dir + project.boundary_basename+'_'+project.basename), 
184      {'verbose': True,
185       'minlat': south,
186       'maxlat': north,
187       'minlon': west,
188       'maxlon': east,
189#       'origin': project.mesh_origin,
190       'origin': domain.geo_reference.get_origin(),
191       'mean_stage': tide,
192       'zscale': 1,                 #Enhance tsunami
193       'fail_on_NaN': False,
194       'inverted_bathymetry': True},
195       evaluate = True,
196       verbose = True,
197      dependencies = source_dir + project.boundary_basename + '.sww')
198
199print 'Available boundary tags', domain.get_boundary_tags()
200
201Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
202                    domain, verbose = True)
203Br = Reflective_boundary(domain)
204Bd = Dirichlet_boundary([tide,0,0])
205domain.set_boundary( {'topright': Bf,'topleft': Bf, 'left':  Bd, 'bottom0': Bd,
206                      'bottom1': Bd, 'bottom2': Bd, 'bottom3': Bd, 
207                        'right': Bd})
208
209#-------------------------------------------------------------------------------                                 
210# Evolve system through time
211#-------------------------------------------------------------------------------
212import time
213t0 = time.time()
214
215for t in domain.evolve(yieldstep = 240, finaltime = 10800): 
216    domain.write_time()
217    domain.write_boundary_statistics(tags = 'topright')     
218
219for t in domain.evolve(yieldstep = 120, finaltime = 16200
220                       ,skip_initial_step = True): 
221    domain.write_time()
222    domain.write_boundary_statistics(tags = 'topright')     
223
224for t in domain.evolve(yieldstep = 60, finaltime = 21600
225                       ,skip_initial_step = True): 
226    domain.write_time()
227    domain.write_boundary_statistics(tags = 'topright')     
228   
229for t in domain.evolve(yieldstep = 120, finaltime = 27000
230                       ,skip_initial_step = True): 
231    domain.write_time()
232    domain.write_boundary_statistics(tags = 'topright')     
233
234for t in domain.evolve(yieldstep = 240, finaltime = 36000
235                       ,skip_initial_step = True): 
236    domain.write_time()
237    domain.write_boundary_statistics(tags = 'topright')   
238 
239print 'That took %.2f seconds' %(time.time()-t0)
240
241print 'finished'
Note: See TracBrowser for help on using the repository browser.