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

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

Work on Shark Bay embayment

File size: 5.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
54   
55   
56#-----------------------------------------------------------------
57# Copy scripts to time stamped output directory and capture screen
58# output to file
59#-----------------------------------------------------------------
60
61copy_code_files(project.output_run_time_dir,
62                __file__, 
63                dirname(project.__file__)+sep+project.__name__+'.py' )
64
65
66start_screen_catcher(project.output_run_time_dir)
67
68
69#--------------------------------------------------------------------------
70# Create the triangular mesh based on overall clipping polygon with a
71# tagged
72# boundary and interior regions defined in project.py along with
73# resolutions (maximal area of per triangle) for each polygon
74#--------------------------------------------------------------------------
75
76print 'Create mesh from regions'
77create_mesh_from_regions(project.small_bounding_polygon,
78                         boundary_tags=project.boundary_tags_small,
79                         maximum_triangle_area=project.res_bounding_polygon,
80                         interior_regions=project.interior_regions,
81                         filename=project.mesh_name+'.msh',
82                         fail_if_polygons_outside=False,
83                         use_cache=False,
84                         verbose=True)
85
86#-------------------------------------------------------------------------
87# Setup computational domain
88#-------------------------------------------------------------------------
89print 'Setup computational domain'
90domain = Domain(project.mesh_name+'.msh',
91                use_cache=False, verbose=True)
92
93print domain.statistics()
94   
95
96#-------------------------------------------------------------------------
97# Setup initial conditions
98#-------------------------------------------------------------------------
99
100print 'Set initial conditions'
101domain.set_quantity('stage', project.tide)
102domain.set_quantity('friction', project.friction)
103domain.set_quantity('elevation', 
104                    filename = project.combined_dir_name + '.pts',
105                    use_cache = True,
106                    verbose = True,
107                    alpha = project.alpha)
108       
109
110#------------------------------------------------------
111# Set domain parameters
112#------------------------------------------------------
113domain.set_name(project.scenario_name)
114domain.set_datadir(project.output_run_time_dir)
115domain.set_default_order(2)              # Apply second order scheme
116domain.set_minimum_storable_height(0.01) # Don't store anything < 1cm
117domain.set_store_vertices_uniquely(False)
118domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
119
120domain.beta_h = 0
121domain.tight_slope_limiters = 1
122   
123
124#-------------------------------------------------------------------------
125# Setup boundary conditions
126#-------------------------------------------------------------------------
127print 'Available boundary tags', domain.get_boundary_tags()
128print 'domain id', id(domain)
129
130def waveform(t):
131    delay = 900
132    A = project.amplitude
133    T = project.period
134    if t < delay:
135        # Let things settle
136        return project.tide
137    else:
138        return project.tide + A*sin(2*pi*(t-delay)/T)
139
140Bf = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
141Br = Reflective_boundary(domain)
142Bd = Dirichlet_boundary([project.tide,0,0])
143
144
145print 'Set boundary'
146domain.set_boundary({'tide': Bd,
147                     'ocean': Bf})
148
149
150#------------------------------
151# Evolve system through time
152#------------------------------
153t0 = time.time()
154for t in domain.evolve(yieldstep = 5,
155                       finaltime = 10000):
156    domain.write_time()
157    domain.write_boundary_statistics(tags='ocean')     
158
159
160#------------------------------
161# Post processing
162#------------------------------
163x, y = domain.get_maximum_inundation_location()
164q = domain.get_maximum_inundation_elevation()
165
166print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
167
168print 'That took %.2f seconds' %(time.time()-t0)
169   
170
171
Note: See TracBrowser for help on using the repository browser.