Changeset 6889
- Timestamp:
- Apr 23, 2009, 5:13:30 PM (16 years ago)
- Location:
- anuga_core
- Files:
-
- 6 edited
-
documentation/user_manual/demos/cairns/project.py (modified) (3 diffs)
-
documentation/user_manual/demos/cairns/runcairns.py (modified) (8 diffs)
-
documentation/user_manual/demos/channel2.py (modified) (1 diff)
-
documentation/user_manual/demos/channel3.py (modified) (1 diff)
-
documentation/user_manual/demos/runup.py (modified) (5 diffs)
-
source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/documentation/user_manual/demos/cairns/project.py
r4932 r6889 1 # -*- coding: cp1252 -*-2 1 """Common filenames and locations for topographic data, meshes and outputs. 3 2 """ 4 3 5 from os import sep, environ, getenv, getcwd6 from os.path import expanduser7 import sys8 from time import localtime, strftime, gmtime9 4 from anuga.utilities.polygon import read_polygon, plot_polygons, \ 10 5 polygon_area, is_inside_polygon 11 6 12 ############################### 7 8 #------------------------------------------------------------------------------ 9 # Define scenario as either slide or fixed_wave. 10 #------------------------------------------------------------------------------ 11 #scenario = 'slide' 12 scenario = 'fixed_wave' 13 14 15 #------------------------------------------------------------------------------ 16 # Filenames 17 #------------------------------------------------------------------------------ 18 demname = 'cairns' 19 meshname = demname + '.msh' 20 21 # Filename for locations where timeseries are to be produced 22 gauge_filename = 'gauges.csv' 23 24 25 #------------------------------------------------------------------------------ 13 26 # Domain definitions 14 # ##############################27 #------------------------------------------------------------------------------ 15 28 16 29 # bounding polygon for study area 17 30 bounding_polygon = read_polygon('extent.csv') 18 31 19 print 'Area of bounding polygon', polygon_area(bounding_polygon)/1000000.0 32 A = polygon_area(bounding_polygon)/1000000.0 33 print 'Area of bounding polygon = %.2f km^2' % A 20 34 21 ############################### 35 36 #------------------------------------------------------------------------------ 22 37 # Interior region definitions 23 # ##############################38 #------------------------------------------------------------------------------ 24 39 25 # interior polygons40 # Read interior polygons 26 41 poly_cairns = read_polygon('cairns.csv') 27 42 poly_island0 = read_polygon('islands.csv') … … 31 46 poly_shallow = read_polygon('shallow.csv') 32 47 33 plot_polygons([bounding_polygon,poly_cairns,poly_island0,poly_island1,\ 34 poly_island2,poly_island3,poly_shallow],\ 35 style='boundingpoly',verbose=False) 48 # Optionally plot points making up these polygons 49 #plot_polygons([bounding_polygon,poly_cairns,poly_island0,poly_island1,\ 50 # poly_island2,poly_island3,poly_shallow],\ 51 # style='boundingpoly',verbose=False) 36 52 37 53 38 ###################################################################39 # Clipping regions for export to asc and regions for clipping data40 ###################################################################41 54 42 # exporting asc grid 55 # Define resolutions (max area per triangle) for each polygon 56 default_res = 10000000 # Background resolution 57 islands_res = 100000 58 cairns_res = 100000 59 shallow_res = 500000 60 61 # Define list of interior regions with associated resolutions 62 interior_regions = [[poly_cairns, cairns_res], 63 [poly_island0, islands_res], 64 [poly_island1, islands_res], 65 [poly_island2, islands_res], 66 [poly_island3, islands_res], 67 [poly_shallow, shallow_res]] 68 69 70 71 #------------------------------------------------------------------------------ 72 # Data for exporting ascii grid 73 #------------------------------------------------------------------------------ 43 74 eastingmin = 363000 44 75 eastingmax = 418000 … … 46 77 northingmax = 8145700 47 78 48 49 slide_origin = [451871, 8128376] # move onto the continental shelf, depth = 500 79 #------------------------------------------------------------------------------ 80 # Data for landslide 81 #------------------------------------------------------------------------------ 82 slide_origin = [451871, 8128376] # Assume to be on continental shelf 50 83 slide_depth = 500. 51 84 52 gauge_filename = 'gauges.csv'53 85 54 86 87 88 -
anuga_core/documentation/user_manual/demos/cairns/runcairns.py
r4869 r6889 9 9 the elevation data and a tsunami wave generated by a submarine mass failure. 10 10 11 Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton and 12 Nick Bartzis, GA - 2006 11 Geoscience Australia, 2004-present 13 12 """ 14 13 … … 23 22 24 23 # Related major packages 25 from anuga. shallow_water import Domain26 from anuga. shallow_waterimport Reflective_boundary27 from anuga. shallow_waterimport Dirichlet_boundary28 from anuga. shallow_waterimport Time_boundary29 from anuga. shallow_waterimport File_boundary30 from anuga. shallow_water import Transmissive_Momentum_Set_Stage_boundary31 from anuga.pmesh.mesh_interface import create_mesh_from_regions 24 from anuga.interface import create_domain_from_regions 25 from anuga.interface import Reflective_boundary 26 from anuga.interface import Dirichlet_boundary 27 from anuga.interface import Time_boundary 28 from anuga.interface import File_boundary 29 from anuga.interface import Transmissive_stage_zero_momentum_boundary 30 32 31 from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf 33 32 from anuga.shallow_water.data_manager import dem2pts 34 33 34 from anuga.shallow_water.smf import slide_tsunami 35 35 36 # Application specific imports 36 37 import project # Definition of file names and polygons 37 38 39 #------------------------------------------------------------------------------40 # Define scenario as either slide or fixed_wave.41 #------------------------------------------------------------------------------42 scenario = 'slide'43 #scenario = 'fixed_wave'44 45 if os.access(scenario, os.F_OK) == 0:46 os.mkdir(scenario)47 basename = scenario + 'source'48 38 49 39 … … 53 43 #------------------------------------------------------------------------------ 54 44 55 # Filenames56 dem_name = 'cairns'57 meshname = 'cairns.msh'58 59 45 # Create DEM from asc data 60 convert_dem_from_ascii2netcdf( dem_name, use_cache=True, verbose=True)46 convert_dem_from_ascii2netcdf(project.demname, use_cache=True, verbose=True) 61 47 62 48 # Create pts file for onshore DEM 63 dem2pts( dem_name, use_cache=True, verbose=True)49 dem2pts(project.demname, use_cache=True, verbose=True) 64 50 65 51 66 52 #------------------------------------------------------------------------------ 67 # Create the triangular mesh based on overall clipping polygon with a tagged68 # boundary and interior regions defined in project.py along with69 # resolutions (maximal area of per triangle) for each polygon53 # Create the triangular mesh and domain based on 54 # overall clipping polygon with a tagged 55 # boundary and interior regions as defined in project.py 70 56 #------------------------------------------------------------------------------ 71 57 72 remainder_res = 10000000 73 islands_res = 100000 74 cairns_res = 100000 75 shallow_res = 500000 76 interior_regions = [[project.poly_cairns, cairns_res], 77 [project.poly_island0, islands_res], 78 [project.poly_island1, islands_res], 79 [project.poly_island2, islands_res], 80 [project.poly_island3, islands_res], 81 [project.poly_shallow, shallow_res]] 82 83 create_mesh_from_regions(project.bounding_polygon, 84 boundary_tags={'top': [0], 85 'ocean_east': [1], 86 'bottom': [2], 87 'onshore': [3]}, 88 maximum_triangle_area=remainder_res, 89 filename=meshname, 90 interior_regions=interior_regions, 91 use_cache=False, 92 verbose=True) 58 domain = create_domain_from_regions(project.bounding_polygon, 59 boundary_tags={'top': [0], 60 'ocean_east': [1], 61 'bottom': [2], 62 'onshore': [3]}, 63 maximum_triangle_area=project.default_res, 64 mesh_filename=project.meshname, 65 interior_regions=project.interior_regions, 66 use_cache=True, 67 verbose=True) 93 68 94 69 95 #------------------------------------------------------------------------------ 96 # Setup computational domain 97 #------------------------------------------------------------------------------ 98 from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance 99 from anuga.caching import cache 100 101 ##domain = cache(Domain(meshname, use_cache=True, verbose=True) 102 103 domain = cache(pmesh_to_domain_instance, 104 (meshname, Domain), 105 dependencies = [meshname]) 106 70 # Print some stats about mesh and domain 107 71 print 'Number of triangles = ', len(domain) 108 72 print 'The extent is ', domain.get_extent() 109 73 print domain.statistics() 110 111 domain.set_name(basename) 112 domain.set_datadir(scenario) 113 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) 114 domain.set_minimum_storable_height(0.01) 115 116 #print 'domain.tight_slope_limiters', domain.tight_slope_limiters 117 domain.tight_slope_limiters = 0 118 print 'domain.tight_slope_limiters', domain.tight_slope_limiters 74 75 #------------------------------------------------------------------------------ 76 # Setup parameters of computational domain 77 #------------------------------------------------------------------------------ 119 78 120 79 80 domain.set_name('cairns_' + project.scenario) # Name of sww file 81 domain.set_datadir('.') # Store sww output here 82 domain.set_minimum_storable_height(0.01) # Store only depth > 1cm 121 83 122 domain.points_file_block_line_size = 50000123 84 124 85 #------------------------------------------------------------------------------ … … 130 91 domain.set_quantity('friction', 0.0) 131 92 domain.set_quantity('elevation', 132 filename= dem_name + '.pts',93 filename=project.demname + '.pts', 133 94 use_cache=True, 134 95 verbose=True, 135 96 alpha=0.1) 97 136 98 137 99 #------------------------------------------------------------------------------ … … 139 101 #------------------------------------------------------------------------------ 140 102 141 if scenario == 'slide':103 if project.scenario == 'slide': 142 104 # Function for submarine slide 143 from anuga.shallow_water.smf import slide_tsunami144 105 tsunami_source = slide_tsunami(length=35000.0, 145 106 depth=project.slide_depth, … … 159 120 print 'Available boundary tags', domain.get_boundary_tags() 160 121 161 Br = Reflective_boundary(domain)162 Bd = Dirichlet_boundary([tide,0,0])163 122 164 # 60 min square wave starting at 1 min, 50m high 165 if scenario == 'fixed_wave': 166 Bw = Transmissive_Momentum_Set_Stage_boundary(domain = domain, 167 function=lambda t: [(60<t<3660)*50, 0, 0]) 123 Bd = Dirichlet_boundary([tide,0,0]) # Mean water level 124 Bs = Transmissive_stage_zero_momentum_boundary(domain) # Neutral boundary 125 126 127 if project.scenario == 'fixed_wave': 128 # Huge 50m wave starting after 60 seconds and lasting 1 hour. 129 Bw = Time_boundary(domain=domain, 130 function=lambda t: [(60<t<3660)*50, 0, 0]) 131 168 132 domain.set_boundary({'ocean_east': Bw, 169 'bottom': B d,133 'bottom': Bs, 170 134 'onshore': Bd, 171 'top': B d})135 'top': Bs}) 172 136 173 # boundary conditions for slide scenario 174 if scenario == 'slide': 137 if project.scenario == 'slide': 138 # Boundary conditions for slide scenario 175 139 domain.set_boundary({'ocean_east': Bd, 176 140 'bottom': Bd, … … 189 153 from anuga.abstract_2d_finite_volumes.quantity import Quantity 190 154 191 if scenario == 'slide':155 if project.scenario == 'slide': 192 156 193 for t in domain.evolve(yieldstep = 10, finaltime =60):194 domain.write_time()195 domain.write_boundary_statistics(tags = 'ocean_east')157 for t in domain.evolve(yieldstep=10, finaltime=60): 158 print domain.timestepping_statistics() 159 print domain.boundary_statistics(tags='ocean_east') 196 160 197 # add slide161 # Add slide 198 162 thisstagestep = domain.get_quantity('stage') 199 163 if allclose(t, 60): … … 202 166 domain.set_quantity('stage', slide + thisstagestep) 203 167 204 for t in domain.evolve(yieldstep = 10, finaltime =5000,168 for t in domain.evolve(yieldstep=10, finaltime=5000, 205 169 skip_initial_step = True): 206 domain.write_time()207 domain.write_boundary_statistics(tags = 'ocean_east')170 print domain.timestepping_statistics() 171 print domain.boundary_statistics(tags='ocean_east') 208 172 209 if scenario == 'fixed_wave': 173 174 if project.scenario == 'fixed_wave': 210 175 211 # save every two mins leading up to wave approaching land212 for t in domain.evolve(yieldstep = 120, finaltime =5000):213 domain.write_time()214 domain.write_boundary_statistics(tags = 'ocean_east')176 # Save every two mins leading up to wave approaching land 177 for t in domain.evolve(yieldstep=120, finaltime=5000): 178 print domain.timestepping_statistics() 179 print domain.boundary_statistics(tags='ocean_east') 215 180 216 # save every 30 secs as wave starts inundating ashore217 for t in domain.evolve(yieldstep = 10, finaltime =10000,218 skip_initial_step =True):219 domain.write_time()220 domain.write_boundary_statistics(tags ='ocean_east')181 # Save every 30 secs as wave starts inundating ashore 182 for t in domain.evolve(yieldstep=10, finaltime=10000, 183 skip_initial_step=True): 184 print domain.timestepping_statistics() 185 print domain.boundary_statistics(tags='ocean_east') 221 186 222 187 print 'That took %.2f seconds' %(time.time()-t0) -
anuga_core/documentation/user_manual/demos/channel2.py
r3982 r6889 53 53 #------------------------------------------------------------------------------ 54 54 for t in domain.evolve(yieldstep = 0.2, finaltime = 40.0): 55 domain.write_time()55 print domain.timestepping_statistics() 56 56 57 57 if domain.get_quantity('stage').\ -
anuga_core/documentation/user_manual/demos/channel3.py
r5315 r6889 19 19 length = 40. 20 20 width = 5. 21 dx = dy = . 05# Resolution: Length of subdivisions on both axes21 dx = dy = .1 # Resolution: Length of subdivisions on both axes 22 22 23 23 points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy), -
anuga_core/documentation/user_manual/demos/runup.py
r6227 r6889 10 10 #------------------------------------------------------------------------------ 11 11 12 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular 12 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 13 13 from anuga.shallow_water import Domain 14 14 from anuga.shallow_water import Reflective_boundary … … 17 17 from anuga.shallow_water import Transmissive_boundary 18 18 19 from math import sin, pi, exp 19 20 20 21 #------------------------------------------------------------------------------ … … 27 28 domain.set_name('runup') # Output to file runup.sww 28 29 domain.set_datadir('.') # Use current directory for output 29 domain.set_quantities_to_be_stored(['stage',# Store all conserved quantities30 'xmomentum',31 'ymomentum'])32 30 33 31 … … 49 47 #------------------------------------------------------------------------------ 50 48 51 from math import sin, pi, exp52 49 Br = Reflective_boundary(domain) # Solid reflective wall 53 50 Bt = Transmissive_boundary(domain) # Continue all values on boundary … … 65 62 66 63 for t in domain.evolve(yieldstep = 0.1, finaltime = 10.0): 67 domain.write_time()64 print domain.timestepping_statistics() 68 65 69 66 -
anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py
r6475 r6889 82 82 """Time dependent boundary returns values for the 83 83 conserved quantities as a function of time. 84 Must specify domain to get access to model time and a function f(t) 85 which must return conserved quantities as a function time 84 Must specify domain to get access to model time and a function of t 85 which must return conserved quantities as a function time. 86 87 Example: 88 B = Time_boundary(domain, 89 function=lambda t: [(60<t<3660)*2, 0, 0]) 90 91 This will produce a boundary condition with is a 2m high square wave 92 starting 60 seconds into the simulation and lasting one hour. 93 Momentum applied will be 0 at all times. 94 86 95 """ 87 96 … … 89 98 # Transmissive_Momentum_Set_Stage_Boundary (cf posting by rrraman) 90 99 def __init__(self, domain=None, 91 f=None, 100 f=None, # Should be removed and replaced by function below 101 function=None, 92 102 default_boundary=None, 93 103 verbose=False): … … 99 109 self.verbose = verbose 100 110 111 112 # FIXME: Temporary code to deal with both f and function 113 if function is not None and f is not None: 114 raise Exception, 'Specify either function or f to Time_boundary' 115 116 if function is None and f is None: 117 raise Exception, 'You must specify a function to Time_boundary' 118 119 if f is None: 120 f = function 121 ##### 122 101 123 try: 102 124 q = f(0.0)
Note: See TracChangeset
for help on using the changeset viewer.
