| 1 | """Script for running a tsunami inundation scenario for Cairns, QLD Australia. |
|---|
| 2 | |
|---|
| 3 | Source data such as elevation and boundary data is assumed to be available in |
|---|
| 4 | directories specified by project.py |
|---|
| 5 | The output sww file is stored in directory named after the scenario, i.e |
|---|
| 6 | slide or fixed_wave. |
|---|
| 7 | |
|---|
| 8 | The scenario is defined by a triangular mesh created from project.polygon, |
|---|
| 9 | the elevation data and a tsunami wave generated by a submarine mass failure. |
|---|
| 10 | |
|---|
| 11 | Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton and |
|---|
| 12 | Nick Bartzis, GA - 2006 |
|---|
| 13 | """ |
|---|
| 14 | |
|---|
| 15 | def convert_arcgis_latlon_list_to_utm(points): |
|---|
| 16 | #Used because arc gis produced csv files put lat lon in |
|---|
| 17 | #reverse order to those accpeted by convert_latlon_to_utm() |
|---|
| 18 | |
|---|
| 19 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
|---|
| 20 | from anuga.coordinate_transforms.redfearn import redfearn |
|---|
| 21 | old_geo = Geo_reference() |
|---|
| 22 | utm_points = [] |
|---|
| 23 | for point in points: |
|---|
| 24 | zone, easting, northing = redfearn(float(point[1]),float(point[0])) |
|---|
| 25 | new_geo = Geo_reference(zone) |
|---|
| 26 | old_geo.reconcile_zones(new_geo) |
|---|
| 27 | utm_points.append([easting,northing]) |
|---|
| 28 | |
|---|
| 29 | return utm_points, old_geo.get_zone() |
|---|
| 30 | |
|---|
| 31 | #------------------------------------------------------------------------------ |
|---|
| 32 | # Import necessary modules |
|---|
| 33 | #------------------------------------------------------------------------------ |
|---|
| 34 | |
|---|
| 35 | # Standard modules |
|---|
| 36 | import os |
|---|
| 37 | import time |
|---|
| 38 | import sys |
|---|
| 39 | |
|---|
| 40 | # Related major packages |
|---|
| 41 | from anuga.shallow_water import Domain |
|---|
| 42 | from anuga.shallow_water import Reflective_boundary |
|---|
| 43 | from anuga.shallow_water import Dirichlet_boundary |
|---|
| 44 | from anuga.shallow_water import Time_boundary |
|---|
| 45 | from anuga.shallow_water import File_boundary |
|---|
| 46 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
|---|
| 47 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, ferret2sww |
|---|
| 48 | from anuga.shallow_water.data_manager import dem2pts |
|---|
| 49 | from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm |
|---|
| 50 | #from anuga.fit_interpolate.fit import fit_to_mesh_file (does not exist) |
|---|
| 51 | |
|---|
| 52 | # Application specific imports |
|---|
| 53 | import project # Definition of file names and polygons |
|---|
| 54 | |
|---|
| 55 | |
|---|
| 56 | #------------------------------------------------------------------------------ |
|---|
| 57 | # Define scenario as either slide or fixed_wave. |
|---|
| 58 | #------------------------------------------------------------------------------ |
|---|
| 59 | |
|---|
| 60 | #scenario = 'coseismic' |
|---|
| 61 | scenario = 'fixed_wave' |
|---|
| 62 | |
|---|
| 63 | if os.access(scenario, os.F_OK) == 0: |
|---|
| 64 | os.mkdir(scenario) |
|---|
| 65 | basename = scenario + 'source' |
|---|
| 66 | |
|---|
| 67 | |
|---|
| 68 | #------------------------------------------------------------------------------ |
|---|
| 69 | # Preparation of topographic data |
|---|
| 70 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
|---|
| 71 | #------------------------------------------------------------------------------ |
|---|
| 72 | |
|---|
| 73 | # Filenames |
|---|
| 74 | dem_name = 'boxingday' |
|---|
| 75 | |
|---|
| 76 | # Create DEM from asc data |
|---|
| 77 | #convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True) |
|---|
| 78 | # Create pts file for onshore DEM |
|---|
| 79 | #dem2pts(dem_name, use_cache=True, verbose=True) |
|---|
| 80 | |
|---|
| 81 | |
|---|
| 82 | #------------------------------------------------------------------------------ |
|---|
| 83 | # Create the triangular mesh based on overall clipping polygon with a tagged |
|---|
| 84 | # boundary and interior regions defined in project.py along with |
|---|
| 85 | # resolutions (maximal area of per triangle) for each polygon |
|---|
| 86 | #------------------------------------------------------------------------------ |
|---|
| 87 | |
|---|
| 88 | remainder_res = 10000000 |
|---|
| 89 | islands_res = 100000 |
|---|
| 90 | cairns_res = 100000 |
|---|
| 91 | shallow_res = 500000 |
|---|
| 92 | #interior_regions = [[project.poly_cairns, cairns_res]#, |
|---|
| 93 | # [project.poly_island0, islands_res], |
|---|
| 94 | # [project.poly_island1, islands_res], |
|---|
| 95 | # [project.poly_island2, islands_res], |
|---|
| 96 | # [project.poly_island3, islands_res], |
|---|
| 97 | # [project.poly_shallow, shallow_res]] |
|---|
| 98 | |
|---|
| 99 | # filenames |
|---|
| 100 | meshname = project.meshname + '.tsh' |
|---|
| 101 | mesh_elevname = project.mesh_elevname + '.tsh' |
|---|
| 102 | |
|---|
| 103 | bounding_polygon, zone = convert_arcgis_latlon_list_to_utm(project.bounding_polygon_latlon) |
|---|
| 104 | |
|---|
| 105 | from caching import cache |
|---|
| 106 | print 'start create mesh from regions' |
|---|
| 107 | |
|---|
| 108 | if scenario == 'coseismic': |
|---|
| 109 | _ = cache(create_mesh_from_regions, |
|---|
| 110 | bounding_polygon, |
|---|
| 111 | {'boundary_tags': {'one': [0], 'two': [1], |
|---|
| 112 | 'three': [2], 'four': [3], |
|---|
| 113 | 'five': [4], 'six': [5], |
|---|
| 114 | 'seven': [6], 'eight': [7], |
|---|
| 115 | 'nine': [8], 'ten': [9], |
|---|
| 116 | 'eleven': [10], 'twelve': [11], |
|---|
| 117 | 'thirteen': [12], 'fourteen': [13]}, |
|---|
| 118 | 'maximum_triangle_area': 5000000, |
|---|
| 119 | 'filename': meshname},#, |
|---|
| 120 | #'interior_regions': interior_regions}, |
|---|
| 121 | verbose = True |
|---|
| 122 | #, evaluate=True |
|---|
| 123 | ) |
|---|
| 124 | |
|---|
| 125 | if scenario == 'fixed_wave': |
|---|
| 126 | south = 7.40911081272 |
|---|
| 127 | north = 8.71484358635 |
|---|
| 128 | east = 99.1683687224 |
|---|
| 129 | west = 97.513856322 |
|---|
| 130 | points = [[south,west],[south,east],[north,east],[north,west]] |
|---|
| 131 | |
|---|
| 132 | bounding_polygon,zone=convert_from_latlon_to_utm(points) |
|---|
| 133 | |
|---|
| 134 | _ = cache(create_mesh_from_regions, |
|---|
| 135 | bounding_polygon, |
|---|
| 136 | {'boundary_tags': {'ocean_west': [0], 'bottom': [1], |
|---|
| 137 | 'onshore': [2], 'top': [3]}, |
|---|
| 138 | 'maximum_triangle_area': 500000000, |
|---|
| 139 | 'filename': meshname},#, |
|---|
| 140 | #'interior_regions': interior_regions}, |
|---|
| 141 | verbose = True |
|---|
| 142 | #, evaluate=True |
|---|
| 143 | ) |
|---|
| 144 | |
|---|
| 145 | ## cache(fit_to_mesh_file,(meshname, |
|---|
| 146 | ## project.combined_dir_name + '.pts', |
|---|
| 147 | ## mesh_elevname), |
|---|
| 148 | ## #{'verbose': True}, |
|---|
| 149 | ## verbose = False |
|---|
| 150 | ## ) |
|---|
| 151 | |
|---|
| 152 | #------------------------------------------------------------------------------ |
|---|
| 153 | # Setup computational domain |
|---|
| 154 | #------------------------------------------------------------------------------ |
|---|
| 155 | |
|---|
| 156 | domain = Domain(meshname, use_cache=False, verbose=True) |
|---|
| 157 | |
|---|
| 158 | print 'Number of triangles = ', len(domain) |
|---|
| 159 | print 'The extent is ', domain.get_extent() |
|---|
| 160 | print domain.statistics() |
|---|
| 161 | |
|---|
| 162 | domain.set_name(basename) |
|---|
| 163 | domain.set_datadir(scenario) |
|---|
| 164 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
|---|
| 165 | domain.set_minimum_storable_height(0.01) |
|---|
| 166 | |
|---|
| 167 | |
|---|
| 168 | #------------------------------------------------------------------------------ |
|---|
| 169 | # Setup initial conditions |
|---|
| 170 | #------------------------------------------------------------------------------ |
|---|
| 171 | |
|---|
| 172 | def elevation(x,y): |
|---|
| 173 | return -10.0 |
|---|
| 174 | |
|---|
| 175 | tide = 0.0 |
|---|
| 176 | domain.set_quantity('stage', tide) |
|---|
| 177 | domain.set_quantity('friction', 0.0) |
|---|
| 178 | |
|---|
| 179 | if scenario == 'fixed_wave': |
|---|
| 180 | # test with coarser bathymetry |
|---|
| 181 | domain.set_quantity('elevation', |
|---|
| 182 | #filename = 'thaicoas_9.pts', |
|---|
| 183 | elevation, |
|---|
| 184 | use_cache=True, |
|---|
| 185 | verbose=True, |
|---|
| 186 | alpha=0.1) |
|---|
| 187 | |
|---|
| 188 | if scenario == 'coseismic': |
|---|
| 189 | domain.set_quantity('elevation', |
|---|
| 190 | filename=project.combined_dir_name + '.pts', |
|---|
| 191 | use_cache=True, |
|---|
| 192 | verbose=True, |
|---|
| 193 | alpha=0.1) |
|---|
| 194 | #------------------------------------------------------------------------------ |
|---|
| 195 | # Setup boundary conditions |
|---|
| 196 | #------------------------------------------------------------------------------ |
|---|
| 197 | |
|---|
| 198 | print 'start ferret2sww' |
|---|
| 199 | south = project.south |
|---|
| 200 | north = project.north |
|---|
| 201 | west = project.west |
|---|
| 202 | east = project.east |
|---|
| 203 | |
|---|
| 204 | """ |
|---|
| 205 | tide = 0.0 |
|---|
| 206 | # Note only need to do when an SWW file for the MOST boundary doesn't exist |
|---|
| 207 | cache(ferret2sww, |
|---|
| 208 | (project.boundary_most_in, |
|---|
| 209 | project.boundary_most_out), |
|---|
| 210 | {'verbose': True, |
|---|
| 211 | 'minlat': south, |
|---|
| 212 | 'maxlat': north, |
|---|
| 213 | 'minlon': west, |
|---|
| 214 | 'maxlon': east, |
|---|
| 215 | # 'origin': project.mesh_origin, |
|---|
| 216 | 'origin': domain.geo_reference.get_origin(), |
|---|
| 217 | 'mean_stage': tide, |
|---|
| 218 | 'zscale': 10, #Enhance tsunami (Does this affect result) |
|---|
| 219 | 'fail_on_NaN': False, |
|---|
| 220 | 'inverted_bathymetry': True}, |
|---|
| 221 | #evaluate = True, |
|---|
| 222 | verbose = True)#, |
|---|
| 223 | #dependencies = 'most_results' + '.sww') |
|---|
| 224 | |
|---|
| 225 | """ |
|---|
| 226 | print 'Available boundary tags', domain.get_boundary_tags() |
|---|
| 227 | |
|---|
| 228 | #Bf = File_boundary('most_results'+'.sww', |
|---|
| 229 | # domain, verbose = True) |
|---|
| 230 | Br = Reflective_boundary(domain) |
|---|
| 231 | Bd = Dirichlet_boundary([tide,0,0]) |
|---|
| 232 | Bw = Time_boundary(domain = domain, |
|---|
| 233 | f=lambda t: [(60<t<3660)*50, 0, 0]) |
|---|
| 234 | |
|---|
| 235 | # 60 min square wave starting at 1 min, 50m high |
|---|
| 236 | if scenario == 'fixed_wave': |
|---|
| 237 | domain.set_boundary({'ocean_west': Br, |
|---|
| 238 | 'bottom': Br, |
|---|
| 239 | 'onshore': Br, |
|---|
| 240 | 'top': Br}) |
|---|
| 241 | |
|---|
| 242 | # boundary conditions for slide scenario |
|---|
| 243 | if scenario == 'coseismic': |
|---|
| 244 | domain.set_boundary( {'one': Bf,'two': Bf, 'three': Bf, 'four': Bf, |
|---|
| 245 | 'five': Bf, 'six': Bd, 'seven': Bd, 'eight': Bf, |
|---|
| 246 | 'nine': Bf, 'ten': Bd, 'eleven': Bd, 'twelve': Bf, |
|---|
| 247 | 'thireteen': Bf, 'fourteen': Bf}) |
|---|
| 248 | |
|---|
| 249 | |
|---|
| 250 | #------------------------------------------------------------------------------ |
|---|
| 251 | # Evolve system through time |
|---|
| 252 | #------------------------------------------------------------------------------ |
|---|
| 253 | #domain.visualise = True |
|---|
| 254 | #domain.visualise_color_stage = True |
|---|
| 255 | |
|---|
| 256 | import time |
|---|
| 257 | t0 = time.time() |
|---|
| 258 | |
|---|
| 259 | from Numeric import allclose |
|---|
| 260 | from anuga.abstract_2d_finite_volumes.quantity import Quantity |
|---|
| 261 | |
|---|
| 262 | if scenario == 'cosesimic': |
|---|
| 263 | |
|---|
| 264 | for t in domain.evolve(yieldstep = 10, finaltime = 60): |
|---|
| 265 | print domain.quantities['elevation'].vertex_values |
|---|
| 266 | domain.write_time() |
|---|
| 267 | |
|---|
| 268 | for t in domain.evolve(yieldstep = 600, finaltime = 21600, |
|---|
| 269 | skip_initial_step = True): |
|---|
| 270 | domain.write_time() |
|---|
| 271 | domain.write_boundary_statistics(tags = 'ocean_west') |
|---|
| 272 | |
|---|
| 273 | if scenario == 'fixed_wave': |
|---|
| 274 | |
|---|
| 275 | # save every two mins leading up to wave approaching land |
|---|
| 276 | for t in domain.evolve(yieldstep = 10, finaltime = 50): |
|---|
| 277 | print domain.quantities['elevation'].vertex_values |
|---|
| 278 | print domain.quantities['stage'].vertex_values |
|---|
| 279 | |
|---|
| 280 | domain.write_time() |
|---|
| 281 | |
|---|
| 282 | #save every 30 secs as wave starts inundating ashore |
|---|
| 283 | #for t in domain.evolve(yieldstep = 10, finaltime = 10000, |
|---|
| 284 | # skip_initial_step = True): |
|---|
| 285 | # domain.write_time() |
|---|
| 286 | |
|---|
| 287 | print 'That took %.2f seconds' %(time.time()-t0) |
|---|
