source: anuga_work/production/MOST_example/run_pt_hedland.py @ 3608

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

write meshes to the output file

File size: 8.3 KB
Line 
1"""Script for running a tsunami inundation scenario for Onslow, WA, Australia.
2
3 Boundary data is assumed to be in this directory.
4Move the SU-AU_clip_blahblah.nc files into this directory, from teh pt_hedland
5data directory.
6
7The file cripple_pts_file creates the elevation file
8pt_hedland_combined_elevation_31204.pts, from
9pt_hedland_combined_elevation.pts - which is created by the other
10run_pt_hedland script.
11 
12directories specified by project.py
13The output sww file is stored in project.outputtimedir
14
15The scenario is defined by a triangular mesh created from project.polygon,
16the elevation data and a simulated submarine landslide.
17
18Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
19"""
20#-------------------------------------------------------------------------------# Import necessary modules
21#-------------------------------------------------------------------------------
22
23# Standard modules
24from os import sep
25from os.path import dirname, basename
26import time
27
28# Related major packages
29from anuga.shallow_water import Domain, Reflective_boundary, \
30                            Dirichlet_boundary, Time_boundary, File_boundary
31from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, \
32     dem2pts, ferret2sww
33from shutil import copy
34from os import mkdir, access, F_OK
35from anuga.geospatial_data.geospatial_data import *
36import sys
37from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
38
39from anuga.fit_interpolate.fit import fit_to_mesh_file
40
41# Application specific imports
42import project                 # Definition of file names and polygons
43
44#-------------------------------------------------------------------------------
45# Copy scripts to time stamped output directory and capture screen
46# output to file
47#-------------------------------------------------------------------------------
48
49# creates copy of code in output dir if dir doesn't exist
50if access(project.outputtimedir,F_OK) == 0 :
51    mkdir (project.outputtimedir)
52copy (dirname(project.__file__) +sep+ project.__name__+'.py', project.outputtimedir + project.__name__+'.py')
53copy (__file__, project.outputtimedir + basename(__file__))
54print 'project.outputtimedir',project.outputtimedir
55
56# normal screen output is stored in
57#screen_output_name = project.outputtimedir + "screen_output.txt"
58#screen_error_name = project.outputtimedir + "screen_error.txt"
59
60# used to catch screen output to file
61#sys.stdout = Screen_Catcher(screen_output_name)
62#sys.stderr = Screen_Catcher(screen_error_name)
63print 'USER:    ', project.user
64
65#-------------------------------------------------------------------------------
66# Preparation of topographic data
67#
68# Convert ASC 2 DEM 2 PTS using source data and store result in source data
69# Do for coarse and fine data
70# Fine pts file to be clipped to area of interest
71#-------------------------------------------------------------------------------
72
73# filenames
74meshname = project.meshname + '.tsh'
75mesh_elevname = project.mesh_elevname  + '.tsh'
76source_dir = project.boundarydir
77
78
79#-------------------------------------------------------------------------------                                 
80# Create the triangular mesh based on overall clipping polygon with a tagged
81# boundary and interior regions defined in project.py along with
82# resolutions (maximal area of per triangle) for each polygon
83#-------------------------------------------------------------------------------
84
85from anuga.pmesh.mesh_interface import create_mesh_from_regions
86
87region_res = 5000000
88coast_res = 500000
89pt_hedland_res = 500000
90interior_regions = [[project.poly_pt_hedland, pt_hedland_res],
91                    [project.poly_region, region_res]]
92
93print 'number of interior regions', len(interior_regions)
94
95print 'start create mesh from regions'
96from caching import cache
97_ = cache(create_mesh_from_regions,
98          project.polyAll,
99          {'boundary_tags': {'topright': [0], 'topleft': [1],
100                             'left': [2], 'bottom0': [3],
101                             'bottom1': [4], 'bottom2': [5],
102                             'bottom3': [6], 'right': [7]},
103           'maximum_triangle_area': 5000000,
104           'filename': meshname,           
105           'interior_regions': interior_regions},
106          verbose = True
107          #, evaluate=True
108          )
109
110cache(fit_to_mesh_file,(meshname,
111                 'pt_hedland_combined_elevation_31204' + '.pts',
112                 mesh_elevname),
113      {'verbose': True}
114      #,evaluate = True     
115      ,verbose = False
116      )
117#-------------------------------------------------------------------------------                                 
118# Setup computational domain
119#-------------------------------------------------------------------------------                                 
120domain = Domain(mesh_elevname, use_cache = False, verbose = True)
121
122print domain.statistics()
123print 'Number of triangles = ', len(domain)
124print 'The extent is ', domain.get_extent()
125print domain.statistics()
126
127domain.set_name(project.basename)
128domain.set_datadir(project.outputtimedir)
129domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
130
131#-------------------------------------------------------------------------------                                 
132# Setup initial conditions
133#-------------------------------------------------------------------------------
134
135tide = 0.
136#high
137#tide = 3.6
138#low
139#tide = -3.9
140
141domain.set_quantity('stage', tide)
142domain.set_quantity('friction', 0.0) 
143print 'hi and file',project.combined_dem_name + '.pts'
144
145#domain.set_quantity('elevation',
146#                    filename = project.combined_dem_name + '.pts',
147#                    use_cache = True,
148#                    verbose = True,
149#                    alpha = 0.1
150 #                   )
151
152#-------------------------------------------------------------------------------                                 
153# Setup boundary conditions (all reflective)
154#-------------------------------------------------------------------------------
155print 'start ferret2sww'
156# skipped as results in file SU-AU_clipped is correct for all WA
157
158
159south = project.south
160north = project.north
161west = project.west
162east = project.east
163
164#note only need to do when an SWW file for the MOST boundary doesn't exist
165cache(ferret2sww,
166      (project.boundary_basename,
167       project.boundary_basename+'_'+project.basename), 
168      {'verbose': True,
169       'minlat': south,
170       'maxlat': north,
171       'minlon': west,
172       'maxlon': east,
173#       'origin': project.mesh_origin,
174       'origin': domain.geo_reference.get_origin(),
175       'mean_stage': tide,
176       'zscale': 10,                 #Enhance tsunami
177       'fail_on_NaN': False,
178       'inverted_bathymetry': True},
179       #evaluate = True,
180       verbose = True,
181      dependencies = source_dir + project.boundary_basename + '.sww')
182
183print 'Available boundary tags', domain.get_boundary_tags()
184
185Bf = File_boundary(project.boundary_basename+'_'+project.basename + '.sww', 
186                    domain, verbose = True)
187Br = Reflective_boundary(domain)
188Bd = Dirichlet_boundary([tide,0,0])
189domain.set_boundary( {'topright': Bf,'topleft': Bf, 'left':  Bd, 'bottom0': Bd,
190                      'bottom1': Bd, 'bottom2': Bd, 'bottom3': Bd, 
191                        'right': Bd})
192
193#-------------------------------------------------------------------------------                                 
194# Evolve system through time
195#-------------------------------------------------------------------------------
196import time
197t0 = time.time()
198
199for t in domain.evolve(yieldstep = 240, finaltime = 10800): 
200    domain.write_time()
201    domain.write_boundary_statistics(tags = 'topright')     
202
203for t in domain.evolve(yieldstep = 120, finaltime = 16200
204                       ,skip_initial_step = True): 
205    domain.write_time()
206    domain.write_boundary_statistics(tags = 'topright')     
207
208for t in domain.evolve(yieldstep = 60, finaltime = 21600
209                       ,skip_initial_step = True): 
210    domain.write_time()
211    domain.write_boundary_statistics(tags = 'topright')     
212   
213for t in domain.evolve(yieldstep = 120, finaltime = 27000
214                       ,skip_initial_step = True): 
215    domain.write_time()
216    domain.write_boundary_statistics(tags = 'topright')     
217
218for t in domain.evolve(yieldstep = 240, finaltime = 36000
219                       ,skip_initial_step = True): 
220    domain.write_time()
221    domain.write_boundary_statistics(tags = 'topright')   
222 
223print 'That took %.2f seconds' %(time.time()-t0)
224
225print 'finished'
Note: See TracBrowser for help on using the repository browser.