| 1 | """Run a tsunami inundation scenario for Busselton, WA, Australia. |
|---|
| 2 | |
|---|
| 3 | The scenario is defined by a triangular mesh created from project.polygon, the |
|---|
| 4 | elevation data is compiled into a pts file through build_elevation.py and a |
|---|
| 5 | simulated tsunami is generated through an sts file from build_boundary.py. |
|---|
| 6 | |
|---|
| 7 | Input: sts file (build_boundary.py for respective event) |
|---|
| 8 | pts file (build_elevation.py) |
|---|
| 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 |
|---|
| 24 | import os.path |
|---|
| 25 | import time |
|---|
| 26 | from time import localtime, strftime, gmtime |
|---|
| 27 | |
|---|
| 28 | # Related major packages |
|---|
| 29 | from Scientific.IO.NetCDF import NetCDFFile |
|---|
| 30 | import Numeric as num |
|---|
| 31 | |
|---|
| 32 | from anuga.interface import create_domain_from_regions |
|---|
| 33 | from anuga.interface import Transmissive_stage_zero_momentum_boundary |
|---|
| 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 |
|---|
| 39 | <<<<<<< .mine |
|---|
| 40 | ======= |
|---|
| 41 | from file_length import file_length |
|---|
| 42 | |
|---|
| 43 | >>>>>>> .r6341 |
|---|
| 44 | from anuga.shallow_water.data_manager import start_screen_catcher |
|---|
| 45 | from anuga.shallow_water.data_manager import copy_code_files |
|---|
| 46 | from anuga.shallow_water.data_manager import urs2sts |
|---|
| 47 | from anuga.utilities.polygon import read_polygon, Polygon_function |
|---|
| 48 | |
|---|
| 49 | # Application specific imports |
|---|
| 50 | from setup_model import project |
|---|
| 51 | |
|---|
| 52 | |
|---|
| 53 | #------------------------------------------------------------------------------- |
|---|
| 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 | #------------------------------------------------------------------------------- |
|---|
| 210 | # Copy scripts to time stamped output directory and capture screen |
|---|
| 211 | # output to file. Copy script must be before screen_catcher |
|---|
| 212 | #------------------------------------------------------------------------------- |
|---|
| 213 | |
|---|
| 214 | copy_code_files(project.output_run, __file__, |
|---|
| 215 | os.path.join(os.path.dirname(project.__file__), |
|---|
| 216 | project.__name__+'.py')) |
|---|
| 217 | start_screen_catcher(project.output_run, 0, 1) |
|---|
| 218 | |
|---|
| 219 | #------------------------------------------------------------------------------- |
|---|
| 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 |
|---|
| 223 | #------------------------------------------------------------------------------- |
|---|
| 224 | |
|---|
| 225 | print 'Create computational domain' |
|---|
| 226 | |
|---|
| 227 | # Create the STS file |
|---|
| 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 | ======= |
|---|
| 234 | print 'project.mux_data_folder=%s' % project.mux_data_folder |
|---|
| 235 | bub.build_urs_boundary(project.mux_input_filename, |
|---|
| 236 | os.path.join(project.event_folder, |
|---|
| 237 | project.scenario_name)) |
|---|
| 238 | >>>>>>> .r6341 |
|---|
| 239 | |
|---|
| 240 | # Read in boundary from ordered sts file |
|---|
| 241 | event_sts = create_sts_boundary(project.event_sts) |
|---|
| 242 | |
|---|
| 243 | # Reading the landward defined points, this incorporates the original clipping |
|---|
| 244 | # polygon minus the 100m contour |
|---|
| 245 | landward_boundary = read_polygon(project.landward_boundary) |
|---|
| 246 | |
|---|
| 247 | # Combine sts polyline with landward points |
|---|
| 248 | bounding_polygon_sts = event_sts + landward_boundary |
|---|
| 249 | |
|---|
| 250 | # Number of boundary segments |
|---|
| 251 | num_ocean_segments = len(event_sts) - 1 |
|---|
| 252 | # Number of landward_boundary points |
|---|
| 253 | num_land_points = file_length(project.landward_boundary) |
|---|
| 254 | |
|---|
| 255 | # Boundary tags refer to project.landward_boundary |
|---|
| 256 | # 4 points equals 5 segments start at N |
|---|
| 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)} |
|---|
| 262 | |
|---|
| 263 | # Build mesh and domain |
|---|
| 264 | domain = create_domain_from_regions(bounding_polygon_sts, |
|---|
| 265 | boundary_tags=boundary_tags, |
|---|
| 266 | maximum_triangle_area=project.bounding_maxarea, |
|---|
| 267 | interior_regions=project.interior_regions, |
|---|
| 268 | mesh_filename=project.meshes, |
|---|
| 269 | use_cache=True, |
|---|
| 270 | verbose=True) |
|---|
| 271 | print domain.statistics() |
|---|
| 272 | |
|---|
| 273 | domain.set_name(project.scenario_name) |
|---|
| 274 | domain.set_datadir(project.output_run) |
|---|
| 275 | domain.set_minimum_storable_height(0.01) # Don't store depth less than 1cm |
|---|
| 276 | |
|---|
| 277 | #------------------------------------------------------------------------------- |
|---|
| 278 | # Setup initial conditions |
|---|
| 279 | #------------------------------------------------------------------------------- |
|---|
| 280 | |
|---|
| 281 | print 'Setup initial conditions' |
|---|
| 282 | |
|---|
| 283 | # Set the initial stage in the offcoast region only |
|---|
| 284 | IC = Polygon_function(project.land_initial_conditions, |
|---|
| 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', |
|---|
| 290 | filename=project.combined_elevation+'.pts', |
|---|
| 291 | use_cache=True, |
|---|
| 292 | verbose=True, |
|---|
| 293 | alpha=project.alpha) |
|---|
| 294 | |
|---|
| 295 | #------------------------------------------------------------------------------- |
|---|
| 296 | # Setup boundary conditions |
|---|
| 297 | #------------------------------------------------------------------------------- |
|---|
| 298 | |
|---|
| 299 | print 'Set boundary - available tags:', domain.get_boundary_tags() |
|---|
| 300 | |
|---|
| 301 | Br = Reflective_boundary(domain) |
|---|
| 302 | Bt = Transmissive_stage_zero_momentum_boundary(domain) |
|---|
| 303 | Bd = Dirichlet_boundary([kwargs['tide'], 0, 0]) |
|---|
| 304 | Bf = Field_boundary(project.event_sts+'.sts', |
|---|
| 305 | domain, mean_stage=project.tide, |
|---|
| 306 | time_thinning=1, |
|---|
| 307 | default_boundary=Bd, |
|---|
| 308 | boundary_polygon=bounding_polygon_sts, |
|---|
| 309 | use_cache=True, |
|---|
| 310 | verbose=True) |
|---|
| 311 | |
|---|
| 312 | domain.set_boundary({'back': Br, |
|---|
| 313 | 'side': Bt, |
|---|
| 314 | 'ocean': Bf}) |
|---|
| 315 | |
|---|
| 316 | #------------------------------------------------------------------------------- |
|---|
| 317 | # Evolve system through time |
|---|
| 318 | #------------------------------------------------------------------------------- |
|---|
| 319 | |
|---|
| 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') |
|---|
| 326 | |
|---|
| 327 | print 'Simulation took %.2f seconds' % (time.time()-t0) |
|---|
