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

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

Played with resolution and timestepping

File size: 7.9 KB
RevLine 
[6190]1"""Script for running a tsunami inundation scenario for Patong Beach, Thailand.
[5971]2
3The scenario is defined by a triangular mesh created from project.polygon,
[6190]4the elevation data is compiled into a pts file through build_patong.py
5and a simulated tsunami for the 2004 event is generated through
6an sts file from build_boundary.py.
[5971]7
8Input: sts file (build_boundary.py)
[6190]9       pts file (build_patong.py)
[5971]10       information from project file
11Outputs: sww file stored in project.output_run_time_dir
12The export_results_all.py and get_timeseries.py is reliant
13on the outputs of this script
14
15"""
16
17#------------------------------------------------------------------------------
18# Import necessary modules
19#------------------------------------------------------------------------------
20
21# Standard modules
22from os import sep
23import os
[6178]24from os.path import dirname
[5971]25import time
26
27# Related major packages
[6190]28from anuga.interface import create_domain_from_regions
29from anuga.interface import Dirichlet_boundary
[6264]30from anuga.interface import Transmissive_stage_zero_momentum_boundary
[6190]31from anuga.interface import Reflective_boundary
32from anuga.interface import Field_boundary
[6219]33from anuga.interface import create_sts_boundary
[6190]34from anuga.interface import csv2building_polygons
[6321]35from file_length import file_length
[6190]36
[6194]37from anuga.caching import cache
[5971]38from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
[6219]39from anuga.utilities.polygon import read_polygon
[5971]40from polygon import Polygon_function
41   
42# Application specific imports
43import project  # Definition of file names and polygons
44numprocs = 1
45myid = 0
46
[6178]47
[5971]48   
[6178]49#------------------------------------------------------------------------------
50# Copy scripts to time stamped output directory and capture screen
51# output to file
52#------------------------------------------------------------------------------
[5971]53
[6178]54#copy script must be before screen_catcher
[5971]55
[6178]56output_dir = project.output_run_time_dir
57print 'output_dir', output_dir
[5971]58
[6323]59copy_code_files(output_dir,__file__, 
60                dirname(project.__file__)+sep+ project.__name__+'.py' )
[5971]61
[6323]62start_screen_catcher(output_dir, myid, numprocs)
[5971]63
[6323]64
[6178]65#-----------------------------------------------------------------------
66# Domain definitions
67#-----------------------------------------------------------------------
[5971]68
[6178]69# Read in boundary from ordered sts file
70urs_boundary_name = os.path.join(project.boundaries_dir,
71                                 project.scenario_name)
72urs_bounding_polygon = create_sts_boundary(urs_boundary_name)
[5971]73
[6178]74# Reading the landward defined points, this incorporates the original clipping
75# polygon minus the 100m contour
76landward_bounding_polygon = read_polygon(project.landward_dir)
[5971]77
[6178]78# Combine sts polyline with landward points
79bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
[5971]80
[6178]81# Counting segments
82N = len(urs_bounding_polygon)-1
[6321]83# Number of landward_boundary points
84M = file_length(project.landward_dir)
[5971]85
[6322]86
[6321]87# Boundary tags refer to project.landward_boundary
88# 4 points equals 5 segments start at N
89boundary_tags={'back': range(N+1,N+M),
90               'side': [N,N+M],
91               'ocean': range(N)}
[5971]92
[6178]93#--------------------------------------------------------------------------
94# Create the computational domain based on overall clipping polygon with
95# a tagged boundary and interior regions defined in project.py along with
96# resolutions (maximal area of per triangle) for each polygon
97#--------------------------------------------------------------------------
[5971]98
[6178]99domain = create_domain_from_regions(bounding_polygon,
100                                    boundary_tags=boundary_tags,
101                                    maximum_triangle_area=project.res_poly_all,
102                                    interior_regions=project.interior_regions,
103                                    mesh_filename=project.meshes_dir_name,
104                                    use_cache=True,
105                                    verbose=True)
[5971]106
[6178]107print domain.statistics()
[5971]108
[6178]109#-------------------------------------------------------------------------
110# Setup initial conditions
111#-------------------------------------------------------------------------
112print 'Setup initial conditions'
[5971]113
[6178]114# sets the initial stage in the offcoast region only
115IC = Polygon_function([(project.poly_mainland, 0)],
116                      default=project.tide,
117                      geo_reference=domain.geo_reference)
[6219]118print 'stage'
119domain.set_quantity('stage', IC,
120                    use_cache=True,
121                    verbose=True)
122print 'friction'
123domain.set_quantity('friction', project.friction)
[6322]124
[6219]125print 'elevation'
[6225]126domain.set_quantity('elevation',
[6178]127                    filename=project.combined_dir_name+'.pts',
[6225]128                    alpha=project.alpha,                   
[6178]129                    use_cache=True,
[6225]130                    verbose=True)
[6133]131
[6264]132if project.use_buildings:
133    # Add buildings from file
134    print 'Reading building polygons'   
135    building_polygons, building_heights = csv2building_polygons(project.building_polygon_file)
136    #clipping_polygons=project.building_area_polygons)
[6225]137
[6264]138    print 'Creating %d building polygons' % len(building_polygons)
139    def create_polygon_function(building_polygons, geo_reference=None):
140        L = []
141        for i, key in enumerate(building_polygons):
142            if i%100==0: print i
143            poly = building_polygons[key]
144            elev = building_heights[key]
145            L.append((poly, elev))
146           
147            buildings = Polygon_function(L, default=0.0,
148                                         geo_reference=geo_reference)
149        return buildings
[6133]150
[6264]151    buildings = cache(create_polygon_function,
152                      building_polygons,
153                      {'geo_reference': domain.geo_reference},
154                      verbose=True)
[5971]155
[6264]156    print 'Adding buildings'
157    domain.add_quantity('elevation',
158                        buildings,
159                        use_cache=True,
160                        verbose=True)
[6194]161
[6219]162
[5971]163
[6178]164#------------------------------------------------------
165# Distribute domain to implement parallelism !!!
166#------------------------------------------------------
[5971]167
[6178]168#  if numprocs > 1:
169#      domain=distribute(domain)
[5971]170
[6178]171#------------------------------------------------------
172# Set domain parameters
173#------------------------------------------------------
174domain.set_name(project.scenario_name)
175domain.set_datadir(output_dir)
176domain.set_default_order(2)                 # Apply second order scheme
177domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
[5971]178
[6178]179#-------------------------------------------------------------------------
180# Setup boundary conditions
181#-------------------------------------------------------------------------
182print 'Set boundary conditions'
[5971]183
[6178]184Br = Reflective_boundary(domain)
[6322]185Bs = Dirichlet_boundary([project.tide,0,0])
[6321]186#Bs = Transmissive_stage_zero_momentum_boundary(domain)
[6178]187Bf = Field_boundary(urs_boundary_name+'.sts', 
188                    domain,
189                    mean_stage= project.tide,
190                    time_thinning=project.time_thinning,
[6322]191                    default_boundary=Bs,
[6178]192                    use_cache=True,
193                    verbose=True,
194                    boundary_polygon=bounding_polygon)
[5971]195
[6178]196domain.set_boundary({'back': Br,
[6322]197                     'side': Bs,
[6178]198                     'ocean': Bf}) 
[5971]199
200
[6178]201#----------------------------------------------------------------------------
202# Evolve system through time
203#--------------------------------------------------------------------
204t0 = time.time()
[5971]205
[6334]206# Skip over the first 6000 seconds
207for t in domain.evolve(yieldstep=2000,
208                       finaltime=6000):
[6264]209    print domain.timestepping_statistics()
210    print domain.boundary_statistics(tags='ocean')
211
212# Start detailed model
[6178]213for t in domain.evolve(yieldstep=project.yieldstep,
[6264]214                       finaltime=project.finaltime,
215                       skip_initial_step=True):
[6178]216    print domain.timestepping_statistics()
217    print domain.boundary_statistics(tags='ocean')
[6192]218   
[6178]219print 'Simulation took %.2f seconds' %(time.time()-t0)
[5971]220     
221   
Note: See TracBrowser for help on using the repository browser.