source: anuga_core/documentation/user_manual/examples/runsydney.py @ 3944

Last change on this file since 3944 was 3944, checked in by sexton, 16 years ago

(i) update sww2timeseries so can handle gauges which don't fall within sww domain (ii) script for generating report for dampier based on sww files from parallel setup

File size: 6.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.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
10"""
11
12
13#------------------------------------------------------------------------------
14# Import necessary modules
15#------------------------------------------------------------------------------
16
17# Standard modules
18import os
19import time
20
21# Related major packages
22from anuga.shallow_water import Domain, Dirichlet_boundary
23from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
24from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files
25from anuga.pmesh.mesh_interface import create_mesh_from_regions
26
27# Application specific imports
28import project                           # Define file names and polygons
29from anuga.shallow_water.smf import slump_tsunami # Function for submarine mudslide
30
31
32
33#------------------------------------------------------------------------------
34# Prepare bathymetric and topographic data
35#------------------------------------------------------------------------------
36
37# filenames
38#coarsedemname = project.coarsedemname + '.pts'
39#finedemname = project.finedemname + '.pts'
40demname = project.demname
41meshname = project.meshname+'.msh'
42
43convert_dem_from_ascii2netcdf(demname, use_cache=True, verbose=True)
44
45# creates pts file
46dem2pts(demname, use_cache=True, verbose=True)
47
48# combining the coarse and fine data
49#combine_rectangular_points_files(finedemname,
50#                                 coarsedemname,
51#                                 combineddemname)
52
53
54
55#------------------------------------------------------------------------------
56# Setup computational domain
57#------------------------------------------------------------------------------
58
59# Interior regions and mesh resolutions
60interior_res = 5000
61shallow_res = 15000
62##interior_regions = [[project.northern_polygon, interior_res],
63##                    [project.harbour_polygon, interior_res],
64##                    [project.southern_polygon, interior_res],
65##                    [project.manly_polygon, interior_res],
66##                    [project.top_polygon, interior_res]]
67interior_regions = [[project.coastal_polygon, interior_res],
68                    [project.shallow_polygon, shallow_res]]
69
70create_mesh_from_regions(project.demopoly,
71                         boundary_tags= {'oceannorth': [0], 
72                                         'onshorenorth1': [1],
73                                         'onshorenorthwest1': [2],
74                                         'onshorenorthwest2': [3],
75                                         'onshorenorth2': [4],
76                                         'onshorewest': [5],
77                                         'onshoresouth': [6],
78                                         'oceansouth': [7],
79                                         'oceaneast': [8]},
80                         maximum_triangle_area=100000,
81                         filename=meshname,           
82                         interior_regions=interior_regions)
83
84
85#Create shallow water domain
86
87domain = Domain(meshname,
88                use_cache = True,
89                verbose = True)                                 
90
91
92print 'Number of triangles = ', len(domain)
93print 'The extent is ', domain.get_extent()
94
95domain.set_name(project.basename)
96domain.set_datadir(project.outputdir)
97domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
98
99
100#------------------------------------------------------------------------------
101# Set up scenario (tsunami_source is a callable object used with set_quantity)
102#------------------------------------------------------------------------------
103
104tsunami_source = slump_tsunami(length=30000.0,
105                               depth=400.0,
106                               slope=6.0,
107                               thickness=176.0, 
108                               radius=3330,
109                               dphi=0.23,
110                               x0=project.slump_origin[0], 
111                               y0=project.slump_origin[1], 
112                               alpha=0.0, 
113                               domain=domain)
114
115
116#------------------------------------------------------------------------------
117# Setup initial conditions
118#------------------------------------------------------------------------------
119
120domain.set_quantity('stage', tsunami_source)
121domain.set_quantity('friction', 0.0) 
122domain.set_quantity('elevation',
123                    filename = demname + '.pts',
124                    use_cache = True,
125                    verbose = True)
126
127#------------------------------------------------------------------------------
128# Setup boundary conditions (all Dirichlet)
129#------------------------------------------------------------------------------
130
131print 'Available boundary tags', domain.get_boundary_tags()
132
133Bd = Dirichlet_boundary([0.0,0.0,0.0])
134domain.set_boundary( {'oceannorth': Bd,
135                      'onshorenorth1': Bd,
136                      'onshorenorthwest1': Bd,
137                      'onshorenorthwest2': Bd,
138                      'onshorenorth2': Bd,
139                      'onshorewest': Bd,
140                      'onshoresouth': Bd,
141                      'oceansouth': Bd,
142                      'oceaneast': Bd} )
143
144#------------------------------------------------------------------------------
145# Evolve system through time
146#------------------------------------------------------------------------------
147
148import time
149t0 = time.time()
150
151for t in domain.evolve(yieldstep = 60, finaltime = 7200):
152    print domain.timestepping_statistics()
153    print domain.boundary_statistics(tags = 'oceaneast')   
154   
155print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.