[3584] | 1 | """Script for running a tsunami inundation scenario for Flagstaff pt, |
---|
| 2 | Wollongong harbour, NSW, Australia. |
---|
| 3 | |
---|
| 4 | Source data such as elevation and boundary data is assumed to be available in |
---|
| 5 | directories specified by project.py |
---|
| 6 | |
---|
| 7 | The scenario is defined by a triangular mesh created from project.polygon, |
---|
| 8 | the elevation data and a hypothetical boundary condition. |
---|
| 9 | |
---|
| 10 | THIS IS THE SEQUENTIAL VERSION |
---|
| 11 | |
---|
| 12 | Ole Nielsen and Duncan Gray, GA - 2005, Nick Bartzis and Jane Sexton, GA - 2006 |
---|
| 13 | """ |
---|
| 14 | |
---|
| 15 | |
---|
| 16 | #------------------------------------------------------------------------------ |
---|
| 17 | # Import necessary modules |
---|
| 18 | #------------------------------------------------------------------------------ |
---|
| 19 | |
---|
| 20 | |
---|
| 21 | # Standard modules |
---|
| 22 | import os, sys, time |
---|
| 23 | from os import sep |
---|
| 24 | from os.path import dirname, basename |
---|
| 25 | from Numeric import zeros, Float |
---|
| 26 | |
---|
| 27 | # Related major packages |
---|
| 28 | from anuga.pyvolution.shallow_water import Domain |
---|
| 29 | from anuga.pyvolution.shallow_water import Dirichlet_boundary |
---|
| 30 | from anuga.pyvolution.shallow_water import Time_boundary |
---|
| 31 | from anuga.pyvolution.shallow_water import Reflective_boundary |
---|
| 32 | from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts |
---|
| 33 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
| 34 | from anuga.pmesh.mesh import importUngenerateFile, Segment |
---|
| 35 | |
---|
| 36 | |
---|
| 37 | # Application specific imports |
---|
| 38 | import project |
---|
| 39 | |
---|
| 40 | |
---|
| 41 | #------------------------------------------------------------------------------ |
---|
| 42 | # Preparation of topographic data |
---|
| 43 | # |
---|
| 44 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
| 45 | #------------------------------------------------------------------------------ |
---|
| 46 | |
---|
| 47 | max_area = project.base_resolution |
---|
| 48 | |
---|
| 49 | |
---|
| 50 | #-------------------------------------------------------------------------- |
---|
| 51 | # Create the triangular mesh based on overall clipping polygon with a |
---|
| 52 | # tagged boundary and interior regions defined in project.py along with |
---|
| 53 | # resolutions (maximal area of per triangle) for each polygon |
---|
| 54 | #-------------------------------------------------------------------------- |
---|
| 55 | |
---|
| 56 | |
---|
| 57 | print 'Generate mesh' |
---|
| 58 | mesh = create_mesh_from_regions(project.bounding_polygon, |
---|
| 59 | boundary_tags=project.boundary_tags, |
---|
| 60 | maximum_triangle_area=max_area, |
---|
| 61 | interior_regions=project.interior_regions) |
---|
| 62 | |
---|
| 63 | # Add buildings |
---|
| 64 | # This should bind to a Reflective boundary |
---|
| 65 | mesh.import_ungenerate_file(project.buildings_filename, tag='wall') |
---|
| 66 | |
---|
| 67 | # Generate and write mesh to file |
---|
| 68 | mesh.generate_mesh(maximum_triangle_area=max_area, |
---|
| 69 | verbose=True) |
---|
| 70 | mesh.export_mesh_file(project.mesh_filename) |
---|
| 71 | |
---|
| 72 | |
---|
| 73 | #-------------------------------------------------------------------------- |
---|
| 74 | # Setup computational domain |
---|
| 75 | #-------------------------------------------------------------------------- |
---|
| 76 | |
---|
| 77 | domain = Domain(project.mesh_filename, use_cache = False, verbose = True) |
---|
| 78 | print domain.statistics() |
---|
| 79 | |
---|
| 80 | domain.set_name(project.basename) |
---|
| 81 | domain.set_datadir(project.outputdir) |
---|
| 82 | |
---|
| 83 | #------------------------------------------------------------------------------ |
---|
| 84 | # Setup initial conditions |
---|
| 85 | #------------------------------------------------------------------------------ |
---|
| 86 | |
---|
| 87 | domain.set_quantity('stage', project.initial_sealevel) |
---|
| 88 | domain.set_quantity('friction', 0.03) |
---|
| 89 | domain.set_quantity('elevation', |
---|
| 90 | filename=project.demname + '.pts', |
---|
| 91 | use_cache=True, |
---|
| 92 | verbose=True) |
---|
| 93 | |
---|
| 94 | #------------------------------------------------------------------------------ |
---|
| 95 | # Setup boundary conditions |
---|
| 96 | #------------------------------------------------------------------------------ |
---|
| 97 | |
---|
| 98 | D = Dirichlet_boundary([project.initial_sealevel, 0, 0]) |
---|
| 99 | R = Reflective_boundary(domain) |
---|
| 100 | W = Time_boundary(domain = domain, |
---|
| 101 | f=lambda t: [project.initial_sealevel + (60<t<480)*6, 0, 0]) |
---|
| 102 | |
---|
| 103 | domain.set_boundary({'exterior': D, |
---|
| 104 | 'side': D, |
---|
| 105 | 'wall': R, |
---|
| 106 | 'ocean': W, |
---|
| 107 | 'ghost': None}) |
---|
| 108 | |
---|
| 109 | |
---|
| 110 | |
---|
| 111 | |
---|
| 112 | #------------------------------------------------------------------------------ |
---|
| 113 | # Evolve system through time |
---|
| 114 | #------------------------------------------------------------------------------ |
---|
| 115 | |
---|
| 116 | t0 = time.time() |
---|
| 117 | for t in domain.evolve(yieldstep = 1, finaltime = 1200): |
---|
| 118 | domain.write_time() |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | print 'That took %.2f seconds' %(time.time()-t0) |
---|