source: anuga_validation/automated_validation_tests/patong_beach_validation/run_model.py @ 7276

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

Merged numpy branch back into the trunk.

In ~/sandpit/anuga/anuga_core/source
svn merge -r 6246:HEAD ../../branches/numpy .

In ~/sandpit/anuga/anuga_validation
svn merge -r 6417:HEAD ../branches/numpy_anuga_validation .

In ~/sandpit/anuga/misc
svn merge -r 6809:HEAD ../branches/numpy_misc .

For all merges, I used numpy version where conflicts existed

The suites test_all.py (in source/anuga) and validate_all.py passed using Python2.5 with numpy on my Ubuntu Linux box.

  • Property svn:executable set to *
File size: 8.5 KB
Line 
1"""Run a tsunami inundation scenario for Busselton, WA, Australia.
2
3The scenario is defined by a triangular mesh created from project.polygon, the
4elevation data is compiled into a pts file through build_elevation.py and a
5simulated tsunami is generated through an sts file from build_boundary.py.
6
7Input: sts file (build_boundary.py for respective event)
8       pts file (build_elevation.py)
9       information from project file
10Outputs: sww file stored in project.output_run_time_dir
11The export_results_all.py and get_timeseries.py is reliant
12on the outputs of this script
13
14Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006
15Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008
16"""
17
18#------------------------------------------------------------------------------
19# Import necessary modules
20#------------------------------------------------------------------------------
21
22# Standard modules
23import os
24import os.path
25import time
26from time import localtime, strftime, gmtime
27
28# Related major packages
29from Scientific.IO.NetCDF import NetCDFFile
30#import numpy as num
31
32from anuga.interface import create_domain_from_regions
33from anuga.interface import Transmissive_stage_zero_momentum_boundary
34from anuga.interface import Dirichlet_boundary
35from anuga.interface import Reflective_boundary
36from anuga.interface import Field_boundary
37from anuga.interface import create_sts_boundary
38from anuga.interface import csv2building_polygons
39from anuga.utilities.system_tools import file_length
40
41from anuga.shallow_water.data_manager import start_screen_catcher
42from anuga.shallow_water.data_manager import copy_code_files
43from anuga.shallow_water.data_manager import urs2sts
44from anuga.utilities.polygon import read_polygon, Polygon_function
45from anuga.caching import cache
46
47import anuga.utilities.log as log
48
49# Application specific imports
50from setup_model import project
51import build_urs_boundary as bub
52
53#-------------------------------------------------------------------------------
54# Copy scripts to time stamped output directory and capture screen
55# output to file. Copy script must be before screen_catcher
56#-------------------------------------------------------------------------------
57
58copy_code_files(project.output_run,
59                [__file__, 
60                 os.path.join(os.path.dirname(project.__file__),
61                              project.__name__+'.py'),
62                 os.path.join(os.path.dirname(project.__file__),
63                              'setup_model.py')],
64                verbose=True
65               )
66#start_screen_catcher(project.output_run, 0, 1)
67
68#-------------------------------------------------------------------------------
69# Create the computational domain based on overall clipping polygon with
70# a tagged boundary and interior regions defined in project.py along with
71# resolutions (maximal area of per triangle) for each polygon
72#-------------------------------------------------------------------------------
73
74log.critical('Create computational domain')
75
76# Create the STS file
77# FIXME (Ole): This is deadly dangerous if buildcode changes (as was the case 24th March 2009)
78# We need to use caching instead!
79
80log.critical( 'project.mux_data_folder=%s' % project.mux_data_folder)
81if not os.path.exists(project.event_sts + '.sts'):
82    bub.build_urs_boundary(project.mux_input_filename, project.event_sts)
83
84# Read in boundary from ordered sts file
85event_sts = create_sts_boundary(project.event_sts)
86
87# Reading the landward defined points, this incorporates the original clipping
88# polygon minus the 100m contour
89landward_boundary = read_polygon(project.landward_boundary)
90
91# Combine sts polyline with landward points
92bounding_polygon_sts = event_sts + landward_boundary
93
94# Number of boundary segments
95num_ocean_segments = len(event_sts) - 1
96# Number of landward_boundary points
97num_land_points = file_length(project.landward_boundary)
98
99# Boundary tags refer to project.landward_boundary
100# 4 points equals 5 segments start at N
101boundary_tags={'back': range(num_ocean_segments+1,
102                             num_ocean_segments+num_land_points),
103               'side': [num_ocean_segments,
104                        num_ocean_segments+num_land_points],
105               'ocean': range(num_ocean_segments)}
106
107# Build mesh and domain
108domain = create_domain_from_regions(bounding_polygon_sts,
109                                    boundary_tags=boundary_tags,
110                                    maximum_triangle_area=project.bounding_maxarea,
111                                    interior_regions=project.interior_regions,
112                                    mesh_filename=project.meshes,
113                                    use_cache=True,
114                                    verbose=True)
115log.critical(domain.statistics())
116
117# FIXME(Ole): How can we make this more automatic?
118domain.geo_reference.zone = project.zone
119
120
121domain.set_name(project.scenario_name)
122domain.set_datadir(project.output_run) 
123domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
124
125#-------------------------------------------------------------------------------
126# Setup initial conditions
127#-------------------------------------------------------------------------------
128
129log.critical('Setup initial conditions')
130
131# Set the initial stage in the offcoast region only
132if project.land_initial_conditions:
133    IC = Polygon_function(project.land_initial_conditions,
134                          default=project.tide,
135                          geo_reference=domain.geo_reference)
136else:
137    IC = 0
138domain.set_quantity('stage', IC, use_cache=True, verbose=True)
139domain.set_quantity('friction', project.friction) 
140domain.set_quantity('elevation', 
141                    filename=project.combined_elevation+'.pts',
142                    use_cache=True,
143                    verbose=True,
144                    alpha=project.alpha)
145
146if project.use_buildings:
147    # Add buildings from file
148    log.critical('Reading building polygons')
149    building_polygons, building_heights = csv2building_polygons(project.building_polygon)
150    #clipping_polygons=project.building_area_polygons)
151
152    log.critical('Creating %d building polygons' % len(building_polygons))
153    def create_polygon_function(building_polygons, geo_reference=None):
154        L = []
155        for i, key in enumerate(building_polygons):
156            if i%100==0: log.critical(i)
157            poly = building_polygons[key]
158            elev = building_heights[key]
159            L.append((poly, elev))
160           
161            buildings = Polygon_function(L, default=0.0,
162                                         geo_reference=geo_reference)
163        return buildings
164
165    log.critical('Creating %d building polygons' % len(building_polygons))
166    buildings = cache(create_polygon_function,
167                      building_polygons,
168                      {'geo_reference': domain.geo_reference},
169                      verbose=True)
170
171    log.critical('Adding buildings')
172    domain.add_quantity('elevation',
173                        buildings,
174                        use_cache=True,
175                        verbose=True)
176
177
178#-------------------------------------------------------------------------------
179# Setup boundary conditions
180#-------------------------------------------------------------------------------
181
182log.critical('Set boundary - available tags:' % domain.get_boundary_tags())
183
184Br = Reflective_boundary(domain)
185Bs = Transmissive_stage_zero_momentum_boundary(domain)
186Bf = Field_boundary(project.event_sts+'.sts',
187                    domain,
188                    mean_stage=project.tide,
189                    time_thinning=1,
190                    default_boundary=Dirichlet_boundary([0, 0, 0]),
191                    boundary_polygon=bounding_polygon_sts,                   
192                    use_cache=True,
193                    verbose=True)
194
195domain.set_boundary({'back': Br,
196                     'side': Bs,
197                     'ocean': Bf}) 
198
199#-------------------------------------------------------------------------------
200# Evolve system through time
201#-------------------------------------------------------------------------------
202
203t0 = time.time()
204
205# Skip over the first 6000 seconds
206for t in domain.evolve(yieldstep=2000,
207                       finaltime=6000):
208    log.critical(domain.timestepping_statistics())
209    log.critical(domain.boundary_statistics(tags='ocean'))
210
211# Start detailed model
212for t in domain.evolve(yieldstep=project.yieldstep,
213                       finaltime=project.finaltime,
214                       skip_initial_step=True):
215    log.critical(domain.timestepping_statistics())
216    log.critical(domain.boundary_statistics(tags='ocean'))
217   
218log.critical('Simulation took %.2f seconds' %(time.time()-t0))
219     
Note: See TracBrowser for help on using the repository browser.