source: anuga_work/production/patong/new_version/run_model.py @ 6418

Last change on this file since 6418 was 6418, checked in by kristy, 14 years ago

new scripts, running with large and small bounding_polygon

File size: 7.9 KB
Line 
1"""Run a tsunami inundation scenario for Busselton, WA, Australia.
2
3The scenario is defined by a triangular mesh created from project.polygon, the
4elevation data is compiled into a pts file through build_elevation.py and a
5simulated tsunami is generated through an sts file from build_boundary.py.
6
7Input: sts file (build_boundary.py for respective event)
8       pts file (build_elevation.py)
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 os.path
25import time
26from time import localtime, strftime, gmtime
27
28# Related major packages
29from Scientific.IO.NetCDF import NetCDFFile
30import Numeric as num
31
32from anuga.interface import create_domain_from_regions
33from anuga.interface import Transmissive_stage_zero_momentum_boundary
34from anuga.interface import Dirichlet_boundary
35from anuga.interface import Reflective_boundary
36from anuga.interface import Field_boundary
37from anuga.interface import create_sts_boundary
38from anuga.interface import csv2building_polygons
39from file_length import file_length
40
41from anuga.shallow_water.data_manager import start_screen_catcher
42from anuga.shallow_water.data_manager import copy_code_files
43from anuga.shallow_water.data_manager import urs2sts
44from anuga.utilities.polygon import read_polygon, Polygon_function
45
46# Application specific imports
47from setup_model import project
48import build_urs_boundary as bub
49
50#-------------------------------------------------------------------------------
51# Copy scripts to time stamped output directory and capture screen
52# output to file. Copy script must be before screen_catcher
53#-------------------------------------------------------------------------------
54
55copy_code_files(project.output_run, __file__, 
56                os.path.join(os.path.dirname(project.__file__),
57                             project.__name__+'.py'))
58start_screen_catcher(project.output_run, 0, 1)
59
60#-------------------------------------------------------------------------------
61# Create the computational domain based on overall clipping polygon with
62# a tagged boundary and interior regions defined in project.py along with
63# resolutions (maximal area of per triangle) for each polygon
64#-------------------------------------------------------------------------------
65
66print 'Create computational domain'
67
68# Create the STS file
69print 'project.mux_data_folder=%s' % project.mux_data_folder
70if not os.path.exists(project.event_sts + '.sts'):
71    bub.build_urs_boundary(project.mux_input_filename, project.event_sts)
72
73# Read in boundary from ordered sts file
74event_sts = create_sts_boundary(project.event_sts)
75
76# Reading the landward defined points, this incorporates the original clipping
77# polygon minus the 100m contour
78landward_boundary = read_polygon(project.landward_boundary)
79
80# Combine sts polyline with landward points
81bounding_polygon_sts = event_sts + landward_boundary
82
83# Number of boundary segments
84num_ocean_segments = len(event_sts) - 1
85# Number of landward_boundary points
86num_land_points = file_length(project.landward_boundary)
87
88# Boundary tags refer to project.landward_boundary
89# 4 points equals 5 segments start at N
90boundary_tags={'back': range(num_ocean_segments+1,
91                             num_ocean_segments+num_land_points),
92               'side': [num_ocean_segments,
93                        num_ocean_segments+num_land_points],
94               'ocean': range(num_ocean_segments)}
95
96# Build mesh and domain
97domain = create_domain_from_regions(bounding_polygon_sts,
98                                    boundary_tags=boundary_tags,
99                                    maximum_triangle_area=project.bounding_maxarea,
100                                    interior_regions=project.interior_regions,
101                                    mesh_filename=project.meshes,
102                                    use_cache=True,
103                                    verbose=True)
104print domain.statistics()
105
106domain.set_name(project.scenario_name)
107domain.set_datadir(project.output_run) 
108domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
109
110#-------------------------------------------------------------------------------
111# Setup initial conditions
112#-------------------------------------------------------------------------------
113
114print 'Setup initial conditions'
115
116# Set the initial stage in the offcoast region only
117if project.land_initial_conditions:
118    IC = Polygon_function(project.land_initial_conditions,
119                          default=project.tide,
120                          geo_reference=domain.geo_reference)
121else:
122    IC = 0
123domain.set_quantity('stage', IC, use_cache=True, verbose=True)
124domain.set_quantity('friction', project.friction) 
125domain.set_quantity('elevation', 
126                    filename=project.combined_elevation+'.pts',
127                    use_cache=True,
128                    verbose=True,
129                    alpha=project.alpha)
130
131if project.use_buildings:
132    # Add buildings from file
133    print 'Reading building polygons'   
134    building_polygons, building_heights = csv2building_polygons(project.building_polygon)
135    #clipping_polygons=project.building_area_polygons)
136
137    print 'Creating %d building polygons' % len(building_polygons)
138    def create_polygon_function(building_polygons, geo_reference=None):
139        L = []
140        for i, key in enumerate(building_polygons):
141            if i%100==0: print i
142            poly = building_polygons[key]
143            elev = building_heights[key]
144            L.append((poly, elev))
145           
146            buildings = Polygon_function(L, default=0.0,
147                                         geo_reference=geo_reference)
148        return buildings
149
150    buildings = cache(create_polygon_function,
151                      building_polygons,
152                      {'geo_reference': domain.geo_reference},
153                      verbose=True)
154
155    print 'Adding buildings'
156    domain.add_quantity('elevation',
157                        buildings,
158                        use_cache=True,
159                        verbose=True)
160
161
162#-------------------------------------------------------------------------------
163# Setup boundary conditions
164#-------------------------------------------------------------------------------
165
166print 'Set boundary - available tags:', domain.get_boundary_tags()
167
168Br = Reflective_boundary(domain)
169Bt = Transmissive_stage_zero_momentum_boundary(domain)
170Bd = Dirichlet_boundary([project.tide, 0, 0])
171Bf = Field_boundary(project.event_sts+'.sts',
172                    domain, mean_stage=project.tide,
173                    time_thinning=1,
174                    default_boundary=Bd,
175                    boundary_polygon=bounding_polygon_sts,                   
176                    use_cache=True,
177                    verbose=True)
178
179domain.set_boundary({'back': Br,
180                     'side': Bd,
181                     'ocean': Bf}) 
182
183#-------------------------------------------------------------------------------
184# Evolve system through time
185#-------------------------------------------------------------------------------
186
187t0 = time.time()
188
189# Skip over the first 6000 seconds
190for t in domain.evolve(yieldstep=2000,
191                       finaltime=6000):
192    print domain.timestepping_statistics()
193    print domain.boundary_statistics(tags='ocean')
194
195# Start detailed model
196for t in domain.evolve(yieldstep=project.yieldstep,
197                       finaltime=project.finaltime,
198                       skip_initial_step=True):
199    print domain.timestepping_statistics()
200    print domain.boundary_statistics(tags='ocean')
201   
202print 'Simulation took %.2f seconds' %(time.time()-t0)
203     
Note: See TracBrowser for help on using the repository browser.