source: anuga_work/production/mandurah_storm_surge_2009/gems_comparison/run_model_GEMS.py @ 7628

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

comparison with GEMS storm surge model for Mandurah

File size: 7.2 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_GEMS import project_GEMS
50##import build_urs_boundary as bub
51
52
53#-------------------------------------------------------------------------------
54# Copy scripts to time stamped output directory and capture screen
55# output to file. Copy script must be before screen_catcher
56#-------------------------------------------------------------------------------
57
58copy_code_files(project_GEMS.output_run, __file__, 
59                os.path.join(os.path.dirname(project_GEMS.__file__),
60                             project_GEMS.__name__+'.py'))
61start_screen_catcher(project_GEMS.output_run, 0, 1)
62
63#-------------------------------------------------------------------------------
64# Create the computational domain based on overall clipping polygon with
65# a tagged boundary and interior regions defined in project_GEMS.py along with
66# resolutions (maximal area of per triangle) for each polygon
67#-------------------------------------------------------------------------------
68
69print 'Create computational domain'
70
71### Create the STS file
72##print 'project.mux_data_folder=%s' % project.mux_data_folder
73##if not os.path.exists(project.event_sts + '.sts'):
74##    bub.build_urs_boundary(project.mux_input_filename, project.event_sts)
75
76### Read in boundary from ordered sts file
77##gems_boundary = create_sts_boundary(project_GEMS.event_sts)
78
79#reading the GEMS boundary points
80gems_boundary = read_polygon(project_GEMS.gems_order)
81##gems_boundary = create_sts_boundary(project_GEMS.event_sts)
82
83# Reading the landward defined points of the original bounding polygon
84landward_boundary = read_polygon(project_GEMS.landward_boundary)
85
86# Combine GEMS input boundary polyline with landward points
87bounding_polygon_gems = gems_boundary + landward_boundary
88
89# Number of boundary segments
90num_ocean_segments = len(gems_boundary) - 1
91# Number of landward_boundary points
92num_land_points = file_length(project_GEMS.landward_boundary)
93
94# Boundary tags refer to project_GEMS.landward_boundary
95# 4 points equals 5 segments start at N
96boundary_tags={'back': range(num_ocean_segments+1,
97                             num_ocean_segments+num_land_points),
98               'side': [num_ocean_segments,
99                        num_ocean_segments+num_land_points],
100               'ocean': range(num_ocean_segments)}
101
102# Build mesh and domain
103domain = create_domain_from_regions(bounding_polygon_gems,
104                                    boundary_tags=boundary_tags,
105                                    maximum_triangle_area=project_GEMS.bounding_maxarea,
106                                    interior_regions=project_GEMS.interior_regions,
107                                    mesh_filename=project_GEMS.meshes,
108                                    use_cache=False, #usually true but changed for testing
109                                    verbose=True)
110print domain.statistics()
111
112domain.starttime(project_GEMS.starttime)
113domain.set_name(project_GEMS.scenario_name)
114domain.set_datadir(project_GEMS.output_run) 
115domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
116
117#-------------------------------------------------------------------------------
118# Setup initial conditions
119#-------------------------------------------------------------------------------
120
121print 'Setup initial conditions'
122
123# Set the initial stage in the offcoast region only
124if project_GEMS.land_initial_conditions:
125    IC = Polygon_function(project_GEMS.land_initial_conditions,
126                          default=project_GEMS.tide,
127                          geo_reference=domain.geo_reference)
128else:
129    IC = 0
130domain.set_quantity('stage', IC, use_cache=True, verbose=True)
131domain.set_quantity('friction', project_GEMS.friction) 
132domain.set_quantity('elevation', 
133                    filename=project_GEMS.combined_elevation+'.pts',
134                    use_cache=True,
135                    verbose=True,
136                    alpha=project_GEMS.alpha)
137
138#-------------------------------------------------------------------------------
139# Setup boundary conditions
140#-------------------------------------------------------------------------------
141
142print 'Set boundary - available tags:', domain.get_boundary_tags()
143
144Br = Reflective_boundary(domain)
145Bt = Transmissive_stage_zero_momentum_boundary(domain)
146Bd = Dirichlet_boundary([project_GEMS.tide, 0, 0])
147Bf = Field_boundary(project_GEMS.event_sts+'.sts',
148                    domain, mean_stage=project_GEMS.tide,
149                    time_thinning=1,
150                    default_boundary=Dirichlet_boundary([0, 0, 0]),
151                    boundary_polygon=bounding_polygon_gems,                   
152                    use_cache=True,
153                    verbose=True)
154
155domain.set_boundary({'back': Br,
156                     'side': Bt,
157                     'ocean': Bf}) 
158
159#-------------------------------------------------------------------------------
160# Evolve system through time
161#-------------------------------------------------------------------------------
162
163t0 = time.time()
164for t in domain.evolve(yieldstep=project_GEMS.yieldstep, 
165                       finaltime=project_GEMS.finaltime,
166                       skip_initial_step=False): 
167    print domain.timestepping_statistics()
168    print domain.boundary_statistics(tags='ocean')
169
170print 'Simulation took %.2f seconds' % (time.time()-t0)
171
Note: See TracBrowser for help on using the repository browser.