source: anuga_work/production/australia_ph2/torres/run_model.py @ 6980

Last change on this file since 6980 was 6980, checked in by Leharne, 15 years ago

updates to torres model

File size: 7.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 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
89print num_ocean_segments
90# Number of landward_boundary points
91num_land_points = file_length(project.landward_boundary)
92print num_land_points
93# Boundary tags refer to project.landward_boundary
94# 4 points equals 5 segments start at N
95boundary_tags={'back': [num_ocean_segments+1],
96               'side': [num_ocean_segments,
97                        num_ocean_segments+2,
98                        num_ocean_segments+3,
99                        num_ocean_segments+4],
100               'ocean': range(num_ocean_segments)}
101##boundary_tags = {'back':[3],
102##                 'side':[2,4],
103##                 'ocean':[0,1,5]}
104
105print 'boundary tags', boundary_tags
106
107# Build mesh and domain
108domain = create_domain_from_regions(bounding_polygon_sts,
109                                    boundary_tags=boundary_tags,
110                                    maximum_triangle_area=project.bounding_maxarea,
111                                    interior_regions=project.interior_regions,
112                                    mesh_filename=project.meshes,
113                                    use_cache=True,
114                                    verbose=True)
115print domain.statistics()
116
117domain.set_name(project.scenario_name)
118domain.set_datadir(project.output_run) 
119domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
120
121#-------------------------------------------------------------------------------
122# Setup initial conditions
123#-------------------------------------------------------------------------------
124
125print 'Setup initial conditions'
126
127# Set the initial stage in the offcoast region only
128if project.land_initial_conditions:
129    IC = Polygon_function(project.land_initial_conditions,
130                          default=project.tide,
131                          geo_reference=domain.geo_reference)
132else:
133    IC = 0
134domain.set_quantity('stage', IC, use_cache=True, verbose=True)
135domain.set_quantity('friction', project.friction) 
136domain.set_quantity('elevation', 
137                    filename=project.combined_elevation+'.pts',
138                    use_cache=True,
139                    verbose=True,
140                    alpha=project.alpha)
141
142#-------------------------------------------------------------------------------
143# Setup boundary conditions
144#-------------------------------------------------------------------------------
145
146print 'Set boundary - available tags:', domain.get_boundary_tags()
147
148# Prepare time boundary
149TB.prepare_timeboundary(project.boundary_csv)
150f = file_function(project.boundary_csv[:-4] + '.tms')
151
152Br = Reflective_boundary(domain)
153Bt = Transmissive_stage_zero_momentum_boundary(domain)
154Bd = Dirichlet_boundary([project.tide, 0, 0])
155
156if project.wave == 'Bf':
157    Bf = Field_boundary(project.event_sts+'.sts',
158                        domain, mean_stage=project.tide,
159                        time_thinning=1,
160                        default_boundary=Dirichlet_boundary([0, 0, 0]),
161                        boundary_polygon=bounding_polygon_sts,                   
162                        use_cache=True,
163                        verbose=True)
164    domain.set_boundary({'back': Br,
165                     'side': Bt,
166                     'ocean': Bf}) 
167   
168elif project.wave == 'Tb':
169    Tb = Time_boundary(domain,f,default_boundary=Dirichlet_boundary([0, 0, 0]) )
170
171    domain.set_boundary({'back': Br,
172                         'side': Bt,
173                         'ocean': Tb})
174else:
175    print 'No wave specified in project script (Bf or Tb)'
176   
177
178#-------------------------------------------------------------------------------
179# Evolve system through time
180#-------------------------------------------------------------------------------
181
182t0 = time.time()
183for t in domain.evolve(yieldstep=project.yieldstep, 
184                       finaltime=project.finaltime,
185                       skip_initial_step=False): 
186    print domain.timestepping_statistics()
187    print domain.boundary_statistics(tags='ocean')
188
189print 'Simulation took %.2f seconds' % (time.time()-t0)
Note: See TracBrowser for help on using the repository browser.