source: anuga_work/production/shark_bay_2007/run_shark_bay_embayment_study.py @ 4884

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

Shark bay embayment study.
Reduced verbose output further (ticket:238)

File size: 6.4 KB
Line 
1"""Script for running embayment study for Steep Point, WA, Australia.
2
3Title: Embayment study for Steep Point
4
5Description:
6
7Source data such as elevation and boundary data is assumed to be available in
8directories specified by project.py
9
10This script is designed to run a series of idealised incident waveforms
11at different frequencies.
12
13
14Author:  Ole Nielsen (Ole.Nielsen@ga.gov.au)
15
16CreationDate: 2007
17
18
19ModifiedBy:
20    $Author: ole $
21    $Date: 2007-11-30 11:47:13 +1100 (Fri, 30 Nov 2007) $
22    $LastChangedRevision: 4867 $   
23
24"""
25
26
27#------------------------------------------------------------------------------
28# Import necessary modules
29#------------------------------------------------------------------------------
30
31# Standard modules
32from os import sep
33from os.path import dirname
34from math import sin, pi
35import time
36
37# Related major packages
38from anuga.shallow_water import Domain
39from anuga.shallow_water import Dirichlet_boundary
40from anuga.shallow_water import File_boundary
41from anuga.shallow_water import Reflective_boundary
42from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
43from Numeric import allclose
44
45from anuga.pmesh.mesh_interface import create_mesh_from_regions
46from anuga.shallow_water.data_manager import start_screen_catcher
47from anuga.shallow_water.data_manager import copy_code_files
48from anuga.shallow_water.data_manager import store_parameters
49
50# Application specific imports
51import project                 # Definition of file names and polygons
52
53
54def run_model(**kwargs):
55   
56    alpha = kwargs['alpha']
57    friction = kwargs['friction']
58    time_thinning = kwargs['time_thinning']
59    scenario_name = kwargs['aa_scenario_name']
60   
61    #-----------------------------------------------------------------
62    # Copy scripts to time stamped output directory and capture screen
63    # output to file
64    #-----------------------------------------------------------------
65
66    #copy_code_files(project.output_run_time_dir,
67    #                __file__,
68    #                dirname(project.__file__)+sep+project.__name__+'.py' )
69
70
71    #start_screen_catcher(project.output_run_time_dir)
72
73
74    #--------------------------------------------------------------------------
75    # Create the triangular mesh based on overall clipping polygon with a
76    # tagged
77    # boundary and interior regions defined in project.py along with
78    # resolutions (maximal area of per triangle) for each polygon
79    #--------------------------------------------------------------------------
80
81    print 'Create mesh from regions'
82    create_mesh_from_regions(project.small_bounding_polygon,
83                             boundary_tags=project.boundary_tags_small,
84                             maximum_triangle_area=project.res_bounding_polygon,
85                             interior_regions=project.interior_regions,
86                             filename=project.mesh_name+'.msh',
87                             fail_if_polygons_outside=False,
88                             use_cache=False,
89                             verbose=True)
90
91    #-------------------------------------------------------------------------
92    # Setup computational domain
93    #-------------------------------------------------------------------------
94    print 'Setup computational domain'
95    domain = Domain(project.mesh_name+'.msh',
96                    use_cache=False, verbose=True)
97       
98    print domain.statistics()
99   
100
101    #-------------------------------------------------------------------------
102    # Setup initial conditions
103    #-------------------------------------------------------------------------
104
105    print 'Set initial conditions'
106    domain.set_quantity('stage', project.tide)
107    domain.set_quantity('friction', friction)
108    domain.set_quantity('elevation', 
109                        filename = project.combined_dir_name + '.pts',
110                        use_cache = True,
111                        verbose = True,
112                        alpha = project.alpha)
113       
114
115    #------------------------------------------------------
116    # Set domain parameters
117    #------------------------------------------------------
118    domain.set_name(scenario_name)
119    domain.set_datadir(project.output_run_time_dir)
120    domain.set_default_order(2)              # Apply second order scheme
121    domain.set_minimum_storable_height(0.01) # Don't store anything < 1cm
122    domain.set_store_vertices_uniquely(False)
123    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
124
125    domain.beta_h = 0
126    domain.tight_slope_limiters = 1
127   
128
129    #-------------------------------------------------------------------------
130    # Setup boundary conditions
131    #-------------------------------------------------------------------------
132    print 'Available boundary tags', domain.get_boundary_tags()
133    print 'domain id', id(domain)
134
135    def waveform(t):
136        delay = 900
137        A = project.amplitude
138        T = project.period
139        if t < delay:
140            # Let things settle
141            return project.tide
142        else:
143            return project.tide + A*sin(2*pi*(t-delay)/T)
144
145    Bf = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
146    Br = Reflective_boundary(domain)
147    Bd = Dirichlet_boundary([project.tide,0,0])
148
149
150    print 'Set boundary'
151    domain.set_boundary({'tide': Bd,
152                         'ocean': Bf})
153
154
155    #------------------------------
156    # Evolve system through time
157    #------------------------------
158    t0 = time.time()
159    for t in domain.evolve(yieldstep = 5,
160                           finaltime = 10000):
161        domain.write_time()
162        domain.write_boundary_statistics(tags = 'ocean')     
163
164
165    #------------------------------
166    # Post processing
167    #------------------------------
168    x, y = domain.get_maximum_inundation_location()
169    q = domain.get_maximum_inundation_elevation()
170
171    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
172
173    print 'That took %.2f seconds' %(time.time()-t0)
174   
175
176#-------------------------------------------------------------
177if __name__ == "__main__":
178
179    run_model(file_name=project.home+'detail.csv',
180              aa_scenario_name=project.scenario_name,
181              ab_time=project.time,
182              res_factor= project.res_factor,
183              tide=project.tide,
184              user=project.user,
185              alpha = project.alpha,
186              friction=project.friction,
187              time_thinning = project.time_thinning,
188              dir_comment=project.dir_comment)
189
190
Note: See TracBrowser for help on using the repository browser.