source: anuga_work/production/wollongong_2006/run_flagstaff_sequential.py @ 3585

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

Work on parallel API

File size: 3.9 KB
Line 
1"""Script for running a tsunami inundation scenario for Flagstaff pt,
2Wollongong harbour, NSW, Australia.
3
4Source data such as elevation and boundary data is assumed to be available in
5directories specified by project.py
6
7The scenario is defined by a triangular mesh created from project.polygon,
8the elevation data and a hypothetical boundary condition.
9
10THIS IS THE SEQUENTIAL VERSION
11
12Ole Nielsen and Duncan Gray, GA - 2005, Nick Bartzis and Jane Sexton, GA - 2006
13"""
14
15
16#------------------------------------------------------------------------------
17# Import necessary modules
18#------------------------------------------------------------------------------
19
20
21# Standard modules
22import os, sys, time 
23from os import sep
24from os.path import dirname, basename
25from Numeric import zeros, Float
26
27# Related major packages
28from anuga.pyvolution.shallow_water import Domain
29from anuga.pyvolution.shallow_water import Dirichlet_boundary
30from anuga.pyvolution.shallow_water import Time_boundary
31from anuga.pyvolution.shallow_water import Reflective_boundary
32from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
33from anuga.pmesh.mesh_interface import create_mesh_from_regions
34from anuga.pmesh.mesh import importUngenerateFile, Segment
35
36
37# Application specific imports
38import project
39
40
41
42
43   
44#--------------------------------------------------------------------------
45# Create the triangular mesh based on overall clipping polygon with a
46# tagged boundary and interior regions defined in project.py along with
47# resolutions (maximal area of per triangle) for each polygon
48#--------------------------------------------------------------------------
49
50
51print 'Generate mesh'
52max_area = project.base_resolution
53mesh = create_mesh_from_regions(project.bounding_polygon,
54                                boundary_tags=project.boundary_tags,
55                                maximum_triangle_area=max_area,
56                                interior_regions=project.interior_regions)
57   
58# Add buildings
59# This should bind to a Reflective boundary
60mesh.import_ungenerate_file(project.buildings_filename, tag='wall')
61
62# Generate and write mesh to file
63mesh.generate_mesh(maximum_triangle_area=max_area,
64                   verbose=True)
65mesh.export_mesh_file(project.mesh_filename)
66
67
68#--------------------------------------------------------------------------
69# Setup computational domain
70#--------------------------------------------------------------------------
71
72domain = Domain(project.mesh_filename, use_cache = False, verbose = True)
73print domain.statistics()
74
75domain.set_name(project.basename)
76domain.set_datadir(project.outputdir)
77
78#------------------------------------------------------------------------------
79# Setup initial conditions
80#------------------------------------------------------------------------------
81
82domain.set_quantity('stage', project.initial_sealevel)
83domain.set_quantity('friction', 0.03)
84domain.set_quantity('elevation', 
85                    filename=project.demname + '.pts',
86                    use_cache=True,
87                    verbose=True)   
88
89#------------------------------------------------------------------------------
90# Setup boundary conditions
91#------------------------------------------------------------------------------
92
93D = Dirichlet_boundary([project.initial_sealevel, 0, 0])
94R = Reflective_boundary(domain)
95W = Time_boundary(domain = domain,
96                  f=lambda t: [project.initial_sealevel + (60<t<480)*6, 0, 0])
97
98domain.set_boundary({'exterior': D,
99                     'side': D,
100                     'wall': R,
101                     'ocean': W,
102                     'ghost': None})
103
104
105
106
107#------------------------------------------------------------------------------
108# Evolve system through time
109#------------------------------------------------------------------------------
110
111t0 = time.time()
112for t in domain.evolve(yieldstep = 1, finaltime = 1200):
113    domain.write_time()
114
115 
116print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.