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

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

Bunbury storm surge modelling

File size: 8.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_name(project.scenario_name)
116domain.set_datadir(project.output_run) 
117domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
118
119#-------------------------------------------------------------------------------
120# Setup initial conditions
121#-------------------------------------------------------------------------------
122
123print 'Setup initial conditions'
124
125# Set the initial stage in the offcoast region only
126if project.land_initial_conditions:
127    IC = Polygon_function(project.land_initial_conditions,
128                          default=project.tide,
129                          geo_reference=domain.geo_reference)
130else:
131    IC = 0
132domain.set_quantity('stage', IC, use_cache=True, verbose=True)
133domain.set_quantity('friction', project.friction) 
134domain.set_quantity('elevation', 
135                    filename=project.combined_elevation+'.pts',
136                    use_cache=True,
137                    verbose=True,
138                    alpha=project.alpha)
139
140#-------------------------------------------------------------------------------
141# Setup boundary conditions
142#-------------------------------------------------------------------------------
143
144print 'Set boundary - available tags:', domain.get_boundary_tags()
145
146Br = Reflective_boundary(domain)
147Bt = Transmissive_stage_zero_momentum_boundary(domain)
148Bd = Dirichlet_boundary([project.tide, 0, 0])
149Bf = Field_boundary(project.event_sts+'.sts',
150                    domain, mean_stage=project.tide,
151                    time_thinning=1,
152                    default_boundary=Dirichlet_boundary([0, 0, 0]),
153                    boundary_polygon=bounding_polygon_gems,                   
154                    use_cache=True,
155                    verbose=True)
156
157domain.set_boundary({'back': Br,
158                     'side': Bt,
159                     'ocean': Bf}) 
160
161#-----------------------------------------------------------------------------
162#  Setup dynamic topography
163#-----------------------------------------------------------------------------
164from project import floodgate_height, floodgate_thickness,floodgate_LHS_x,floodgate_LHS_y,floodgate_RHS_x,floodgate_RHS_y,floodgate_start,floodgate_duration
165
166floodgate_end = floodgate_start + floodgate_duration #the time in seconds at which the floodgates are fully closed
167
168def floodgate(x,y):
169    z = 0*x
170    N = len(x)
171    for i in range(N):
172        ymin = ((floodgate_RHS_y - floodgate_LHS_y)/(floodgate_RHS_x - floodgate_LHS_x))*(x[i]-floodgate_LHS_x)+ floodgate_LHS_y - 0.5*floodgate_thickness
173        ymax = ((floodgate_RHS_y - floodgate_LHS_y)/(floodgate_RHS_x - floodgate_LHS_x))*(x[i]-floodgate_LHS_x)+ floodgate_LHS_y + 0.5*floodgate_thickness
174        xmin = floodgate_LHS_x
175        xmax = floodgate_RHS_x
176        if ymin < y[i] < ymax and xmin < x[i] < xmax:
177            z[i] += (floodgate_height/((floodgate_duration)/project.yieldstep))
178
179    return z
180   
181#-------------------------------------------------------------------------------
182# Evolve system through time
183#-------------------------------------------------------------------------------
184
185raising_floodgates = False
186lowering_floodgates = False
187
188t0 = time.time()
189for t in domain.evolve(yieldstep=project.yieldstep, 
190                       finaltime=project.finaltime,
191                       skip_initial_step=False): 
192    print domain.timestepping_statistics()
193    print domain.boundary_statistics(tags='ocean')
194
195#raise floodgate at a given time interval (to make floodgate lower again repeat this code below
196#for a later time interval and change rising to false and falling to true
197    if floodgate_start < t < floodgate_end:
198        rising = True
199        falling = False
200
201    if rising:
202        domain.add_quantity('elevation',floodgate)
203
204    if falling:
205        domain.add_quantity('elevation', lambda x,y: -1*floodgate(x,y))
206
207print 'Simulation took %.2f seconds' % (time.time()-t0)
208
Note: See TracBrowser for help on using the repository browser.