source: anuga_work/production/shark_bay_2007/run_shark_bay_frequency_sweep.py @ 4689

Last change on this file since 4689 was 4689, checked in by ole, 17 years ago

Steep point back log

File size: 9.0 KB
RevLine 
[4625]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,umask
19from os.path import dirname, basename
20from os import mkdir, access, F_OK
21from shutil import copy
22from math import sin, pi
23import time
24import sys
25
26# Related major packages
27from anuga.shallow_water import Domain
28from anuga.shallow_water import Dirichlet_boundary
29from anuga.shallow_water import File_boundary
30from anuga.shallow_water import Reflective_boundary
31from anuga.shallow_water import Field_boundary
32from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
33from Numeric import allclose
34
35from anuga.pmesh.mesh_interface import create_mesh_from_regions
36from anuga.shallow_water.data_manager import start_screen_catcher
37from anuga.shallow_water.data_manager import copy_code_files
38from anuga.shallow_water.data_manager import store_parameters
39from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
40from anuga_parallel.parallel_abstraction import get_processor_name
41from anuga.caching import myhash
42
43# Application specific imports
44import project                 # Definition of file names and polygons
45
46def run_model(**kwargs):
47   
48    tide = kwargs['tide']
49    alpha = kwargs['alpha']
50    friction = kwargs['friction']
51    time_thinning = kwargs['time_thinning']
52    scenario_name = kwargs['aa_scenario_name']
53   
54    #------------------------------------------------------------------------------
55    # Copy scripts to time stamped output directory and capture screen
56    # output to file
57    #------------------------------------------------------------------------------
58
59    #copy script must be before screen_catcher
60    print 'tide',tide
61    kwargs['est_num_trigs']=project.trigs_min
62    kwargs['num_cpu']=numprocs
63    kwargs['host']=project.host
64    kwargs['res_factor']=project.res_factor
65    kwargs['starttime']=project.starttime
66    kwargs['yieldstep']=project.yieldstep
67    kwargs['finaltime']=project.finaltime
68   
69    kwargs['output_dir']=project.output_run_time_dir
70    kwargs['bathy_file']=project.combined_dir_name + '.pts'
71#    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
72    kwargs['boundary_file']=project.boundary_name + '.sww'
73#    kwargs['Completed']=''
74
75    print 'output_dir',kwargs['output_dir']
76    if myid == 0:
77        copy_code_files(kwargs['output_dir'],__file__, 
78                        dirname(project.__file__)+sep+project.__name__+'.py' )
79
80        store_parameters(**kwargs)
81
82    barrier()
83
84    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
85
86    print 'Processor Name:', get_processor_name()
87
88    #--------------------------------------------------------------------------
89    # Create the triangular mesh based on overall clipping polygon with a
90    # tagged
91    # boundary and interior regions defined in project.py along with
92    # resolutions (maximal area of per triangle) for each polygon
93    #--------------------------------------------------------------------------
94
95    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
96    # causes problems with the ability to cache set quantity which takes alot of times
97
98    # FIXME (Ole): Move to build_shark_bay.py
99    if myid == 0:
100   
101        print 'start create mesh from regions'
102       
103        create_mesh_from_regions(project.small_bounding_polygon,
104                                 boundary_tags=project.boundary_tags_small,
105                                 maximum_triangle_area=project.res_bounding_polygon,
106                                 interior_regions=project.interior_regions,
107                                 filename=project.mesh_name+'.msh',
[4689]108                                 fail_if_polygons_outside=False,
[4625]109                                 use_cache=False,
110                                 verbose=True)
111
112    #-------------------------------------------------------------------------
113    # Setup computational domain
114    #-------------------------------------------------------------------------
115    print 'Setup computational domain'
116
117    #domain = cache(Domain, (mesh_name), {'use_cache':True, 'verbose':True}, verbose=True)
118    #above don't work
119    domain = Domain(project.mesh_name+'.msh',
120                    use_cache=False, verbose=True)
121       
122    print domain.statistics()
123    print 'triangles',len(domain)
124   
125    kwargs['act_num_trigs']=len(domain)
126
127    #-------------------------------------------------------------------------
128    # Setup initial conditions
129    #-------------------------------------------------------------------------
130    if myid == 0:
131
132        IC = project.tide
133
134        print 'Set initial conditions'
135        domain.set_quantity('stage', IC)
136        domain.set_quantity('friction', friction)
137        print 'Set elevation'
138        domain.set_quantity('elevation', 
139                            filename = kwargs['bathy_file'],
140                            use_cache = True,
141                            verbose = True,
142                            alpha = alpha)
143       
144
145    #------------------------------------------------------
146    # Distribute domain to implement parallelism !!!
147    #------------------------------------------------------
148    barrier()
149    if numprocs > 1:
150        domain=distribute(domain)
151
152    #------------------------------------------------------
153    # Set domain parameters
154    #------------------------------------------------------
155    print 'domain id', id(domain)
156    domain.set_name(scenario_name)
157    domain.set_datadir(kwargs['output_dir'])
158    domain.set_default_order(2) # Apply second order scheme
159    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
160    domain.set_store_vertices_uniquely(False)
161    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
162    domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
163    print 'domain id', id(domain)
164    domain.beta_h = 0
[4631]165    #domain.tight_slope_limiters = 1
[4625]166   
167
168    #-------------------------------------------------------------------------
169    # Setup boundary conditions
170    #-------------------------------------------------------------------------
171    print 'Available boundary tags', domain.get_boundary_tags()
172    print 'domain id', id(domain)
173    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
174
175    if project.boundary_event == 'experimental':
176        # Artificial       
177        #Bf = Dirichlet_boundary([2,0,0])
178        def waveform(t):
179            delay = 900
180            A = project.amplitude
181            T = project.period
182            if t < delay:
183                # Let things settle
184                return project.tide
185            else:
186                return project.tide + A*sin(2*pi*(t-delay)/T)
187
188        Bf = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
189    else:   
190        Bf = Field_boundary(kwargs['boundary_file'],
191                            domain,
192                            time_thinning=time_thinning,
193                            mean_stage=tide, 
194                            use_cache=False,
195                            verbose=True)
196                   
197    kwargs['input_start_time']=domain.starttime
198
199    print 'finished reading boundary file'
200
201    Br = Reflective_boundary(domain)
202    Bd = Dirichlet_boundary([tide,0,0])
203
204
205    print'set_boundary'
206
207    domain.set_boundary({'tide': Bd,
208                         'ocean': Bf})     
209    print'finish set boundary'
210
211    #----------------------------------------------------------------------------
212    # Evolve system through time
213    #----------------------------------------------------------------------------
214
215    t0 = time.time()
216
217    for t in domain.evolve(yieldstep = 5,
[4689]218                           finaltime = 10000): 
[4625]219        domain.write_time()
220        domain.write_boundary_statistics(tags = 'ocean')     
221
222
223    x, y = domain.get_maximum_inundation_location()
224    q = domain.get_maximum_inundation_elevation()
225
226    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
227
228    print 'That took %.2f seconds' %(time.time()-t0)
229   
230    kwargs['completed']=str(time.time()-t0)
231   
232    store_parameters(**kwargs)
233
234#-------------------------------------------------------------
235if __name__ == "__main__":
236
237    run_model(file_name=project.home+'detail.csv', aa_scenario_name=project.scenario_name,
238              ab_time=project.time, res_factor= project.res_factor, tide=project.tide, user=project.user,
239              alpha = project.alpha, friction=project.friction,
240              time_thinning = project.time_thinning,
241              dir_comment=project.dir_comment)
242
243
Note: See TracBrowser for help on using the repository browser.