source: anuga_work/production/wollongong_2006/run_flagstaff_parallel_api.py @ 3584

Last change on this file since 3584 was 3584, checked in by ole, 18 years ago

Parallel work

File size: 7.4 KB
Line 
1"""Script for running a tsunami inundation scenario for Flagstaff pt,
2Wollongong harbour, NSW, Australia.
3
4Source data such as elevation and boundary data is assumed to be available in
5directories specified by project.py
6
7The scenario is defined by a triangular mesh created from project.polygon,
8the elevation data and a hypothetical boundary condition.
9
10Ole Nielsen and Duncan Gray, GA - 2005, Nick Bartzis and Jane Sexton, GA - 2006
11"""
12
13
14#------------------------------------------------------------------------------
15# Import necessary modules
16#------------------------------------------------------------------------------
17
18
19# Standard modules
20import os, sys, time 
21from os import sep
22from os.path import dirname, basename
23from Numeric import zeros, Float
24
25# Related major packages
26from anuga.shallow_water import Domain
27from anuga.shallow_water import Dirichlet_boundary
28from anuga.shallow_water import Time_boundary
29from anuga.shallow_water import Reflective_boundary
30from anuga.pmesh.mesh_interface import create_mesh_from_regions
31from anuga.pmesh.mesh import importUngenerateFile, Segment
32
33# Parallelism
34import pypar   # The Python-MPI interface
35from anuga_parallel.pmesh_divide  import pmesh_divide_metis
36from anuga_parallel.build_submesh import build_submesh
37from anuga_parallel.build_local   import build_local_mesh
38from anuga_parallel.build_commun  import send_submesh, rec_submesh, extract_hostmesh
39from anuga_parallel.parallel_shallow_water import Parallel_Domain
40
41
42# Application specific imports
43import project
44
45
46#------------------------------------------------------------------------------
47# Read in processor information
48#------------------------------------------------------------------------------
49
50numprocs = pypar.size()
51myid = pypar.rank()
52processor_name = pypar.Get_processor_name()
53print 'I am processor %d of %d on node %s' %(myid, numprocs, processor_name)
54
55
56#------------------------------------------------------------------------------
57# Preparation of topographic data
58#
59# Convert ASC 2 DEM 2 PTS using source data and store result in source data
60#------------------------------------------------------------------------------
61
62max_area = project.base_resolution
63if myid == 0:
64   
65    #--------------------------------------------------------------------------
66    # Create the triangular mesh based on overall clipping polygon with a
67    # tagged boundary and interior regions defined in project.py along with
68    # resolutions (maximal area of per triangle) for each polygon
69    #--------------------------------------------------------------------------
70
71
72    print 'Generate mesh'
73    # Generate basic mesh
74    mesh = create_mesh_from_regions(project.bounding_polygon,
75                                    boundary_tags=project.boundary_tags,
76                                    maximum_triangle_area=max_area,
77                                    interior_regions=project.interior_regions)
78   
79    # Add buildings that will bind to a Reflective boundary
80    mesh.import_ungenerate_file(project.buildings_filename, tag='wall')
81
82    # Generate and write mesh to file
83    mesh.generate_mesh(maximum_triangle_area=max_area,
84                       verbose=True)
85   
86    mesh.export_mesh_file(project.mesh_filename)
87
88
89    #--------------------------------------------------------------------------
90    # Setup computational domain
91    #--------------------------------------------------------------------------
92
93    domain = Domain(project.mesh_filename, use_cache = False, verbose = True)
94    print domain.statistics()
95
96    domain.set_quantity('elevation', 
97                        filename=project.demname + '.pts',
98                        use_cache=True,
99                        verbose=True)   
100
101
102    # Subdivide the mesh
103    print 'Subdivide mesh'
104    nodes, triangles, boundary, triangles_per_proc, quantities = \
105           pmesh_divide_metis(domain, numprocs)
106
107    # Build the mesh that should be assigned to each processor,
108    # this includes ghost nodes and the communicaiton pattern
109    print 'Build submeshes'   
110    submesh = build_submesh(nodes, triangles, boundary,\
111                            quantities, triangles_per_proc)
112
113    # Send the mesh partition to the appropriate processor
114    print 'Distribute submeshes'       
115    for p in range(1, numprocs):
116      send_submesh(submesh, triangles_per_proc, p)
117
118    # Build the local mesh for processor 0
119    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
120              extract_hostmesh(submesh, triangles_per_proc)
121
122    print 'Communication done'       
123   
124else:
125    # Read in the mesh partition that belongs to this
126    # processor (note that the information is in the
127    # correct form for the GA data structure)
128
129    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict \
130            = rec_submesh(0)
131
132
133
134
135#------------------------------------------------------------------------------
136# Start the computations on each subpartion
137#------------------------------------------------------------------------------
138
139
140# Build the domain for this processor
141domain = Parallel_Domain(points, vertices, boundary,
142                         full_send_dict  = full_send_dict,
143                         ghost_recv_dict = ghost_recv_dict)
144
145# Name etc currently has to be set here as they are not transferred from the
146# original domain
147domain.set_name(project.basename)
148domain.set_datadir(project.outputdir)
149
150
151#------------------------------------------------------------------------------
152# Setup initial conditions
153#------------------------------------------------------------------------------
154
155
156domain.set_quantity('elevation', quantities['elevation']) # Distribute elevation
157domain.set_quantity('stage', project.initial_sealevel)
158domain.set_quantity('friction', 0.03)
159
160#
161# FIXME (Ole): This one segfaults which is bad, because set_quantity is
162# time consuming and should be done here rather than on processor 0
163# It did not segfault today 2 Aug 2006 !!!
164# But values are zero ??....
165#
166#domain.set_quantity('elevation',
167#                    filename=project.demname + '.pts',
168#                    use_cache=False,
169#                    verbose=True)
170
171
172#------------------------------------------------------------------------------
173# Setup boundary conditions
174#------------------------------------------------------------------------------
175
176D = Dirichlet_boundary([project.initial_sealevel, 0, 0])
177R = Reflective_boundary(domain)
178W = Time_boundary(domain = domain,
179                  f=lambda t: [project.initial_sealevel + (60<t<480)*6, 0, 0])
180
181domain.set_boundary({'exterior': D,
182                     'side': D,
183                     'wall': R,
184                     'ocean': W,
185                     'ghost': None})
186
187
188
189print 'P%d: Ready to evolve. Value of store is %s' %(myid, str(domain.store))
190
191
192#------------------------------------------------------------------------------
193# Evolve system through time
194#------------------------------------------------------------------------------
195
196t0 = time.time()
197for t in domain.evolve(yieldstep = 1, finaltime = 1200):
198    if myid == 0:
199        domain.write_time()
200        #domain.write_boundary_statistics(tags = 'ocean')     
201 
202
203if myid == 0:
204    print 'That took %.2f seconds' %(time.time()-t0)
205    print 'Communication time %.2f seconds'%domain.communication_time
206    print 'Reduction Communication time %.2f seconds'\
207          %domain.communication_reduce_time
208    print 'Broadcast time %.2f seconds'\
209          %domain.communication_broadcast_time
210
211pypar.finalize()
Note: See TracBrowser for help on using the repository browser.