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

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

Work on Patong and some diagnostic output

File size: 6.8 KB
Line 
1"""Script for running a tsunami inundation scenario for Patong Beach, Thailand.
2
3The scenario is defined by a triangular mesh created from project.polygon,
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.
7
8Input: sts file (build_boundary.py)
9       pts file (build_patong.py)
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
24from os.path import dirname
25import time
26
27# Related major packages
28from anuga.interface import create_domain_from_regions
29from anuga.interface import Domain
30from anuga.interface import Dirichlet_boundary
31from anuga.interface import File_boundary
32from anuga.interface import Reflective_boundary
33from anuga.interface import Field_boundary
34from anuga.interface import export_grid, create_sts_boundary
35from anuga.interface import csv2building_polygons
36
37from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
38from anuga.fit_interpolate.benchmark_least_squares import mem_usage
39from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
40from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
41from polygon import Polygon_function
42   
43# Application specific imports
44import project  # Definition of file names and polygons
45numprocs = 1
46myid = 0
47
48
49   
50#------------------------------------------------------------------------------
51# Copy scripts to time stamped output directory and capture screen
52# output to file
53#------------------------------------------------------------------------------
54
55#copy script must be before screen_catcher
56
57output_dir = project.output_run_time_dir
58print 'output_dir', output_dir
59
60#copy_code_files(output_dir,__file__,
61#                dirname(project.__file__)+sep+ project.__name__+'.py' )
62
63#start_screen_catcher(output_dir, myid, numprocs)
64
65
66#-----------------------------------------------------------------------
67# Domain definitions
68#-----------------------------------------------------------------------
69
70# Read in boundary from ordered sts file
71urs_boundary_name = os.path.join(project.boundaries_dir,
72                                 project.scenario_name)
73urs_bounding_polygon = create_sts_boundary(urs_boundary_name)
74
75# Reading the landward defined points, this incorporates the original clipping
76# polygon minus the 100m contour
77landward_bounding_polygon = read_polygon(project.landward_dir)
78
79# Combine sts polyline with landward points
80bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
81
82# Counting segments
83N = len(urs_bounding_polygon)-1
84
85# boundary tags refer to project.landward 4 points equals 5 segments start at N
86boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}
87
88#--------------------------------------------------------------------------
89# Create the computational domain based on overall clipping polygon with
90# a tagged boundary and interior regions defined in project.py along with
91# resolutions (maximal area of per triangle) for each polygon
92#--------------------------------------------------------------------------
93
94domain = create_domain_from_regions(bounding_polygon,
95                                    boundary_tags=boundary_tags,
96                                    maximum_triangle_area=project.res_poly_all,
97                                    interior_regions=project.interior_regions,
98                                    mesh_filename=project.meshes_dir_name,
99                                    use_cache=True,
100                                    verbose=True)
101
102print 'memory usage after creation of domain', mem_usage()
103print domain.statistics()
104
105
106#-------------------------------------------------------------------------
107# Setup initial conditions
108#-------------------------------------------------------------------------
109print 'Setup initial conditions'
110
111# sets the initial stage in the offcoast region only
112IC = Polygon_function([(project.poly_mainland, 0)],
113                      default=project.tide,
114                      geo_reference=domain.geo_reference)
115domain.set_quantity('stage', IC)
116domain.set_quantity('friction', project.friction) 
117domain.set_quantity('elevation', 
118                    filename=project.combined_dir_name+'.pts',
119                    use_cache=True,
120                    verbose=True,
121                    alpha=project.alpha)
122
123# Add buildings from file
124print 'Reading building polygons'   
125building_polygons, building_heights = csv2building_polygons(project.building_polygon_file)
126
127print 'Creating building polygons'   
128L = []
129for key in building_polygons:
130    poly = building_polygons[key]
131    elev = building_heights[key]
132    L.append((poly, elev))
133
134buildings = Polygon_function(L, default=0.0)
135domain.add_quantity('elevation', buildings)
136
137
138#------------------------------------------------------
139# Distribute domain to implement parallelism !!!
140#------------------------------------------------------
141
142#  if numprocs > 1:
143#      domain=distribute(domain)
144
145#------------------------------------------------------
146# Set domain parameters
147#------------------------------------------------------
148domain.set_name(project.scenario_name)
149domain.set_datadir(output_dir)
150domain.set_default_order(2)                 # Apply second order scheme
151domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
152
153#-------------------------------------------------------------------------
154# Setup boundary conditions
155#-------------------------------------------------------------------------
156print 'Set boundary conditions'
157
158Br = Reflective_boundary(domain)
159Bd = Dirichlet_boundary([project.tide,0,0])
160Bf = Field_boundary(urs_boundary_name+'.sts', 
161                    domain,
162                    mean_stage= project.tide,
163                    time_thinning=project.time_thinning,
164                    default_boundary=Bd,
165                    use_cache=True,
166                    verbose=True,
167                    boundary_polygon=bounding_polygon)
168
169domain.set_boundary({'back': Br,
170                     'side': Bd,
171                     'ocean': Bf}) 
172
173
174#----------------------------------------------------------------------------
175# Evolve system through time
176#--------------------------------------------------------------------
177t0 = time.time()
178
179for t in domain.evolve(yieldstep=project.yieldstep,
180                       finaltime=project.finaltime):
181    print domain.timestepping_statistics()
182    print domain.boundary_statistics(tags='ocean')
183   
184print 'Simulation took %.2f seconds' %(time.time()-t0)
185     
186   
Note: See TracBrowser for help on using the repository browser.