Changeset 6382 for anuga_work/production/australia_ph2/broome/run_model.py
- Timestamp:
- Feb 23, 2009, 4:16:51 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/production/australia_ph2/broome/run_model.py
r6340 r6382 22 22 # Standard modules 23 23 import os 24 import os.path 24 25 import time 26 from time import localtime, strftime, gmtime 25 27 26 28 # Related major packages 29 from Scientific.IO.NetCDF import NetCDFFile 30 import Numeric as num 31 27 32 from anuga.interface import create_domain_from_regions 28 33 from anuga.interface import Transmissive_stage_zero_momentum_boundary … … 36 41 from anuga.shallow_water.data_manager import start_screen_catcher 37 42 from anuga.shallow_water.data_manager import copy_code_files 43 from anuga.shallow_water.data_manager import urs2sts 38 44 from anuga.utilities.polygon import read_polygon, Polygon_function 39 45 40 46 # Application specific imports 41 import project # Definition of file names and polygons 47 from setup_model import project 48 import build_urs_boundary as bub 42 49 43 44 #------------------------------------------------------------------------------ 50 #------------------------------------------------------------------------------- 45 51 # Copy scripts to time stamped output directory and capture screen 46 52 # output to file. Copy script must be before screen_catcher 47 #------------------------------------------------------------------------------ 53 #------------------------------------------------------------------------------- 54 48 55 copy_code_files(project.output_run, __file__, 49 os.path. dirname(project.__file__)+os.sep+\50 project.__name__+'.py')56 os.path.join(os.path.dirname(project.__file__), 57 project.__name__+'.py')) 51 58 start_screen_catcher(project.output_run, 0, 1) 52 59 53 54 #------------------------------------------------------------------------------ 60 #------------------------------------------------------------------------------- 55 61 # Create the computational domain based on overall clipping polygon with 56 62 # a tagged boundary and interior regions defined in project.py along with 57 63 # resolutions (maximal area of per triangle) for each polygon 58 #------------------------------------------------------------------------------ 64 #------------------------------------------------------------------------------- 65 59 66 print 'Create computational domain' 67 68 # Create the STS file 69 print 'project.mux_data_folder=%s' % project.mux_data_folder 70 if not os.path.exists(project.event_sts + '.sts'): 71 bub.build_urs_boundary(project.mux_input_filename, project.event_sts) 60 72 61 73 # Read in boundary from ordered sts file … … 70 82 71 83 # Number of boundary segments 72 N = len(event_sts)-1 84 num_ocean_segments = len(event_sts) - 1 85 # Number of landward_boundary points 86 num_land_points = file_length(project.landward_boundary) 73 87 74 # Number of landward_boundary points75 M = file_length(project.landward_boundary)76 77 88 # Boundary tags refer to project.landward_boundary 78 89 # 4 points equals 5 segments start at N 79 boundary_tags={'back': range(N+1, N+M), 80 'side': [N, N+M], 81 'ocean': range(N)} 90 boundary_tags={'back': range(num_ocean_segments+1, 91 num_ocean_segments+num_land_points), 92 'side': [num_ocean_segments, 93 num_ocean_segments+num_land_points], 94 'ocean': range(num_ocean_segments)} 82 95 83 96 # Build mesh and domain … … 95 108 domain.set_minimum_storable_height(0.01) # Don't store depth less than 1cm 96 109 110 #------------------------------------------------------------------------------- 111 # Setup initial conditions 112 #------------------------------------------------------------------------------- 97 113 98 #------------------------------------------------------------------------------99 # Setup initial conditions100 #------------------------------------------------------------------------------101 114 print 'Setup initial conditions' 102 115 103 116 # Set the initial stage in the offcoast region only 104 ##IC = Polygon_function(project.land_initial_conditions, 105 ## default=project.tide, 106 ## geo_reference=domain.geo_reference) 107 domain.set_quantity('stage', 0, use_cache=True, verbose=True) 117 if project.land_initial_conditions: 118 IC = Polygon_function(project.land_initial_conditions, 119 default=project.tide, 120 geo_reference=domain.geo_reference) 121 else: 122 IC = 0 123 domain.set_quantity('stage', IC, use_cache=True, verbose=True) 108 124 domain.set_quantity('friction', project.friction) 109 125 domain.set_quantity('elevation', … … 113 129 alpha=project.alpha) 114 130 131 #------------------------------------------------------------------------------- 132 # Setup boundary conditions 133 #------------------------------------------------------------------------------- 115 134 116 #------------------------------------------------------------------------------117 # Setup boundary conditions118 #------------------------------------------------------------------------------119 135 print 'Set boundary - available tags:', domain.get_boundary_tags() 120 136 121 137 Br = Reflective_boundary(domain) 122 138 Bt = Transmissive_stage_zero_momentum_boundary(domain) 139 Bd = Dirichlet_boundary([project.tide, 0, 0]) 123 140 Bf = Field_boundary(project.event_sts+'.sts', 124 141 domain, mean_stage=project.tide, … … 130 147 131 148 domain.set_boundary({'back': Br, 132 'side': B t,149 'side': Bd, 133 150 'ocean': Bf}) 134 151 152 #------------------------------------------------------------------------------- 153 # Evolve system through time 154 #------------------------------------------------------------------------------- 135 155 136 #------------------------------------------------------------------------------137 # Evolve system through time138 #------------------------------------------------------------------------------139 156 t0 = time.time() 140 157 for t in domain.evolve(yieldstep=project.yieldstep, … … 144 161 print domain.boundary_statistics(tags='ocean') 145 162 146 print 'Simulation took %.2f seconds' % (time.time()-t0)163 print 'Simulation took %.2f seconds' % (time.time()-t0)
Note: See TracChangeset
for help on using the changeset viewer.