[5393] | 1 | """Script for running the 2004 boxing Day tsunami inundation scenario for |
---|
| 2 | Phuket, Thailand. |
---|
| 3 | |
---|
| 4 | Source data such as elevation and boundary data is assumed to be available in |
---|
| 5 | directories specified by project.py |
---|
| 6 | The output sww file is stored in directory named after the scenario, i.e |
---|
| 7 | slide or fixed_wave. |
---|
| 8 | |
---|
| 9 | The scenario is defined by a triangular mesh created from project.polygon, |
---|
| 10 | the elevation data and a tsunami wave generated by a submarine mass failure. |
---|
| 11 | |
---|
| 12 | Author: John Jakeman, The Australian National University (2008) |
---|
| 13 | """ |
---|
| 14 | |
---|
| 15 | #------------------------------------------------------------------------------ |
---|
| 16 | # Import necessary modules |
---|
| 17 | #------------------------------------------------------------------------------ |
---|
| 18 | |
---|
| 19 | # Standard modules |
---|
| 20 | import os |
---|
| 21 | import time |
---|
| 22 | import sys |
---|
| 23 | from time import localtime, strftime |
---|
| 24 | |
---|
| 25 | # Related major packages |
---|
| 26 | from anuga.shallow_water import Domain |
---|
| 27 | from anuga.shallow_water import Reflective_boundary |
---|
| 28 | from anuga.shallow_water import Dirichlet_boundary |
---|
| 29 | from anuga.shallow_water import Time_boundary |
---|
| 30 | from anuga.shallow_water import File_boundary |
---|
| 31 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
| 32 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, ferret2sww |
---|
| 33 | from anuga.shallow_water.data_manager import dem2pts |
---|
| 34 | from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm |
---|
| 35 | from anuga.utilities.polygon import number_mesh_triangles |
---|
| 36 | from anuga.fit_interpolate.fit import fit_to_mesh_file |
---|
| 37 | from anuga.caching import cache |
---|
| 38 | from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance |
---|
| 39 | |
---|
| 40 | # Application specific imports |
---|
| 41 | import project # Definition of file names and polygons |
---|
| 42 | |
---|
| 43 | |
---|
| 44 | #------------------------------------------------------------------------------ |
---|
| 45 | # Define scenario as either slide or fixed_wave. |
---|
| 46 | #------------------------------------------------------------------------------ |
---|
| 47 | |
---|
| 48 | #scenario = 'poor_simulation' |
---|
| 49 | scenario = 'good_simulation' |
---|
| 50 | |
---|
| 51 | if os.access(scenario, os.F_OK) == 0: |
---|
| 52 | os.mkdir(scenario) |
---|
| 53 | |
---|
| 54 | timestamp = strftime('%Y%m%d_%H%M%S',localtime()) |
---|
[5401] | 55 | basename = scenario[:4] + '_polyline' |
---|
[5393] | 56 | |
---|
| 57 | |
---|
| 58 | #------------------------------------------------------------------------------ |
---|
| 59 | # Preparation of topographic data |
---|
| 60 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
| 61 | #------------------------------------------------------------------------------ |
---|
| 62 | |
---|
| 63 | if scenario == 'good_simualation': |
---|
| 64 | msg = 'Must use combine_good_data to create bathymetry .pts file' |
---|
| 65 | assert os.path.exists(project.good_combined_dir_name+'.pts'), msg |
---|
| 66 | |
---|
| 67 | #------------------------------------------------------------------------------ |
---|
| 68 | # Create the triangular mesh based on overall clipping polygon with a tagged |
---|
| 69 | # boundary and interior regions defined in project.py along with |
---|
| 70 | # resolutions (maximal area of per triangle) for each polygon |
---|
| 71 | #------------------------------------------------------------------------------ |
---|
[5401] | 72 | extent_res = 1000000.0 |
---|
[5393] | 73 | contour20m_res = 50000.0 |
---|
| 74 | island_res = 5000.0 |
---|
| 75 | bay_res = 2000.0 |
---|
[5401] | 76 | patong_res = 400.0 |
---|
[5393] | 77 | |
---|
| 78 | |
---|
| 79 | # make polygon that contains land that does not affect result. |
---|
| 80 | |
---|
| 81 | interior_regions = [[project.patong, patong_res], |
---|
| 82 | [project.bay, bay_res], |
---|
[5401] | 83 | [project.contour20m, contour20m_res]]#, |
---|
| 84 | #[project.island_north, island_res], |
---|
| 85 | #[project.island_south, island_res], |
---|
| 86 | #[project.island_south2, island_res]] |
---|
[5393] | 87 | |
---|
| 88 | |
---|
| 89 | #for coarse run to test gauges |
---|
[5401] | 90 | #extent_res = 10000000.0 |
---|
| 91 | #interior_regions=None |
---|
[5393] | 92 | |
---|
| 93 | from Numeric import arange,allclose |
---|
[5401] | 94 | boundary_tags={'ocean': arange(0,41).tolist(), 'otherocean': [41,44], 'land': [42,43]} |
---|
[5393] | 95 | |
---|
[5401] | 96 | #trigs_min = number_mesh_triangles(interior_regions, project.bounding_polygon,extent_res) |
---|
[5393] | 97 | #print 'Minimum number of traingles ', trigs_min |
---|
| 98 | |
---|
| 99 | # filenames |
---|
| 100 | meshname = project.meshname + '_polyline.tsh' |
---|
| 101 | mesh_elevname = project.mesh_elevname + '_polyline.tsh' |
---|
| 102 | |
---|
| 103 | print 'start create mesh from regions' |
---|
| 104 | |
---|
| 105 | _ = cache(create_mesh_from_regions, |
---|
[5401] | 106 | project.bounding_polygon, |
---|
[5393] | 107 | {'boundary_tags': boundary_tags, |
---|
| 108 | 'maximum_triangle_area': extent_res, |
---|
| 109 | 'filename': meshname, |
---|
| 110 | 'interior_regions': interior_regions}, |
---|
| 111 | verbose = True, |
---|
[5401] | 112 | dependencies = ['project.py'] |
---|
[5393] | 113 | #, evaluate=True |
---|
| 114 | ) |
---|
| 115 | |
---|
| 116 | |
---|
| 117 | #------------------------------------------------------------------------------ |
---|
| 118 | # Setup computational domain |
---|
| 119 | #------------------------------------------------------------------------------ |
---|
| 120 | print 'Converting mesh to domain' |
---|
| 121 | |
---|
| 122 | #domain = Domain(meshname, use_cache=False, verbose=True) |
---|
| 123 | |
---|
| 124 | domain = cache(pmesh_to_domain_instance, |
---|
| 125 | (meshname, Domain), |
---|
| 126 | dependencies = [meshname]) |
---|
| 127 | |
---|
| 128 | print 'The extent is ', domain.get_extent() |
---|
| 129 | print domain.statistics() |
---|
| 130 | |
---|
[5401] | 131 | domain.set_name(basename+timestamp) |
---|
[5393] | 132 | #domain.set_name(basename) |
---|
| 133 | domain.set_datadir(scenario) |
---|
| 134 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
| 135 | domain.set_minimum_storable_height(0.01) |
---|
| 136 | |
---|
| 137 | domain.beta_h = 0 |
---|
| 138 | domain.tight_slope_limiters = 1 |
---|
| 139 | domain.set_default_order(2) |
---|
| 140 | print 'domain.tight_slope_limiters', domain.tight_slope_limiters |
---|
| 141 | |
---|
| 142 | domain.points_file_block_line_size = 50000 |
---|
| 143 | |
---|
| 144 | #------------------------------------------------------------------------------ |
---|
| 145 | # Setup initial conditions |
---|
| 146 | #------------------------------------------------------------------------------ |
---|
| 147 | |
---|
| 148 | tide = 0.0 |
---|
| 149 | domain.set_quantity('stage', tide) |
---|
| 150 | domain.set_quantity('friction', 0.01) |
---|
| 151 | |
---|
| 152 | if scenario == 'poor_simulation': |
---|
| 153 | domain.set_quantity('elevation', |
---|
| 154 | filename=project.poor_combined_dir_name + '.pts', |
---|
| 155 | use_cache=True, |
---|
| 156 | verbose=True, |
---|
| 157 | alpha=0.1) |
---|
| 158 | |
---|
| 159 | if scenario == 'good_simulation': |
---|
| 160 | domain.set_quantity('elevation', |
---|
| 161 | filename=project.good_combined_dir_name + '.pts', |
---|
| 162 | use_cache=True, |
---|
| 163 | verbose=True, |
---|
| 164 | alpha=0.1) |
---|
| 165 | #------------------------------------------------------------------------------ |
---|
| 166 | # Setup boundary conditions |
---|
| 167 | #------------------------------------------------------------------------------ |
---|
| 168 | boundary_urs_in='data/boxing' |
---|
[5401] | 169 | boundary_urs_out=project.base_name |
---|
[5393] | 170 | |
---|
| 171 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
| 172 | Bf = File_boundary(boundary_urs_out+'.sts', |
---|
| 173 | domain, time_thinning=1, |
---|
[5401] | 174 | use_cache=True,verbose = True,boundary_polygon=project.bounding_polygon) |
---|
[5393] | 175 | |
---|
| 176 | Br = Reflective_boundary(domain) |
---|
| 177 | Bd = Dirichlet_boundary([tide,0.0,0.0]) |
---|
| 178 | |
---|
| 179 | domain.set_boundary({'ocean': Bf, |
---|
| 180 | 'otherocean': Bd, |
---|
| 181 | 'land': Br, |
---|
| 182 | 'both': Bd}) |
---|
| 183 | |
---|
| 184 | #------------------------------------------------------------------------------ |
---|
| 185 | # Evolve system through time |
---|
| 186 | #------------------------------------------------------------------------------ |
---|
| 187 | import time |
---|
| 188 | t0 = time.time() |
---|
| 189 | |
---|
| 190 | #from Numeric import allclose |
---|
| 191 | #from anuga.abstract_2d_finite_volumes.quantity import Quantity |
---|
| 192 | |
---|
| 193 | # Add new loop that uses larger yieldstep until wave first reaches a point of |
---|
| 194 | # the ANUGA boundary. Or find a way to clip MOST sww boundary file to only |
---|
| 195 | # start when boundary stage first becomes non-zero. |
---|
| 196 | |
---|
[5401] | 197 | for t in domain.evolve(yieldstep = 20.0, finaltime = 18000.0, |
---|
[5393] | 198 | skip_initial_step = False): |
---|
| 199 | if allclose(t,10800.0): |
---|
| 200 | print 'Changing urs file boundary to dirichlet. Urs gauges only have 3 hours of data' |
---|
| 201 | domain.set_boundary({'ocean': Bd, |
---|
| 202 | 'otherocean': Bd, |
---|
| 203 | 'land': Br, |
---|
| 204 | 'both': Bd}) |
---|
| 205 | |
---|
| 206 | domain.write_time() |
---|
| 207 | |
---|
| 208 | print 'That took %.2f seconds' %(time.time()-t0) |
---|