[2461] | 1 | """Script for running a tsunami inundation scenario for Sydney, NSW, 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 project.outputdir |
---|
| 6 | |
---|
| 7 | The scenario is defined by a triangular mesh created from project.polygon, |
---|
| 8 | the elevation data and a simulated submarine landslide. |
---|
| 9 | |
---|
| 10 | """ |
---|
| 11 | |
---|
| 12 | |
---|
[2887] | 13 | #------------------------------------------------------------------------------ |
---|
[2478] | 14 | # Import necessary modules |
---|
[2887] | 15 | #------------------------------------------------------------------------------ |
---|
[2461] | 16 | |
---|
| 17 | # Standard modules |
---|
| 18 | import os |
---|
| 19 | import time |
---|
| 20 | |
---|
| 21 | # Related major packages |
---|
[7808] | 22 | import anuga |
---|
[2461] | 23 | |
---|
| 24 | # Application specific imports |
---|
[2887] | 25 | import project # Define file names and polygons |
---|
[2461] | 26 | |
---|
| 27 | |
---|
| 28 | |
---|
[2887] | 29 | #------------------------------------------------------------------------------ |
---|
[2599] | 30 | # Prepare bathymetric and topographic data |
---|
[2887] | 31 | #------------------------------------------------------------------------------ |
---|
[2478] | 32 | |
---|
[2461] | 33 | # filenames |
---|
[3869] | 34 | #coarsedemname = project.coarsedemname + '.pts' |
---|
| 35 | #finedemname = project.finedemname + '.pts' |
---|
| 36 | demname = project.demname |
---|
[2461] | 37 | meshname = project.meshname+'.msh' |
---|
| 38 | |
---|
[7808] | 39 | anuga.asc2dem(demname, use_cache=True, verbose=True) |
---|
[2461] | 40 | |
---|
[3869] | 41 | # creates pts file |
---|
[7808] | 42 | anuga.dem2pts(demname, use_cache=True, verbose=True) |
---|
[3869] | 43 | |
---|
[2461] | 44 | |
---|
[2887] | 45 | #------------------------------------------------------------------------------ |
---|
[2478] | 46 | # Setup computational domain |
---|
[2887] | 47 | #------------------------------------------------------------------------------ |
---|
[2478] | 48 | |
---|
| 49 | # Interior regions and mesh resolutions |
---|
[2461] | 50 | interior_res = 5000 |
---|
[3944] | 51 | shallow_res = 15000 |
---|
| 52 | ##interior_regions = [[project.northern_polygon, interior_res], |
---|
| 53 | ## [project.harbour_polygon, interior_res], |
---|
| 54 | ## [project.southern_polygon, interior_res], |
---|
| 55 | ## [project.manly_polygon, interior_res], |
---|
| 56 | ## [project.top_polygon, interior_res]] |
---|
| 57 | interior_regions = [[project.coastal_polygon, interior_res], |
---|
| 58 | [project.shallow_polygon, shallow_res]] |
---|
[2461] | 59 | |
---|
[7808] | 60 | anuga.create_mesh_from_regions(project.demopoly, |
---|
[3136] | 61 | boundary_tags= {'oceannorth': [0], |
---|
[3190] | 62 | 'onshorenorth1': [1], |
---|
| 63 | 'onshorenorthwest1': [2], |
---|
| 64 | 'onshorenorthwest2': [3], |
---|
| 65 | 'onshorenorth2': [4], |
---|
| 66 | 'onshorewest': [5], |
---|
| 67 | 'onshoresouth': [6], |
---|
| 68 | 'oceansouth': [7], |
---|
| 69 | 'oceaneast': [8]}, |
---|
[2478] | 70 | maximum_triangle_area=100000, |
---|
| 71 | filename=meshname, |
---|
| 72 | interior_regions=interior_regions) |
---|
[2461] | 73 | |
---|
| 74 | |
---|
[2805] | 75 | #Create shallow water domain |
---|
[2461] | 76 | |
---|
[7808] | 77 | domain = anuga.Domain(meshname, |
---|
[2866] | 78 | use_cache = True, |
---|
| 79 | verbose = True) |
---|
[2805] | 80 | |
---|
[2866] | 81 | |
---|
[2461] | 82 | print 'Number of triangles = ', len(domain) |
---|
| 83 | print 'The extent is ', domain.get_extent() |
---|
| 84 | |
---|
| 85 | domain.set_name(project.basename) |
---|
| 86 | domain.set_datadir(project.outputdir) |
---|
| 87 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
| 88 | |
---|
| 89 | |
---|
[2887] | 90 | #------------------------------------------------------------------------------ |
---|
[2461] | 91 | # Set up scenario (tsunami_source is a callable object used with set_quantity) |
---|
[2887] | 92 | #------------------------------------------------------------------------------ |
---|
[2461] | 93 | |
---|
[7808] | 94 | tsunami_source = anuga.slump_tsunami(length=30000.0, |
---|
[2461] | 95 | depth=400.0, |
---|
| 96 | slope=6.0, |
---|
| 97 | thickness=176.0, |
---|
| 98 | radius=3330, |
---|
| 99 | dphi=0.23, |
---|
| 100 | x0=project.slump_origin[0], |
---|
| 101 | y0=project.slump_origin[1], |
---|
| 102 | alpha=0.0, |
---|
| 103 | domain=domain) |
---|
| 104 | |
---|
| 105 | |
---|
[2887] | 106 | #------------------------------------------------------------------------------ |
---|
[2461] | 107 | # Setup initial conditions |
---|
[2887] | 108 | #------------------------------------------------------------------------------ |
---|
[2461] | 109 | |
---|
| 110 | domain.set_quantity('stage', tsunami_source) |
---|
| 111 | domain.set_quantity('friction', 0.0) |
---|
| 112 | domain.set_quantity('elevation', |
---|
[3869] | 113 | filename = demname + '.pts', |
---|
[2461] | 114 | use_cache = True, |
---|
| 115 | verbose = True) |
---|
| 116 | |
---|
[2887] | 117 | #------------------------------------------------------------------------------ |
---|
[3136] | 118 | # Setup boundary conditions (all Dirichlet) |
---|
[2887] | 119 | #------------------------------------------------------------------------------ |
---|
[2461] | 120 | |
---|
| 121 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
| 122 | |
---|
[7808] | 123 | Bd = anuga.Dirichlet_boundary([0.0,0.0,0.0]) |
---|
[3190] | 124 | domain.set_boundary( {'oceannorth': Bd, |
---|
| 125 | 'onshorenorth1': Bd, |
---|
| 126 | 'onshorenorthwest1': Bd, |
---|
| 127 | 'onshorenorthwest2': Bd, |
---|
| 128 | 'onshorenorth2': Bd, |
---|
| 129 | 'onshorewest': Bd, |
---|
| 130 | 'onshoresouth': Bd, |
---|
| 131 | 'oceansouth': Bd, |
---|
| 132 | 'oceaneast': Bd} ) |
---|
[2461] | 133 | |
---|
[2887] | 134 | #------------------------------------------------------------------------------ |
---|
[2461] | 135 | # Evolve system through time |
---|
[2887] | 136 | #------------------------------------------------------------------------------ |
---|
[2461] | 137 | |
---|
| 138 | import time |
---|
| 139 | t0 = time.time() |
---|
| 140 | |
---|
[3190] | 141 | for t in domain.evolve(yieldstep = 60, finaltime = 7200): |
---|
[2497] | 142 | print domain.timestepping_statistics() |
---|
[3136] | 143 | print domain.boundary_statistics(tags = 'oceaneast') |
---|
[2461] | 144 | |
---|
| 145 | print 'That took %.2f seconds' %(time.time()-t0) |
---|