Changeset 6178
- Timestamp:
- Jan 15, 2009, 5:45:46 PM (16 years ago)
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
r6175 r6178 348 348 349 349 if time_limit is not None: 350 #if verbose is True: 351 # print '****** Time limit', time_limit 352 # print '****** Start time', starttime 353 # print '****** Time in ', time[0], time[-1] 354 350 355 # Adjust given time limit to given start time 351 356 time_limit = time_limit - starttime 357 352 358 353 359 # Find limit point -
anuga_core/source/anuga/pmesh/mesh_interface.py
r6177 r6178 90 90 'minimum_triangle_angle': minimum_triangle_angle, 91 91 'fail_if_polygons_outside': fail_if_polygons_outside, 92 'verbose': verbose} # FIXME (Ole): Should be bypassed one day 93 # What should be bypassed? Verbose? 94 95 #print 'kwargs', kwargs 92 'verbose': verbose} # FIXME (Ole): Should be bypassed one day. See ticket:14 96 93 97 94 # Call underlying engine with or without caching -
anuga_core/source/anuga/shallow_water/__init__.py
r6045 r6178 8 8 # Make selected classes available directly 9 9 from shallow_water_domain import Domain,\ 10 Transmissive_boundary, Reflective_boundary,\ 11 Dirichlet_boundary, Time_boundary, File_boundary,\ 12 Transmissive_momentum_set_stage_boundary,\ 13 Dirichlet_discharge_boundary,\ 14 Field_boundary 10 create_domain_from_regions,\ 11 Transmissive_boundary, Reflective_boundary,\ 12 Dirichlet_boundary, Time_boundary, File_boundary,\ 13 Transmissive_momentum_set_stage_boundary,\ 14 Dirichlet_discharge_boundary,\ 15 Field_boundary 15 16 16 17 -
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r6173 r6178 94 94 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ 95 95 import Transmissive_boundary 96 96 from anuga.pmesh.mesh_interface import create_mesh_from_regions 97 97 from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric 98 98 from anuga.geospatial_data.geospatial_data import ensure_geospatial … … 117 117 118 118 119 #--------------------- 119 120 #--------------------------- 121 # Create domain from regions 122 #--------------------------- 123 124 def create_domain_from_regions(bounding_polygon, 125 boundary_tags, 126 maximum_triangle_area=None, 127 mesh_filename=None, 128 interior_regions=None, 129 interior_holes=None, 130 poly_geo_reference=None, 131 mesh_geo_reference=None, 132 minimum_triangle_angle=28.0, 133 fail_if_polygons_outside=True, 134 use_cache=False, 135 verbose=True): 136 """Create domain from bounding polygons and resolutions. 137 138 bounding_polygon is a list of points in Eastings and Northings, 139 relative to the poly_geo_reference. 140 141 Boundary tags is a dictionary of symbolic tags. For every tag there 142 is a list of indices referring to segments associated with that tag. 143 If a segment is omitted it will be assigned the default tag ''. 144 145 maximum_triangle_area is the maximal area per triangle 146 for the bounding polygon, excluding the interior regions. 147 148 Interior_regions is a list of tuples consisting of (polygon, 149 resolution) for each region to be separately refined. Do not have 150 polygon lines cross or be on-top of each other. Also do not have 151 polygon close to each other. 152 153 NOTE: If a interior_region is outside the bounding_polygon it should 154 throw an error 155 156 Interior_holes is a list of ploygons for each hole. 157 158 This function does not allow segments to share points - use underlying 159 pmesh functionality for that 160 161 poly_geo_reference is the geo_reference of the bounding polygon and 162 the interior polygons. 163 If none, assume absolute. Please pass one though, since absolute 164 references have a zone. 165 166 mesh_geo_reference is the geo_reference of the mesh to be created. 167 If none is given one will be automatically generated. It was use 168 the lower left hand corner of bounding_polygon (absolute) 169 as the x and y values for the geo_ref. 170 171 Returns the shallow water domain instance 172 173 Note, interior regions should be fully nested, as overlaps may cause 174 unintended resolutions. 175 176 fail_if_polygons_outside: If True (the default) Exception in thrown 177 where interior polygons fall outside bounding polygon. If False, these 178 will be ignored and execution continued. 179 180 181 """ 182 183 184 # Build arguments and keyword arguments for use with caching or apply. 185 args = (bounding_polygon, 186 boundary_tags) 187 188 kwargs = {'maximum_triangle_area': maximum_triangle_area, 189 'mesh_filename': mesh_filename, 190 'interior_regions': interior_regions, 191 'interior_holes': interior_holes, 192 'poly_geo_reference': poly_geo_reference, 193 'mesh_geo_reference': mesh_geo_reference, 194 'minimum_triangle_angle': minimum_triangle_angle, 195 'fail_if_polygons_outside': fail_if_polygons_outside, 196 'verbose': verbose} #FIXME (Ole): See ticket:14 197 198 # Call underlying engine with or without caching 199 if use_cache is True: 200 try: 201 from anuga.caching import cache 202 except: 203 msg = 'Caching was requested, but caching module'+\ 204 'could not be imported' 205 raise msg 206 207 208 domain = cache(_create_domain_from_regions, 209 args, kwargs, 210 verbose=verbose, 211 compression=False) 212 else: 213 domain = apply(_create_domain_from_regions, 214 args, kwargs) 215 216 return domain 217 218 219 def _create_domain_from_regions(bounding_polygon, 220 boundary_tags, 221 maximum_triangle_area=None, 222 mesh_filename=None, 223 interior_regions=None, 224 interior_holes=None, 225 poly_geo_reference=None, 226 mesh_geo_reference=None, 227 minimum_triangle_angle=28.0, 228 fail_if_polygons_outside=True, 229 verbose=True): 230 """_create_domain_from_regions - internal function. 231 232 See create_domain_from_regions for documentation. 233 """ 234 235 create_mesh_from_regions(bounding_polygon, 236 boundary_tags, 237 maximum_triangle_area=maximum_triangle_area, 238 interior_regions=interior_regions, 239 filename=mesh_filename, 240 interior_holes=interior_holes, 241 poly_geo_reference=poly_geo_reference, 242 mesh_geo_reference=mesh_geo_reference, 243 minimum_triangle_angle=minimum_triangle_angle, 244 fail_if_polygons_outside=fail_if_polygons_outside, 245 use_cache=False, 246 verbose=verbose) 247 248 domain = Domain(mesh_filename, use_cache=False, verbose=verbose) 249 250 251 return domain 252 253 254 255 256 257 258 # 120 259 # Shallow water domain 121 260 #--------------------- -
anuga_work/production/patong/project.py
r6168 r6178 48 48 finaltime=15000 # final time for simulation 15000 49 49 50 setup=' final' # Final can be replaced with trial or basic.50 setup='trial' # Final can be replaced with trial or basic. 51 51 # Either will result in a coarser mesh that will allow a 52 52 # faster, but less accurate, simulation. -
anuga_work/production/patong/run_patong.py
r6175 r6178 23 23 from os import sep 24 24 import os 25 from os.path import dirname, basename 26 from os import mkdir, access, F_OK 27 from shutil import copy 25 from os.path import dirname 28 26 import time 29 import sys30 27 31 28 # Related major packages … … 35 32 from anuga.shallow_water import Reflective_boundary 36 33 from anuga.shallow_water import Field_boundary 37 from Numeric import allclose38 34 from anuga.shallow_water.data_manager import export_grid, create_sts_boundary, csv2building_polygons 39 from anuga. pmesh.mesh_interface import create_mesh_from_regions35 from anuga.shallow_water import create_domain_from_regions 40 36 from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters 41 37 from anuga.caching import myhash … … 51 47 myid = 0 52 48 53 def run_model(**kwargs):54 55 #------------------------------------------------------------------------------56 # Copy scripts to time stamped output directory and capture screen57 # output to file58 #------------------------------------------------------------------------------59 60 #copy script must be before screen_catcher61 62 print 'output_dir',kwargs['output_dir']63 64 copy_code_files(kwargs['output_dir'],__file__,65 dirname(project.__file__)+sep+ project.__name__+'.py' )66 67 store_parameters(**kwargs)68 69 start_screen_catcher(kwargs['output_dir'], myid, numprocs)70 49 71 50 72 #----------------------------------------------------------------------- 73 # Domain definitions 74 #----------------------------------------------------------------------- 51 #------------------------------------------------------------------------------ 52 # Copy scripts to time stamped output directory and capture screen 53 # output to file 54 #------------------------------------------------------------------------------ 75 55 76 # Read in boundary from ordered sts file 77 urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name)) 56 #copy script must be before screen_catcher 78 57 79 # Reading the landward defined points, this incorporates the original clipping 80 # polygon minus the 100m contour 81 landward_bounding_polygon = read_polygon(project.landward_dir) 58 output_dir = project.output_run_time_dir 59 print 'output_dir', output_dir 82 60 83 # Combine sts polyline with landward points 84 bounding_polygon = urs_bounding_polygon + landward_bounding_polygon 85 86 # counting segments 87 N = len(urs_bounding_polygon)-1 61 #copy_code_files(output_dir,__file__, 62 # dirname(project.__file__)+sep+ project.__name__+'.py' ) 88 63 89 # boundary tags refer to project.landward 4 points equals 5 segments start at N 90 boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)} 91 92 #-------------------------------------------------------------------------- 93 # Create the triangular mesh based on overall clipping polygon with a tagged 94 # boundary and interior regions defined in project.py along with 95 # resolutions (maximal area of per triangle) for each polygon 96 #-------------------------------------------------------------------------- 97 98 # IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it 99 # causes problems with the ability to cache set quantity which takes alot of times 100 101 print 'start create mesh from regions' 102 103 create_mesh_from_regions(bounding_polygon, 104 boundary_tags=boundary_tags, 105 maximum_triangle_area=project.res_poly_all, 106 interior_regions=project.interior_regions, 107 filename=project.meshes_dir_name, 108 use_cache=True, 109 verbose=True) 110 111 #------------------------------------------------------------------------- 112 # Setup computational domain 113 #------------------------------------------------------------------------- 114 print 'Setup computational domain' 115 116 domain = Domain(project.meshes_dir_name, use_cache=False, verbose=True) 117 print 'memory usage before del domain',mem_usage() 118 119 print domain.statistics() 120 print 'triangles',len(domain) 121 122 kwargs['act_num_trigs']=len(domain) 64 #start_screen_catcher(output_dir, myid, numprocs) 123 65 124 66 125 #------------------------------------------------------------------------- 126 # Setup initial conditions 127 #------------------------------------------------------------------------- 128 print 'Setup initial conditions' 67 #----------------------------------------------------------------------- 68 # Domain definitions 69 #----------------------------------------------------------------------- 129 70 130 # sets the initial stage in the offcoast region only 131 IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'], 132 geo_reference = domain.geo_reference) 133 domain.set_quantity('stage', IC) 134 #domain.set_quantity('stage',kwargs['tide'] ) 135 domain.set_quantity('friction', kwargs['friction']) 136 137 print 'Start Set quantity',kwargs['elevation_file'] 71 # Read in boundary from ordered sts file 72 urs_boundary_name = os.path.join(project.boundaries_dir, 73 project.scenario_name) 74 urs_bounding_polygon = create_sts_boundary(urs_boundary_name) 138 75 139 domain.set_quantity('elevation', 140 filename = kwargs['elevation_file'], 141 use_cache = False, 142 verbose = True, 143 alpha = kwargs['alpha']) 76 # Reading the landward defined points, this incorporates the original clipping 77 # polygon minus the 100m contour 78 landward_bounding_polygon = read_polygon(project.landward_dir) 144 79 145 # Add buildings from file 146 building_polygons, building_heights = csv2building_polygons(project.building_polygon_file) 80 # Combine sts polyline with landward points 81 bounding_polygon = urs_bounding_polygon + landward_bounding_polygon 147 82 148 L = [] 149 for key in building_polygons: 150 poly = building_polygons[key] 151 elev = building_heights[key] 152 L.append((poly, elev)) 153 domain.add_quantity('elevation', Polygon_function(L, default=0.0)) 154 155 print 'Finished Set quantity' 83 # Counting segments 84 N = len(urs_bounding_polygon)-1 156 85 157 ## #------------------------------------------------------ 158 ## # Distribute domain to implement parallelism !!! 159 ## #------------------------------------------------------ 160 ## 161 ## if numprocs > 1: 162 ## domain=distribute(domain) 86 # boundary tags refer to project.landward 4 points equals 5 segments start at N 87 boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)} 163 88 164 #------------------------------------------------------ 165 # Set domain parameters 166 #------------------------------------------------------ 167 print 'domain id', id(domain) 168 domain.set_name(kwargs['scenario_name']) 169 domain.set_datadir(kwargs['output_dir']) 170 domain.set_default_order(2) # Apply second order scheme 171 domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm 172 domain.set_store_vertices_uniquely(False) 173 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) 174 domain.tight_slope_limiters = 1 175 print 'domain id', id(domain) 89 #-------------------------------------------------------------------------- 90 # Create the computational domain based on overall clipping polygon with 91 # a tagged boundary and interior regions defined in project.py along with 92 # resolutions (maximal area of per triangle) for each polygon 93 #-------------------------------------------------------------------------- 176 94 177 #------------------------------------------------------------------------- 178 # Setup boundary conditions179 #-------------------------------------------------------------------------180 print 'Available boundary tags', domain.get_boundary_tags()181 print 'domain id', id(domain)182 183 boundary_urs_out=project.boundaries_dir + sep + project.scenario_name95 domain = create_domain_from_regions(bounding_polygon, 96 boundary_tags=boundary_tags, 97 maximum_triangle_area=project.res_poly_all, 98 interior_regions=project.interior_regions, 99 mesh_filename=project.meshes_dir_name, 100 use_cache=True, 101 verbose=True) 184 102 185 Br = Reflective_boundary(domain) 186 Bd = Dirichlet_boundary([kwargs['tide'],0,0]) 187 188 print 'Available boundary tags', domain.get_boundary_tags() 189 Bf = Field_boundary(boundary_urs_out+'.sts', # Change from file_boundary 190 domain, mean_stage= project.tide, 191 time_thinning=1, 192 time_limit=project.finaltime, 193 default_boundary=Bd, 194 use_cache=True, 195 verbose=True, 196 boundary_polygon=bounding_polygon) 103 print 'memory usage after creation of domain', mem_usage() 104 print domain.statistics() 197 105 198 domain.set_boundary({'back': Br,199 'side': Bd,200 'ocean': Bf})201 106 202 kwargs['input_start_time']=domain.starttime 107 #------------------------------------------------------------------------- 108 # Setup initial conditions 109 #------------------------------------------------------------------------- 110 print 'Setup initial conditions' 203 111 204 print'finish set boundary' 112 # sets the initial stage in the offcoast region only 113 IC = Polygon_function([(project.poly_mainland, 0)], 114 default=project.tide, 115 geo_reference=domain.geo_reference) 116 domain.set_quantity('stage', IC) 117 domain.set_quantity('friction', project.friction) 118 domain.set_quantity('elevation', 119 filename=project.combined_dir_name+'.pts', 120 use_cache=True, 121 verbose=True, 122 alpha=project.alpha) 205 123 206 #---------------------------------------------------------------------------- 207 # Evolve system through time 208 #-------------------------------------------------------------------- 209 t0 = time.time() 124 # Add buildings from file 125 print 'Creating building polygons' 126 building_polygons, building_heights = csv2building_polygons(project.building_polygon_file) 210 127 211 for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime'] 212 ,skip_initial_step = False): 213 domain.write_time() 214 domain.write_boundary_statistics(tags = 'ocean') 128 L = [] 129 for key in building_polygons: 130 poly = building_polygons[key] 131 elev = building_heights[key] 132 L.append((poly, elev)) 215 133 216 # these outputs should be checked with the resultant inundation map 217 x, y = domain.get_maximum_inundation_location() 218 q = domain.get_maximum_inundation_elevation() 219 print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q) 134 buildings = Polygon_function(L, default=0.0) 135 domain.add_quantity('elevation', buildings) 220 136 221 print 'Simulation took %.2f seconds' %(time.time()-t0)222 137 223 #kwargs 'completed' must be added to write the final parameters to file 224 kwargs['completed']=str(time.time()-t0) 225 226 store_parameters(**kwargs) 227 228 229 230 #------------------------------------------------------------- 231 if __name__ == "__main__": 232 233 kwargs={} 234 kwargs['file_name']=project.dir_comment 235 kwargs['finaltime']=project.finaltime 236 kwargs['output_dir']=project.output_run_time_dir 237 kwargs['elevation_file']=project.combined_dir_name+'.pts' 238 kwargs['scenario_name']=project.scenario_name 239 kwargs['tide']=project.tide 240 kwargs['alpha'] = project.alpha 241 kwargs['friction']=project.friction 242 #kwargs['num_cpu']=numprocs 243 #kwargs['host']=project.host 244 #kwargs['starttime']=project.starttime 245 #kwargs['yieldstep']=project.yieldstep 246 #kwargs['user']=project.user 247 #kwargs['time_thinning'] = project.time_thinning 248 249 run_model(**kwargs) 138 #------------------------------------------------------ 139 # Distribute domain to implement parallelism !!! 140 #------------------------------------------------------ 141 142 # if numprocs > 1: 143 # domain=distribute(domain) 144 145 #------------------------------------------------------ 146 # Set domain parameters 147 #------------------------------------------------------ 148 domain.set_name(project.scenario_name) 149 domain.set_datadir(output_dir) 150 domain.set_default_order(2) # Apply second order scheme 151 domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm 152 153 #------------------------------------------------------------------------- 154 # Setup boundary conditions 155 #------------------------------------------------------------------------- 156 print 'Set boundary conditions' 157 158 Br = Reflective_boundary(domain) 159 Bd = Dirichlet_boundary([project.tide,0,0]) 160 Bf = Field_boundary(urs_boundary_name+'.sts', 161 domain, 162 mean_stage= project.tide, 163 time_thinning=project.time_thinning, 164 default_boundary=Bd, 165 use_cache=True, 166 verbose=True, 167 boundary_polygon=bounding_polygon) 168 169 domain.set_boundary({'back': Br, 170 'side': Bd, 171 'ocean': Bf}) 172 173 174 #---------------------------------------------------------------------------- 175 # Evolve system through time 176 #-------------------------------------------------------------------- 177 t0 = time.time() 178 179 for t in domain.evolve(yieldstep=project.yieldstep, 180 finaltime=project.finaltime) 181 print domain.timestepping_statistics() 182 print domain.boundary_statistics(tags='ocean') 183 184 print 'Simulation took %.2f seconds' %(time.time()-t0) 250 185 251 186
Note: See TracChangeset
for help on using the changeset viewer.