source: anuga_work/production/shark_bay_2007/run_shark_bay_embayment_study.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: 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.tight_slope_limiters = True
121   
122
123#-------------------------------------------------------------------------
124# Setup boundary conditions
125#-------------------------------------------------------------------------
126print 'Available boundary tags', domain.get_boundary_tags()
127print 'domain id', id(domain)
128
129def waveform(t):
130    delay = 900
131    A = project.amplitude
132    T = project.period
133    if t < delay:
134        # Let things settle
135        return project.tide
136    else:
137        return project.tide + A*sin(2*pi*(t-delay)/T)
138
139Bf = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
140Br = Reflective_boundary(domain)
141Bd = Dirichlet_boundary([project.tide,0,0])
142
143
144print 'Set boundary'
145domain.set_boundary({'tide': Bd,
146                     'ocean': Bf})
147
148
149#------------------------------
150# Evolve system through time
151#------------------------------
152t0 = time.time()
153for t in domain.evolve(yieldstep = 5,
154                       finaltime = 10000):
155    domain.write_time()
156    domain.write_boundary_statistics(tags='ocean')     
157
158
159#------------------------------
160# Post processing
161#------------------------------
162x, y = domain.get_maximum_inundation_location()
163q = domain.get_maximum_inundation_elevation()
164
165print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
166
167print 'That took %.2f seconds' %(time.time()-t0)
168   
169
170
Note: See TracBrowser for help on using the repository browser.