source: production/onslow_2006/run_onslow.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.7 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
13
14#-------------------------------------------------------------------------------# Import necessary modules
15#-------------------------------------------------------------------------------
16
17# Standard modules
18import os
19import time
20from shutil import copy
21from os import mkdir, access, F_OK
22import sys
23
24# Related major packages
25from anuga.pyvolution.shallow_water import Domain, Reflective_boundary, \
26                            Dirichlet_boundary, Time_boundary, File_boundary
27from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
28from anuga.pyvolution.combine_pts import combine_rectangular_points_files
29from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
30from anuga.geospatial_data.geospatial_data import *
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 (project.codedirname, project.outputtimedir + project.codename)
45copy (project.codedir + 'run_onslow.py', project.outputtimedir + 'run_onslow.py')
46print'output dir', 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)
54#sys.stderr = Screen_Catcher(screen_output_name)
55sys.stderr = Screen_Catcher(screen_error_name)
56
57print 'USER:    ', project.user
58
59#-------------------------------------------------------------------------------
60# Preparation of topographic data
61#
62# Convert ASC 2 DEM 2 PTS using source data and store result in source data
63# Do for coarse and fine data
64# Fine pts file to be clipped to area of interest
65#-------------------------------------------------------------------------------
66
67# filenames
68onshore_dem_name = project.onshore_dem_name
69islands_dem_name = project.islands_dem_name
70coast_points = project.coast_dem_name
71offshore_points = project.offshore_dem_name
72meshname = project.meshname+'.msh'
73source_dir = project.boundarydir
74
75copied_files = False
76
77# files to be used
78files_used = [onshore_dem_name, offshore_points, coast_points,]
79
80# fine data (clipping the points file to smaller area)
81# creates DEM from asc data
82convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
83
84#creates pts file for onshore DEM
85dem2pts(onshore_dem_name,
86        easting_min=project.eastingmin,
87        easting_max=project.eastingmax,
88        northing_min=project.northingmin,
89        northing_max= project.northingmax,
90        use_cache=True, 
91        verbose=True)
92
93convert_dem_from_ascii2netcdf(islands_dem_name, use_cache=True, verbose=True)
94
95#creates pts file for islands DEM
96dem2pts(islands_dem_name, use_cache=True, verbose=True)
97
98print'create G1'
99G1 = Geospatial_data(file_name = project.offshore_dem_name + '.xya')
100print'create G2'
101G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
102print'create G3'
103G3 = Geospatial_data(file_name = project.coast_dem_name + '.xya')
104print'create G4'
105G4 = Geospatial_data(file_name = project.islands_dem_name + '.pts')
106print'add G1+G2+G3+G4'
107G = G1 + G2 + G3 + G4
108print'export G'
109G.export_points_file(project.combined_dem_name + '.pts')
110
111#-------------------------------------------------------------------------------                                 
112# Create the triangular mesh based on overall clipping polygon with a tagged
113# boundary and interior regions defined in project.py along with
114# resolutions (maximal area of per triangle) for each polygon
115#-------------------------------------------------------------------------------
116
117from pmesh.mesh_interface import create_mesh_from_regions
118
119#new
120region_res = 200000
121coast_res = 25000
122onslow_res = 5000
123interior_regions = [[project.poly_onslow, onslow_res],
124                    [project.poly_coast, coast_res],
125                    [project.poly_region, region_res]]
126
127print 'number of interior regions', len(interior_regions)
128
129from caching import cache
130_ = cache(create_mesh_from_regions,
131          project.polyAll,
132          {'boundary_tags': {'top': [0], 'topleft': [1],
133                             'topleft1': [2], 'bottomleft': [3],
134                             'bottom': [4], 'bottomright': [5],
135                             'topright':[6]},
136           'maximum_triangle_area': 100000,
137           'filename': meshname,           
138           'interior_regions': interior_regions},
139          verbose = True, evaluate=True)
140
141
142#-------------------------------------------------------------------------------                                 
143# Setup computational domain
144#-------------------------------------------------------------------------------                                 
145
146domain = pmesh_to_domain_instance(meshname, Domain,
147                                  use_cache = False,
148                                  verbose = True)
149
150print 'Number of triangles = ', len(domain)
151print 'The extent is ', domain.get_extent()
152print domain.statistics()
153
154domain.set_name(project.basename)
155domain.set_datadir(project.outputtimedir)
156domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
157
158
159#-------------------------------------------------------------------------------                                 
160# Setup initial conditions
161#-------------------------------------------------------------------------------
162
163tide = 0.0
164
165domain.set_quantity('stage', tide)
166domain.set_quantity('friction', 0.0) 
167print 'hi and file',project.combined_dem_name + '.pts'
168
169domain.set_quantity('elevation', 
170#                    0.
171#                    filename = project.onshore_dem_name + '.pts',
172                    filename = project.combined_dem_name + '.pts',
173#                    filename = project.offshore_dem_name + '.pts',
174                    use_cache = True,
175                    verbose = True,
176                    alpha = 0.1
177                    )
178
179print 'hi1'
180
181#-------------------------------------------------------------------------------                                 
182# Setup boundary conditions (all reflective)
183#-------------------------------------------------------------------------------
184print 'start ferret2sww'
185from anuga.pyvolution.data_manager import ferret2sww
186
187south = project.south
188north = project.north
189west = project.west
190east = project.east
191
192#note only need to do when an SWW file for the MOST boundary doesn't exist
193cache(ferret2sww,
194      (source_dir + project.boundary_basename,
195       source_dir + project.boundary_basename), 
196#      (project.MOST_dir + project.boundary_basename,
197#       source_dir + project.boundary_basename),
198      {'verbose': True,
199# note didn't work with the below
200#       'minlat': south - 1,
201#       'maxlat': north + 1,
202#       'minlon': west - 1,
203#       'maxlon': east + 1,
204       'minlat': south,
205       'maxlat': north,
206       'minlon': west,
207       'maxlon': east,
208#       'origin': project.mesh_origin,
209       'origin': domain.geo_reference.get_origin(),
210       'mean_stage': tide,
211       'zscale': 1,                 #Enhance tsunami
212       'fail_on_NaN': False,
213       'inverted_bathymetry': True},
214      #evaluate = True,
215       verbose = True,
216       dependencies = source_dir + project.boundary_basename + '.sww')
217
218
219print 'Available boundary tags', domain.get_boundary_tags()
220
221Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
222                    domain, verbose = True)
223Br = Reflective_boundary(domain)
224Bd = Dirichlet_boundary([tide,0,0])
225
226
227# 7 min square wave starting at 1 min, 6m high
228Bw = Time_boundary(domain = domain,
229                   f=lambda t: [(60<t<480)*6, 0, 0])
230
231domain.set_boundary( {'top': Bf, 'topleft': Bf,
232                      'topleft1': Bf, 'bottomleft': Bd,
233                      'bottom': Br, 'bottomright': Br, 'topright': Bd} )
234
235#-------------------------------------------------------------------------------                                 
236# Evolve system through time
237#-------------------------------------------------------------------------------
238import time
239t0 = time.time()
240
241for t in domain.evolve(yieldstep = 240, finaltime = 7200): 
242    domain.write_time()
243    domain.write_boundary_statistics(tags = 'top')     
244
245for t in domain.evolve(yieldstep = 120, finaltime = 12600
246                       ,skip_initial_step = True): 
247    domain.write_time()
248    domain.write_boundary_statistics(tags = 'top')     
249
250for t in domain.evolve(yieldstep = 60, finaltime = 19800
251                       ,skip_initial_step = True): 
252    domain.write_time()
253    domain.write_boundary_statistics(tags = 'top')     
254   
255for t in domain.evolve(yieldstep = 120, finaltime = 25200
256                       ,skip_initial_step = True): 
257    domain.write_time()
258    domain.write_boundary_statistics(tags = 'top')     
259
260for t in domain.evolve(yieldstep = 240, finaltime = 36000
261                       ,skip_initial_step = True): 
262    domain.write_time()
263    domain.write_boundary_statistics(tags = 'top')
264   
265print 'That took %.2f seconds' %(time.time()-t0)
266
267print 'finished'
Note: See TracBrowser for help on using the repository browser.