source: anuga_work/production/australia_ph2/sydney/run_model.py @ 6307

Last change on this file since 6307 was 6307, checked in by jgriffin, 15 years ago
File size: 5.8 KB
RevLine 
[6298]1"""Run a tsunami inundation scenario for sydney region, NSW, Australia.
[6287]2
3The scenario is defined by a triangular mesh created from project.polygon, the
[6298]4elevation data is compiled into a pts file through build_sydney.py and a
[6287]5simulated tsunami is generated through an sts file from build_boundary.py.
6
7Input: sts file (build_boundary.py for respective event)
[6298]8       pts file (build_sydney.py)
[6287]9       information from project file
10Outputs: sww file stored in project.output_run_time_dir
11The export_results_all.py and get_timeseries.py is reliant
12on the outputs of this script
13
14Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006
15Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008
16"""
17
18#------------------------------------------------------------------------------
19# Import necessary modules
20#------------------------------------------------------------------------------
21
22# Standard modules
23import os
24import time
25
26# Related major packages
27from anuga.interface import create_domain_from_regions
28from anuga.interface import Transmissive_stage_zero_momentum_boundary
29from anuga.interface import Dirichlet_boundary
30from anuga.interface import Reflective_boundary
31from anuga.interface import Field_boundary
32from anuga.interface import create_sts_boundary
33from anuga.interface import csv2building_polygons
[6298]34from file_length import file_length
[6287]35
36from anuga.shallow_water.data_manager import start_screen_catcher
37from anuga.shallow_water.data_manager import copy_code_files
38from anuga.utilities.polygon import read_polygon, Polygon_function
39   
40# Application specific imports
41import project  # Definition of file names and polygons
42
43
44#------------------------------------------------------------------------------
45# Copy scripts to time stamped output directory and capture screen
46# output to file. Copy script must be before screen_catcher.
47#------------------------------------------------------------------------------
48copy_code_files(project.output_run, __file__, 
49                os.path.dirname(project.__file__)+os.sep+\
50                project.__name__+'.py' )
51start_screen_catcher(project.output_run, 0, 1)
52
53
54#------------------------------------------------------------------------------
55# Create the computational domain based on overall clipping polygon with
56# a tagged boundary and interior regions defined in project.py along with
57# resolutions (maximal area of per triangle) for each polygon.
58#------------------------------------------------------------------------------
59print 'Create computational domain'
60
61# Read in boundary from ordered sts file
62event_sts = create_sts_boundary(project.event_sts)
63
64# Reading the landward defined points, this incorporates the original clipping
65# polygon minus the 100m contour.
66landward_boundary = read_polygon(project.landward_boundary)
67
68# Combine sts polyline with landward points
69bounding_polygon_sts = event_sts + landward_boundary
70
71# Number of boundary segments
72N = len(event_sts)-1
[6298]73# Number of landward_boundary points
74M = file_length(project.landward_boundary)
[6287]75
76# Boundary tags refer to project.landward_boundary
77# 4 points equals 5 segments start at N
[6298]78boundary_tags={'back': range(N+1,N+M),
79               'side': [N,N+M],
[6287]80               'ocean': range(N)}
81
82# Build mesh and domain
83domain = create_domain_from_regions(bounding_polygon_sts,
84                                    boundary_tags=boundary_tags,
[6298]85                                    maximum_triangle_area=project.bounding_maxarea,
[6287]86                                    interior_regions=project.interior_regions,
87                                    mesh_filename=project.meshes,
88                                    use_cache=True,
89                                    verbose=True)
90print domain.statistics()
91
92domain.set_name(project.scenario_name)
93domain.set_datadir(project.output_run) 
94domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
95
96
97#------------------------------------------------------------------------------
98# Setup initial conditions
99#------------------------------------------------------------------------------
100print 'Setup initial conditions'
101
102# Set the initial stage in the offcoast region only
103IC = Polygon_function(project.land_initial_conditions,
104                      default=project.tide,
105                      geo_reference=domain.geo_reference)
106domain.set_quantity('stage', IC, use_cache=True, verbose=True)
107domain.set_quantity('friction', project.friction) 
108domain.set_quantity('elevation', 
109                    filename=project.combined_elevation+'.pts',
[6307]110                    use_cache=False,
[6287]111                    verbose=True,
112                    alpha=project.alpha)
113
114
115#------------------------------------------------------------------------------
116# Setup boundary conditions
117#------------------------------------------------------------------------------
118print 'Set boundary - available tags:', domain.get_boundary_tags()
119
120Br = Reflective_boundary(domain)
121Bt = Transmissive_stage_zero_momentum_boundary(domain)
122Bf = Field_boundary(project.event_sts+'.sts',
123                    domain, mean_stage=project.tide,
124                    time_thinning=1,
125                    default_boundary=Bd,
126                    boundary_polygon=bounding_polygon_sts,                   
127                    use_cache=True,
128                    verbose=True)
129
130domain.set_boundary({'back': Br,
131                     'side': Bt,
132                     'ocean': Bf}) 
133
134
135#------------------------------------------------------------------------------
136# Evolve system through time
137#------------------------------------------------------------------------------
138t0 = time.time()
139for t in domain.evolve(yieldstep=project.yieldstep, 
140                       finaltime=project.finaltime,
141                       skip_initial_step=False): 
142    print domain.timestepping_statistics()
143    print domain.boundary_statistics(tags='ocean')
144
145print 'Simulation took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.