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

Last change on this file since 3944 was 3944, checked in by sexton, 17 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
RevLine 
[2461]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
[2887]13#------------------------------------------------------------------------------
[2478]14# Import necessary modules
[2887]15#------------------------------------------------------------------------------
[2461]16
17# Standard modules
18import os
19import time
20
21# Related major packages
[3563]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
[3534]25from anuga.pmesh.mesh_interface import create_mesh_from_regions
[2461]26
27# Application specific imports
[2887]28import project                           # Define file names and polygons
[3563]29from anuga.shallow_water.smf import slump_tsunami # Function for submarine mudslide
[2461]30
31
32
[2887]33#------------------------------------------------------------------------------
[2599]34# Prepare bathymetric and topographic data
[2887]35#------------------------------------------------------------------------------
[2478]36
[2461]37# filenames
[3869]38#coarsedemname = project.coarsedemname + '.pts'
39#finedemname = project.finedemname + '.pts'
40demname = project.demname
[2461]41meshname = project.meshname+'.msh'
42
[3869]43convert_dem_from_ascii2netcdf(demname, use_cache=True, verbose=True)
[2461]44
[3869]45# creates pts file
46dem2pts(demname, use_cache=True, verbose=True)
47
[2461]48# combining the coarse and fine data
[3869]49#combine_rectangular_points_files(finedemname,
50#                                 coarsedemname,
51#                                 combineddemname)
[2461]52
53
[2478]54
[2887]55#------------------------------------------------------------------------------
[2478]56# Setup computational domain
[2887]57#------------------------------------------------------------------------------
[2478]58
59# Interior regions and mesh resolutions
[2461]60interior_res = 5000
[3944]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]]
[2461]69
[3136]70create_mesh_from_regions(project.demopoly,
71                         boundary_tags= {'oceannorth': [0], 
[3190]72                                         'onshorenorth1': [1],
73                                         'onshorenorthwest1': [2],
74                                         'onshorenorthwest2': [3],
75                                         'onshorenorth2': [4],
76                                         'onshorewest': [5],
77                                         'onshoresouth': [6],
78                                         'oceansouth': [7],
79                                         'oceaneast': [8]},
[2478]80                         maximum_triangle_area=100000,
81                         filename=meshname,           
82                         interior_regions=interior_regions)
[2461]83
84
[2805]85#Create shallow water domain
[2461]86
[2866]87domain = Domain(meshname,
88                use_cache = True,
89                verbose = True)                                 
[2805]90
[2866]91
[2461]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
[2887]100#------------------------------------------------------------------------------
[2461]101# Set up scenario (tsunami_source is a callable object used with set_quantity)
[2887]102#------------------------------------------------------------------------------
[2461]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
[2887]116#------------------------------------------------------------------------------
[2461]117# Setup initial conditions
[2887]118#------------------------------------------------------------------------------
[2461]119
120domain.set_quantity('stage', tsunami_source)
121domain.set_quantity('friction', 0.0) 
122domain.set_quantity('elevation',
[3869]123                    filename = demname + '.pts',
[2461]124                    use_cache = True,
125                    verbose = True)
126
[2887]127#------------------------------------------------------------------------------
[3136]128# Setup boundary conditions (all Dirichlet)
[2887]129#------------------------------------------------------------------------------
[2461]130
131print 'Available boundary tags', domain.get_boundary_tags()
132
[3136]133Bd = Dirichlet_boundary([0.0,0.0,0.0])
[3190]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} )
[2461]143
[2887]144#------------------------------------------------------------------------------
[2461]145# Evolve system through time
[2887]146#------------------------------------------------------------------------------
[2461]147
148import time
149t0 = time.time()
150
[3190]151for t in domain.evolve(yieldstep = 60, finaltime = 7200):
[2497]152    print domain.timestepping_statistics()
[3136]153    print domain.boundary_statistics(tags = 'oceaneast')   
[2461]154   
155print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.