source: anuga_work/production/onslow_2006/run_onslow_old.py @ 3650

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

change imports so reflect the new structure

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.outputdir
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#------------------------------------------------------------------------------
15# Import necessary modules
16#------------------------------------------------------------------------------
17
18# Standard modules
19import os
20import time
21
22# Related major packages
23from anuga.pyvolution.shallow_water import Domain, Reflective_boundary, \
24                            Dirichlet_boundary, Time_boundary, File_boundary
25from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
26from anuga.pyvolution.combine_pts import combine_rectangular_points_files
27from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
28#from anuga.geospatial_data.geospatial_data import *
29
30# Application specific imports
31import project                 # Definition of file names and polygons
32from smf import slump_tsunami  # Function for submarine mudslide
33
34from shutil import copy
35from os import mkdir, access, F_OK
36
37#------------------------------------------------------------------------------
38# Preparation of topographic data
39#
40# Convert ASC 2 DEM 2 PTS using source data and store result in source data
41# Do for coarse and fine data
42# Fine pts file to be clipped to area of interest
43#------------------------------------------------------------------------------
44
45# filenames
46coarsedemname = project.coarsedemname
47
48onshore_dem_name = project.onshore_dem_name
49
50meshname = project.meshname+'.msh'
51
52#############################################################
53# making a copy of project and run to a time stamped directory
54# with the ouput
55source_dir = project.boundarydir
56
57#if dir doesn't exists then makes dir
58if access(project.outputdir,F_OK) == 0 :
59    mkdir (project.outputdir)
60# creates copy of code in output dir
61copy (project.codedirname, project.outputdir + project.codename)
62copy (project.codedir + 'run_onslow.py', project.outputdir + 'run_onlsow.py')
63#print "copied to"+ project.outputdir + project.codename + 'and run_onlsow.py'
64#################################################################
65'''
66# coarse data
67convert_dem_from_ascii2netcdf(coarsedemname, use_cache=True, verbose=True)
68dem2pts(coarsedemname, use_cache=True, verbose=True)
69
70
71# fine data (clipping the points file to smaller area)
72convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
73dem2pts(onshore_dem_name,
74        easting_min=project.eastingmin,
75        easting_max=project.eastingmax,
76        northing_min=project.northingmin,
77        northing_max= project.northingmax,
78        use_cache=True,
79        verbose=True)
80
81print 'before off xya to object'
82offshore_pts = Geospatial_data(project.offshore_dem_name + '.xya')
83print 'before offshore to dict'
84offshore_dict = geospatial_data2points_dictionary(offshore_pts)
85
86print 'offshore to pts file'
87export_points_file(project.offshore_dem_name + '.pts', offshore_dict)
88'''
89# combining the coarse and fine data
90# NOTE MUST HAVE FINE FIRST!
91'''
92combine_rectangular_points_files(
93                                 project.onshore_dem_name + '.pts',
94                                 project.coarsedemname + '.pts',
95                                 project.combined_dem_name + '.pts')
96
97print 'create G1'
98G1 = Geospatial_data(project.onshore_dem_name + '.pts')
99print 'G1 dict'
100G1_points_dict = geospatial_data2points_dictionary(G1)
101print 'G1 xya file'
102export_points_file(project.datadir + 'offshore_dem.xya', G1_points_dict)
103print 'create G2'
104G2 = Geospatial_data(project.offshore_dem_name + '.xya')
105print 'G1+g2'
106G = G1 + G2
107print 'G dict'
108G_points_dict = geospatial_data2points_dictionary(G)
109print 'G to xya'
110export_points_file(project.combined_dem_name + '.xya', G_points_dict)
111
112
113add_points_files(project.onshore_dem_name + '.pts',
114                  project.offshore_dem_name + '.xya',
115                  project.combined_dem_name + '.pts')
116
117print 'please'
118
119add_points_files(project.onshore_dem_name + '.pts',
120                  project.offshore_dem_name + '.pts',
121                  project.combined_dem_name + '.pts')
122
123print ' finished with points'
124'''
125#-------------------------------------------------------------------------------                                 
126# Create the triangular mesh based on overall clipping polygon with a tagged
127# boundary and interior regions defined in project.py along with
128# resolutions (maximal area of per triangle) for each polygon
129#------------------------------------------------------------------------------
130from anuga.pmesh.mesh_interface import create_mesh_from_regions
131
132# original
133interior_res = 50000
134interior_regions = [[project.poly_onslow, interior_res],
135                    [project.poly_thevenard, interior_res],
136                    [project.poly_direction, interior_res]]
137                    #[project.testpoly, interior_res]]
138print 'number of interior regions', len(interior_regions)
139
140from caching import cache
141_ = cache(create_mesh_from_regions,
142          project.polyAll,
143          {'boundary_tags': {'top': [0], 'topleft': [1],
144                             'left': [2], 'bottom': [3],
145                             'bottomright': [4], 'topright': [5]},
146           'maximum_triangle_area': 1000000,
147           'filename': meshname,           
148           'interior_regions': interior_regions},
149          verbose = True)
150
151
152#------------------------------------------------------------------------------
153# Setup computational domain
154#------------------------------------------------------------------------------
155
156domain = pmesh_to_domain_instance(meshname, Domain,
157                                  use_cache = True,
158                                  verbose = True)
159
160print 'Number of triangles = ', len(domain)
161print 'The extent is ', domain.get_extent()
162print domain.statistics()
163
164domain.set_name(project.basename)
165domain.set_datadir(project.outputdir)
166domain.set_quantities_to_be_stored(['stage'])
167
168print 'hi'
169#------------------------------------------------------------------------------
170# Set up scenario (tsunami_source is a callable object used with set_quantity)
171#------------------------------------------------------------------------------
172'''
173tsunami_source = slump_tsunami(length=30000.0,
174                               depth=400.0,
175                               slope=6.0,
176                               thickness=176.0,
177                               radius=3330,
178                               dphi=0.23,
179                               x0=project.slump_origin[0],
180                               y0=project.slump_origin[1],
181                               alpha=0.0,
182                               domain=domain)
183
184'''
185#------------------------------------------------------------------------------
186# Setup initial conditions
187#------------------------------------------------------------------------------
188
189tide = 0.
190
191domain.set_quantity('stage', tide)
192domain.set_quantity('friction', 0.0) 
193print 'hi1', project.combined_dem_name + '.pts'
194domain.set_quantity('elevation', 
195#                    0.
196#                    filename = project.onshore_dem_name + '.pts',
197                    filename = project.combined_dem_name + '.pts',
198#                    filename = project.coarsedemname + '.pts',
199                    use_cache = True,
200                    verbose = True
201                    )
202
203print 'hi2'
204#------------------------------------------------------------------------------
205# Setup boundary conditions (all reflective)
206#------------------------------------------------------------------------------
207'''
208from anuga.pyvolution.data_manager import ferret2sww
209
210south = project.south
211north = project.north
212west = project.west
213east = project.east
214
215cache(ferret2sww,
216      (source_dir + project.boundary_basename,
217       source_dir + project.boundary_basename),
218      {'verbose': True,
219# note didn't work with the below
220#       'minlat': south - 1,
221#       'maxlat': north + 1,
222#       'minlon': west - 1,
223#       'maxlon': east + 1,
224       'minlat': south,
225       'maxlat': north,
226       'minlon': west,
227       'maxlon': east,
228#       'origin': project.mesh_origin,
229       'origin': domain.geo_reference.get_origin(),
230       'mean_stage': tide,
231       'zscale': 1,                 #Enhance tsunami
232       'fail_on_NaN': False,
233       'inverted_bathymetry': True},
234      #evaluate = True,
235       verbose = True)
236'''
237
238print 'Available boundary tags', domain.get_boundary_tags()
239
240#Bf = File_boundary(source_dir + project.boundary_basename + '.sww',
241#                    domain, verbose = True)
242Br = Reflective_boundary(domain)
243Bd = Dirichlet_boundary([tide,0,0])
244
245print 'hi3'
246# 7 min square wave starting at 1 min, 6m high
247Bw = Time_boundary(domain = domain,
248                   f=lambda t: [(60<t<480)*6, 0, 0])
249
250domain.set_boundary( {'top': Bw, 'topleft': Bw,
251                             'left': Br, 'bottom': Br,
252                             'bottomright': Br, 'topright': Br} )
253print 'hi4'
254
255#------------------------------------------------------------------------------
256# Evolve system through time
257#------------------------------------------------------------------------------
258
259import time
260t0 = time.time()
261
262for t in domain.evolve(yieldstep = 50, finaltime = 50): 
263    domain.write_time()
264    domain.write_boundary_statistics(tags = 'top')     
265   
266print 'That took %.2f seconds' %(time.time()-t0)
267print 'hi5'
Note: See TracBrowser for help on using the repository browser.