source: anuga_work/production/carnarvon/run_carnarvon_250m.py @ 5793

Last change on this file since 5793 was 5793, checked in by kristy, 16 years ago
File size: 9.5 KB
Line 
1"""Script for running a tsunami inundation scenario for carnarvon, WA, Australia.
2
3The scenario is defined by a triangular mesh created from project_250m.polygon,
4the elevation data is compiled into a pts file through build_carnarvon.py
5and a simulated 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_carnarvon.py)
9       information from project file
10Outputs: sww file stored in project_250m.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
23from os import sep
24import os
25from os.path import dirname, basename
26from os import mkdir, access, F_OK
27from shutil import copy
28import time
29import sys
30
31# Related major packages
32from anuga.shallow_water import Domain
33from anuga.shallow_water import Dirichlet_boundary
34from anuga.shallow_water import File_boundary
35from anuga.shallow_water import Reflective_boundary
36from anuga.shallow_water import Field_boundary
37from Numeric import allclose
38from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
39from anuga.pmesh.mesh_interface import create_mesh_from_regions
40from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
41from anuga_parallel.parallel_abstraction import get_processor_name
42from anuga.caching import myhash
43from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
44from anuga.fit_interpolate.benchmark_least_squares import mem_usage
45from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
46from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
47from polygon import Polygon_function
48   
49# Application specific imports
50import project_250m  # Definition of file names and polygons
51numprocs = 1
52myid = 0
53
54def run_model(**kwargs):
55   
56    #------------------------------------------------------------------------------
57    # Copy scripts to time stamped output directory and capture screen
58    # output to file
59    #------------------------------------------------------------------------------
60    print "Processor Name:",get_processor_name()
61
62    #copy script must be before screen_catcher
63
64    print 'output_dir',kwargs['output_dir']
65   
66    copy_code_files(kwargs['output_dir'],__file__, 
67             dirname(project_250m.__file__)+sep+ project_250m.__name__+'.py' )
68
69    store_parameters(**kwargs)
70
71    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
72
73    print "Processor Name:",get_processor_name()
74   
75    #-----------------------------------------------------------------------
76    # Domain definitions
77    #-----------------------------------------------------------------------
78
79    # Read in boundary from ordered sts file
80    urs_bounding_polygon=create_sts_boundary(os.path.join(project_250m.boundaries_dir_event,project_250m.scenario_name))
81
82    # Reading the landward defined points, this incorporates the original clipping
83    # polygon minus the 100m contour
84    landward_bounding_polygon = read_polygon(project_250m.landward_dir)
85
86    # Combine sts polyline with landward points
87    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
88   
89    # counting segments
90    N = len(urs_bounding_polygon)-1
91
92    # boundary tags refer to project_250m.landward 4 points equals 5 segments start at N
93    boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}
94
95
96    #--------------------------------------------------------------------------
97    # Create the triangular mesh based on overall clipping polygon with a tagged
98    # boundary and interior regions defined in project_250m.py along with
99    # resolutions (maximal area of per triangle) for each polygon
100    #--------------------------------------------------------------------------
101
102    # IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
103    # causes problems with the ability to cache set quantity which takes alot of times
104       
105    print 'start create mesh from regions'
106
107    create_mesh_from_regions(bounding_polygon,
108                         boundary_tags=boundary_tags,
109                         maximum_triangle_area=project_250m.res_poly_all,
110                         interior_regions=project_250m.interior_regions,
111                         filename=project_250m.meshes_dir_name,
112                         use_cache=True,
113                         verbose=True)
114   
115    #-------------------------------------------------------------------------
116    # Setup computational domain
117    #-------------------------------------------------------------------------
118    print 'Setup computational domain'
119
120    domain = Domain(project_250m.meshes_dir_name, use_cache=False, verbose=True)
121    print 'memory usage before del domain',mem_usage()
122       
123    print domain.statistics()
124    print 'triangles',len(domain)
125   
126    kwargs['act_num_trigs']=len(domain)
127
128
129    #-------------------------------------------------------------------------
130    # Setup initial conditions
131    #-------------------------------------------------------------------------
132    print 'Setup initial conditions'
133
134    # sets the initial stage in the offcoast region only
135    IC = Polygon_function( [(project_250m.poly_mainland, 0),(project_250m.poly_island1, 0)
136                            ,(project_250m.poly_island2, 0)], default = kwargs['tide'],
137                             geo_reference = domain.geo_reference)
138    domain.set_quantity('stage', IC)
139    #domain.set_quantity('stage',kwargs['tide'] )
140    domain.set_quantity('friction', kwargs['friction']) 
141   
142    print 'Start Set quantity',kwargs['elevation_file']
143
144    domain.set_quantity('elevation', 
145                        filename = kwargs['elevation_file'],
146                        use_cache = False,
147                        verbose = True,
148                        alpha = kwargs['alpha'])
149    print 'Finished Set quantity'
150
151##   #------------------------------------------------------
152##    # Distribute domain to implement parallelism !!!
153##    #------------------------------------------------------
154##
155##    if numprocs > 1:
156##        domain=distribute(domain)
157
158    #------------------------------------------------------
159    # Set domain parameters
160    #------------------------------------------------------
161    print 'domain id', id(domain)
162    domain.set_name(kwargs['scenario_name'])
163    domain.set_datadir(kwargs['output_dir'])
164    domain.set_default_order(2)                 # Apply second order scheme
165    domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
166    domain.set_store_vertices_uniquely(False)
167    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
168    domain.tight_slope_limiters = 1
169    print 'domain id', id(domain)
170
171    #-------------------------------------------------------------------------
172    # Setup boundary conditions
173    #-------------------------------------------------------------------------
174    print 'Available boundary tags', domain.get_boundary_tags()
175    print 'domain id', id(domain)
176   
177    boundary_urs_out=project_250m.boundaries_dir_event + sep + project_250m.scenario_name
178
179    Br = Reflective_boundary(domain)
180    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
181   
182    print 'Available boundary tags', domain.get_boundary_tags()
183    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
184                   domain, mean_stage= project_250m.tide,
185                   time_thinning=1,
186                   default_boundary=Bd,
187                   use_cache=True,
188                   verbose = True,
189                   boundary_polygon=bounding_polygon)
190
191    domain.set_boundary({'back': Br,
192                         'side': Bd,
193                         'ocean': Bf}) 
194
195    kwargs['input_start_time']=domain.starttime
196
197    print'finish set boundary'
198
199    #----------------------------------------------------------------------------
200    # Evolve system through time
201    #--------------------------------------------------------------------
202    t0 = time.time()
203
204    for t in domain.evolve(yieldstep = project_250m.yieldstep, finaltime = kwargs['finaltime']
205                       ,skip_initial_step = False): 
206        domain.write_time()
207        domain.write_boundary_statistics(tags = 'ocean')
208
209    # these outputs should be checked with the resultant inundation map
210    x, y = domain.get_maximum_inundation_location()
211    q = domain.get_maximum_inundation_elevation()
212    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
213
214    print 'Simulation took %.2f seconds' %(time.time()-t0)
215
216    #kwargs 'completed' must be added to write the final parameters to file
217    kwargs['completed']=str(time.time()-t0)
218     
219    store_parameters(**kwargs)
220     
221    print 'memory usage before del domain1',mem_usage()
222   
223   
224#-------------------------------------------------------------
225if __name__ == "__main__":
226   
227    kwargs={}
228    kwargs['finaltime']=project_250m.finaltime
229    kwargs['output_dir']=project_250m.output_run_time_dir
230    kwargs['elevation_file']=project_250m.combined_dir_name+'.pts'
231    kwargs['scenario_name']=project_250m.scenario_name
232    kwargs['tide']=project_250m.tide
233    kwargs['alpha'] = project_250m.alpha
234    kwargs['friction']=project_250m.friction
235     
236    run_model(**kwargs)
237     
238   
Note: See TracBrowser for help on using the repository browser.