source: anuga_work/development/joaqium_luis/run.py @ 7840

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

Moved data_audit to utilities

Caching work on Joaquim Luis example

File size: 6.0 KB
Line 
1"""Script for running a tsunami inundation scenario for Boca do Rio.
2Adopted from cairns files
3
4"""
5
6#------------------------------------------------------------------------------
7# Import necessary modules
8#------------------------------------------------------------------------------
9
10# Standard modules
11import os
12import time
13import sys
14
15# Related major packages
16from anuga.shallow_water import Domain
17from anuga.shallow_water import Transmissive_boundary
18from anuga.shallow_water import Reflective_boundary
19from anuga.shallow_water import Dirichlet_boundary
20from anuga.shallow_water import Time_boundary
21from anuga.shallow_water import File_boundary
22from anuga.shallow_water import Field_boundary
23
24from anuga.pmesh.mesh_interface import create_mesh_from_regions
25from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
26from anuga.shallow_water.data_manager import dem2pts
27from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
28from anuga.abstract_2d_finite_volumes.util import file_function
29
30from anuga.caching import cache
31from anuga.utilities.polygon import read_polygon, plot_polygons, \
32                                    polygon_area, is_inside_polygon
33
34#------------------------------------------------------------------------------
35# Define scenario as either ...
36#------------------------------------------------------------------------------
37scenario = 'br'
38
39if os.access(scenario, os.F_OK) == 0:
40    os.mkdir(scenario)
41basename = scenario + 'source'
42
43#------------------
44# Initial condition
45#------------------
46tide = 0.0
47
48
49#------------------------------------------------------------------------------
50# Preparation of topographic data
51# Convert ASC 2 DEM 2 PTS using source data and store result in source data
52#------------------------------------------------------------------------------
53
54# Filenames
55#dem_name = 'sub_region_bat'
56#meshname = 'sub_region.msh'
57dem_name = 'bocarra' 
58meshname = 'bocarra.msh'
59
60
61#---------
62# Polygons
63#---------
64
65# bounding polygon for study area
66bounding_polygon = read_polygon('out_bocarra.csv')
67
68print 'Area of bounding polygon in km?', polygon_area(bounding_polygon)/1000000.0
69
70# interior polygons
71#poly_midle = read_polygon('midle_rect.csv')
72poly_shallow = read_polygon('shallow.csv')
73
74
75#------------
76# Resolutions
77#------------
78remainder_res = 2000
79midle_res = 3500
80shallow_res = 100
81
82interior_regions = [[poly_shallow, shallow_res]]
83
84def setup_domain(tide,
85                 bounding_polygon,
86                 remainder_res,                 
87                 interior_regions,
88                 mesh_filename,
89                 points_filename):
90    """Set up domain except boundary conditions
91    This function is defined in order to cache it.
92    """
93
94
95
96
97
98    create_mesh_from_regions(bounding_polygon,
99                             boundary_tags={'west': [0],
100                                            'north': [1],
101                                            'east': [2],
102                                            'south': [3]},
103                             maximum_triangle_area=remainder_res,
104                             filename=mesh_filename,
105                             interior_regions=interior_regions,
106                             use_cache=False,
107                             verbose=True)
108
109    #--------------------------------------------------------------------------
110    # Setup computational domain
111    #--------------------------------------------------------------------------
112
113    domain = Domain(mesh_filename,
114                    use_cache=False,
115                    verbose=True)
116
117    print 'Number of triangles = ', len(domain)
118    print 'The extent is ', domain.get_extent()
119    print domain.statistics()
120
121    domain.set_name(basename) 
122    domain.set_datadir(scenario)
123    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
124    domain.set_minimum_storable_height(0.05)
125
126    #--------------------------------------------------------------------------
127    # Setup initial conditions
128    #--------------------------------------------------------------------------
129
130    domain.set_quantity('stage', tide)
131    domain.set_quantity('friction', 0.025) 
132    domain.set_quantity('elevation',
133                        filename=points_filename,
134                        use_cache=False,
135                        verbose=True)
136
137    return domain
138   
139
140#--------------------------------   
141# Call (and cache) setup function
142#--------------------------------   
143
144# Run set_up domain without caching
145#domain = setup_domain(tide, bounding_polygon, poly_shallow,
146#                      remainder_res, midle_res, shallow_res,
147#                      meshname,
148#                      dem_name + '.pts')
149
150# Run set_up domain with caching
151domain = cache(setup_domain,
152               (tide,
153                bounding_polygon,
154                remainder_res, 
155                interior_regions, 
156                meshname, 
157                dem_name + '.pts'),
158               #dependencies =
159               verbose=True)           
160
161
162
163#------------------------------------------------------------------------------
164# Setup boundary conditions
165#------------------------------------------------------------------------------
166
167print 'Available boundary tags', domain.get_boundary_tags()
168Br = Reflective_boundary(domain)
169Bd = Dirichlet_boundary([tide,0,0])
170Bf = Field_boundary('swan_tsu_stageLand.sww', domain, time_thinning=1, mean_stage=tide,
171                    use_cache=True, verbose=True)
172   
173# Boundary condition for sww feed at the east boundary
174domain.set_boundary({'west': Bf,'south':Bf,'east': Br,'north': Br})
175
176   
177#------------------------------------------------------------------------------
178# Evolve system through time
179#------------------------------------------------------------------------------
180
181import time
182t0 = time.time()
183
184from Numeric import allclose
185from anuga.abstract_2d_finite_volumes.quantity import Quantity
186
187# Save every 1 sec leading up to wave approaching land
188for t in domain.evolve(yieldstep = 10, finaltime = 90):
189    domain.write_time()
190    domain.write_boundary_statistics(tags = ['west','south'])
191
192           
193print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.