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

Last change on this file since 6701 was 6701, checked in by rwilson, 15 years ago

Put 'do not evolve' option code in.

File size: 9.3 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 sys
24import os
25import os.path
26import time
27from time import localtime, strftime, gmtime
28
29# Related major packages
30from Scientific.IO.NetCDF import NetCDFFile
31import Numeric as num
32
33from anuga.interface import create_domain_from_regions
34from anuga.interface import Transmissive_stage_zero_momentum_boundary
35from anuga.interface import Dirichlet_boundary
36from anuga.interface import Reflective_boundary
37from anuga.interface import Field_boundary
38from anuga.interface import create_sts_boundary
39from anuga.interface import csv2building_polygons
40from file_length import file_length
41
42from anuga.shallow_water.data_manager import start_screen_catcher
43from anuga.shallow_water.data_manager import copy_code_files
44from anuga.shallow_water.data_manager import urs2sts
45from anuga.utilities.polygon import read_polygon, Polygon_function
46from anuga.caching import cache
47
48# Application specific imports
49from setup_model import project
50import build_urs_boundary as bub
51
52
53#-------------------------------------------------------------------------------
54# Examine args - look for switches to control execution
55#
56# We expect:
57#   -n --no-evolve  Stop after writing cached data - don't enter evolve loop
58#-------------------------------------------------------------------------------
59
60##
61# @brief Return a usage string.
62def usage():
63    '''Return a 'usage' string.'''
64
65    return '''Usage: run_model.py [<options>]
66where <options> is one or more of:
67                  --help       print this help
68                  --no-evolve  quit after creating cache data
69'''
70
71
72try:
73    try:
74        opts, args = getopt.getopt(sys.argv[1:], 'hn', ['help', 'no-evolve'])
75    except getopt.error, msg:
76        raise RuntimeError, Usage(msg)
77except Usage, err:
78    print >>sys.stderr, err.msg
79    print >>sys.stderr, "for help use --help"
80    sys.exit(2)
81
82# process options
83DoEvolve = True
84for opt, arg in opts:
85    if opt in ('-h', '--help'):
86        print usage()
87        sys.exit(0)
88    elif opt in ('-n', '--no-evolve'):
89        DoEvolve = False
90
91#-------------------------------------------------------------------------------
92# Copy scripts to time stamped output directory and capture screen
93# output to file. Copy script must be before screen_catcher
94#-------------------------------------------------------------------------------
95
96copy_code_files(project.output_run, __file__, 
97                os.path.join(os.path.dirname(project.__file__),
98                             project.__name__+'.py'))
99start_screen_catcher(project.output_run, 0, 1)
100
101#-------------------------------------------------------------------------------
102# Create the computational domain based on overall clipping polygon with
103# a tagged boundary and interior regions defined in project.py along with
104# resolutions (maximal area of per triangle) for each polygon
105#-------------------------------------------------------------------------------
106
107print 'Create computational domain'
108
109# Create the STS file
110# FIXME (Ole): This is deadly dangerous if buildcode changes (as was the case 24th March 2009)
111# We need to use caching instead!
112
113print 'project.mux_data_folder=%s' % project.mux_data_folder
114if not os.path.exists(project.event_sts + '.sts'):
115    bub.build_urs_boundary(project.mux_input_filename, project.event_sts)
116
117# Read in boundary from ordered sts file
118event_sts = create_sts_boundary(project.event_sts)
119
120# Reading the landward defined points, this incorporates the original clipping
121# polygon minus the 100m contour
122landward_boundary = read_polygon(project.landward_boundary)
123
124# Combine sts polyline with landward points
125bounding_polygon_sts = event_sts + landward_boundary
126
127# Number of boundary segments
128num_ocean_segments = len(event_sts) - 1
129# Number of landward_boundary points
130num_land_points = file_length(project.landward_boundary)
131
132# Boundary tags refer to project.landward_boundary
133# 4 points equals 5 segments start at N
134boundary_tags={'back': range(num_ocean_segments+1,
135                             num_ocean_segments+num_land_points),
136               'side': [num_ocean_segments,
137                        num_ocean_segments+num_land_points],
138               'ocean': range(num_ocean_segments)}
139
140# Build mesh and domain
141domain = create_domain_from_regions(bounding_polygon_sts,
142                                    boundary_tags=boundary_tags,
143                                    maximum_triangle_area=project.bounding_maxarea,
144                                    interior_regions=project.interior_regions,
145                                    mesh_filename=project.meshes,
146                                    use_cache=True,
147                                    verbose=True)
148print domain.statistics()
149
150# FIXME(Ole): How can we make this more automatic?
151domain.geo_reference.zone = project.zone
152
153
154domain.set_name(project.scenario_name)
155domain.set_datadir(project.output_run) 
156domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
157
158#-------------------------------------------------------------------------------
159# Setup initial conditions
160#-------------------------------------------------------------------------------
161
162print 'Setup initial conditions'
163
164# Set the initial stage in the offcoast region only
165if project.land_initial_conditions:
166    IC = Polygon_function(project.land_initial_conditions,
167                          default=project.tide,
168                          geo_reference=domain.geo_reference)
169else:
170    IC = 0
171domain.set_quantity('stage', IC, use_cache=True, verbose=True)
172domain.set_quantity('friction', project.friction) 
173domain.set_quantity('elevation', 
174                    filename=project.combined_elevation+'.pts',
175                    use_cache=True,
176                    verbose=True,
177                    alpha=project.alpha)
178
179if project.use_buildings:
180    # Add buildings from file
181    print 'Reading building polygons'   
182    building_polygons, building_heights = csv2building_polygons(project.building_polygon)
183    #clipping_polygons=project.building_area_polygons)
184
185    def create_polygon_function(building_polygons, geo_reference=None):
186        L = []
187        for i, key in enumerate(building_polygons):
188            if i%100==0: print i
189            poly = building_polygons[key]
190            elev = building_heights[key]
191            L.append((poly, elev))
192           
193            buildings = Polygon_function(L, default=0.0,
194                                         geo_reference=geo_reference)
195        return buildings
196
197    print 'Creating %d building polygons' % len(building_polygons)
198    buildings = cache(create_polygon_function,
199                      building_polygons,
200                      {'geo_reference': domain.geo_reference},
201                      verbose=True)
202
203    print 'Adding buildings'
204    domain.add_quantity('elevation',
205                        buildings,
206                        use_cache=True,
207                        verbose=True)
208
209
210#-------------------------------------------------------------------------------
211# Setup boundary conditions
212#-------------------------------------------------------------------------------
213
214print 'Set boundary - available tags:', domain.get_boundary_tags()
215
216Br = Reflective_boundary(domain)
217Bs = Transmissive_stage_zero_momentum_boundary(domain)
218Bf = Field_boundary(project.event_sts+'.sts',
219                    domain,
220                    mean_stage=project.tide,
221                    time_thinning=1,
222                    default_boundary=Dirichlet_boundary([0, 0, 0]),
223                    boundary_polygon=bounding_polygon_sts,                   
224                    use_cache=True,
225                    verbose=True)
226
227domain.set_boundary({'back': Br,
228                     'side': Bs,
229                     'ocean': Bf}) 
230
231#-------------------------------------------------------------------------------
232# Evolve system through time
233#-------------------------------------------------------------------------------
234
235if DoEvolve:
236    t0 = time.time()
237
238    # Skip over the first 6000 seconds
239    for t in domain.evolve(yieldstep=2000,
240                           finaltime=6000):
241        print domain.timestepping_statistics()
242        print domain.boundary_statistics(tags='ocean')
243
244    # Start detailed model
245    for t in domain.evolve(yieldstep=project.yieldstep,
246                           finaltime=project.finaltime,
247                           skip_initial_step=True):
248        print domain.timestepping_statistics()
249        print domain.boundary_statistics(tags='ocean')
250else:
251    print 'Model not evolved, as requested'
252
253print 'Simulation took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.