source: anuga_work/production/patong/run_patong.py @ 6133

Last change on this file since 6133 was 6133, checked in by ole, 15 years ago

Work towards buildings in Patong simulation

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