[5411] | 1 | """Anuga convergence study using a true-scale version of the Okushiri |
---|
| 2 | Island tsunami wavetank experiment |
---|
[5343] | 3 | |
---|
[5411] | 4 | This script runs a modified version of the Okushiri Island benchmark |
---|
[5343] | 5 | |
---|
| 6 | as published at |
---|
| 7 | |
---|
| 8 | THE THIRD INTERNATIONAL WORKSHOP ON LONG-WAVE RUNUP MODELS |
---|
| 9 | June 17-18 2004 |
---|
| 10 | Wrigley Marine Science Center |
---|
| 11 | Catalina Island, California |
---|
| 12 | http://www.cee.cornell.edu/longwave/ |
---|
| 13 | |
---|
| 14 | This version up-scales the original 1:400 scale wave-tank experiment, to |
---|
[5411] | 15 | "true-scale" with mesh defined using interior polygons. |
---|
[5343] | 16 | |
---|
[5411] | 17 | The original validation data is available at |
---|
[5343] | 18 | http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2 |
---|
| 19 | where a detailed description of the problem is also available. |
---|
| 20 | |
---|
[5411] | 21 | Run create_okushiri_truescale.py to convert bathymetry and input boundary |
---|
| 22 | condition into NetCDF format. |
---|
[5343] | 23 | |
---|
| 24 | """ |
---|
[5411] | 25 | #---------------------------------------- |
---|
| 26 | # Import necessary modules |
---|
| 27 | #---------------------------------------- |
---|
[5343] | 28 | |
---|
[5411] | 29 | # Standard modules |
---|
| 30 | from os import sep,umask |
---|
| 31 | from os.path import dirname, basename |
---|
| 32 | from os import mkdir, access, F_OK |
---|
| 33 | from shutil import copy |
---|
| 34 | import time |
---|
| 35 | import sys |
---|
| 36 | |
---|
| 37 | # Related major packages |
---|
[5343] | 38 | from anuga.shallow_water import Domain |
---|
| 39 | from anuga.shallow_water import Reflective_boundary |
---|
| 40 | from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary |
---|
| 41 | from anuga.abstract_2d_finite_volumes.util import file_function |
---|
[5411] | 42 | from anuga.shallow_water.data_manager import copy_code_files, start_screen_catcher |
---|
| 43 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
[5343] | 44 | |
---|
| 45 | |
---|
[5411] | 46 | import project_truescale # Definition of filenames and interior polygons |
---|
[5343] | 47 | |
---|
[5411] | 48 | copy_code_files(project_truescale.output_run_time_dir,'run_okushiri_truescale.py', |
---|
| 49 | 'project_truescale.py' ) |
---|
| 50 | myid = 0 |
---|
| 51 | numprocs = 1 |
---|
| 52 | start_screen_catcher(project_truescale.output_run_time_dir, myid, numprocs) |
---|
| 53 | |
---|
| 54 | #-------------------------------------------------------------------------- |
---|
| 55 | # Create the triangular mesh based on overall bounding polygon with a |
---|
| 56 | # tagged boundary and interior regions defined in project_truescale.py |
---|
| 57 | # along with resolutions (maximal area of per triangle) for each polygon |
---|
| 58 | #-------------------------------------------------------------------------- |
---|
| 59 | |
---|
| 60 | create_mesh_from_regions(project_truescale.poly_all, |
---|
| 61 | boundary_tags={'wall': [0, 1, 3],'wave': [2]}, |
---|
| 62 | maximum_triangle_area=project_truescale.res_poly_all, |
---|
| 63 | #interior_regions=project_truescale.interior_regions, |
---|
| 64 | filename=project_truescale.mesh_name+'.msh', |
---|
| 65 | verbose=True) |
---|
| 66 | |
---|
| 67 | #------------------------------ |
---|
[5343] | 68 | # Create Domain from mesh |
---|
[5411] | 69 | #------------------------------ |
---|
| 70 | domain = Domain(project_truescale.mesh_name+'.msh', use_cache=True, verbose=True) |
---|
[5343] | 71 | print domain.statistics() |
---|
| 72 | |
---|
[5396] | 73 | z = domain.get_vertex_coordinates() |
---|
[5343] | 74 | |
---|
[5396] | 75 | # Write vertex coordinates to file |
---|
[5411] | 76 | filename=project.vertex_filename |
---|
[5396] | 77 | fid=open(filename,'w') |
---|
| 78 | fid.write('x (m), y (m)\n') |
---|
| 79 | for i in range(len(z)): |
---|
| 80 | pt=z[i] |
---|
| 81 | x=pt[0] |
---|
| 82 | y=pt[1] |
---|
| 83 | fid.write('%.6f, %.6f\n' %(x, y)) |
---|
| 84 | |
---|
[5343] | 85 | #------------------------- |
---|
| 86 | # Initial Conditions |
---|
| 87 | #------------------------- |
---|
| 88 | domain.set_quantity('friction', 0.0) |
---|
| 89 | domain.set_quantity('stage', 0.0) |
---|
| 90 | domain.set_quantity('elevation', |
---|
| 91 | filename=project_truescale.bathymetry_filename, |
---|
| 92 | alpha=0.02, |
---|
| 93 | verbose=True, |
---|
| 94 | use_cache=True) |
---|
| 95 | |
---|
| 96 | |
---|
| 97 | #------------------------- |
---|
| 98 | # Set simulation parameters |
---|
| 99 | #------------------------- |
---|
| 100 | domain.set_name(project_truescale.output_filename) # Name of output sww file |
---|
[5411] | 101 | domain.set_datadir(project_truescale.output_run_time_dir) # Name of output directory |
---|
[5343] | 102 | domain.set_default_order(2) # Apply second order scheme |
---|
| 103 | domain.set_all_limiters(0.9) # Max second order scheme (old lim) |
---|
[5411] | 104 | domain.set_minimum_storable_height(0.1) # Don't store h < 0.1m |
---|
[5343] | 105 | |
---|
| 106 | |
---|
| 107 | #------------------------- |
---|
| 108 | # Boundary Conditions |
---|
| 109 | #------------------------- |
---|
| 110 | |
---|
| 111 | # Create boundary function from timeseries provided in file |
---|
| 112 | function = file_function(project_truescale.boundary_filename, |
---|
| 113 | domain, verbose=True) |
---|
| 114 | |
---|
| 115 | # Create and assign boundary objects |
---|
| 116 | Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) |
---|
| 117 | Br = Reflective_boundary(domain) |
---|
| 118 | domain.set_boundary({'wave': Bts, 'wall': Br}) |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | # Select triangle containing ch5 for diagnostic output |
---|
| 122 | # around known gauge |
---|
| 123 | triangle_id = domain.get_triangle_containing_point([1808.4, 478.4]) |
---|
| 124 | # This should get triangle id 32833 with centroid (4.5244, 1.1972) |
---|
| 125 | |
---|
| 126 | |
---|
| 127 | #------------------------- |
---|
| 128 | # Evolve through time |
---|
| 129 | #------------------------- |
---|
| 130 | import time |
---|
| 131 | t0 = time.time() |
---|
| 132 | |
---|
[5411] | 133 | for t in domain.evolve(yieldstep=project_truescale.yieldstep, finaltime = 450): |
---|
[5343] | 134 | print domain.timestepping_statistics(track_speeds=False, |
---|
| 135 | triangle_id=triangle_id) |
---|
| 136 | print domain.boundary_statistics() |
---|
| 137 | |
---|
| 138 | print 'That took %.2f seconds' %(time.time()-t0) |
---|