source: anuga_work/production/sydney_2006/run_sydney_slide.py @ 4276

Last change on this file since 4276 was 4276, checked in by sexton, 17 years ago

update project files for slide modelling

File size: 9.0 KB
Line 
1"""Script for running a tsunami inundation scenario for Sydney, NSW, Australia.
2
3Source data such as elevation and boundary data is assumed to be available in
4directories specified by project_slide.py
5The output sww file is stored in project_slide.outputtimedir
6
7The scenario is defined by a triangular mesh created from project_slide.polygon,
8the elevation data and a tsunami wave generated by s submarine mass failure.
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.path import dirname, basename
22from os import mkdir, access, F_OK, sep
23import sys
24
25# Related major packages
26from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary
27from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
28from anuga.geospatial_data.geospatial_data import *
29from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
30
31# Application specific imports
32import project_slide              # Definition of file names and polygons
33
34#-------------------------------------------------------------------------------
35# Copy scripts to time stamped output directory and capture screen
36# output to file
37#-------------------------------------------------------------------------------
38
39# creates copy of code in output dir
40copy_code_files(project_slide.outputtimedir,__file__,dirname(project_slide.__file__)+sep+ project_slide.__name__+'.py' )
41myid = 0
42numprocs = 1
43start_screen_catcher(project_slide.outputtimedir, myid, numprocs)
44
45print 'USER:    ', project_slide.user
46
47#-------------------------------------------------------------------------------
48# Preparation of topographic data
49#
50# Convert ASC 2 DEM 2 PTS using source data and store result in source data
51#-------------------------------------------------------------------------------
52
53# filenames
54on_offshore25_dem_name = project_slide.on_offshore25_dem_name
55on_offshore100_dem_name = project_slide.on_offshore100_dem_name
56nsw_dem_name = project_slide.nsw_dem_name
57meshname = project_slide.meshname+'.msh'
58
59# creates DEM from asc data
60convert_dem_from_ascii2netcdf(on_offshore25_dem_name, use_cache=True, verbose=True)
61convert_dem_from_ascii2netcdf(on_offshore100_dem_name, use_cache=True, verbose=True)
62convert_dem_from_ascii2netcdf(nsw_dem_name, use_cache=True, verbose=True)
63
64#creates pts file for onshore DEM
65dem2pts(on_offshore25_dem_name,
66        easting_min=project_slide.eastingmin25,
67        easting_max=project_slide.eastingmax25,
68        northing_min=project_slide.northingmin25,
69        northing_max= project_slide.northingmax25,
70        use_cache=True, verbose=True)
71dem2pts(on_offshore100_dem_name,
72        easting_min=project_slide.eastingmin100,
73        easting_max=project_slide.eastingmax100,
74        northing_min=project_slide.northingmin100,
75        northing_max= project_slide.northingmax100,
76        use_cache=True, verbose=True)
77dem2pts(nsw_dem_name,
78        easting_min=project_slide.eastingmin_nsw,
79        easting_max=project_slide.eastingmax_nsw,
80        northing_min=project_slide.northingmin_nsw,
81        northing_max= project_slide.northingmax_nsw,
82        use_cache=True, verbose=True)
83
84print 'create offshore'
85G11 = Geospatial_data(file_name = project_slide.offshore_dem_name1 + '.xya')+\
86      Geospatial_data(file_name = project_slide.offshore_dem_name2 + '.xya')+\
87      Geospatial_data(file_name = project_slide.offshore_dem_name3 + '.xya')
88G12 = Geospatial_data(file_name = project_slide.offshore_dem_name4 + '.xya')+\
89      Geospatial_data(file_name = project_slide.offshore_dem_name5 + '.xya')+\
90      Geospatial_data(file_name = project_slide.offshore_dem_name6 + '.xya')+\
91      Geospatial_data(file_name = project_slide.offshore_dem_name7 + '.xya')+\
92      Geospatial_data(file_name = project_slide.offshore_dem_name8 + '.xya')+\
93      Geospatial_data(file_name = project_slide.offshore_dem_name9 + '.xya')
94print 'create onshore'
95G2 = Geospatial_data(file_name = project_slide.on_offshore25_dem_name + '.pts')
96G3 = Geospatial_data(file_name = project_slide.on_offshore100_dem_name + '.pts')
97G4 = Geospatial_data(file_name = project_slide.nsw_dem_name + '.pts')
98print 'add'
99G = G11.clip(Geospatial_data(project_slide.poly_surveyclip)) +\
100    G12.clip(Geospatial_data(project_slide.polyAll)) +\
101    G2.clip(Geospatial_data(project_slide.poly_25mclip)) +\
102    G3.clip(Geospatial_data(project_slide.poly_origsyd)) +\
103    (G4.clip(Geospatial_data(project_slide.polyAll))).clip_outside(Geospatial_data(project_slide.poly_surveyclip)).clip_outside(Geospatial_data(project_slide.poly_origsyd))
104print 'export points'
105G.export_points_file(project_slide.combined_dem_name + '.pts')
106#G.export_points_file(project_slide.combined_dem_name + '.xya')
107
108#----------------------------------------------------------------------------
109# Create the triangular mesh based on overall clipping polygon with a tagged
110# boundary and interior regions defined in project_slide.py along with
111# resolutions (maximal area of per triangle) for each polygon
112#-------------------------------------------------------------------------------
113
114from anuga.pmesh.mesh_interface import create_mesh_from_regions
115remainder_res = 250000.
116local_res = 50000.
117coast_res = 500.
118interior_regions = [[project_slide.poly_syd1, local_res],
119                    [project_slide.poly_coast, coast_res]]
120
121from caching import cache
122_ = cache(create_mesh_from_regions,
123          project_slide.polyAll,
124           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
125                              'e3': [3], 'e4':[4]},
126           'maximum_triangle_area': remainder_res,
127           'filename': meshname,
128           'interior_regions': interior_regions},
129          verbose = True, evaluate=False)
130print 'created mesh'
131
132#-------------------------------------------------------------------------------                                 
133# Setup computational domain
134#-------------------------------------------------------------------------------                                 
135domain = Domain(meshname, use_cache = True, verbose = True)
136
137print 'Number of triangles = ', len(domain)
138print 'The extent is ', domain.get_extent()
139print domain.statistics()
140
141domain.set_name(project_slide.basename)
142domain.set_datadir(project_slide.outputtimedir)
143domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
144domain.set_minimum_storable_height(0.01)
145
146#-------------------------------------------------------------------------------                                 
147# Setup initial conditions
148#-------------------------------------------------------------------------------
149
150tide = 0.0
151domain.set_quantity('stage', tide)
152domain.set_quantity('friction', 0.0) 
153domain.set_quantity('elevation', 
154                    filename = project_slide.combined_dem_name + '.pts',
155                    use_cache = True,
156                    verbose = True,
157                    alpha = 0.1
158                    )
159
160#-------------------------------------------------------------------------------
161# Set up scenario (tsunami_source is a callable object used with set_quantity)
162#-------------------------------------------------------------------------------
163from smf import slide_tsunami
164
165tsunami_source = slide_tsunami(length=project_slide.yacaaba_length,
166                               width=project_slide.yacaaba_width,
167                               depth=project_slide.yacaaba_depth,
168                               slope=project_slide.yacaaba_slope,
169                               thickness=project_slide.yacaaba_thickness, 
170                               x0=project_slide.slide_origin_yacaaba_b[0], 
171                               y0=project_slide.slide_origin_yacaaba_b[1], 
172                               alpha=project_slide.yacaaba_alpha, 
173                               domain=domain)
174
175#-------------------------------------------------------------------------------                                 
176# Setup boundary conditions
177#-------------------------------------------------------------------------------
178print 'Available boundary tags', domain.get_boundary_tags()
179
180Br = Reflective_boundary(domain)
181Bd = Dirichlet_boundary([tide,0,0])
182
183domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd} )
184
185
186#-------------------------------------------------------------------------------                                 
187# Evolve system through time
188#-------------------------------------------------------------------------------
189import time
190from Numeric import allclose
191from anuga.abstract_2d_finite_volumes.quantity import Quantity
192
193t0 = time.time()
194
195for t in domain.evolve(yieldstep = 30, finaltime = 5000): 
196    domain.write_time()
197    domain.write_boundary_statistics(tags = 'e2')
198    stagestep = domain.get_quantity('stage') 
199
200    if allclose(t, 30):
201        slide = Quantity(domain)
202        slide.set_values(tsunami_source)
203        domain.set_quantity('stage', slide + stagestep)
204   
205print 'That took %.2f seconds' %(time.time()-t0)
206
207print 'finished'
Note: See TracBrowser for help on using the repository browser.