source: anuga_work/production/newcastle_2006/run_newcastle_slide.py @ 4036

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

preparing for slide scenarios

File size: 8.3 KB
Line 
1"""Script for running a tsunami inundation scenario for Newcastle, NSW, 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 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.outputtimedir,__file__,dirname(project.__file__)+sep+ project.__name__+'.py' )
41myid = 0
42numprocs = 1
43start_screen_catcher(project.outputtimedir, myid, numprocs)
44
45print 'USER:    ', project.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_offshore_dem_name = project.on_offshore_dem_name
55meshname = project.meshname+'.msh'
56
57# creates DEM from asc data
58convert_dem_from_ascii2netcdf(on_offshore_dem_name, use_cache=True, verbose=True)
59
60#creates pts file for onshore DEM
61dem2pts(on_offshore_dem_name, use_cache=True, verbose=True)
62
63print 'create offshore'
64G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')+\
65     Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')+\
66     Geospatial_data(file_name = project.offshore_dem_name3 + '.xya')+\
67     Geospatial_data(file_name = project.offshore_dem_name4 + '.xya')+\
68     Geospatial_data(file_name = project.offshore_dem_name5 + '.xya')+\
69     Geospatial_data(file_name = project.offshore_dem_name6 + '.xya')+\
70     Geospatial_data(file_name = project.offshore_dem_name7 + '.xya')+\
71     Geospatial_data(file_name = project.offshore_dem_name8 + '.xya')+\
72     Geospatial_data(file_name = project.offshore_dem_name9 + '.xya')+\
73     Geospatial_data(file_name = project.offshore_dem_name10 + '.xya')+\
74     Geospatial_data(file_name = project.offshore_dem_name11 + '.xya')+\
75     Geospatial_data(file_name = project.offshore_dem_name12 + '.xya')+\
76     Geospatial_data(file_name = project.offshore_dem_name13 + '.xya')+\
77     Geospatial_data(file_name = project.offshore_dem_name14 + '.xya')+\
78     Geospatial_data(file_name = project.offshore_dem_name15 + '.xya')+\
79     Geospatial_data(file_name = project.offshore_dem_name16 + '.xya')+\
80     Geospatial_data(file_name = project.offshore_dem_name17 + '.xya')+\
81     Geospatial_data(file_name = project.offshore_dem_name18 + '.xya')+\
82     Geospatial_data(file_name = project.offshore_dem_name19 + '.xya')+\
83     Geospatial_data(file_name = project.offshore_dem_name20 + '.xya')+\
84     Geospatial_data(file_name = project.offshore_dem_name21 + '.xya')+\
85     Geospatial_data(file_name = project.offshore_dem_name22 + '.xya')+\
86     Geospatial_data(file_name = project.offshore_interp_dem_name + '.pts')
87print 'create onshore'
88G2 = Geospatial_data(file_name = project.on_offshore_dem_name + '.pts')
89print 'add'
90G = G1 + G2 + G3
91print 'export points'
92G.export_points_file(project.combined_dem_name + '.pts')
93G.export_points_file(project.combined_dem_name + '.xya')
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 anuga.pmesh.mesh_interface import create_mesh_from_regions
102remainder_res = 500000
103local_res = 25000
104newcastle_res = 5000
105coast_res = 500
106interior_regions = [[project.poly_newcastle1, local_res],
107                    [project.poly_newcastle2, newcastle_res],
108                    [project.poly_newcastle3, coast_res]]
109
110from caching import cache
111_ = cache(create_mesh_from_regions,
112          project.polyAll,
113           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
114                              'e3': [3], 'e4':[4], 'e5': [5],
115                              'e6': [6]},
116           'maximum_triangle_area': remainder_res,
117           'filename': meshname,
118           'interior_regions': interior_regions},
119          verbose = True, evaluate=False)
120print 'created mesh'
121
122#-------------------------------------------------------------------------------                                 
123# Setup computational domain
124#-------------------------------------------------------------------------------                                 
125domain = Domain(meshname, use_cache = True, verbose = True)
126
127print 'Number of triangles = ', len(domain)
128print 'The extent is ', domain.get_extent()
129print domain.statistics()
130
131domain.set_name(project.basename)
132domain.set_datadir(project.outputtimedir)
133domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
134domain.set_minimum_storable_height(0.01)
135
136#-------------------------------------------------------------------------------                                 
137# Setup initial conditions
138#-------------------------------------------------------------------------------
139
140tide = 0.0
141domain.set_quantity('stage', tide)
142domain.set_quantity('friction', 0.0) 
143domain.set_quantity('elevation', 
144                    filename = project.combined_dem_name + '.pts',
145                    use_cache = True,
146                    verbose = True,
147                    alpha = 0.1
148                    )
149
150#-------------------------------------------------------------------------------
151# Set up scenario (tsunami_source is a callable object used with set_quantity)
152#-------------------------------------------------------------------------------
153from smf import slide_tsunami
154
155tsunami_source = slide_tsunami(length=30000.0,
156                               depth=400.0,
157                               slope=6.0,
158                               thickness=176.0, 
159                               radius=3330,
160                               dphi=0.23,
161                               x0=project.slump_origin[0], 
162                               y0=project.slump_origin[1], 
163                               alpha=0.0, 
164                               domain=domain)
165
166#-------------------------------------------------------------------------------                                 
167# Setup boundary conditions
168#-------------------------------------------------------------------------------
169print 'Available boundary tags', domain.get_boundary_tags()
170
171Br = Reflective_boundary(domain)
172Bd = Dirichlet_boundary([tide,0,0])
173
174domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd,
175                      'e5': Bd,  'e6': Bd} )
176
177
178#-------------------------------------------------------------------------------                                 
179# Evolve system through time
180#-------------------------------------------------------------------------------
181import time
182t0 = time.time()
183
184for t in domain.evolve(yieldstep = 30, finaltime = 480): 
185    domain.write_time()
186    domain.write_boundary_statistics(tags = 'e14')
187    stagestep = domain.get_quantity('stage') 
188
189    if allclose(t, 30):
190        slide = Quantity(domain)
191        slide.set_values(tsunami_source)
192        domain.set_quantity('stage', slide + stagestep)
193   
194print 'That took %.2f seconds' %(time.time()-t0)
195
196print 'finished'
Note: See TracBrowser for help on using the repository browser.