source: anuga_validation/automated_validation_tests/patong_beach_validation/run_model.py @ 7750

Last change on this file since 7750 was 7750, checked in by hudson, 14 years ago

Updated validation suite to work with new API - confirmed that Patong passes.

  • Property svn:executable set to *
File size: 8.6 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
30
31from anuga.interface import create_domain_from_regions
32from anuga.interface import Transmissive_stage_zero_momentum_boundary
33from anuga.interface import Dirichlet_boundary
34from anuga.interface import Reflective_boundary
35from anuga.interface import Field_boundary
36from anuga.interface import create_sts_boundary
37from anuga.interface import csv2building_polygons
38from anuga.utilities.system_tools import file_length
39
40from anuga.utilities.file_utils import copy_code_files
41from anuga.shallow_water.data_manager import urs2sts
42from anuga.utilities.polygon import read_polygon, Polygon_function
43from anuga.caching import cache
44
45import project
46import anuga.utilities.log as log
47
48#-------------------------------------------------------------------------------
49# Copy scripts to time stamped output directory and capture screen
50# output to file. Copy script must be before screen_catcher
51#-------------------------------------------------------------------------------
52
53# Make output dir and set log filename before anything is logged
54os.mkdir(project.output_run)
55# Tell log module to store log file in output dir
56log.log_filename = os.path.join(project.output_run, 'anuga.log')
57log.critical('Log filename: %s' % log.log_filename)
58
59
60output_dirname = os.path.dirname(project.__file__)
61copy_code_files(project.output_run,
62                [__file__, 
63                 os.path.join(output_dirname, project.__name__+'.py'),
64                 os.path.join(output_dirname, 'setup_model.py')],
65                verbose=True
66                )
67
68# Application specific imports
69from setup_model import project
70import build_urs_boundary as bub
71
72
73#-------------------------------------------------------------------------------
74# Create the computational domain based on overall clipping polygon with
75# a tagged boundary and interior regions defined in project.py along with
76# resolutions (maximal area of per triangle) for each polygon
77#-------------------------------------------------------------------------------
78
79log.critical('Create computational domain')
80
81# Create the STS file
82# FIXME (Ole): This is deadly dangerous if buildcode changes (as was the case 24th March 2009)
83# We need to use caching instead!
84
85log.critical( 'project.mux_data_folder=%s' % project.mux_data_folder)
86if not os.path.exists(project.event_sts + '.sts'):
87    bub.build_urs_boundary(project.mux_input_filename, project.event_sts)
88
89# Read in boundary from ordered sts file
90event_sts = create_sts_boundary(project.event_sts)
91
92# Reading the landward defined points, this incorporates the original clipping
93# polygon minus the 100m contour
94landward_boundary = read_polygon(project.landward_boundary)
95
96# Combine sts polyline with landward points
97bounding_polygon_sts = event_sts + landward_boundary
98
99# Number of boundary segments
100num_ocean_segments = len(event_sts) - 1
101# Number of landward_boundary points
102num_land_points = file_length(project.landward_boundary)
103
104# Boundary tags refer to project.landward_boundary
105# 4 points equals 5 segments start at N
106boundary_tags={'back': range(num_ocean_segments+1,
107                             num_ocean_segments+num_land_points),
108               'side': [num_ocean_segments,
109                        num_ocean_segments+num_land_points],
110               'ocean': range(num_ocean_segments)}
111
112# Build mesh and domain
113domain = create_domain_from_regions(bounding_polygon_sts,
114                                    boundary_tags=boundary_tags,
115                                    maximum_triangle_area=project.bounding_maxarea,
116                                    interior_regions=project.interior_regions,
117                                    mesh_filename=project.meshes,
118                                    use_cache=True,
119                                    verbose=True)
120log.critical(domain.statistics())
121
122# FIXME(Ole): How can we make this more automatic?
123domain.geo_reference.zone = project.zone
124
125
126domain.set_name(project.scenario_name)
127domain.set_datadir(project.output_run) 
128domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
129
130#-------------------------------------------------------------------------------
131# Setup initial conditions
132#-------------------------------------------------------------------------------
133
134log.critical('Setup initial conditions')
135
136# Set the initial stage in the offcoast region only
137if project.land_initial_conditions:
138    IC = Polygon_function(project.land_initial_conditions,
139                          default=project.tide,
140                          geo_reference=domain.geo_reference)
141else:
142    IC = 0
143domain.set_quantity('stage', IC, use_cache=True, verbose=True)
144domain.set_quantity('friction', project.friction) 
145domain.set_quantity('elevation', 
146                    filename=project.combined_elevation+'.pts',
147                    use_cache=True,
148                    verbose=True,
149                    alpha=project.alpha)
150
151if project.use_buildings:
152    # Add buildings from file
153    log.critical('Reading building polygons')
154    building_polygons, building_heights = csv2building_polygons(project.building_polygon)
155    #clipping_polygons=project.building_area_polygons)
156
157    log.critical('Creating %d building polygons' % len(building_polygons))
158    def create_polygon_function(building_polygons, geo_reference=None):
159        L = []
160        for i, key in enumerate(building_polygons):
161            if i%100==0: log.critical(i)
162            poly = building_polygons[key]
163            elev = building_heights[key]
164            L.append((poly, elev))
165           
166            buildings = Polygon_function(L, default=0.0,
167                                         geo_reference=geo_reference)
168        return buildings
169
170    log.critical('Creating %d building polygons' % len(building_polygons))
171    buildings = cache(create_polygon_function,
172                      building_polygons,
173                      {'geo_reference': domain.geo_reference},
174                      verbose=True)
175
176    log.critical('Adding buildings')
177    domain.add_quantity('elevation',
178                        buildings,
179                        use_cache=True,
180                        verbose=True)
181
182
183#-------------------------------------------------------------------------------
184# Setup boundary conditions
185#-------------------------------------------------------------------------------
186
187log.critical('Set boundary - available tags:' % domain.get_boundary_tags())
188
189Br = Reflective_boundary(domain)
190Bs = Transmissive_stage_zero_momentum_boundary(domain)
191Bf = Field_boundary(project.event_sts+'.sts',
192                    domain,
193                    mean_stage=project.tide,
194                    time_thinning=1,
195                    default_boundary=Dirichlet_boundary([0, 0, 0]),
196                    boundary_polygon=bounding_polygon_sts,                   
197                    use_cache=True,
198                    verbose=True)
199
200domain.set_boundary({'back': Br,
201                     'side': Bs,
202                     'ocean': Bf}) 
203
204#-------------------------------------------------------------------------------
205# Evolve system through time
206#-------------------------------------------------------------------------------
207
208t0 = time.time()
209
210# Skip over the first 6000 seconds
211for t in domain.evolve(yieldstep=2000,
212                       finaltime=6000):
213    log.critical(domain.timestepping_statistics())
214    log.critical(domain.boundary_statistics(tags='ocean'))
215
216# Start detailed model
217for t in domain.evolve(yieldstep=project.yieldstep,
218                       finaltime=project.finaltime,
219                       skip_initial_step=True):
220    log.critical(domain.timestepping_statistics())
221    log.critical(domain.boundary_statistics(tags='ocean'))
222   
223log.critical('Simulation took %.2f seconds' %(time.time()-t0))
224     
Note: See TracBrowser for help on using the repository browser.