Changeset 5084
- Timestamp:
- Feb 27, 2008, 10:02:55 AM (17 years ago)
- Location:
- anuga_work/debug/bad_time_step
- Files:
-
- 3 edited
- 2 moved
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/debug/bad_time_step/create_mesh.py
r5073 r5084 14 14 15 15 """ 16 16 17 #Basic geometry 17 18 18 xright = 19.0 19 19 ybottom = 0 -
anuga_work/debug/bad_time_step/project.py
r5075 r5084 12 12 from anuga.utilities.system_tools import get_user_name 13 13 14 boundary_file = 'boundary_mellow.tsm' 14 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ 15 import File_boundary_time 16 boundary_file = 'boundary.tms' 17 #boundary_file = 'boundary_mellow.tms' 15 18 16 19 class Project: -
anuga_work/debug/bad_time_step/run_dam.py
r5079 r5084 6 6 7 7 Duncan Gray, GA - 2007 8 9 10 11 8 """ 12 9 … … 16 13 #---------------------------------------------------------------------------- 17 14 18 # Standard modules19 import time 20 from time import localtime, strftime21 import sys 22 from shutil import copy23 from os import path, sep24 from os.path import dirname #, basename15 # Related major packages 16 from anuga.shallow_water import Domain 17 from anuga.shallow_water import Reflective_boundary 18 from anuga.shallow_water import Dirichlet_boundary 19 from anuga.shallow_water import Time_boundary 20 from anuga.shallow_water import File_boundary 21 from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary 25 22 26 23 27 # Related major packages 28 from anuga.shallow_water import Domain, Reflective_boundary, \ 29 Dirichlet_boundary, Time_boundary, \ 30 File_boundary, \ 31 Transmissive_Momentum_Set_Stage_boundary 32 from anuga.fit_interpolate.interpolate import interpolate_sww2csv 33 from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \ 34 copy_code_files, file_function 35 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ 36 import File_boundary_time 24 from anuga.abstract_2d_finite_volumes.util import file_function 25 #from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 26 from anuga.pmesh.mesh_interface import create_mesh_from_regions 27 28 37 29 38 30 # Scenario specific imports 39 import project # Definition of file names and polygons40 31 import create_mesh 32 41 33 42 34 def elevation_function(x,y): … … 52 44 return z 53 45 54 def main(friction=0.01, outputdir_name=None, is_trial_run=False): 55 basename = 'zz' + str(friction) 56 if is_trial_run is True: 57 outputdir_name += '_test' 58 yieldstep = 0.1 59 finaltime = 15. 60 maximum_triangle_area=0.01 61 #maximum_triangle_area=0.1 #OK 62 else: 63 yieldstep = 0.02 64 finaltime = 15.1 65 maximum_triangle_area=0.0001 66 67 maximum_triangle_area=0.001 68 69 pro_instance = project.Project([], 70 outputdir_name=outputdir_name, 71 home='.') 72 print 'The output dir is', pro_instance.outputdir 73 copy_code_files(pro_instance.outputdir,__file__, 74 dirname(project.__file__) \ 75 + sep + project.__name__+'.py') 76 copy (pro_instance.codedir + 'run_dam.py', 77 pro_instance.outputdir + 'run_dam.py') 78 copy (pro_instance.codedir + 'create_mesh.py', 79 pro_instance.outputdir + 'create_mesh.py') 46 47 #------------------------------------------------------------------------- 48 # Model parameters 49 #------------------------------------------------------------------------- 50 friction=0.01 51 maximum_triangle_area=0.01 52 #maximum_triangle_area=0.1 #OK 53 54 basename = 'debug_boundary' 55 mesh_filename = basename + '.msh' 56 57 boundary_file = 'boundary.tms' 58 #boundary_file = 'boundary_mellow.tms' 59 60 61 62 #------------------------------------------------------------------------- 63 # Create the triangular mesh 64 #------------------------------------------------------------------------- 65 66 # Original unstructured mesh () 67 #create_mesh.generate(mesh_filename, 68 # maximum_triangle_area=maximum_triangle_area) 69 70 71 # Simple mesh with same dimensions 72 73 xright = 19.0 74 ybottom = 0 75 ytop = 0.45 76 xleft = 0.0 77 #xslope = 4.0 78 79 #Outline 80 point_sw = [xleft, ybottom] 81 point_se = [xright, ybottom] 82 point_nw = [xleft, ytop] 83 point_ne = [xright, ytop] 84 85 86 create_mesh_from_regions([point_sw, 87 point_se, 88 point_ne, 89 point_nw], 90 {'wall': [0,1,2], 'wave': [3]}, 91 maximum_triangle_area=maximum_triangle_area, 92 filename=mesh_filename) 93 94 95 #TODO 96 #points, elements, boundary = rectangular_cross(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)): 97 98 99 100 #------------------------------------------------------------------------- 101 # Setup computational domain 102 #------------------------------------------------------------------------- 103 domain = Domain(mesh_filename, use_cache = False, verbose = True) 104 print domain.statistics() 105 106 domain.set_name(basename) 107 domain.set_datadir('.') 108 109 110 #------------------------------------------------------------------------- 111 # Setup initial conditions 112 #------------------------------------------------------------------------- 113 114 domain.set_quantity('stage', 0.06) 115 domain.set_quantity('friction', friction) 116 domain.set_quantity('elevation', elevation_function) 117 118 119 print 'Available boundary tags', domain.get_boundary_tags() 120 121 # Create boundary function from timeseries provided in file 122 function = file_function(boundary_file, domain, verbose=True) 123 124 125 Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) 126 Br = Reflective_boundary(domain) 127 domain.set_boundary( {'wall': Br, 'wave': Bts} ) 128 129 130 #------------------------------------------------------------------------- 131 # Evolve system through time 132 #------------------------------------------------------------------------- 133 134 for t in domain.evolve(0.1, 15): 135 print domain.timestepping_statistics(track_speeds=False) 136 print domain.boundary_statistics(tags=['wave']) 80 137 81 mesh_filename = pro_instance.meshdir + basename + '.msh'138 print 'finished' 82 139 83 #--------------------------------------------------------------------------84 # Copy scripts to output directory and capture screen85 # output to file86 #--------------------------------------------------------------------------87 88 # creates copy of code in output dir89 if is_trial_run is False:90 start_screen_catcher(pro_instance.outputdir, rank, pypar.size())91 92 print 'USER: ', pro_instance.user93 94 95 #-------------------------------------------------------------------------96 # Create the triangular mesh97 #-------------------------------------------------------------------------98 99 # this creates the mesh100 #gate_position = 12.0101 create_mesh.generate(mesh_filename,102 maximum_triangle_area=maximum_triangle_area)103 104 head,tail = path.split(mesh_filename)105 copy (mesh_filename,106 pro_instance.outputdir + tail )107 #-------------------------------------------------------------------------108 # Setup computational domain109 #-------------------------------------------------------------------------110 domain = Domain(mesh_filename, use_cache = False, verbose = True)111 112 113 print 'Number of triangles = ', len(domain)114 print 'The extent is ', domain.get_extent()115 print domain.statistics()116 117 118 domain.set_name(basename)119 domain.set_datadir(pro_instance.outputdir)120 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])121 domain.set_minimum_storable_height(0.001)122 #domain.set_store_vertices_uniquely(True) # for writting to sww123 124 #-------------------------------------------------------------------------125 # Setup initial conditions126 #-------------------------------------------------------------------------127 128 domain.set_quantity('stage', 0.06)129 domain.set_quantity('friction', friction)130 domain.set_quantity('elevation', elevation_function)131 132 133 print 'Available boundary tags', domain.get_boundary_tags()134 135 # Create boundary function from timeseries provided in file136 function = file_function(project.boundary_file, domain, verbose=True)137 138 #print dir(function)139 #print function.time, function.quantities_range140 #import sys; sys.exit()141 142 Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)143 144 Br = Reflective_boundary(domain)145 #Bd = Dirichlet_boundary([0.3,0,0])146 #Bts = Time_boundary(domain, function)147 domain.set_boundary( {'wall': Br, 'wave': Bts} )148 149 #-------------------------------------------------------------------------150 # Evolve system through time151 #-------------------------------------------------------------------------152 t0 = time.time()153 154 for t in domain.evolve(yieldstep, finaltime):155 print domain.timestepping_statistics(track_speeds=False)156 print domain.boundary_statistics(tags=['wave'])157 158 print 'That took %.2f seconds' %(time.time()-t0)159 print 'finished'160 161 points = [[2.8,0.225], #-1.8m from SWL162 [5.1,0.225], #0.5m from SWL163 [6.6,0.225], #2m from SWL164 [6.95,0.255], #2.35m from SWL165 [7.6,0.255], #3m from SWL166 [8.2,0.255], #3.5m from SWL167 [9.2,0.255] #4.5m from SWL168 ]169 170 # points = [[gate_position - 0.65,0.2],171 # [gate_position - 0.55,0.2],172 # [gate_position - 0.45,0.2],173 # [gate_position - 0.35,0.2],174 # [gate_position - 0.25,0.2]175 # ]176 177 #-------------------------------------------------------------------------178 # Calculate gauge info179 #-------------------------------------------------------------------------180 181 if False:182 interpolate_sww2csv(pro_instance.outputdir + basename +".sww",183 points,184 pro_instance.outputdir + "depth_manning_"+str(friction)+".csv",185 pro_instance.outputdir + "velocity_x.csv",186 pro_instance.outputdir + "velocity_y.csv")187 188 return pro_instance189 190 #-------------------------------------------------------------191 if __name__ == "__main__":192 main( is_trial_run = True,193 outputdir_name='Hinwood_draft')
Note: See TracChangeset
for help on using the changeset viewer.