[6263] | 1 | """Run a tsunami inundation scenario for Busselton, WA, Australia. |
---|
[6253] | 2 | |
---|
[6263] | 3 | The scenario is defined by a triangular mesh created from project.polygon, the |
---|
[6284] | 4 | elevation data is compiled into a pts file through build_elevation.py and a |
---|
[6263] | 5 | simulated tsunami is generated through an sts file from build_boundary.py. |
---|
[6253] | 6 | |
---|
| 7 | Input: sts file (build_boundary.py for respective event) |
---|
[6284] | 8 | pts file (build_elevation.py) |
---|
[6253] | 9 | information from project file |
---|
| 10 | Outputs: sww file stored in project.output_run_time_dir |
---|
| 11 | The export_results_all.py and get_timeseries.py is reliant |
---|
| 12 | on the outputs of this script |
---|
| 13 | |
---|
| 14 | Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006 |
---|
| 15 | Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008 |
---|
| 16 | """ |
---|
| 17 | |
---|
| 18 | #------------------------------------------------------------------------------ |
---|
| 19 | # Import necessary modules |
---|
| 20 | #------------------------------------------------------------------------------ |
---|
| 21 | |
---|
| 22 | # Standard modules |
---|
| 23 | import os |
---|
[6343] | 24 | import os.path |
---|
[6253] | 25 | import time |
---|
[6343] | 26 | from time import localtime, strftime, gmtime |
---|
[6253] | 27 | |
---|
| 28 | # Related major packages |
---|
[6343] | 29 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 30 | import Numeric as num |
---|
| 31 | |
---|
[6259] | 32 | from anuga.interface import create_domain_from_regions |
---|
[6279] | 33 | from anuga.interface import Transmissive_stage_zero_momentum_boundary |
---|
[6259] | 34 | from anuga.interface import Dirichlet_boundary |
---|
| 35 | from anuga.interface import Reflective_boundary |
---|
| 36 | from anuga.interface import Field_boundary |
---|
| 37 | from anuga.interface import create_sts_boundary |
---|
| 38 | from anuga.interface import csv2building_polygons |
---|
[6343] | 39 | <<<<<<< .mine |
---|
| 40 | ======= |
---|
[6328] | 41 | from file_length import file_length |
---|
[6259] | 42 | |
---|
[6343] | 43 | >>>>>>> .r6341 |
---|
[6259] | 44 | from anuga.shallow_water.data_manager import start_screen_catcher |
---|
| 45 | from anuga.shallow_water.data_manager import copy_code_files |
---|
[6343] | 46 | from anuga.shallow_water.data_manager import urs2sts |
---|
[6259] | 47 | from anuga.utilities.polygon import read_polygon, Polygon_function |
---|
[6343] | 48 | |
---|
[6253] | 49 | # Application specific imports |
---|
[6324] | 50 | from setup_model import project |
---|
[6253] | 51 | |
---|
[6261] | 52 | |
---|
[6312] | 53 | #------------------------------------------------------------------------------- |
---|
[6343] | 54 | # Get gauges (timeseries of index points) |
---|
| 55 | #------------------------------------------------------------------------------- |
---|
| 56 | def 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. |
---|
| 126 | def 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] | 214 | copy_code_files(project.output_run, __file__, |
---|
[6312] | 215 | os.path.join(os.path.dirname(project.__file__), |
---|
| 216 | project.__name__+'.py')) |
---|
[6279] | 217 | start_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] | 225 | print 'Create computational domain' |
---|
[6253] | 226 | |
---|
[6312] | 227 | # Create the STS file |
---|
[6343] | 228 | <<<<<<< .mine |
---|
| 229 | build_urs_boundary(project.mux_data_folder, |
---|
| 230 | project.mux_input, |
---|
| 231 | os.path.join(project.event_folder, |
---|
| 232 | project.scenario_name)) |
---|
| 233 | ======= |
---|
[6329] | 234 | print 'project.mux_data_folder=%s' % project.mux_data_folder |
---|
| 235 | bub.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] | 241 | event_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] | 245 | landward_boundary = read_polygon(project.landward_boundary) |
---|
[6253] | 246 | |
---|
[6259] | 247 | # Combine sts polyline with landward points |
---|
[6279] | 248 | bounding_polygon_sts = event_sts + landward_boundary |
---|
[6253] | 249 | |
---|
[6259] | 250 | # Number of boundary segments |
---|
[6329] | 251 | num_ocean_segments = len(event_sts) - 1 |
---|
[6328] | 252 | # Number of landward_boundary points |
---|
[6329] | 253 | num_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] | 257 | boundary_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] | 264 | domain = 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) |
---|
| 271 | print domain.statistics() |
---|
[6253] | 272 | |
---|
[6259] | 273 | domain.set_name(project.scenario_name) |
---|
[6279] | 274 | domain.set_datadir(project.output_run) |
---|
[6259] | 275 | domain.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] | 281 | print 'Setup initial conditions' |
---|
[6253] | 282 | |
---|
[6259] | 283 | # Set the initial stage in the offcoast region only |
---|
[6279] | 284 | IC = Polygon_function(project.land_initial_conditions, |
---|
[6259] | 285 | default=project.tide, |
---|
| 286 | geo_reference=domain.geo_reference) |
---|
| 287 | domain.set_quantity('stage', IC, use_cache=True, verbose=True) |
---|
| 288 | domain.set_quantity('friction', project.friction) |
---|
| 289 | domain.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] | 299 | print 'Set boundary - available tags:', domain.get_boundary_tags() |
---|
[6253] | 300 | |
---|
[6259] | 301 | Br = Reflective_boundary(domain) |
---|
[6279] | 302 | Bt = Transmissive_stage_zero_momentum_boundary(domain) |
---|
[6312] | 303 | Bd = Dirichlet_boundary([kwargs['tide'], 0, 0]) |
---|
[6279] | 304 | Bf = 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] | 312 | domain.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] | 320 | t0 = time.time() |
---|
| 321 | for 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] | 327 | print 'Simulation took %.2f seconds' % (time.time()-t0) |
---|