source: anuga_work/production/bunbury_storm_surge_2009/run_model.py @ 7657

Last change on this file since 7657 was 7657, checked in by fountain, 14 years ago

updates to bunbury model to refine elevation around leschenault inlet and estuary and include storm gates

File size: 7.6 KB
Line 
1"""Run a storm surge inundation scenario for Mandurah, WA, Australia using
2input from the third party consultant GEMS.
3
4The scenario is defined by a triangular mesh created from project.polygon, the
5elevation data is compiled into a pts file through build_elevation.py and a
6simulated storm surge is generated from GEMS data.
7
8Input: sts files (csv format time series at ocean boundary for respective event)
9       pts file (build_elevation.py)
10       information from project file
11Outputs: sww file stored in project.output_run_time_dir
12The export_results_all.py and get_timeseries.py is reliant
13on the outputs of this script
14
15Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006
16Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008
17Nariman Habili, Leharne Fountain - 2009
18"""
19
20#------------------------------------------------------------------------------
21# Import necessary modules
22#------------------------------------------------------------------------------
23
24# Standard modules
25import os
26##import os.path
27import time
28##from time import localtime, strftime, gmtime
29
30# Related major packages
31##from Scientific.IO.NetCDF import NetCDFFile
32import numpy as num
33
34from anuga.interface import create_domain_from_regions
35from anuga.interface import Transmissive_stage_zero_momentum_boundary
36from anuga.interface import Dirichlet_boundary
37from anuga.interface import Reflective_boundary
38from anuga.interface import Field_boundary
39from anuga.interface import create_sts_boundary
40##from anuga.interface import csv2building_polygons
41from file_length import file_length
42
43from anuga.shallow_water.data_manager import start_screen_catcher
44from anuga.shallow_water.data_manager import copy_code_files
45##from anuga.shallow_water.data_manager import urs2sts
46from anuga.utilities.polygon import read_polygon, Polygon_function
47
48# Application specific imports
49from setup_model import project
50##import build_urs_boundary as bub
51
52# Imports for floodgates
53from project import floodgate_boundary, floodgate_resolution
54
55
56#-------------------------------------------------------------------------------
57# Copy scripts to time stamped output directory and capture screen
58# output to file. Copy script must be before screen_catcher
59#-------------------------------------------------------------------------------
60
61copy_code_files(project.output_run, __file__, 
62                os.path.join(os.path.dirname(project.__file__),
63                             project.__name__+'.py'))
64start_screen_catcher(project.output_run, 0, 1)
65
66#-------------------------------------------------------------------------------
67# Create the computational domain based on overall clipping polygon with
68# a tagged boundary and interior regions defined in project.py along with
69# resolutions (maximal area of per triangle) for each polygon
70#-------------------------------------------------------------------------------
71
72print 'Create computational domain'
73
74### Create the STS file
75##print 'project.mux_data_folder=%s' % project.mux_data_folder
76##if not os.path.exists(project.event_sts + '.sts'):
77##    bub.build_urs_boundary(project.mux_input_filename, project.event_sts)
78
79### Read in boundary from ordered sts file
80##gems_boundary = create_sts_boundary(project.event_sts)
81
82#reading the GEMS boundary points
83gems_boundary = read_polygon(project.gems_order)
84##gems_boundary = create_sts_boundary(project.event_sts)
85
86# Reading the landward defined points of the original bounding polygon
87landward_boundary = read_polygon(project.landward_boundary)
88
89# Combine GEMS input boundary polyline with landward points
90bounding_polygon_gems = gems_boundary + landward_boundary
91
92# Number of boundary segments
93num_ocean_segments = len(gems_boundary) - 1
94# Number of landward_boundary points
95num_land_points = file_length(project.landward_boundary)
96
97# Boundary tags refer to project.landward_boundary
98# 4 points equals 5 segments start at N
99boundary_tags={'back': range(num_ocean_segments+1,
100                             num_ocean_segments+num_land_points),
101               'side': [num_ocean_segments,
102                        num_ocean_segments+num_land_points],
103               'ocean': range(num_ocean_segments)}
104
105# Build mesh and domain
106domain = create_domain_from_regions(bounding_polygon_gems,
107                                    boundary_tags=boundary_tags,
108                                    maximum_triangle_area=project.bounding_maxarea,
109                                    interior_regions=project.interior_regions,
110                                    mesh_filename=project.meshes,
111                                    use_cache=False, #usually true but changed for testing
112                                    verbose=True)
113print domain.statistics()
114
115domain.set_starttime(project.starttime)
116domain.set_name(project.scenario_name)
117domain.set_datadir(project.output_run) 
118domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
119
120#-------------------------------------------------------------------------------
121# Setup initial conditions
122#-------------------------------------------------------------------------------
123
124print 'Setup initial conditions'
125
126# Set the initial stage in the offcoast region only
127if project.land_initial_conditions:
128    IC = Polygon_function(project.land_initial_conditions,
129                          default=project.tide,
130                          geo_reference=domain.geo_reference)
131else:
132    IC = 0
133domain.set_quantity('stage', IC, use_cache=True, verbose=True)
134domain.set_quantity('friction', project.friction) 
135domain.set_quantity('elevation', 
136                    filename=project.combined_elevation+'.pts',
137                    use_cache=True,
138                    verbose=True,
139                    alpha=project.alpha)
140
141#--------------------------
142# Add storm surge gate
143#--------------------------
144storm_gate_polygon, storm_gate_height = csv2building_polygons(project.storm_gate, floor_height=4)
145
146gate = []
147for key in storm_gate_polygon:
148    poly = storm_gate_polygon[key]
149    elev = storm_gate_height[key]
150    gate.append((poly, elev))
151
152domain.add_quantity('elevation', 
153                    Polygon_function(gate, geo_reference=domain.geo_reference),
154                    use_cache=True,
155                    verbose=True)
156
157#-------------------------------------------------------------------------------
158# Setup boundary conditions
159#-------------------------------------------------------------------------------
160
161print 'Set boundary - available tags:', domain.get_boundary_tags()
162
163Br = Reflective_boundary(domain)
164Bt = Transmissive_stage_zero_momentum_boundary(domain)
165Bd = Dirichlet_boundary([project.tide, 0, 0])
166Bf = Field_boundary(project.event_sts+'.sts',
167                    domain, mean_stage=project.tide,
168                    time_thinning=1,
169                    default_boundary=Dirichlet_boundary([0, 0, 0]),
170                    boundary_polygon=bounding_polygon_gems,                   
171                    use_cache=True,
172                    verbose=True)
173
174domain.set_boundary({'back': Br,
175                     'side': Bt,
176                     'ocean': Bf}) 
177
178#-------------------------------------------------------------------------------
179# Evolve system through time
180#-------------------------------------------------------------------------------
181
182t0 = time.time()
183for t in domain.evolve(yieldstep=project.yieldstep, 
184                       finaltime=project.finaltime,
185                       skip_initial_step=False): 
186    print domain.timestepping_statistics()
187    print domain.boundary_statistics(tags='ocean')
188
189
190print 'Simulation took %.2f seconds' % (time.time()-t0)
191
Note: See TracBrowser for help on using the repository browser.