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

Last change on this file since 7567 was 7567, checked in by ole, 14 years ago

Got rid of environment variables and move log file to output dir.
The logging module will need some work, though.

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