source: anuga_work/production/australia_ph2/derby/run_model.py @ 6856

Last change on this file since 6856 was 6767, checked in by myall, 15 years ago

newer project scripts with shared caching and loop over multiple events possible. added scripts for new model, Derby

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