source: anuga_work/production/carnarvon/run_carnarvon.py @ 6022

Last change on this file since 6022 was 6022, checked in by kristy, 15 years ago

added file_name to kwarg

File size: 9.7 KB
RevLine 
[5790]1"""Script for running a tsunami inundation scenario for carnarvon, WA, Australia.
[5005]2
3The scenario is defined by a triangular mesh created from project.polygon,
[5790]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.
[5005]6
[5790]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.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
[5005]16"""
17
18#------------------------------------------------------------------------------
19# Import necessary modules
20#------------------------------------------------------------------------------
21
22# Standard modules
23from os import sep
[5755]24import os
[5005]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
[5755]38from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
[5005]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
[5755]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
[5790]47from polygon import Polygon_function
48   
[5005]49# Application specific imports
[5790]50import project  # Definition of file names and polygons
[5755]51numprocs = 1
52myid = 0
[5005]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']
[5790]65   
66    copy_code_files(kwargs['output_dir'],__file__, 
67             dirname(project.__file__)+sep+ project.__name__+'.py' )
[5005]68
[5790]69    store_parameters(**kwargs)
[5005]70
71    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
72
73    print "Processor Name:",get_processor_name()
[5755]74   
75    #-----------------------------------------------------------------------
76    # Domain definitions
77    #-----------------------------------------------------------------------
[5005]78
[5755]79    # Read in boundary from ordered sts file
[5790]80    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir_event,project.scenario_name))
[5005]81
[5755]82    # Reading the landward defined points, this incorporates the original clipping
83    # polygon minus the 100m contour
[5790]84    landward_bounding_polygon = read_polygon(project.landward_dir)
[5755]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
[5790]92    # boundary tags refer to project.landward 4 points equals 5 segments start at N
[5793]93    boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}
[5790]94
[5005]95    #--------------------------------------------------------------------------
[5790]96    # Create the triangular mesh based on overall clipping polygon with a tagged
[5005]97    # boundary and interior regions defined in project.py along with
98    # resolutions (maximal area of per triangle) for each polygon
99    #--------------------------------------------------------------------------
100
[5790]101    # IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
[5005]102    # causes problems with the ability to cache set quantity which takes alot of times
[5790]103       
104    print 'start create mesh from regions'
[5005]105
[5790]106    create_mesh_from_regions(bounding_polygon,
107                         boundary_tags=boundary_tags,
108                         maximum_triangle_area=project.res_poly_all,
109                         interior_regions=project.interior_regions,
110                         filename=project.meshes_dir_name,
111                         use_cache=True,
112                         verbose=True)
113   
[5005]114    #-------------------------------------------------------------------------
115    # Setup computational domain
116    #-------------------------------------------------------------------------
117    print 'Setup computational domain'
118
[5790]119    domain = Domain(project.meshes_dir_name, use_cache=False, verbose=True)
[5005]120    print 'memory usage before del domain',mem_usage()
121       
122    print domain.statistics()
123    print 'triangles',len(domain)
124   
125    kwargs['act_num_trigs']=len(domain)
126
127
128    #-------------------------------------------------------------------------
129    # Setup initial conditions
130    #-------------------------------------------------------------------------
[5790]131    print 'Setup initial conditions'
[5005]132
[5790]133    # sets the initial stage in the offcoast region only
134    IC = Polygon_function( [(project.poly_mainland, 0),(project.poly_island1, 0)
135                            ,(project.poly_island2, 0)], default = kwargs['tide'],
136                             geo_reference = domain.geo_reference)
137    domain.set_quantity('stage', IC)
138    #domain.set_quantity('stage',kwargs['tide'] )
139    domain.set_quantity('friction', kwargs['friction']) 
140   
141    print 'Start Set quantity',kwargs['elevation_file']
[5005]142
[5790]143    domain.set_quantity('elevation', 
144                        filename = kwargs['elevation_file'],
145                        use_cache = False,
146                        verbose = True,
147                        alpha = kwargs['alpha'])
148    print 'Finished Set quantity'
[5005]149
[5790]150##   #------------------------------------------------------
151##    # Distribute domain to implement parallelism !!!
[5755]152##    #------------------------------------------------------
153##
[5790]154##    if numprocs > 1:
155##        domain=distribute(domain)
[5005]156
157    #------------------------------------------------------
158    # Set domain parameters
159    #------------------------------------------------------
160    print 'domain id', id(domain)
[5790]161    domain.set_name(kwargs['scenario_name'])
[5005]162    domain.set_datadir(kwargs['output_dir'])
[5790]163    domain.set_default_order(2)                 # Apply second order scheme
164    domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
[5005]165    domain.set_store_vertices_uniquely(False)
166    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
[5755]167    domain.tight_slope_limiters = 1
[5005]168    print 'domain id', id(domain)
169
170    #-------------------------------------------------------------------------
171    # Setup boundary conditions
172    #-------------------------------------------------------------------------
173    print 'Available boundary tags', domain.get_boundary_tags()
174    print 'domain id', id(domain)
[5755]175   
[5790]176    boundary_urs_out=project.boundaries_dir_event + sep + project.scenario_name
[5005]177
178    Br = Reflective_boundary(domain)
179    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
[5755]180   
181    print 'Available boundary tags', domain.get_boundary_tags()
182    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
183                   domain, mean_stage= project.tide,
184                   time_thinning=1,
185                   default_boundary=Bd,
186                   use_cache=True,
187                   verbose = True,
188                   boundary_polygon=bounding_polygon)
[5005]189
190    domain.set_boundary({'back': Br,
191                         'side': Bd,
192                         'ocean': Bf}) 
[5755]193
194    kwargs['input_start_time']=domain.starttime
195
[5005]196    print'finish set boundary'
197
198    #----------------------------------------------------------------------------
199    # Evolve system through time
200    #--------------------------------------------------------------------
201    t0 = time.time()
202
[5755]203    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
204                       ,skip_initial_step = False): 
[5005]205        domain.write_time()
[5755]206        domain.write_boundary_statistics(tags = 'ocean')
[5005]207
[5790]208    # these outputs should be checked with the resultant inundation map
[5005]209    x, y = domain.get_maximum_inundation_location()
210    q = domain.get_maximum_inundation_elevation()
211    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
212
[5790]213    print 'Simulation took %.2f seconds' %(time.time()-t0)
[5005]214
215    #kwargs 'completed' must be added to write the final parameters to file
216    kwargs['completed']=str(time.time()-t0)
[5790]217     
218    store_parameters(**kwargs)
219     
[5005]220    print 'memory usage before del domain1',mem_usage()
221   
222   
223#-------------------------------------------------------------
224if __name__ == "__main__":
225   
226    kwargs={}
[6022]227    kwargs['file_name']=project.dir_comment
[5005]228    kwargs['finaltime']=project.finaltime
229    kwargs['output_dir']=project.output_run_time_dir
[5790]230    kwargs['elevation_file']=project.combined_dir_name+'.pts'
231    kwargs['scenario_name']=project.scenario_name
[5005]232    kwargs['tide']=project.tide
233    kwargs['alpha'] = project.alpha
234    kwargs['friction']=project.friction
[5790]235    #kwargs['num_cpu']=numprocs
236    #kwargs['host']=project.host
237    #kwargs['starttime']=project.starttime
238    #kwargs['yieldstep']=project.yieldstep
239    #kwargs['user']=project.user
240    #kwargs['time_thinning'] = project.time_thinning
241     
[5005]242    run_model(**kwargs)
243     
[5790]244   
Note: See TracBrowser for help on using the repository browser.