source: anuga_work/production/onslow_2006/ticket160_run.py @ 5442

Last change on this file since 5442 was 5442, checked in by ole, 16 years ago

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.

File size: 8.6 KB
Line 
1"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
2
3Source data such as elevation and boundary data is assumed to be available in
4directories specified by project_urs.py
5The output sww file is stored in project_urs.output_run_time_dir
6
7The scenario is defined by a triangular mesh created from project_urs.polygon,
8the elevation data and a simulated tsunami generated with URS code.
9
10Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
11"""
12
13#------------------------------------------------------------------------------
14# Import necessary modules
15#------------------------------------------------------------------------------
16
17# Standard modules
18from os import sep
19from os.path import dirname, basename
20from os import mkdir, access, F_OK
21from shutil import copy
22import time
23import sys
24
25# Related major packages
26from anuga.shallow_water import Domain
27from anuga.shallow_water import Dirichlet_boundary
28from anuga.shallow_water import File_boundary
29from anuga.shallow_water import Reflective_boundary
30from anuga.shallow_water import Field_boundary
31from Numeric import allclose
32
33from anuga.pmesh.mesh_interface import create_mesh_from_regions
34from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files,store_parameters
35from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
36from anuga_parallel.parallel_abstraction import get_processor_name
37from anuga.caching import myhash
38# Application specific imports
39import aticket160_project as project_urs
40
41def run_model(**kwargs):
42   
43    tide = kwargs['tide']
44    alpha = kwargs['alpha']
45    friction = kwargs['friction']
46    time_thinning = kwargs['time_thinning']
47    scenario_name = kwargs['aa_scenario_name']
48   
49   
50   
51    #------------------------------------------------------------------------------
52    # Copy scripts to time stamped output directory and capture screen
53    # output to file
54    #------------------------------------------------------------------------------
55
56    #copy script must be before screen_catcher
57    print 'tide',tide
58    kwargs['est_num_trigs']=project_urs.trigs_min
59    kwargs['num_cpu']=numprocs
60    kwargs['host']=project_urs.host
61    kwargs['res_factor']=project_urs.res_factor
62    kwargs['starttime']=project_urs.starttime   
63    kwargs['yieldstep']=project_urs.yieldstep
64    kwargs['finaltime']=project_urs.finaltime
65   
66    kwargs['output_dir']=project_urs.output_run_time_dir
67    kwargs['bathy_file']=project_urs.combined_dir_name + '.pts'
68    kwargs['boundary_file']=project_urs.boundaries_dir_name + '.sww'
69#    kwargs['Completed']=''
70
71    print 'output_dir',kwargs['output_dir']
72   
73    print 'myid',myid
74    if myid == 0:
75        copy_code_files(kwargs['output_dir'],__file__, 
76                 dirname(project_urs.__file__)+sep+ project_urs.__name__+'.py' )
77
78        store_parameters(**kwargs)
79
80    barrier()
81
82    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
83
84    print "Processor Name:",get_processor_name()
85#    barrier()
86
87    # filenames
88    #boundaries_name = project_urs.boundaries_name
89    meshes_dir_name = project_urs.meshes_dir_name+'.msh'
90    #boundaries_dir_name = project_urs.boundaries_dir_name
91
92#    tide = project_urs.tide
93
94    # creates copy of code in output dir
95    print 'min triangles', project_urs.trigs_min,
96    print 'Note: This is generally about 20% less than the final amount'
97
98    #--------------------------------------------------------------------------
99    # Create the triangular mesh based on overall clipping polygon with a
100    # tagged
101    # boundary and interior regions defined in project_urs.py along with
102    # resolutions (maximal area of per triangle) for each polygon
103    #--------------------------------------------------------------------------
104
105    if myid == 0:
106   
107        print 'start create mesh from regions'
108
109        create_mesh_from_regions(project_urs.poly_all,
110                             boundary_tags={'back': [4, 5], 'side': [3, 6],
111                                            'ocean': [0, 1, 2]},
112                             maximum_triangle_area=project_urs.res_poly_all,
113                             interior_regions=project_urs.interior_regions,
114                             filename=meshes_dir_name,
115                             use_cache=False,
116                             verbose=True)
117    barrier()
118   
119
120    #-------------------------------------------------------------------------
121    # Setup computational domain
122    #-------------------------------------------------------------------------
123    print 'Setup computational domain'
124    #from caching import cache
125
126    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
127    #above don't work
128    domain = Domain(meshes_dir_name, use_cache=False, verbose=True)
129#    print 'domain id', id(domain)
130#    print 'myhash', myhash(domain)     
131       
132    print domain.statistics()
133#    domain.
134    print 'triangles',len(domain)
135   
136    kwargs['act_num_trigs']=len(domain)
137#    print len(domain.mesh)
138
139
140    #-------------------------------------------------------------------------
141    # Setup initial conditions
142    #-------------------------------------------------------------------------
143    if myid == 0:
144
145        print 'Setup initial conditions'
146
147        from polygon import Polygon_function
148        #following sets the stage/water to be offcoast only
149        IC = Polygon_function( [(project_urs.poly_mainland, -1.0)], default = tide,
150                                 geo_reference = domain.geo_reference)
151        domain.set_quantity('stage', IC)
152        domain.set_quantity('friction', friction) 
153       
154        print 'Start Set quantity'
155
156        domain.set_quantity('elevation', 
157                            filename = kwargs['bathy_file'],
158                            use_cache = True,
159                            verbose = True,
160                            alpha = alpha)
161        print 'Finished Set quantity'
162    barrier()
163
164
165    #------------------------------------------------------
166    # Distribute domain to implement parallelism !!!
167    #------------------------------------------------------
168
169    if numprocs > 1:
170        domain=distribute(domain)
171
172    #------------------------------------------------------
173    # Set domain parameters
174    #------------------------------------------------------
175    print 'domain id', id(domain)
176    domain.set_name(scenario_name)
177    domain.set_datadir(kwargs['output_dir'])
178    domain.set_default_order(2) # Apply second order scheme
179    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
180    domain.set_store_vertices_uniquely(False)
181    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
182    domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
183    domain.H0=0.01
184
185    #-------------------------------------------------------------------------
186    # Setup boundary conditions
187    #-------------------------------------------------------------------------
188    print 'Available boundary tags', domain.get_boundary_tags()
189    print 'domain id', id(domain)
190    #print 'Reading Boundary file',project_urs.boundaries_dir_namea + '.sww'
191
192    #Bf = Field_boundary(kwargs['boundary_file'],
193     #               domain, time_thinning=time_thinning, mean_stage=tide,
194      #              use_cache=True, verbose=True)
195                   
196    kwargs['input_start_time']=domain.starttime
197
198    print 'finished reading boundary file'
199
200    Br = Reflective_boundary(domain)
201    Bd = Dirichlet_boundary([tide,0,0])
202    Bf = Dirichlet_boundary([5.,0,0])
203
204    print'set_boundary'
205
206    domain.set_boundary({'back': Br,
207                         'side': Bd,
208                         'ocean': Bf}) 
209    print'finish set boundary'
210
211    #----------------------------------------------------------------------------
212    # Evolve system through time
213    #----------------------------------------------------------------------------
214
215    t0 = time.time()
216    print "in time"
217    for t in domain.evolve(yieldstep = 10, finaltime = 20):
218        print "yeah", 
219        domain.write_time()
220
221    if myid == 0:
222        store_parameters(**kwargs)
223    barrier()
224
225
226#-------------------------------------------------------------
227if __name__ == "__main__":
228
229    run_model(file_name=project_urs.home+'detail.csv',
230              aa_scenario_name=project_urs.scenario_name,
231              ab_time=project_urs.time,
232              res_factor= project_urs.res_factor,
233              tide=project_urs.tide, user=project_urs.user,
234              alpha = project_urs.alpha, friction=project_urs.friction,
235              time_thinning = project_urs.time_thinning,
236              dir_comment=project_urs.dir_comment)
237
238
Note: See TracBrowser for help on using the repository browser.