[3125] | 1 | """Script for running a tsunami inundation scenario for Onslow, WA, 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.outputtimedir |
---|
| 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 | Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006 |
---|
| 11 | """ |
---|
| 12 | |
---|
| 13 | |
---|
[3178] | 14 | #---------------------------------------------------------------------------- |
---|
| 15 | # Import necessary modules |
---|
| 16 | #---------------------------------------------------------------------------- |
---|
[3125] | 17 | |
---|
| 18 | # Standard modules |
---|
| 19 | import time |
---|
| 20 | import sys |
---|
| 21 | from shutil import copy |
---|
[3226] | 22 | from os import mkdir, access, F_OK, path |
---|
[3125] | 23 | |
---|
| 24 | # Related major packages |
---|
| 25 | from pyvolution.shallow_water import Domain, Reflective_boundary, \ |
---|
| 26 | Dirichlet_boundary, Time_boundary, File_boundary |
---|
| 27 | from pyvolution.util import Screen_Catcher |
---|
| 28 | from pyvolution.region import Set_region |
---|
| 29 | |
---|
| 30 | # Application specific imports |
---|
| 31 | import project # Definition of file names and polygons |
---|
| 32 | import create_mesh |
---|
| 33 | |
---|
[3145] | 34 | |
---|
[3178] | 35 | #----------------------------------------------------------------------------- |
---|
[3145] | 36 | # Setup archiving of simulations |
---|
[3178] | 37 | #----------------------------------------------------------------------------- |
---|
[3145] | 38 | |
---|
[3125] | 39 | copy (project.codedirname, project.outputtimedir + 'project.py') |
---|
| 40 | copy (project.codedir + 'run_dam.py', project.outputtimedir + 'run_dam.py') |
---|
[3178] | 41 | copy (project.codedir + 'create_mesh.py', |
---|
| 42 | project.outputtimedir + 'create_mesh.py') |
---|
[3125] | 43 | print'output dir', project.outputtimedir |
---|
| 44 | |
---|
| 45 | #normal screen output is stored in |
---|
| 46 | screen_output_name = project.outputtimedir + "screen_output.txt" |
---|
| 47 | screen_error_name = project.outputtimedir + "screen_error.txt" |
---|
| 48 | |
---|
[3178] | 49 | #----------------------------------------------------------------------------- |
---|
[3125] | 50 | # Create the triangular mesh |
---|
[3178] | 51 | #----------------------------------------------------------------------------- |
---|
[3125] | 52 | |
---|
| 53 | create_mesh.generate() # this creates the mesh |
---|
| 54 | |
---|
[3226] | 55 | head,tail = path.split(project.mesh_filename) |
---|
| 56 | copy (project.mesh_filename, |
---|
| 57 | project.outputtimedir + tail ) |
---|
[3178] | 58 | #----------------------------------------------------------------------------- |
---|
[3125] | 59 | # Setup computational domain |
---|
[3178] | 60 | #----------------------------------------------------------------------------- |
---|
[3125] | 61 | domain = Domain(project.mesh_filename, use_cache = False, verbose = True) |
---|
| 62 | |
---|
| 63 | |
---|
| 64 | print 'Number of triangles = ', len(domain) |
---|
| 65 | print 'The extent is ', domain.get_extent() |
---|
| 66 | print domain.statistics() |
---|
| 67 | |
---|
| 68 | domain.set_name(project.basename) |
---|
| 69 | domain.set_datadir(project.outputtimedir) |
---|
| 70 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
| 71 | |
---|
[3178] | 72 | #----------------------------------------------------------------------------- |
---|
[3125] | 73 | # Setup initial conditions |
---|
[3178] | 74 | #----------------------------------------------------------------------------- |
---|
[3125] | 75 | |
---|
| 76 | tide = 0.0 |
---|
| 77 | |
---|
| 78 | def elevation_tilt(x, y): |
---|
| 79 | return -x*0.15 |
---|
| 80 | |
---|
| 81 | domain.set_quantity('stage', elevation_tilt) |
---|
| 82 | domain.set_quantity('friction', 0.03) |
---|
| 83 | domain.set_quantity('elevation',elevation_tilt) |
---|
| 84 | |
---|
| 85 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
| 86 | |
---|
| 87 | domain.set_region(Set_region('dam','stage',0.4,location = 'unique vertices')) |
---|
| 88 | |
---|
| 89 | Br = Reflective_boundary(domain) |
---|
| 90 | Bd = Dirichlet_boundary([tide,0,0]) |
---|
| 91 | |
---|
| 92 | #domain.set_boundary( {'wall': Br, 'wave': Bw} ) |
---|
| 93 | domain.set_boundary( {'wall': Br, 'wave': Br} ) |
---|
| 94 | |
---|
[3178] | 95 | #----------------------------------------------------------------------------- |
---|
[3125] | 96 | # Evolve system through time |
---|
[3178] | 97 | #----------------------------------------------------------------------------- |
---|
[3125] | 98 | import time |
---|
| 99 | t0 = time.time() |
---|
| 100 | |
---|
[3145] | 101 | for t in domain.evolve(yieldstep = 0.1, finaltime = 5): |
---|
[3125] | 102 | domain.write_time() |
---|
| 103 | |
---|
| 104 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
| 105 | |
---|
| 106 | print 'finished' |
---|