Changeset 4153 for anuga_work/production/broome_2006
- Timestamp:
- Jan 9, 2007, 3:55:58 PM (18 years ago)
- Location:
- anuga_work/production/broome_2006
- Files:
-
- 1 added
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/production/broome_2006/project.py
r3992 r4153 7 7 import sys 8 8 from time import localtime, strftime, gmtime 9 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon9 from anuga.utilities.polygon import read_polygon, plot_polygons, is_inside_polygon, number_mesh_triangles 10 10 #from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm 11 from anuga.utilities.system_tools import get_user_name 11 12 12 if sys.platform == 'win32': 13 home = getenv('INUNDATIONHOME') 14 user = getenv('USERPROFILE') 13 # file and system info 14 #--------------------------------- 15 codename = 'project.py' 15 16 16 else: 17 home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation') 18 user = getenv('LOGNAME') 19 print 'USER:', user 17 home = getenv('INUNDATIONHOME') #Sandpit's parent dir 18 user = get_user_name() 20 19 21 20 # INUNDATIONHOME is the inundation directory, not the data directory. 22 21 home += sep +'data' 23 22 23 #time stuff 24 time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir 25 gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir 26 build_time = time+'_build' 27 run_time = time+'_run' 28 print 'gtime: ', gtime 29 30 tide = 0 31 24 32 #Making assumptions about the location of scenario data 25 33 state = 'western_australia' 26 scenario_dir_name = 'broome_tsunami_scenario_2006' 34 scenario_name = 'broome' 35 scenario = 'broome_tsunami_scenario_2006' 27 36 28 37 # onshore data provided by WA DLI 29 38 onshore_name = 'dted2_z51' # original 30 39 31 # AHO + DPI data 40 #island 41 #island_name = 'rott_dli_ext' # original 42 43 # offshore 44 coast_name = 'coastline' 32 45 offshore_name1 = 'XY100011610' 33 46 offshore_name2 = 'XY100011611' … … 53 66 offshore_name22 = 'XYWADPI' 54 67 55 # developed by NM&I 56 coast_name = 'coastline' 68 #final topo name 69 combined_name ='perth_combined_elevation' 70 combined_smaller_name = 'perth_combined_elevation_smaller' 57 71 58 # TIN model to fill in data gap59 bathy_interp = 'nearest_neighbour' #'interpolate'60 72 61 boundary_basename = 'SU-AU' # Mw ? 73 topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep 74 topographies_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep 75 topographies_time_dir = topographies_dir+build_time+sep 62 76 63 #swollen/ all data output 64 basename = 'source' 65 codename = 'project.py' 77 #input topo file location 78 onshore_in_dir_name = topographies_in_dir + onshore_name 79 island_in_dir_name = topographies_in_dir + island_name 80 island_in_dir_name1 = topographies_in_dir + island_name1 81 island_in_dir_name2 = topographies_in_dir + island_name2 82 island_in_dir_name3 = topographies_in_dir + island_name3 66 83 67 #Derive subdirectories and filenames 68 local_time = strftime('%Y%m%d_%H%M%S',gmtime()) 69 meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep 70 datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep 71 gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep 72 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep 73 boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep 74 outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep 75 outputtimedir = outputdir + local_time + sep 76 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep 84 coast_in_dir_name = topographies_in_dir + coast_name 85 offshore_in_dir_name = topographies_in_dir + offshore_name 86 offshore1_in_dir_name = topographies_in_dir + offshore1_name 77 87 78 gauge_filename = gaugedir + 'broome_gauges.csv' 79 community_filename = gaugedir + 'CHINS_v2.csv' 80 community_broome = gaugedir + 'community_broome.csv' 81 codedir = getcwd()+sep 82 codedirname = codedir + 'project.py' 83 meshname = outputtimedir + 'mesh_' + basename 88 onshore_dir_name = topographies_dir + onshore_name 89 island_dir_name = topographies_dir + island_name 90 island_dir_name1 = topographies_dir + island_name1 91 island_dir_name2 = topographies_dir + island_name2 92 island_dir_name3 = topographies_dir + island_name3 84 93 85 # Necessary if using point datasets, rather than grid 86 onshore_dem_name = datadir + onshore_name 87 offshore_interp_dem_name = datadir + bathy_interp 88 offshore_dem_name1 = datadir + offshore_name1 89 offshore_dem_name2 = datadir + offshore_name2 90 offshore_dem_name3 = datadir + offshore_name3 91 offshore_dem_name4 = datadir + offshore_name4 92 offshore_dem_name5 = datadir + offshore_name5 93 offshore_dem_name6 = datadir + offshore_name6 94 offshore_dem_name7 = datadir + offshore_name7 95 offshore_dem_name8 = datadir + offshore_name8 96 offshore_dem_name9 = datadir + offshore_name9 97 offshore_dem_name10 = datadir + offshore_name10 98 offshore_dem_name11 = datadir + offshore_name11 99 offshore_dem_name12 = datadir + offshore_name12 100 offshore_dem_name13 = datadir + offshore_name13 101 offshore_dem_name14 = datadir + offshore_name14 102 offshore_dem_name15 = datadir + offshore_name15 103 offshore_dem_name16 = datadir + offshore_name16 104 offshore_dem_name17 = datadir + offshore_name17 105 offshore_dem_name18 = datadir + offshore_name18 106 offshore_dem_name19 = datadir + offshore_name19 107 offshore_dem_name20 = datadir + offshore_name20 108 offshore_dem_name21 = datadir + offshore_name21 109 offshore_dem_name22 = datadir + offshore_name22 94 coast_dir_name = topographies_dir + coast_name 95 offshore_dir_name = topographies_dir + offshore_name 110 96 111 coast_dem_name = datadir + coast_name 97 #final topo files 98 combined_dir_name = topographies_dir + combined_name 99 combined_time_dir_name = topographies_time_dir + combined_name 100 combined_smaller_name_dir = topographies_dir + combined_smaller_name 101 #combined_time_dir_final_name = topographies_time_dir + combined_final_name 112 102 113 combined_dem_name = datadir + 'broome_combined_elevation' 103 meshes_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'meshes'+sep 104 meshes_dir_name = meshes_dir + scenario_name 105 106 polygons_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'polygons'+sep 107 tide_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'tide_data'+sep 108 109 110 boundaries_source = '????' 111 #boundaries locations 112 boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep 113 boundaries_in_dir_name = boundaries_in_dir + scenario_name 114 boundaries_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep 115 boundaries_dir_name = boundaries_dir + scenario_name 116 #boundaries_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+build_time+sep 117 #boundaries_time_dir_name = boundaries_time_dir + boundaries_name #Used by post processing 118 119 #output locations 120 output_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep 121 output_build_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+build_time+sep 122 output_run_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+run_time+sep 123 output_run_time_dir_name = output_run_time_dir + scenario_name #Used by post processing 124 125 #gauges 126 gauge_name = 'broome_gauges.csv' 127 gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep 128 gauges_dir_name = gauges_dir + gauge_name 129 130 community_filename = gauges_dir + 'CHINS_v2.csv' 131 community_broome = gauges_dir + 'community_broome.csv' 132 133 114 134 115 135 ############################### 116 136 # Domain definitions 117 137 ############################### 138 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon 118 139 119 # bounding box for clipping MOST/URS output (much bigger than study area) 120 ##south = degminsec2decimal_degrees(-19,0,0) 121 ##north = degminsec2decimal_degrees(-17,15,0) 122 ##west = degminsec2decimal_degrees(120,0,0) 123 ##east = degminsec2decimal_degrees(122,30,0) 124 ## 125 ##d0 = [south, west] 126 ##d1 = [south, east] 127 ##d2 = [north, east] 128 ##d3 = [north, west] 129 ##poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3]) 130 ##refzone = zone 140 poly_all = read_polygon(polygons_dir+'extent.csv') 141 res_poly_all = 100000 131 142 132 # bounding polygon for study area 133 polyAll = read_polygon(polygondir+'extent.csv') 143 ############################### 144 # Interior region definitions 145 ############################### 134 146 135 # plot bounding polygon and make sure BC info surrounds it 136 #plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False) 137 print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0 147 poly_1 = read_polygon(polygons_dir+'Broome_Local_Polygon_update.csv') 148 res_1 = 20000 149 150 poly_2 = read_polygon(polygons_dir+'Broome_Close2_update.csv') 151 res_2 = 1000 152 153 poly_3 = read_polygon(polygons_dir+'Broome_Coast_update2.csv') 154 res_3 = 1000 155 #assert zone == refzone 156 157 interior_regions = [[poly_1,res_1],[poly_2,res_2] 158 ,[poly_3,res_3]] 159 160 trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all) 161 162 print 'min number triangles', trigs_min 138 163 139 164 ################################################################### … … 142 167 143 168 # exporting asc grid 144 eastingmin = 412036.43 145 eastingmax = 427184.37 146 northingmin = 8007984 147 northingmax = 8029830 148 149 150 169 eastingmin = 406215.87 170 eastingmax = 440208.78 171 northingmin = 7983427.73 172 northingmax = 8032834.52 151 173 152 174 153 ###############################154 # Interior region definitions155 ###############################156 175 157 # broome digitized polygons158 poly_broome1 = read_polygon(polygondir+'Broome_Local_Polygon_update.csv')159 poly_broome2 = read_polygon(polygondir+'Broome_Close2_update.csv')160 poly_broome3 = read_polygon(polygondir+'Broome_Coast_update2.csv')161 #poly_broome4 = read_polygon(polygondir+'Cable_Beach_revised.csv')162 163 #plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],'boundingpoly2',verbose=False)164 print 'Area of local polygon', polygon_area(poly_broome1)/1000000.0165 print 'Area of close polygon', polygon_area(poly_broome2)/1000000.0166 print 'Area of coastal polygon', polygon_area(poly_broome3)/1000000.0167 #print 'Area of cable beach polygon', polygon_area(poly_broome4)/1000000.0168 169 for i in poly_broome3:170 v = is_inside_polygon(i,poly_broome1, verbose=False)171 if v == False: print v172 173 def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):174 from anuga.utilities.polygon import polygon_area175 176 # TO DO check if any of the regions fall inside one another177 no_triangles = 0.0178 area = polygon_area(bounding_poly)179 for i,j in interior_regions:180 this_area = polygon_area(i)181 no_triangles += this_area/j182 area -= this_area183 print j, this_area/1000000., area/1000000.184 no_triangles += area/remainder_res185 return int(no_triangles/0.7) -
anuga_work/production/broome_2006/run_broome.py
r4063 r4153 1 """Script for running a tsunami inundation scenario for Broome, WA, Australia.1 """Script for running tsunami inundation scenario for Dampier, WA, Australia. 2 2 3 3 Source data such as elevation and boundary data is assumed to be available in 4 4 directories specified by project.py 5 The output sww file is stored in project.output timedir5 The output sww file is stored in project.output_run_time_dir 6 6 7 7 The scenario is defined by a triangular mesh created from project.polygon, 8 the elevation data and a tsunami wave generated by MOST.9 10 Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 20068 the elevation data and a simulated submarine landslide. 9 10 Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006 11 11 """ 12 12 13 #------------------------------------------------------------------------------ -13 #------------------------------------------------------------------------------ 14 14 # Import necessary modules 15 #------------------------------------------------------------------------------ -15 #------------------------------------------------------------------------------ 16 16 17 17 # Standard modules 18 import os 18 from os import sep, mkdir, access, F_OK 19 from os.path import dirname, basename 20 from shutil import copy 19 21 import time 20 from shutil import copy21 from os.path import dirname, basename22 from os import mkdir, access, F_OK, sep23 22 import sys 24 23 25 24 # Related major packages 26 from anuga.shallow_water import Domain , Reflective_boundary, \27 Dirichlet_boundary, Time_boundary, File_boundary28 from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts29 from anuga. abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files25 from anuga.shallow_water import Domain 26 from anuga.shallow_water import Dirichlet_boundary, File_boundary, Reflective_boundary 27 from Numeric import allclose 28 from anuga.pmesh.mesh_interface import create_mesh_from_regions 30 29 from anuga.geospatial_data.geospatial_data import * 31 30 from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files 31 from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier 32 from anuga.utilities.polygon import plot_polygons, polygon_area 32 33 33 34 # Application specific imports 34 35 import project # Definition of file names and polygons 35 36 36 #------------------------------------------------------------------------------ -37 #------------------------------------------------------------------------------ 37 38 # Copy scripts to time stamped output directory and capture screen 38 39 # output to file 39 #------------------------------------------------------------------------------- 40 41 # creates copy of code in output dir 42 copy_code_files(project.outputtimedir,__file__,dirname(project.__file__)+sep+ project.__name__+'.py' ) 43 myid = 0 44 numprocs = 1 45 #start_screen_catcher(project.outputtimedir, myid, numprocs) 46 47 print 'USER: ', project.user 48 49 #------------------------------------------------------------------------------- 50 # Preparation of topographic data 51 # 52 # Convert ASC 2 DEM 2 PTS using source data and store result in source data 53 #------------------------------------------------------------------------------- 40 #------------------------------------------------------------------------------ 41 42 start_screen_catcher(project.output_run_time_dir, myid, numprocs) 54 43 55 44 # filenames 56 onshore_dem_name = project.onshore_dem_name 57 offshore_interp_dem_name = project.offshore_interp_dem_name 58 coast_points = project.coast_dem_name 59 meshname = project.meshname+'.msh' 60 61 # creates DEM from asc data 62 convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True) 63 64 #creates pts file for onshore DEM 65 dem2pts(onshore_dem_name, use_cache=True, verbose=True) 66 67 # creates DEM from asc data 68 convert_dem_from_ascii2netcdf(offshore_interp_dem_name, use_cache=True, verbose=True) 69 70 #creates pts file for offshore interpolated DEM 71 dem2pts(offshore_interp_dem_name, use_cache=True, verbose=True) 72 73 print 'create offshore' 74 G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')+\ 75 Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')+\ 76 Geospatial_data(file_name = project.offshore_dem_name3 + '.xya')+\ 77 Geospatial_data(file_name = project.offshore_dem_name4 + '.xya')+\ 78 Geospatial_data(file_name = project.offshore_dem_name5 + '.xya')+\ 79 Geospatial_data(file_name = project.offshore_dem_name6 + '.xya')+\ 80 Geospatial_data(file_name = project.offshore_dem_name7 + '.xya')+\ 81 Geospatial_data(file_name = project.offshore_dem_name8 + '.xya')+\ 82 Geospatial_data(file_name = project.offshore_dem_name9 + '.xya')+\ 83 Geospatial_data(file_name = project.offshore_dem_name10 + '.xya')+\ 84 Geospatial_data(file_name = project.offshore_dem_name11 + '.xya')+\ 85 Geospatial_data(file_name = project.offshore_dem_name12 + '.xya')+\ 86 Geospatial_data(file_name = project.offshore_dem_name13 + '.xya')+\ 87 Geospatial_data(file_name = project.offshore_dem_name14 + '.xya')+\ 88 Geospatial_data(file_name = project.offshore_dem_name15 + '.xya')+\ 89 Geospatial_data(file_name = project.offshore_dem_name16 + '.xya')+\ 90 Geospatial_data(file_name = project.offshore_dem_name17 + '.xya')+\ 91 Geospatial_data(file_name = project.offshore_dem_name18 + '.xya')+\ 92 Geospatial_data(file_name = project.offshore_dem_name19 + '.xya')+\ 93 Geospatial_data(file_name = project.offshore_dem_name20 + '.xya')+\ 94 Geospatial_data(file_name = project.offshore_dem_name21 + '.xya')+\ 95 Geospatial_data(file_name = project.offshore_dem_name22 + '.xya')+\ 96 Geospatial_data(file_name = project.offshore_interp_dem_name + '.pts') 97 print 'create onshore' 98 G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts') 99 print 'create coast' 100 G3 = Geospatial_data(file_name = project.coast_dem_name + '.xya') 101 print 'add' 102 G = G1 + G2 + G3 103 print 'export points' 104 G.export_points_file(project.combined_dem_name + '.pts') 105 G.export_points_file(project.combined_dem_name + '.xya') 106 107 #---------------------------------------------------------------------------- 108 # Create the triangular mesh based on overall clipping polygon with a tagged 45 boundaries_name = project.scenario 46 meshes_dir_name = project.meshes_dir_name+'.msh' 47 boundaries_dir_name = project.boundaries_dir_name 48 49 tide = project.tide 50 51 # creates copy of code in output dir 52 if myid == 0: 53 copy_code_files(project.output_run_time_dir,__file__, 54 dirname(project.__file__)+sep+ project.__name__+'.py' ) 55 barrier() 56 57 print 'USER: ', project.user 58 print 'min triangles', project.trigs_min, 59 print 'Note: This is generally about 20% less than the final amount' 60 61 #-------------------------------------------------------------------------- 62 # Create the triangular mesh based on overall clipping polygon with a 63 # tagged 109 64 # boundary and interior regions defined in project.py along with 110 65 # resolutions (maximal area of per triangle) for each polygon 111 #------------------------------------------------------------------------------- 112 113 from anuga.pmesh.mesh_interface import create_mesh_from_regions 114 remainder_res = 500000 115 local_res = 25000 116 broome_res = 5000 117 coast_res = 500 118 interior_regions = [[project.poly_broome1, local_res], 119 [project.poly_broome2, broome_res], 120 [project.poly_broome3, coast_res]] 121 122 from project import number_mesh_triangles 123 print 'estimate of number of triangles', \ 124 number_mesh_triangles(interior_regions, project.polyAll, remainder_res) 125 66 #-------------------------------------------------------------------------- 67 68 if myid == 0: 69 70 print 'start create mesh from regions' 71 create_mesh_from_regions(project.poly_all, 72 boundary_tags={'back': [1, 2], 'side': [0,3], 73 'ocean': [4, 5, 6]}, 74 maximum_triangle_area=project.res_poly_all, 75 interior_regions=project.interior_regions, 76 filename=meshes_dir_name, 77 use_cache=True, 78 verbose=True) 79 80 # to sync all processors are ready 81 barrier() 82 83 #------------------------------------------------------------------------- 84 # Setup computational domain 85 #------------------------------------------------------------------------- 86 print 'Setup computational domain' 87 #domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True) 88 domain = Domain(meshes_dir_name, use_cache=True, verbose=True) 89 print domain.statistics() 90 ''' 91 print 'starting to create boundary conditions' 92 boundaries_in_dir_name = project.boundaries_in_dir_name 93 94 from anuga.shallow_water.data_manager import urs2sww, ferret2sww 95 96 # put above distribute 97 print 'boundary file is: ',boundaries_dir_name 126 98 from caching import cache 127 _ = cache(create_mesh_from_regions, 128 project.polyAll, 129 {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2], 130 'e3': [3], 'e4':[4], 'e5': [5], 131 'e6': [6]}, 132 'maximum_triangle_area': remainder_res, 133 'filename': meshname, 134 'interior_regions': interior_regions}, 135 verbose = True, evaluate=False) 136 print 'created mesh' 137 138 #------------------------------------------------------------------------------- 139 # Setup computational domain 140 #------------------------------------------------------------------------------- 141 domain = Domain(meshname, use_cache = True, verbose = True) 142 143 print 'Number of triangles = ', len(domain) 144 print 'The extent is ', domain.get_extent() 145 print domain.statistics() 146 147 domain.set_name(project.basename) 148 domain.set_datadir(project.outputtimedir) 149 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) 150 domain.set_minimum_storable_height(0.01) 151 152 #------------------------------------------------------------------------------- 99 if myid == 0: 100 cache(ferret2sww, 101 (boundaries_in_dir_name, 102 # boundaries_time_dir_name), 103 boundaries_dir_name), 104 {'verbose': True, 105 'minlat': project.south_boundary, 106 'maxlat': project.north_boundary, 107 'minlon': project.west_boundary, 108 'maxlon': project.east_boundary, 109 # 'minlat': project.south, 110 # 'maxlat': project.north, 111 # 'minlon': project.west, 112 # 'maxlon': project.east, 113 'mint': 0, 'maxt': 35100, 114 'origin': domain.geo_reference.get_origin(), 115 'mean_stage': project.tide, 116 # 'zscale': 1, #Enhance tsunami 117 'fail_on_NaN': False}, 118 verbose = True, 119 ) 120 barrier() 121 ''' 122 123 #------------------------------------------------------------------------- 153 124 # Setup initial conditions 154 #------------------------------------------------------------------------------- 155 156 tide = 0.0 157 domain.set_quantity('stage', tide) 158 domain.set_quantity('friction', 0.0) 159 domain.set_quantity('elevation', 160 filename = project.combined_dem_name + '.pts', 125 #------------------------------------------------------------------------- 126 if myid == 0: 127 128 print 'Setup initial conditions' 129 130 domain.set_quantity('stage', tide) 131 domain.set_quantity('friction', 0.01) 132 #combined_time_dir_name = project.topographies_dir+build_time+project.combined_name 133 print 'Start Set quantity' 134 135 domain.set_quantity('elevation', 136 filename = project.combined_dir_name+'.txt', 137 # filename = project.combined_smaller_name_dir+'.xya', 161 138 use_cache = True, 162 139 verbose = True, 163 alpha = 0.1 164 ) 165 166 #------------------------------------------------------------------------------- 167 # Setup boundary conditions 168 #------------------------------------------------------------------------------- 169 ''' 170 print 'start urs2sww' 171 print '', project.boundary_basename 172 from anuga.shallow_water.data_manager import urs2sww 173 174 south = project.south 175 north = project.north 176 west = project.west 177 east = project.east 178 179 #note only need to do when an SWW file for the MOST boundary doesn't exist 180 cache(urs2sww, 181 (source_dir + project.boundary_basename, 182 source_dir + project.boundary_basename), 183 {'verbose': True, 184 'minlat': south, 185 'maxlat': north, 186 'minlon': west, 187 'maxlon': east, 188 #'origin': domain.geo_reference.get_origin(), 189 'mean_stage': tide, 190 'zscale': 1, #Enhance tsunami 191 'fail_on_NaN': False, 192 'inverted_bathymetry': True}, 193 #evaluate = True, 194 verbose = True, 195 dependencies = source_dir + project.boundary_basename + '.sww') 196 197 ''' 140 alpha = 0.1) 141 print 'Finished Set quantity' 142 barrier() 143 144 #------------------------------------------------------ 145 # Distribute domain to implement parallelism !!! 146 #------------------------------------------------------ 147 148 if numprocs > 1: 149 domain=distribute(domain) 150 151 #------------------------------------------------------ 152 # Set domain parameters 153 #------------------------------------------------------ 154 155 domain.set_name(project.scenario_name) 156 domain.set_datadir(project.output_run_time_dir) 157 domain.set_default_order(2) # Apply second order scheme 158 domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm 159 domain.set_store_vertices_uniquely(False) 160 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) 161 domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK) 162 163 #------------------------------------------------------------------------- 164 # Setup boundary conditions 165 #------------------------------------------------------------------------- 198 166 print 'Available boundary tags', domain.get_boundary_tags() 199 167 200 #Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 201 # domain, verbose = True) 168 print 'Reading Boundary file' 169 print 'domain id', id(domain) 170 #Bf = File_boundary(boundaries_dir_name + '.sww', 171 # domain, time_thinning=5, use_cache=True, verbose=True) 172 173 print 'finished reading boundary file' 202 174 Br = Reflective_boundary(domain) 203 175 Bd = Dirichlet_boundary([tide,0,0]) 204 205 # 7 min square wave starting at 1 min, 6m high 206 Bw = Time_boundary(domain = domain, 207 f=lambda t: [(60<t<480)*10, 0, 0]) 208 209 domain.set_boundary( {'e0': Bd, 'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd, 210 'e5': Bd, 'e6': Bd} ) 211 212 213 #------------------------------------------------------------------------------- 176 Bo = Dirichlet_boundary([tide+5.0,0,0]) 177 178 print'set_boundary' 179 domain.set_boundary({'back': Br, 180 'side': Bd, 181 'ocean': Bd}) 182 print'finish set boundary' 183 184 #---------------------------------------------------------------------------- 214 185 # Evolve system through time 215 #---------------------------------------------------------------------------- ---216 import time 186 #---------------------------------------------------------------------------- 187 217 188 t0 = time.time() 218 189 219 for t in domain.evolve(yieldstep = 240, finaltime = 480):190 for t in domain.evolve(yieldstep = 60, finaltime = 10000): 220 191 domain.write_time() 221 domain.write_boundary_statistics(tags = 'e14') 192 domain.write_boundary_statistics(tags = 'ocean') 193 if allclose(t, 120): 194 domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo}) 195 196 if allclose(t, 1020): 197 domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd}) 198 222 199 223 200 print 'That took %.2f seconds' %(time.time()-t0) 224 201 225 print 'finished'
Note: See TracChangeset
for help on using the changeset viewer.