source: anuga_work/production/busselton/standardised_version/run_model.py @ 6343

Last change on this file since 6343 was 6343, checked in by kristy, 15 years ago

Added build_boundary

File size: 13.1 KB
RevLine 
[6263]1"""Run a tsunami inundation scenario for Busselton, WA, Australia.
[6253]2
[6263]3The scenario is defined by a triangular mesh created from project.polygon, the
[6284]4elevation data is compiled into a pts file through build_elevation.py and a
[6263]5simulated tsunami is generated through an sts file from build_boundary.py.
[6253]6
7Input: sts file (build_boundary.py for respective event)
[6284]8       pts file (build_elevation.py)
[6253]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
[6343]24import os.path
[6253]25import time
[6343]26from time import localtime, strftime, gmtime
[6253]27
28# Related major packages
[6343]29from Scientific.IO.NetCDF import NetCDFFile
30import Numeric as num
31
[6259]32from anuga.interface import create_domain_from_regions
[6279]33from anuga.interface import Transmissive_stage_zero_momentum_boundary
[6259]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
[6343]39<<<<<<< .mine
40=======
[6328]41from file_length import file_length
[6259]42
[6343]43>>>>>>> .r6341
[6259]44from anuga.shallow_water.data_manager import start_screen_catcher
45from anuga.shallow_water.data_manager import copy_code_files
[6343]46from anuga.shallow_water.data_manager import urs2sts
[6259]47from anuga.utilities.polygon import read_polygon, Polygon_function
[6343]48
[6253]49# Application specific imports
[6324]50from setup_model import project
[6253]51
[6261]52
[6312]53#-------------------------------------------------------------------------------
[6343]54# Get gauges (timeseries of index points)
55#-------------------------------------------------------------------------------
56def get_sts_gauge_data(filename, verbose=False):
57    fid = NetCDFFile(filename+'.sts', 'r')      #Open existing file for read
58    permutation = fid.variables['permutation'][:]
59    x = fid.variables['x'][:] + fid.xllcorner   #x-coordinates of vertices
60    y = fid.variables['y'][:] + fid.yllcorner   #y-coordinates of vertices
61    points = num.transpose(num.asarray([x.tolist(), y.tolist()]))
62    time = fid.variables['time'][:] + fid.starttime
63    elevation = fid.variables['elevation'][:]
64       
65    basename = 'sts_gauge'
66    quantity_names = ['stage', 'xmomentum', 'ymomentum']
67    quantities = {}
68    for i, name in enumerate(quantity_names):
69        quantities[name] = fid.variables[name][:]
70
71    #---------------------------------------------------------------------------
72    # Get maximum wave height throughout timeseries at each index point
73    #---------------------------------------------------------------------------
74
75    maxname = 'max_sts_stage.csv'
76    fid_max = open(os.path.join(project.event_sts, maxname), 'w')
77    fid_max.write('index, x, y, max_stage \n')   
78    for j in range(len(x)):
79        index = permutation[j]
80        stage = quantities['stage'][:,j]
81        xmomentum = quantities['xmomentum'][:,j]
82        ymomentum = quantities['ymomentum'][:,j]
83
84        fid_max.write('%d, %.6f, %.6f, %.6f\n' % (index, x[j], y[j], max(stage)))
85     
86    #---------------------------------------------------------------------------
87    # Get minimum wave height throughout timeseries at each index point
88    #---------------------------------------------------------------------------
89
90    minname = 'min_sts_stage.csv'
91    fid_min = open(os.path.join(project.event_sts, minname), 'w')
92    fid_min.write('index, x, y, max_stage \n')   
93    for j in range(len(x)):
94        index = permutation[j]
95        stage = quantities['stage'][:,j]
96        xmomentum = quantities['xmomentum'][:,j]
97        ymomentum = quantities['ymomentum'][:,j]
98
99        fid_min.write('%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage)))
100
101        out_file = os.path.join(project.event_sts,
102                                basename+'_'+str(index)+'.csv')
103        fid_sts = open(out_file, 'w')
104        fid_sts.write('time, stage, xmomentum, ymomentum \n')
105
106        #-----------------------------------------------------------------------
107        # End of the get gauges
108        #-----------------------------------------------------------------------
109        for k in range(len(time)-1):
110            fid_sts.write('%.6f, %.6f, %.6f, %.6f\n'
111                          % (time[k], stage[k], xmomentum[k], ymomentum[k]))
112
113        fid_sts.close()     
114
115    fid.close()
116
117    return quantities,elevation,time
118
119
120##
121# @brief Build boundary STS files from one or more MUX files.
122# @param mux_dir Directory containing the MUX files.
123# @param event_file Path to meta-file containing mux files+weight data.
124# @param output_dir Directory to write STS data to.
125# @note 'event_file' is produced by EventSelection.
126def build_urs_boundary(mux_dir, event_file, output_dir):
127    '''Build a boundary STS file from a set of MUX files.'''
128
129    print 'build_urs_boundary: mux_dir=%s' % mux_dir
130    print 'build_urs_boundary: event_file=%s' % event_file
131    print 'build_urs_boundary: output_dir=%s' % output_dir
132
133    # if we are using an EventSelection multi-mux file
134    if project.multi_mux:   
135        # get the mux+weight data from the event file
136        mux_event_file = os.path.join(mux_dir, event_file)
137        try:
138            fd = open(mux_event_file, 'r')
139            mux_data = fd.readlines()
140            fd.close()
141        except IOError, e:
142            msg = 'File %s cannot be read: %s' % (mux_event_file, str(e))
143            raise Exception, msg
144        except:
145            raise
146
147        # first line of file is # filenames+weight in rest of file
148        num_lines = int(mux_data[0].strip())
149        mux_data = mux_data[1:]
150        print 'number of sources %d' % num_lines
151
152        # quick sanity check on input mux meta-file
153        if num_lines != len(mux_data):
154            msg = ('Bad file %s: %d data lines, but line 1 count is %d'
155                   % (event_file, len(mux_data), num_lines))
156            raise Exception, msg
157
158        # Create filename and weights lists.
159        # Must chop GRD filename just after '*.grd'.
160        mux_filenames = []
161        for line in mux_data:
162            muxname = line.strip().split()[0]
163            split_index = muxname.index('.grd')
164            muxname = muxname[:split_index+len('.grd')]
165            muxname = os.path.join(mux_dir, muxname)
166            mux_filenames.append(muxname)
167       
168        mux_weights = [float(line.strip().split()[1]) for line in mux_data]
169
170        # Call legacy function to create STS file.
171        # This should be replaced in the future.
172        print 'reading', project.urs_order
173        print 'creating STS file'
174        print 'mux_filenames=%s' % str(mux_filenames)
175        print 'basename_out=%s' % str(output_dir)
176        print 'ordering_filename=%s' % str(project.urs_order)
177        print 'weights=%s' % str(mux_weights)
178        print 'mean_stage=%s' % str(project.tide)
179        urs2sts(mux_filenames,
180                basename_out=output_dir,
181                ordering_filename=project.urs_order,
182                weights=mux_weights,
183                mean_stage=project.tide,
184                verbose=True)
185    else:                           # a single mux stem file
186        urs_filenames = [os.path.join(mux_dir, event_file)]
187
188        weight_factor = 1.0
189        weights = weight_factor*num.ones(len(urs_filenames), num.float)
190           
191        order_filename = os.path.join(project.order_filename_dir)
192
193        print 'reading', order_filename
194        # Create ordered sts file
195        print 'creating sts file'
196        urs2sts(urs_filenames,
197                basename_out=os.path.join(project.boundaries_dir,
198                                          project.scenario_name),
199                ordering_filename=order_filename,
200                weights=weights,
201                mean_stage=project.tide,
202                verbose=True)
203
204    # report on stuff so far
205    quantities, elevation, time = get_sts_gauge_data(project.event_folder,
206                                                     verbose=False)
207    print len(elevation), len(quantities['stage'][0,:])
208
209#-------------------------------------------------------------------------------
[6259]210# Copy scripts to time stamped output directory and capture screen
[6263]211# output to file. Copy script must be before screen_catcher
[6312]212#-------------------------------------------------------------------------------
213
[6279]214copy_code_files(project.output_run, __file__, 
[6312]215                os.path.join(os.path.dirname(project.__file__),
216                             project.__name__+'.py'))
[6279]217start_screen_catcher(project.output_run, 0, 1)
[6253]218
[6312]219#-------------------------------------------------------------------------------
[6259]220# Create the computational domain based on overall clipping polygon with
221# a tagged boundary and interior regions defined in project.py along with
222# resolutions (maximal area of per triangle) for each polygon
[6312]223#-------------------------------------------------------------------------------
224
[6259]225print 'Create computational domain'
[6253]226
[6312]227# Create the STS file
[6343]228<<<<<<< .mine
229build_urs_boundary(project.mux_data_folder,
230                   project.mux_input,
231                   os.path.join(project.event_folder,
232                                project.scenario_name))
233=======
[6329]234print 'project.mux_data_folder=%s' % project.mux_data_folder
235bub.build_urs_boundary(project.mux_input_filename,
[6312]236                       os.path.join(project.event_folder,
237                                    project.scenario_name))
[6343]238>>>>>>> .r6341
[6312]239
[6259]240# Read in boundary from ordered sts file
[6279]241event_sts = create_sts_boundary(project.event_sts)
[6253]242
[6259]243# Reading the landward defined points, this incorporates the original clipping
244# polygon minus the 100m contour
[6279]245landward_boundary = read_polygon(project.landward_boundary)
[6253]246
[6259]247# Combine sts polyline with landward points
[6279]248bounding_polygon_sts = event_sts + landward_boundary
[6253]249
[6259]250# Number of boundary segments
[6329]251num_ocean_segments = len(event_sts) - 1
[6328]252# Number of landward_boundary points
[6329]253num_land_points = file_length(project.landward_boundary)
[6253]254
[6281]255# Boundary tags refer to project.landward_boundary
256# 4 points equals 5 segments start at N
[6329]257boundary_tags={'back': range(num_ocean_segments+1,
258                             num_ocean_segments+num_land_points),
259               'side': [num_ocean_segments,
260                        num_ocean_segments+num_land_points],
261               'ocean': range(num_ocean_segments)}
[6253]262
[6259]263# Build mesh and domain
[6279]264domain = create_domain_from_regions(bounding_polygon_sts,
[6259]265                                    boundary_tags=boundary_tags,
[6284]266                                    maximum_triangle_area=project.bounding_maxarea,
[6259]267                                    interior_regions=project.interior_regions,
[6279]268                                    mesh_filename=project.meshes,
[6259]269                                    use_cache=True,
270                                    verbose=True)
271print domain.statistics()
[6253]272
[6259]273domain.set_name(project.scenario_name)
[6279]274domain.set_datadir(project.output_run) 
[6259]275domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
[6254]276
[6312]277#-------------------------------------------------------------------------------
278# Setup initial conditions
279#-------------------------------------------------------------------------------
[6263]280
[6259]281print 'Setup initial conditions'
[6253]282
[6259]283# Set the initial stage in the offcoast region only
[6279]284IC = Polygon_function(project.land_initial_conditions,
[6259]285                      default=project.tide,
286                      geo_reference=domain.geo_reference)
287domain.set_quantity('stage', IC, use_cache=True, verbose=True)
288domain.set_quantity('friction', project.friction) 
289domain.set_quantity('elevation', 
[6279]290                    filename=project.combined_elevation+'.pts',
[6259]291                    use_cache=True,
292                    verbose=True,
[6267]293                    alpha=project.alpha)
[6253]294
[6312]295#-------------------------------------------------------------------------------
296# Setup boundary conditions
297#-------------------------------------------------------------------------------
[6253]298
[6259]299print 'Set boundary - available tags:', domain.get_boundary_tags()
[6253]300
[6259]301Br = Reflective_boundary(domain)
[6279]302Bt = Transmissive_stage_zero_momentum_boundary(domain)
[6312]303Bd = Dirichlet_boundary([kwargs['tide'], 0, 0])
[6279]304Bf = Field_boundary(project.event_sts+'.sts',
[6259]305                    domain, mean_stage=project.tide,
306                    time_thinning=1,
307                    default_boundary=Bd,
[6279]308                    boundary_polygon=bounding_polygon_sts,                   
[6259]309                    use_cache=True,
310                    verbose=True)
[6253]311
[6259]312domain.set_boundary({'back': Br,
[6279]313                     'side': Bt,
[6259]314                     'ocean': Bf}) 
[6253]315
[6312]316#-------------------------------------------------------------------------------
317# Evolve system through time
318#-------------------------------------------------------------------------------
[6262]319
[6259]320t0 = time.time()
321for t in domain.evolve(yieldstep=project.yieldstep, 
322                       finaltime=project.finaltime,
323                       skip_initial_step=False): 
324    print domain.timestepping_statistics()
325    print domain.boundary_statistics(tags='ocean')
[6253]326
[6312]327print 'Simulation took %.2f seconds' % (time.time()-t0)
Note: See TracBrowser for help on using the repository browser.